Pegasus: sound continuous invariant generation

Continuous invariants are an important component in deductive verification of hybrid and continuous systems. Just like discrete invariants are used to reason about correctness in discrete systems without having to unroll their loops, continuous invariants are used to reason about differential equations without having to solve them. Automatic generation of continuous invariants remains one of the biggest practical challenges to the automation of formal proofs of safety for hybrid systems. There are at present many disparate methods available for generating continuous invariants; however, this wealth of diverse techniques presents a number of challenges, with different methods having different strengths and weaknesses. To address some of these challenges, we develop Pegasus: an automatic continuous invariant generator which allows for combinations of various methods, and integrate it with the KeYmaera X theorem prover for hybrid systems. We describe some of the architectural aspects of this integration, comment on its methods and challenges, and present an experimental evaluation on a suite of benchmarks.


Introduction
Safety verification problems for ordinary differential equations (ODEs) are continuous analogues to Hoare triples: the objective is to show that an ODE cannot evolve out of a designated set of safe states from any of its designated initial states.The role of continuous invariants is broadly analogous to that of inductive invariants for discrete program verification.A continuous invariant is a set of states that can never be left when following the ODE from that set; such an invariant implies safety when it contains all of the initial states and is also a subset of the safe states.The problem of automatically generating invariants (also known as invariant synthesis) is one of the greatest practical challenges in deductive verification of both continuous and discrete systems.In theory, it is even the only challenge for hybrid systems safety [53].
The proliferation of published techniques [5,36,41,57,63,65,75,83,85] for continuous invariant generation-targeting various classes of systems, and having different strengths and weaknesses-presents a challenge: ideally, one does not want to be restricted by the limitations of one particular generation technique (or a small family of techniques).Instead, it is far more desirable to have a framework that accommodates existing generation methods, allows for their combination, and is extensible with new methods as they become available.In this article we (partially) meet the above challenge by developing a single framework which allows us to combine invariant generation methods into novel invariant generation strategies.In our work, we are guided by the following considerations: 1. Specialized invariant generation methods are effective only when the problem falls within their domain; their use must therefore be targeted.2. A combination of invariant generation methods can be more practical than any of the methods considered in isolation.A flexible and reconfigurable mechanism for combining these methods is thus highly desirable.3. Reasoning with automatically generated invariants needs to be done in a sound fashion: any deficiencies in the generation procedure must not compromise the final verification result.
Our interest in automatic invariant generation is motivated by the pressing need to enhance the level of proof automation in deductive verification tools for hybrid systems.In this work we target the KeYmaera X theorem prover [24].
the principles behind this coupling, the techniques used to generate invariants, and the mechanism used for combining them into more powerful invariant generation strategies.An evaluation of this integration on a set of verification benchmarks is presented-with very promising results.The present article extends our previous work [78] with: 1. Extensive coverage of the methods for generating continuous invariants employed by Pegasus (Section 3.3), including extended descriptions of several invariant generation methods, as well as new material on conic abstractions [6] and on the theory and practice of generating rational first integrals for non-linear and linear systems [20,21,28,44,45,72].The extended article also includes a detailed account of the pitfalls and caveats associated with the various invariant generation and checking methods (Sections 3-6).2. New insights on invariant generation strategies based on combining various invariant generation methods (Section 5), including various configuration options for the differential saturation [57] strategy and a new strategy based on differential divide-and-conquer [75].3.An extended benchmark suite with 51 new problems on top of the 90 existing ones (Section 6), together with extended experimental evaluation and analysis of various invariant generation strategy configurations.
Structure of this article.Mathematical preliminaries and definitions are reviewed in 2. Section 3 recalls the problem of continuous invariant checking and describes our architecture for sound invariant checking and generation.Sections 3.3 and 5 describe some of the methods employed by Pegasus for generating continuous invariants, along with mechanisms for their combination.Section 6 presents an empirical evaluation of our integration with KeYmaera X on a suite of verification benchmarks.Section 7 reviews related work and Section 8 discusses the outlook and possible further extensions.Section 9 ends with a summary and concluding remarks.

Preliminaries
Ordinary Differential Equations.An n-dimensional autonomous system of firstorder ODEs has the form: x = f (x), where x = (x 1 , . . ., x n ) ∈ R n is a vector of state variables, x = (x 1 , . . ., x n ) denotes their time-derivatives, i.e. dxi dt for each i = 1, . . ., n, and f (x) = (f 1 (x), . . ., f n (x)) specify the RHS of the equations that these time-derivatives must obey along solutions to the ODEs.Geometrically, such a system of ODEs defines a vector field f : R n → R n , associating to each point x ∈ R n the vector f (x) = (f 1 (x), . . ., f n (x)) ∈ R n specifying in which direction the continuous system evolves at x. Whenever the state of the system is required to be confined within some prescribed set of states Q ⊆ R n , called evolution constraint2 , we will write x = f (x) & Q.
If no evolution constraint is specified, Q is all of R n .A solution to the initial value problem for the system of ODEs x = f (x) with initial value x 0 ∈ R n is a differentiable function x(x 0 , t) : (a, b) → R n defined for all times t ∈ (a, b) ⊆ R ∪ {∞, −∞} where a < 0 < b, and such that x(x 0 , 0) = x 0 and d dt x(x 0 , t) = f (x(x 0 , t)) for all t ∈ (a, b).The Lie derivative of a continuously differentiable function p : R n → R with respect to vector field f is defined as p ≡ n i=1 ∂p ∂xi f i and provably [56,60] equals the time-derivative of p evaluated along the solutions to the system x = f (x).
Semi-algebraic Sets.A set S ⊆ R n is semi-algebraic iff it is characterized by a finite Boolean combination of polynomial equations and inequalities: where p ij ∈ R[x 1 , . . ., x n ] are polynomials.By quantifier elimination, every first-order formula of real arithmetic characterizes a semi-algebraic set and can be put into form (1) (see e.g.Mishra [46,§8.6]).By an abuse of notation, this article uses formulas and the sets they characterize interchangeably.
Continuous Invariants in Verification.Safety specifications for ODEs and hybrid systems can be rigorously verified in formal logics, such as differential dynamic logic (dL ) [52,55,56] as implemented in the KeYmaera X proof assistant [24] and hybrid Hoare logic [40] as implemented in the HHL prover [86].The use of appropriate continuous invariants is key to these verification approaches as they allow the complexities of the continuous dynamics to be handled rigorously even for ODEs without closed-form solutions.For example, the dL formula Init → [x = f (x) & Q] Safe states that the safety property Safe is satisfied throughout the continuous evolution of the system x = f (x) & Q whenever the system begins its evolution from a state satisfying Init.The invariant reasoning principle for verifying such a safety property is given by the following sound rule of inference in dL , with three premisses above the bar and the conclusion below: (Safety) In this rule, the first and third premiss respectively state that the initial set Init is contained within the set I, and that I lies entirely inside the safe set of states Safe.The second premiss states that I is a continuous invariant, i.e.I is maintained throughout the continuous evolution of the system whenever it starts inside I, that is, the following dL formula is true in all states: Thus, the problem of verifying safety properties of ODEs reduces to finding an invariant I that can be proved to satisfy all three premisses.Semantically, a continuous invariant can also be defined as follows: Definition 1 (Continuous invariant) Given a system x = f (x) & Q, the set I ⊆ R n is a continuous invariant iff the following statement holds: For any given set of initial states Init ⊆ R n , a continuous invariant I such that Init ⊆ I provides a sound over-approximation of the states reachable by the system from Init by following the solutions to the ODEs within the domain constraint Q.Indeed, the exact set of states reachable by a continuous system from Init provides the smallest such invariant. 3While Def. 1 above features the solution x(x 0 , t), which may not be available explicitly, a crucial advantage afforded by continuous invariants is the possibility of checking whether a given set is a continuous invariant without computing the solution, i.e. by working directly with the ODEs.

Sound Invariant Checking and Generation
The problem of checking whether a semi-algebraic set I ⊆ R n is a continuous invariant of a polynomial system of ODEs x = f (x) & Q was shown to be decidable by Liu, Zhan, and Zhao [41].This decision procedure, henceforth referred to as LZZ, provides a way of automatically checking continuous invariants (2) by exploiting facts about higher-order Lie derivatives of multivariate polynomials appearing in the syntactic description of I and the Noetherian property of the ring R[x] [27,41]; its implementation requires an algorithm for constructing Gröbner bases [14], as well as a decision procedure for the universal fragment of real arithmetic [68].A logical alternative for invariant checking is provided by the complete dL axiomatization for differential equation invariants [60].Whereas using LZZ results in a yes/no answer to an invariance question (2), dL makes it possible to construct a formal proof of invariance from a small set of ODE axioms [60] whenever the property holds (or a refutation whenever it does not).

Invariant Generation with Template Enumeration
Given a means to perform invariant checking with real arithmetic, an obvious solution to the invariant generation problem (which has been suggested by numerous authors [41,57,80]) involves the method of template enumeration, which yields a theoretically complete semi-algorithm, in the sense that it terminates with a positive answer iff that is possible with the given templates.A template is a parametric formula, such as e.g.
composed from polynomials in the state variables (in this example x, y) with symbolic coefficients (here a 0 ,a 1 ,a 2 ,a 3 ,a 4 ,a 5 and b 0 ,b 1 ,b 2 ), which are interpreted over the reals.All it takes in theory is to exhaustively enumerate parametric templates matching all real arithmetic formulas describing all semialgebraic sets, and use a quantifier elimination algorithm (such as CAD [13]) to identify whether choices for the template parameters exist that meet the required arithmetic constraints.While templates make this British Museum Algorithm-like approach more successful than, e.g.exhaustively enumerating all proofs [31], the method is nevertheless quite impractical for the resulting real arithmetic [54].To appreciate why, let us only remark that quantifier elimination algorithms for real arithmetic used in practice have doubly-exponential time complexity in the number of variables [64].Template enumeration treats every monomial coefficient in the template as a fresh variable, leading to exponentially many real arithmetic variables, which makes this approach highly unscalable.In practice, invariant generation is achieved by using incompletebut more efficient-generation methods.These methods are numerous and vary considerably in their strengths and limitations, creating a wide spectrum of possible trade-offs in performance, the quality, and the form of invariants that one can generate.Effectively navigating this spectrum is an important practical challenge that this article seeks to address.

Soundness: Proof Assistants and Invariant Generation
There are a number of design decisions that can be exercised in how reasoning with continuous invariants is performed within a deductive verification framework.A fundamental design decision is how tightly (i) continuous invariance checking and (ii) continuous invariant generation are to be coupled with the implementation of the prover.This space of design choices is exemplified by the HHL prover and the KeYmaera X prover.
The HHL prover [11,86] implements (i) the LZZ decision procedure for invariant checking and (ii) the method of template enumeration for invariant generation based on real quantifier elimination and Gröbner bases.From the perspective of the HHL prover, these are trusted external oracles for checking the validity of statements about continuous invariance; trusting the output of the HHL prover includes trusting the implementation of its LZZ procedure and the invariant generator (and any arithmetic tool either of them use).
In contrast, KeYmaera X [24] pursues an LCF-style approach, seeking to minimize the soundness-critical code that needs to be trusted in its output.For continuous invariants, it achieves this by (i) checking invariance within the axiomatic framework of dL (rather than trusting external checking procedures) and (ii) accepting conjectured invariants generated from a variety of sources but separately checking the result.Invariant checking in KeYmaera X is automatic [60], which is made possible by the use of specialized proof tactics [23]; these additionally allow it to use a variety of other (incomplete, but computationally inexpensive) methods for proving continuous invariance [27].
Fig. 1: Alternative prover architectures for checking conjectured continuous invariants, i.e. formulas for the form The difference between these two approaches (Fig. 1) is broadly analogous to the use of trusted decision procedures in PVS [17] and oracles in HOL [7,88] on the one hand, and LCF-style proof reconstruction (e.g. in Isabelle [87]) on the other.
Remark 2 KeYmaera X also supports witness checking for the universal fragment of real arithmetic [59] resulting from ODE invariance checking [60].In theory, this leads to a complete LCF-style approach, but in practice, the performance of real arithmetic witness generation is only competitive with secondtier quantifier elimination [59].

Syntactic Representation of Invariants
A subtle issue that arises when interfacing with provers like KeYmaera X or the HHL prover is which terms can be syntactically represented in the prover.The choice of representation limits the kinds of invariants that can be described (or generated), but it is an important consideration for computational efficiency and soundness purposes.For example, Noetherian functions support a sound and complete axiomatization of invariants [60] but can lead to undecidable arithmetic.Rational functions and roots could also be supported [8] but would increase the complexity of the required symbolic computations.For decidability of the invariance and arithmetic questions, this article only considers semi-algebraic invariants, i.e., those built from polynomials as in (1).
A similar issue arises even when restricted to polynomial terms.Naïvely, for maximum flexibility, one would like to describe invariants using polynomials p ∈ R[x] that have arbitrary real-valued coefficients.In practice though, only computable subfields K of R can be effectively represented and used on a computer.Thus, any computational tool must necessarily work with polynomials p ∈ K[x] over some choice of representation for the field of coefficients K. Real algebraic numbers K = Q would work as coefficients, but they increase the complexity of symbolic computations due to the added need to work with polynomial ideal arithmetic for coefficients and can also lead to some subtleties with the nondifferentiability of the resulting root function itself [8].On the other extreme, floating point numbers are computationally efficient but they do not form a field, and would also cause numerical errors that make it harder to obtain sound and exact answers in the end.For these reasons, KeYmaera X works with polynomials p ∈ Q[x] that have rational coefficients. 4his results in fast evaluations and symbolic computations, and a reasonable (although nontrivial) complexity for the resulting real arithmetic validity decision problem.Many invariant generation techniques described in this article implements are fairly general and agnostic to the precise choice of field K. Thus, the rest of this article elides this subtlety and describes the invariant generation algorithms over p ∈ R[x], i.e., with R as the coefficient field.Pegasus is a continuous invariant generator implemented in the Wolfram Language with an interface accessible through both Mathematica and KeYmaera X. 5When KeYmaera X is faced with a continuous safety verification problem that it is unable to prove directly, it automatically invokes Pegasus to help find an appropriate invariant (if possible).KeYmaera X checks all the invariants it is supplied with-including those provided by Pegasus (see Fig. 2).This de-sign ensures that correctness of Pegasus is not integral to the soundness of KeYmaera X.It also presents implementation opportunities for Pegasus:

KeYmaera X
1. Pegasus can freely integrate numerical procedures and heuristic methods while providing best-effort guarantees of correctness.Final correctness checks for the generated invariants are left to the purview of KeYmaera X.6 2. Pegasus records proof hints corresponding to various methods that were used to generate continuous invariants.These hints enable KeYmaera X to build more efficient shortcut proofs of continuous invariance [27].
Pegasus currently implements an array of powerful invariant generation methods, which we describe below, beginning with a large family of related methods that are based on qualitative analysis, which can be best explained using the machinery of discrete abstraction of continuous systems.We first briefly recall the main idea behind this approach.

Exact Discrete Abstraction
Discrete abstraction is the subject of numerous works [1,82,84].Briefly, the steps are: (i) discretize the continuous state space of a system by defining predicates that correspond to discrete states, (ii) compute a (local) transition relation between the discrete states obtained from the previous step, yielding a discrete transition system which abstracts the behaviour of the original continuous system, and finally (iii) compute reachable sets in the discrete abstraction to obtain an over-approximation of the reachable sets of the continuous system.
A discrete abstraction is sound iff the relation computed in step (ii) has a transition between two discrete states whenever there is a corresponding continuous trajectory of the original system between the two neighboring sets corresponding to those discrete states.The abstraction is exact iff these are the only transitions computed in step (ii).Soundness of the discrete abstraction guarantees that any invariant extracted from the discretization corresponds to an invariant for the original system.Fig. 3 illustrates a discretization of a system of ODEs (Fig. 3a), which results in 9 discrete states in a sound and exact abstraction (Fig. 3b).The state space is discretized using predicates built from sign conditions on polynomials, . The discrete states of the abstraction are given by formulas such as S 1 ≡ p 1 < 0 ∧ p 2 = 0, S 2 ≡ p 1 < 0 ∧ p 2 > 0, and so on.
The ability to construct sound and exact discrete abstractions [75] has an important consequence: if an appropriate semi-algebraic continuous invariant I exists at all, it can always be extracted from a discrete abstraction built from discretizing the state space using sign conditions on the polynomials describing I.The problem of (semi-algebraic) invariant generation therefore reduces to finding appropriate polynomials whose sign conditions can yield suitable discrete abstractions and computing reachable states in these abstractions.
Remark 3 Reachable sets (from the initial states) in discrete abstractions are the smallest invariants with respect to ⊆ (set inclusion) that are representable in that abstraction.The smallest invariant is the most informative because it allows one to prove the most safety properties, but it may not be the most useful invariant in practice.In particular, one often wants to work with invariants that have low descriptive complexity and are easy to prove in the formal proof calculus.This leads naturally to consider alternative ways of extracting invariants.Pegasus is able to extract reachable sets of discrete abstractions, but favours less costly techniques, such as differential saturation [57], which often succeed in quickly extracting more conservative invariants.
Finding "good" polynomials that can abstract the system in useful ways and allow proving properties of interest is generally difficult.While abstraction using predicates that are extracted from the verification problem itself can be surprisingly effective, in certain cases useful predicates may not be syntactically extracted from the problem statement.In order to improve the quality of discrete abstractions, Pegasus employs a separate classifier, which extracts features from the verification problem which can then be used to suggest polynomials that are more tailored to the problem at hand.Certain systems have structure that, to a human expert, might suggest an "obvious" choice of good predicates.Below we sketch some basic examples of what is currently possible.

Targeted Qualitative Analysis
As a motivating example, consider the class of one-dimensional ODEs x = f (x), where f ∈ R[x].A standard way of studying qualitative behaviour in these systems is to inspect the graph of the function f (x) [79].Fig. 4 illustrates such a graph of f (x), along with a vector field induced by such a system on the real line.The ODE x = f (x) is at an equilibrium without any motion at points where f (x) = 0.By computing the real roots of the polynomial in the Fig. 4: Qualitative analysis of one-dimensional ODEs x = f (x) right-hand side, i.e the real roots r 1 , . . ., r k ∈ R of f (x), we may form a list of polynomials x − r 1 , . . ., x − r k that can be used for an algebraic decomposition of R into invariant subregions corresponding to real intervals from which an over-approximation of the reachable set can be constructed.Such an algebraic decomposition can be further refined by augmenting the list of polynomials with x − b 1 , . . ., x − b l , where b 1 , . . ., b l ∈ R are the boundary points of the initial set in the safety specification.From this augmented list, one can exactly construct the reachable set of the system by computing the reachable set of the corresponding exact abstraction.
is one-dimensional one can exploit another useful fact: every one-dimensional system is a gradient system, i.e. its motion is generated by a potential function F (x) which can be computed directly by integrating In higher dimensions, the behaviour of linear systems x = Ax with a constant coefficient matrix A can be studied qualitatively by examining the eigenvalues and eigenvectors7 of the matrix A [2]. Pegasus implements methods targeted at linear systems that take advantage of facts such as these to suggest useful abstractions from which invariants can be extracted.The current strategy is similar in spirit to the abstraction methods proposed in the work of Tiwari [81], and works by computing linear forms describing the invariant half-spaces in the state space of linear systems.Briefly, whenever the system matrix A has a real eigenvalue λ ∈ R, by considering an eigenvector v of the transpose matrix A T , which is associated with the eigenvalue λ (recall that the eigenvalues of square matrices A and A T are the same), one may construct the linear form p = v T x, which has the property that [81, §2]: Such linear forms correspond to a special case of so-called Darboux polynomials, which will be described in more detail in Section 4.4.2 and have the property that p > 0, p = 0, and p < 0 define invariant regions in state space (the fact that λ is a real number in the equation p = λp also allows us to construct invariants p ≤ k, where k is an appropriately chosen offset depending on the sign of λ).
Additionally, when all the eigenvalues of the system matrix A have strictly negative real parts, the origin 0 is asymptotically stable and one may construct a Lyapunov function (see [35,Ch. 3]) for the linear system by solving the Lyapunov equation A T P + P A = Q where Q is some given negative definite matrix, and the solution P is positive definite; the quadratic Lyapunov function V for the stable system is given by V (x) = x T P x.Every sub-level set V ≤ k defines a continuous invariant of the system; Fig. 5 (right) illustrates the kind of invariants that can be obtained by using Lyapunov functions together with invariant half-planes to perform abstraction of linear systems.
Example 1 The linear systems in Fig. 5 exhibit different qualitative behaviours.The invariants (shown in blue), demonstrate unreachability of the unsafe states (shown in red) from the initial states (shown as green discs in Fig. 5).
Fig. 5: Automatically generated invariants for linear systems In the leftmost system, all eigenvalues of the system matrix A are purely imaginary.Pegasus generates annular invariants containing the green discs because trajectories of such systems are always elliptical.For the middle system, the (asymptotic) behaviour of its trajectories is determined by the eigenvectors of its system matrix (eigenvalues are real and of opposite sign [2]).Pegasus uses these eigenvectors to generate two invariant half-planes, one for each green disc.Invariant half-planes are also generated for the rightmost system which is asymptotically stable (all real parts of eigenvalues are negative [2]).Pegasus further refines these half-planes with elliptical regions containing the green discs because elliptical regions are invariants for such systems.
In textbook examples of linear systems, one usually finds matrices with eigenvalues and eigenvectors that can be described using rational numbers.However, the situation is not always that nice in practice: eigenvectors of matrices will often feature irrational components, which in the case of the example above leads to invariant half-planes described by linear polynomials with irrational coefficients.It is therefore important to have the means of working with irrational real numbers in the invariant generator and the prover.
In special cases when the verification problem features a purely algebraic initial set, the strongest algebraic invariants for linear systems (i.e. the smallest continuous invariants that can be described by polynomial equalities p = 0) can be computed following the method of Rodríguez-Carbonell & Tiwari [65], which we implement in Pegasus.
Remark 5 Bogomolov et al. [6] introduced a technique called conic abstractions that combines discrete abstraction of affine systems with an associated reachability analysis method.It is particularly powerful for diagonalizable systems, where the authors' experiments suggest it outperforms other tools for linear reachability analysis, like SpaceEx [22].The eponymous idea behind the method is to partition state space into a number of regions (i.e., cones), so that within each cone the change in angle of the vector field (i.e., the twisting) is bounded by a tunable parameter θ.Given any point in the vector field, then, this construction gives a known range of possible slopes for the vector at that point.This is useful information for the subsequent reachability analysis-instead of simply computing the transition relation between neighboring cones, as in Section 4.1, their algorithm uses the twisting information to determine what portions of each cone is potentially reachable from an initial set.We experimented with the conic abstraction method in a limited setting: bounded linear 2-dimensional systems.The major obstacle inhibiting a complete implementation is that Mathematica's native support for polyhedra computations does not quite meet the demands of the algorithm.Our limited implementation is not able to return an exact invariant region-instead, we produce promising visualizations of the invariant generated for two examples from Fig. 5 (see Fig. 6). 8With better support for polyhedra computations, this could be an exciting direction for future implementation by interfacing Pegasus with the Parma Polyhedra Library.

Qualitative Analysis for Non-Linear Systems
General non-linear polynomial systems present a hard class of problems for invariant generation.A number of useful heuristics can be applied to partition the continuous state space, in the hope that the resulting abstraction exhibits a suitable invariant.For example, factorizing the RHS of a differential equation x i = f i (x) yields a set of irreducible polynomial factors p 1 , . . ., p k such that f i = k j=1 p j , which implies that the flow along the curves p j = 0 vanish in the x i direction.This information can be used to cheaply approximate the transition relation in the discrete abstraction and to efficiently extract invariant candidates.For the non-linear ODE in Fig. 3, the discretization polynomials p 1 , p 2 are chosen such that x 2 = 0 and x 1 = 0 on their respective level curves.
Fig. 6: A visualization of our implementation of the conic abstractions method (each example is shown row-wise).The left figures show the generated conic partition into 20 cones (alternating red and blue colours).The right figures show the reachable set computation (in blue) from the same green initial sets as in Fig. 5.These reachable sets, which are also invariant sets, also suffice to show that the ODE never reaches any unsafe states (in red).The method automatically produces finer partitions of the state space (using more cones) when the direction of the vector field changes more drastically.For example, the top partition concentrates several cones around its unstable manifold [12,79] (the line y = 1 6 (1 + √ 13)x), while the bottom partition has more evenly spaced out cones.This yields a useful discrete abstraction e.g. S 4 is an invariant for the resulting abstraction (Fig. 3b).Other useful sources of polynomials for qualitative analysis of non-linear systems are found in e.g. the summands and irreducible factors of the right-hand sides of the ODEs, the Lie derivatives of the factors, and physically meaningful quantities such as the divergence of the system's vector field.
Locally orthogonal linear forms A particularly simple geometric idea can sometimes be profitably applied to generate linear polynomials for abstraction.It may be concisely described as follows: for a system of ODEs x = f (x), which may be non-linear, and given a regular point x 0 ∈ R n with f (x 0 ) = 0, one may use the linear form f (x 0 )•(x−x 0 ) which has the property that its zero set is locally orthogonal to the direction of the vector f (x 0 ).These linear forms only yield local invariants, but sufficiently many of them put together have a good chance of describing invariant regions.In problems where the evolution domain constraint describes a bounded set, it is possible to obtain useful abstractions by choosing a finite number of sample points x 0 within the set and Fig. 7: Abstraction using locally orthogonal linear forms from a set of points partitioning the constraint with the corresponding locally orthogonal linear forms (as illustrated in Fig. 7).Of course, choosing "good" points is the main problem in this method; one possible heuristic is to use evenly-spaced points forming a grid covering the domain constraint.

General-Purpose Methods
Beyond qualitative analysis, Pegasus implements several general-purpose invariant generation techniques which represent restricted, but tractable fragments of the general method of template enumeration.The search for symbolic parameters in these methods is not performed using real quantifier elimination, but instead takes place in more tractable theories.

Polynomial First Integrals
its Lie derivative p with respect to the vector field f is the zero polynomial.First integrals are also known as conserved quantities because they have an important property: their value never changes along the solutions to ODEs; that is to say, for any k ∈ R, p = k is an invariant of the system.For a single first integral p, if one were to use (the signs of) the polynomial p−k to build an abstraction, the abstract state space would not feature any transitions between its states (illustrated in Fig. 8).Thus, one has the freedom to choose values k for which the resulting discrete abstraction suitably partitions the state space.For example, if the initial states lie entirely within p < k and the unsafe ones within p > k, then p < k is an invariant separating those sets.
Pegasus can search for all polynomial first integrals up to a configurable degree bound by solving a system of linear equations whose solutions provide the coefficients of the bounded degree polynomial template for the first integral.This is known as the method of undetermined coefficients and we illustrate the main steps of the method in the following example.
Example 2 (Kasner's equations) Consider the non-linear system of ODEs describing a special case of Einstein's gravitational equations [34] and a polynomial template of maximum degree 2 in the state variables x, y, z: Computing the Lie derivative of this template with respect to the system, i.e. (p a,2 ) = ∂pa,2 ∂x x + ∂pa,2 ∂y y + ∂pa,2 ∂z z gives a parametric polynomial of degree 3: In order to find a first integral, one is required to solve the equation (p a,2 ) = 0, but a polynomial is 0 precisely when all of its coefficients are 0. Thus, by equating all coefficients of the Lie derivative to 0, finding a first integral reduces to solving a linear system of equations over the symbolic coefficients a 0 , . . .,a 9 : Solutions are efficiently found using linear algebra [29,71].In this example, a non-trivial solution yields the polynomial first integral xy + xz + yz.Moreover, all first integrals of degree (up to) two provide concrete instances of the coefficients a and so must correspond to a solution of these equations.
When a polynomial first integral p is computed, one has the freedom of choosing its initial value, which is guaranteed to remain invariant throughout the evolution of the system.In the above example, one may choose any real number k and partition the state space into invariant regions defined by the sign conditions on the polynomial xy +xz +yz −k.Generally, to obtain a tight over-approximation of the reachable set from the initial set of states given in the verification problem, one may choose k by attempting to maximize and minimize the value of the first integral p on the initial set of states within the domain constraint, i.e., one may search for the real values (if they exist): If finite values k max and k min can be obtained, one may generate a continuous invariant Maximizing/minimizing multivariate polynomials subject to semi-algebraic constraints often leads to irrational and real algebraic numbers as exact maxima/minima.Numerical algorithms will yield values that are near-optimal, which may require them to be increased/decreased by some amount before a genuine invariant is constructed as described above.
The set Init ∩ Q may have multiple connected components, and tighter invariants may be obtained from first integrals when the value k is optimized subject to each connected component separately.A cheap way to approximate the connected components is to normalize Init ∩ Q to disjunctive normal form and consider each disjunct as a separate component.
If more than one independent first integral for a system is found, one may construct finer abstractions and generate tighter invariants over-approximating the reachable set.A particularly interesting case is when an n-dimensional system of ODEs has n − 1 functionally independent algebraic first integrals: such a system is said to be algebraically integrable [29,48].In such a system, given any state x 0 ∈ R n , one may evaluate the first integrals p 1 , p 2 , . . ., p n−1 at that state to obtain a continuous invariant given by: If the first integrals are functionally independent, i.e. when the matrix whose columns are formed by the gradients ∇p i ≡ ∂pi ∂x1 , ∂pi ∂x2 , . . .∂pi ∂xn T has full rank at x 0 (i.e. when the vectors ∇p i evaluated at x 0 are linearly independent, see e.g.[48]), the resulting conjunctive formula (locally) describes a 1-dimensional invariant curve in n-dimensional state space and provides the tightest possible algebraic invariant containing x 0 .Local invariance is a necessary criterion, because only local invariants can be global invariants.

Darboux Polynomials
Darboux polynomials were first introduced in 1878 [16] to study integrability of polynomial ODEs.A polynomial p ∈ R[x] is said to be a Darboux polynomial for the system x = f (x) if and only if p = αp for some polynomial α ∈ R[x], which is known as the cofactor of p. Like first integrals, discrete abstractions produced with Darboux polynomials result in three states with no transitions between them (as illustrated in Fig. 8, but with k = 0).Unlike first integrals, only p = 0 is guaranteed to be an invariant of the system.Darboux polynomials have been used for predicate abstraction of continuous systems by Zaki et al. [91], who successfully applied them to verify electrical circuit designs.The problem of generating Darboux polynomials is generally far more difficult than that of generating polynomial first integrals (which represent the special case of Darboux polynomials where the cofactor α is 0 in the equation p = αp).A modification of the method of undetermined coefficients described in the previous section can likewise be applied to search for Darboux polynomials.However, in order to apply this method, one is required to provide a polynomial template for both the Darboux polynomial and for its cofactor.Whenever one has a polynomial system of ODEs x = f (x) in which the maximum polynomial degree of the components f 1 , f 2 , . . ., f n of f is some r ≥ 0, then the maximum possible degree of the Lie derivative of a polynomial p of maximum degree d is given by d + r − 1.Consequently, to search for a Darboux polynomial of maximum degree d, the maximum degree of the cofactor α in the equation p = αp that one needs to consider is given by r − 1.To apply the method of undetermined coefficients, one requires a template p a,d for the Darboux polynomial and a separate template α b,r−1 for the cofactor.The equation to be solved is the following: By expanding the resulting polynomial on the left-hand side and equating each of its monomial coefficients to 0, one obtains a system of equations in the symbolic parameters a, b; however, while this system is linear in the parameter variables a and b considered separately, it is a non-linear system of equations in a, b simultaneously.In practice, solving such a non-linear system is far more computationally expensive than solving a linear system in the case of polynomial first integrals; the naïve method of undetermined coefficients does not provide a practically appealing solution to the problem of Darboux polynomial generation.
Fortunately, automatic generation of Darboux polynomials is an active area of research, owing largely to their importance as a crucial component in the Prelle-Singer method [62] for computing elementary closed-form solutions to ODEs.In order to implement the Prelle-Singer method, more sophisticated algorithms for Darboux polynomial generation have been developed in the computer algebra community, e.g. two algorithms were reported by Man [44].Indeed, in our experiments we have found the algorithms ps 1 and new ps 1 [44] to be much more practical and implement them in Pegasus.

Remark 6
We remark also that several algorithms for generating (what are essentially) Darboux polynomials have more recently been developed within the verification community [36,63,71].However, our experience with some of these procedures has been less positive.The method by Rebiha et al. [63] was in practice found to be very inefficient and incomplete, i.e. unable in general to find all the Darboux polynomials matching a given polynomial template; the technique by Kong et al. [36] is significantly faster but is likewise incomplete.
Determining whether an arbitrary system of polynomial ODEs possesses a Darboux polynomial (and finding a bound on its degree if it does) remains an open problem [92, §4.1].

Rational First Integrals
Beyond polynomial functions, a much larger class of algebraic conserved quantities is that of rational first integrals; these are first integrals represented by rational functions, i.e. functions of the form a b , where a, b are polynomials and b = 0. Searching for this kind of first integral is (unsurprisingly) more difficult than is the case with polynomials; however, it is made possible by exploiting an idea from the seminal work of Darboux (see e.g.Schlomiuk [72]).
Theorem 1 Let p 1 , p 2 , . . ., p k be independent Darboux polynomials for the system x = f (x), with p i = α i p i , where α i is some polynomial cofactor for each i = 1, . . ., k.If the equation has a non-trivial integer solution, i.e. λ = (λ 1 , λ 2 , . . ., λ k ) ∈ Z k \ {0}, then the system has a rational first integral r λ ∈ R(x) given by the product Proof By applying the product rule to compute the Lie derivative r λ , we get 3) it follows that r λ = 0 and r λ is therefore a first integral.
Remark 7 Obviously, if the solution to (3) is such that λ ∈ Z k ≥0 , then the first integral is polynomial; at least one negative component in λ is therefore required in order to construct a non-polynomial rational first integral.We also note that one may search for rational solutions to Eq. (3), i.e. λ ∈ Q k , which will in general result in first integrals featuring radicals.Any such first integral can be turned into a rational first integral by raising it to an integer power corresponding to the least common multiple of the denominators of the rational numbers λ 1 , . . ., λ k .In general, λ 1 , . . ., λ k need not be rational or even real numbers in order for the construction given in Theorem 1 to work; however, irrational solutions lead to first integrals that are not rational functions.
In light of the above theorem, a straightforward procedure for generating rational first integrals (which has previously been suggested by Man [45]) involves (i) generating Darboux polynomials p 1 , p 2 . . ., p k for the system x = f (x), e.g. using an implementation of Man's algorithms [44], and (ii) finding integer (or rational) solutions to the linear system of equations (3) in Theorem 1.If the coefficients of the cofactors α 1 , α 2 , . . ., α k in equation ( 3) are all rational numbers, the problem reduces to solving a system of linear Diophantine equations, for which there exist polynomial-time algorithms.If a rational first integral r λ = a b is found, then a b = l defines an invariant hypersurface for any choice of l ∈ R, assuming b = 0; rewriting this, we get that a − lb = 0 is invariant for any l ∈ R (when b = 0).Example 3 Consider the following non-linear system of ODEs [21]: Using our implementation of Man's algorithm [44], we obtain the following list of Darboux polynomials in under one second of computation time: Solving Eq. 3 in Theorem 1, we obtain the solution (λ 1 , λ 2 , λ 3 ) = (2, 1, −1), from which we obtain the rational first integral (illustrated in Fig. 9) .
Remark 8 Before attempting to search for algebraic first integrals (whether polynomials or rational functions) it is helpful to have static criteria that determine whether such first integrals can arise in a given system of ODEs.Criteria for non-existence of various kinds of first integrals have been studied by numerous authors (notably by Poincaré [92, §7.2]) and typically make use of the linearization x = Ax of the system x = f (x) around a point of equilibrium (i.e. a point x * where f (x * ) = 0).In particular, a sufficient criterion for non-existence of rational first integrals in non-linear systems of ODEs is given by Shi [73, Theorem 1]; it requires that the eigenvalues λ 1 , . . ., λ n of the matrix A are such that k 1 λ 1 + • • • + k n λ n = 0 does not have a non-trivial integer solution (k 1 , . . ., k n ) ∈ Z n \ 0. A similar criterion, which furthermore accounts for repeated eigenvalues, is given by Goriely [29, Ch. 5, Prop.5.5].
Combining Darboux Polynomials and Rational First Integrals.As a first hint of its flexibility for combining invariant generation methods, Pegasus implements rational first integral generation by combining several ideas described thus far in Section 3.3 as follows.This flexibility is further exploited in the discussion of strategies in Section 5.
1. Compute a list of Darboux polynomials p 1 , . . ., p n of some maximum polynomial degree d using generation methods from Section 4.4.2. 2. Abstract the state space into sign invariant cells using those polynomials, e.g., S 1 ≡ p 1 < 0 ∧ p 2 = 0, S 2 ≡ p 1 < 0 ∧ p 2 > 0, S 3 ≡ p 1 < 0 ∧ p 2 < 0, etc., as described in Section 4.1.Notably, the resulting abstraction has no transitions between its discrete states, as illustrated in Fig. 8. 3. Prune away those invariant cells that do not intersect the initial set of states, e.g., delete S 1 if Init ∩S 1 = ∅ since S 1 is then unreachable.Similarly, prune away cells that do not intersect the unsafe set, e.g., delete S 2 if Unsafe ∩ S 2 = ∅ because no initial states in S 2 can reach the unsafe set.4. The remaining unpruned conflict cells, say S 3 , define new invariant generation subproblems, where the original domain constraint Q is restricted to Q ∧ S 3 .Each of the Darboux polynomials are sign-invariant in these cells; moreover, those Darboux polynomials that are sign-definite (either strictly positive or negative) in each cell, e.g.p 1 , p 2 with domain constraint p 1 < 0 ∧ p 2 > 0 for S 3 , can be used to compute rational first integrals r λ (following Theorem 1).The denominator of r λ is guaranteed to be a product of (powers of) sign-definite polynomials so these rational functions are always defined within each conflict cell. 5. Using their respective rational first integrals r λ , refine each conflict cell by maximizing and minimizing the values of r λ to obtain invariant sub-level sets k min ≤ r λ ∧ r λ ≤ k max over the initial set (restricted to that cell), as described in Section 4.4.1.6.If conflict cells remain, increase the polynomial degree d and go to step 1.
Rational First Integrals of Linear Systems.In the case of linear systems of ODEs x = Ax, more efficient methods exist that allow us to directly construct rational first integrals from the eigenvalues and eigenvectors of the system matrix A. These explicit constructions are described, e.g. in the work of Gorbuzov & Pranevich [28] and Falconi & Llibre [20]; in Pegasus, we implement and deploy techniques described in the former.
It is instructive to compare the results obtained by Lafferriere, Pappas and Yovine [38] (which state that semi-algebraic reachable sets of linear ODEs x = Ax can be constructed from semi-algebraic initial sets in cases when A is diagonalizable and all of its eigenvalues are rational) to analogous results independently obtained in the study of integrability of linear systems, e.g.[28,Property 1.1], which states that a linear system x = Ax has a basis of rational first integrals (i.e. is algebraically integrable) whenever the eigenvalues of A are rational and of multiplicity 1.Indeed, such a basis of rational first integrals enables one to construct reachable sets described by polynomials.

Barrier Certificates
The method of barrier certificates is a popular Lyapunov-like technique for safety verification of continuous and hybrid systems [61].Barrier certificates are differentiable functions p that define an invariant region p ≤ 0 which separates the initial states (wholly contained within p ≤ 0) from the unsafe states (wholly contained within p > 0).In order to ensure continuous invariance of the region defined by p ≤ 0, the Lie derivative p of the barrier certificate needs to satisfy certain criteria; differences in these criteria give rise to a number of variations of barrier certificates in the literature.The original work by Prajna and Jadbabaie [61] introduced convex barrier certificates, which employ the differential inequality p ≤ 0 to guarantee invariance of p ≤ 0 under the flow of the system.Later work by Kong et al. [37] introduced so-called exponential-type barrier certificates, which provide a generalization employing the differential inequality p ≤ λp, where λ ∈ R; this was generalized further yet in the work of Dai et al. [15], who introduced general barrier certificates employing the differential inequality p ≤ ω(p), where ω is a specifically crafted scalar function to guarantee invariance of p ≤ 0. All of the above developments are fundamentally based on the classical notion of comparison systems [66, Ch II, §3, Ch.IX] in the theory of ODEs.A unified understanding of these generalizations is described in our earlier work [77], in which we introduced a further generalization of the barrier certificate framework: vector barrier certificates, employing multidimensional comparison systems in a way analogous to vector Lyapunov functions introduced by Bellman [4].
Barrier certificates are practically interesting because one may apply the method of undetermined coefficients to automatically search for them using tractable techniques: either sum-of-squares programming (SOS) [61] or linear programming (LP) [89].Pegasus is able to search for convex [61], exponentialtype [37], and vector barrier certificates [77] using both SOS and LP techniques.However, the resulting barrier certificates often suffer from numerical inaccuracies arising from the use of semi-definite solvers and interior point methods [67].Pegasus currently uses a simple rounding heuristic on the numerical result and explicitly checks invariance for the resulting (exact) barrier certificate candidates using real quantifier elimination.An example barrier certificate generation technique implemented in Pegasus, and an illustration of its numerical issues is given next.
Example 4 Consider the safety verification problem illustrated in Fig. 10 (left).The task is to generate an invariant showing that ODE solutions starting within the initial set Init (in green) do not enter the unsafe set Unsafe (in red).A candidate continuous invariant p ≤ 0 (shown in blue in Fig. 10, left) is found using numerical barrier certificate generation techniques.
The (exponential-type) barrier certificate p is generated from a polynomial template p a,d of degree d over variables x, y, by solving (and then substituting) for appropriate concrete values of the template coefficients a.For clarity below, the notation p a,d is used in steps where the generation algorithm produces constraints on the coefficients a, while p always refers to the final, generated barrier certificate.Logically, it suffices to find real values for a so that the following formulas are simultaneously valid: Constraints ( 4) and ( 5) ensure that the generated barrier separates the initial set from the unsafe set, e.g., in Fig. 10 (left) the green initial region is wholly contained in the blue candidate invariant region p ≤ 0, while the red unsafe region lies entirely outside.Constraint (6) ensures that the sub-level set p ≤ 0 is a continuous invariant, intuitively, the vector field points "inwards" along the boundary of p ≤ 0 (blue region in Fig. 10), so it is impossible to flow from within p ≤ 0 to p > 0. A more general version of these constraints, and a soundness proof, is available elsewhere [37].
Sum-of-squares (SOS) programming [49] provides a tractable way of solving for coefficients a. Suppose that Init, Unsafe are described with polynomial inequalities Init ≡ a i=1 I i ≥ 0, Unsafe ≡ b i=1 U i ≥ 0. The inequalities (4)-( 6) are respectively implied by the following SOS inequalities, where ε > 0 is a small positive constant and σ Ii , σ Ui are template SOS polynomials [49]: Sum-of-squares solvers, such as SOSTOOLS [49], witness the inequalities ( 7)-( 9) by finding an SOS representation for their LHS.For example, a set of polynomials g 1 , . . ., g n satisfying the polynomial identity −p a,d − a i=1 σ Ii I i = n i=1 g 2 i proves (7) because the RHS of this inequality is a sum-of-squares , which is non-negative.These polynomial identities are found efficiently by semidefinite programming [51], which is also where numerical solvers are used.In practice, Pegasus loops through a range of values for the parameters d, λ, ε as well as the degrees of the SOS polynomials σ Ii , σ Ui and attempts to solve these constraints for each concrete choice of parameters.
The efficiency of numerical solvers is also a drawback because the generated coefficients a need not strictly satisfy all the required constraints.This is why Pegasus (and KeYmaera X) treats the generated barrier certificate p only as a candidate invariant and performs additional arithmetical checks to ensure that the constraints are truly met.As a cautionary example, Fig. 10 (left) rather misleadingly suggests that p ≤ 0 is an invariant within its small plot domain.Indeed, Fig. 10 (right) is a zoomed out version of the same plot which shows that the constraint (6) fails to hold for large values of x, y.
Linear Programming (LP) was employed as an alternative to sum-of-squares programming by Sankaranarayanan et al. [70] to generate Lyapunov functions, and later applied by Yang et al. [89] to similarly generate barrier certificates.The main idea behind this approach is to employ a linear relaxation, whereby non-negativity of a polynomial p is witnessed, subject to non-negativity of (basis) polynomials p 1 , p 2 , . . ., p k , i.e.
In cases when the evolution constraint Q is described using a conjunction of polynomial inequalities Q ≡ q 1 ≥ 0 ∧ • • • ∧ q l ≥ 0 (e.g. in the case of hyperboxes or polyhedra), one may form all products p i = q α1i 1 • • • q α li l of maximum degree less than or equal to that of the parametric template p a,d and use them to solve the linear relaxation for p 1 ≥ 0 ∧ • • • ∧ p k ≥ 0 → p a,d ≥ 0 using linear programming, obtaining a polynomial which is non-negative on Q.The conditions for barrier certificates are encoded in an obvious way.
In using SOS or LP to search for barrier certificates, one is not concerned with optimizing the value of any particular objective function (which can be set to be the zero function); one is rather interested in finding a feasible solution to a set of constraints.For LP, it is possible to use an SMT solver which supports the theory of linear real arithmetic (LRA, e.g., as supported by Z3) to search for models of formulas describing the constraints to obtain instantiations of the parameter variables in the template; however, in our experience, implementations of linear programming solvers (especially employing interior point algorithms) in Mathematica and MATLAB offer considerably better performance compared to Z3 (which implements the Dual Simplex algorithm [19]).

Strategies for Invariant Generation
The implementation of primitive invariant generation methods from Section 3.3 in a single framework is a significant undertaking in itself.The overall goal be- Fig. 11: Invariant synthesis using the differential saturation loop in Pegasus hind Pegasus, however, is to enable these heterogeneous methods to be effectively deployed and fruitfully combined into strategies for invariant generation that are tailored to specific classes of verification problems.Different invariant generation strategies are invoked in Pegasus, depending on the classification of the input problem it receives from the problem classifier.In this section, and for the evaluation, we focus on the most challenging and general class of non-linear systems in which no further structure is known or assumed beyond the fact that the right-hand sides of the ODEs are polynomials.

Differential Saturation
The main invariant generation strategy Pegasus uses for general non-linear systems is based on a differential saturation procedure [57].Briefly, the procedure loops through a prescribed sequence of invariant generation methods and successively attempts to strengthen the domain constraint using invariants found by those methods until the desired safety condition is proved. 9otably, this loop allows Pegasus to exploit the strengths of different invariant generation methods, even if it is a priori unclear whether one is better than the other.The precise sequencing of invariant generation methods is also important in this strategy to avoid redundancy.Pegasus orders the methods by computational efficiency, e.g. it first searches for first integrals, followed by Darboux polynomials and barrier certificates.This sequencing allows slower methods to exploit invariants that are quickly generated by earlier methods.

Example 5
The synergy between individual methods exploited by differential saturation is illustrated in Fig. 11 for an example from our benchmarks.Initially (leftmost plot), the entire plane (in blue) is under consideration and Pegasus wants to show the safety property that trajectories from the initial states (in green) never reach the unsafe states (in red).In the second plot, Pegasus confines its search to the domain x 1 > 0 using the generated Darboux polynomial x 1 .In the third plot, using x 1 > 0, qualitative analysis finds the invariant x 2 > 0 (whose invariance depends on x 1 > 0) which further confines the evolution domain.Finally (rightmost plot), Pegasus finds a barrier certificate (of polynomial degree 2) that suffices to show the safety property within the strengthened domain (which, by construction, is invariant).The final invariant region contains several sharp corners and thus cannot be directly obtained as the sub-level set of a single polynomial barrier certificate.Instead, it incorporates a conjunction of invariants discovered earlier by other means.
Remark 9 Pegasus extracts proof hints from the internal reasoning sequence used in its differential saturation strategy, e.g., it tracks the order of construction of the invariants x 1 > 0, x 2 > 0, . . .from Example 5 and how they were individually proved.These hints are useful for deductive tools like KeYmaera X because they can be used to guide its proofs for the generated invariants in a corresponding, step-by-step manner, with the most appropriate verification technique for the invariant used at each step.
Given an input safety verification problem, it is a priori unknown which of the invariant generation methods used for differential saturation would succeed; and even for those that do succeed, it is difficult to predict the precise duration required.The overall strategy in Pegasus imposes carefully balanced timeouts, where each method called by differential saturation attempts to: detect their applicability efficiently to conserve time budgets for other methods if they are not applicable, -keep track of intermediate results and report partial results (if applicable) when their individual timeouts are hit, -efficiently check when they are done.
Pegasus uses configuration parameters to adjust timeouts and method behaviour, e.g., maximum degree of barrier certificate templates.In addition, Pegasus supports configuration of the overall strategy behaviour in terms of combining method results, how it handles method timeouts, and how it detects when the methods succeeded.In the current implementation, and in Section 6, we explore the following strategy configuration options: C1 Auto-Reduction: whether or not to filter redundant invariants when combining results C2 Heuristic Search: whether or not to apply qualitative analysis and other heuristic search methods C3 Budget Redistribution: strict method timeouts or redistribution of unused time budget to later methods C4 Subsystem Splitting: whether or not to analyse subsystems separately Option C1 allows Pegasus to find invariants of lower descriptive complexity, which may be more insightful for users and easier to prove in KeYmaera X. Options C2-C4 allow expert users finer control over how Pegasus searches for invariants.For example, C4 is useful when the input problem is known to consist of many subsystems of ODEs [57] that can be tackled separately.The tradeoff between these options is qualitatively evaluated in Section 6.
-2 -1 0 1 The differential saturation strategy uses a pot of primitive invariant generation methods without (directly) adding more logical or mathematical considerations.The differential divide-and-conquer (DDC) proof rule [75] is an example logical technique that also fits well into the Pegasus framework.Briefly, the rule says that if p = 0 is an invariant for both the forwards ODE x = f (x) and the backwards ODE x = −f (x), then the state space partitions into three invariant subspaces p < 0, p = 0, p > 0, and it suffices to consider the invariant generation subproblems (restricted to each subspace) separately.All Darboux polynomials p (Section 4.4.2) meet the forwards-and-backwards invariance criteria and can be used to partition the state space.Indeed, this DDC strategy is already implicitly used in the invariant generation method for rational first integrals in Section 4.4.3, which partitions the state space using Darboux polynomials, and then generates rational first integrals on the resulting subproblems.Pegasus generalizes this by looking for invariants on each subproblem instead, i.e., by replacing steps 4 and 5 from the method described in Section 4.4.2 as follows: 4* For each unpruned conflict cell S, define a new invariant generation subproblem, with the original domain constraint Q restricted to Q ∧ S. 5* Call the differential saturation strategy (Section 5.1) to find an invariant on all newly generated subproblems.

Example 6
The differential divide-and-conquer strategy is illustrated in Fig. 12 for a tweaked version Example 5 with larger initial set and smaller unsafe set.As before, initially (leftmost plot), the entire plane (in blue) is under consideration and Pegasus wants to show the safety property that trajectories from the initial states (in green) never reach the unsafe states (in red).Pegasus partitions the problem into three subproblems, shown in the subsequent plots, using the Darboux polynomial x 1 ; in those plots, only the part of the plane relevant to each subproblem is drawn.Note that in the third plot, the domain constraint x 1 = 0 is slightly (but soundly) enlarged to −0.2 ≤ x 1 ≤ 0.2 for visibility in the illustration as it would otherwise be an infinitesimal strip.In the second (domain constraint x 1 < 0) and third (domain constraint x 1 = 0, enlarged) plots, the subproblems are proved trivially because they contain no

Evaluation
This section presents a qualitative evaluation of the invariant generation capabilities of Pegasus and its interaction with the ODE proving tactics of KeYmaera X.The insights obtained from these benchmarks provide useful default configuration options for Pegasus, e.g., those described in Section 5.

Benchmark Suite
The benchmark suite consists of 141 continuous safety verification problems, with 90 earlier problems [78] and 51 new ones, all drawn from the literature [5,15,18,21,26,28,30,32,33,36,41,42,50,69,76,89,90,91].Some of the problems are drawn from papers that present and discuss properties of a system of ODEs without explicitly providing initial and safe conditions; in such cases, we design our own initial and safe sets based on the provided discussion.
The suite consists of linear, affine, multi-affine, and polynomial problems, see Fig. 13: 70 two-dimensional systems, 27 three-dimensional systems, 30 higher-dimensional (≥4, ≤16) systems, and 14 product systems that were   formed by randomly combining pairs of two-and three-dimensional systems.
The experiment was run on commodity hardware 10 .We briefly analyse how the invariant generation methods individually perform on the benchmark set when run with a timeout of 120s.Fig. 14 illustrates that none of the methods outperforms the others on a given problem class, with the exception of Darboux polynomials which works very well on homogeneous problems.Overall, Qualitative Analysis performs well except on linear systems (with some variation in duration), Darboux Polynomials and Barrier certificates perform consistently well across classes, but on average require significant computation time (Darboux Polynomials vary widely in duration except for homogeneous polynomial or affine problems), and First Integrals are inexpensive when successful.

Differential Saturation Performance
We analyse the differential saturation strategy compared to each invariant generation method in isolation, measuring the duration of invariant generation, duration of checking the generated invariants, and the total proof duration.
We analyse the effect of exposing proof hints with the generated invariants, and the effect of strategy configuration options C1-C4 from Section 5. Class P ------------PMPMMP P L ---L M -----M L -L P --P A P P A P PM L -L P P Empty columns are unsolved.Legend: the combined Differential Saturation (DS) strategy against Qualitative Analysis (QA), First Integrals (FI), Darboux Polynomials (DP), and Barrier Certificates (BC), on total proof duration (T), generation duration (G), and checking duration (C).Earlier reported results [78] are also shown for comparison (DS '19).ODE classification is annotated at the top: homogeneous polynomial (H), polynomial (P), linear (L), affine (A), multi-affine (M), dashes indicate same class as the enclosing labels.

Differential Saturation versus Individual Generation Methods
The results comparing differential saturation against individual methods for each benchmark problem are shown in Fig. 15.Several experimental insights can be drawn from these results: (i) different invariant generation methods generally solve different subsets of the problems, (ii) invariant generation almost always dominates overall proof duration although invariant checking becomes more expensive as problem dimension increases, (iii) when multiple methods solve a problem, qualitative analysis and first integrals are often quickest, followed by Darboux polynomials and then barrier certificates, (iv) the differential saturation strategy effectively combines invariant generation methods; it solves 11 additional problems (of which 6 are product systems) that no individual method solves by itself.Differential saturation is especially effective on product systems because each part of the product may be only solvable using a specific method, (v) Pegasus has improved compared to its earlier version [78] and now solves 7 previously unsolved benchmarks.
To further evaluate the effectiveness of combining methods by differential saturation, Fig. 16 plots the accumulated duration for solving the fastest n problems.The main insights here are: (i) differential saturation solves the largest number of problems per accumulated time, which means that, despite sequential execution, it often succeeds in trying out the most efficient method first and fails fast when earlier methods fail to apply, (ii) first integrals are inexpensive (especially in terms of checking) when they solve problems, (iii) checking barrier certificates and Darboux polynomials is much faster than generating them, and (iv) qualitative analysis is less expensive for generation than other methods, but is most expensive for checking due to missing proof hints.

Differential Saturation Configuration Options
We explored the effect of configuration options on the invariant generation and subsequent checking duration by disabling features of the differential saturation procedure.Specifically, we executed differential saturation: C1ÄR No Auto-Reduction, which is expected to speed up generation but may cause redundant cuts.C2ḦS No Heuristic Search, which is expected to produce more principled invariants and more specific proof hints but solve fewer problems.C3BR No Budget Redistribution, which is expected to result in a more predictable generation duration but solve fewer problems.

C4 & &
SS No Subsystem Splitting, which is expected to result in faster performance on problems without clear subsystems, but solve fewer problems overall (e.g., the product problems should benefit from C4).

PH
No Proof Hints, which is expected to slow down invariant checking but have no effect on invariant generation.
Figure 17 shows the benefits and drawbacks of each configuration option, separated in invariant generation and checking, while Fig. 18 summarizes the cumulative effect of configuration options. 11xcept for Auto-Reduction (C1), which is a post-processing step after generating invariants, disabling features mostly results in faster generation on some problems at the expense of slowing down or not solving others at all (see Fig. 17b).Overall, the configuration options have little effect on most examples for proof checking (see Fig. 17c), but can make a difference on some select single examples: -  Techniques developed for qualitative simulation have been applied to prove temporal properties of continuous systems in the work of Shults and Kuipers [74], as well as Loeser, Iwasaki and Fikes [43].Zhao [93] developed a tool, MAPS, to automatically identify significant features of dynamical systems, such as stability regions, equilibria, and limit cycles.Since our ultimate goal is sound invariant generation, we are less interested in a full qualitative analysis of the state space.In the verification community, discrete abstraction of hybrid systems was studied by Alur et al. [1].The case of systems whose continuous motion is governed by non-linear ODEs was studied in the work of Tiwari and Khanna [82,84].Tiwari further studied reachability of linear systems [81], using information from real eigenvectors and ideas from qualitative abstraction to generate invariants.Zaki et al. [91] were the first to apply Darboux polynomials to verification of continuous systems using discrete abstraction.Numerous works employ barrier certificates for verification [15,37,61,77,89].Since we implement many of the above techniques as methods for invariant generation in our framework, our work draws heavily upon ideas developed previously in the verification and hybrid systems communities.Previously [75], we introduced a construction of exact abstractions and applied rudimentary methods from qualitative analysis to compute invariants; in certain ways, our present work also builds on this experience, incorporating some of the techniques as special methods in a more general framework.The coupling between KeYmaera X and Pegasus that we pursue in our work is quite distinct from the use of trusted oracles in the work of Wang et al. [86] (for the HHL prover) and provides a sound framework for reasoning with continuous invariants that is significantly less exposed to soundness issues in external tools.
A complete semi-algorithm for computing algebraic invariants (described by zero sets of polynomial functions) for polynomial systems of ODEs was developed by Ghorbal and Platzer [26].An interesting development along very similar lines was also recently pursued by Boreale [10], whose method makes use of the algebraic nature of the precondition (initial set) in the verification problem in order to speed up the algebraic invariant generation.Both of these (semi-)algorithms involve enumeration of polynomial templates; the biggest practical difficulty stems from the computational cost of minimizing the rank of symbolic matrices in [26], and computing the generators of real radical ideals in [10], both of which are difficult problems with the latter having few algorithms with robust implementations currently in existence.In the future, we hope to extend Pegasus with an implementation of these techniques, thereby extending our current capabilities.

Outlook and Challenges
The improvements in continuous invariant generation have a significant impact on the overall proof automation capabilities of KeYmaera X and serve to increase overall system usability and improve user experience.Better proof automation will certainly also be useful in future applications of provably correct runtime monitoring frameworks, such as ModelPlex [47], as well as frameworks for generating verified controller executables, such as VeriPhy [9].Some interesting directions for extending our work include implementation of reachable set computation algorithms for all classes of problems where this is possible.For instance, semi-algebraic reachable sets for diagonalizable classes of linear systems with tame eigenvalues [25,38].The complexity of invariants obtained using these methods may not always make them practical, but they would provide a valuable fallback in cases where simpler invariants cannot be obtained using our currently implemented methods.
A more pressing challenge lies in expanding the collection of safety verification problems for continuous systems.While we have done our best to find compelling examples from the literature, a larger corpus of problems would allow for a more comprehensive empirical evaluation of invariant generation strategies and could reveal interesting new insights that can suggest more effective strategies.
Correctness of decision procedures for real arithmetic is another important challenge.KeYmaera X currently uses Mathematica's implementation of real quantifier elimination to close first-order real arithmetic goals, primarily due to the impressive performance afforded by this implementation (compared to currently existing alternatives).Removing this reliance by efficiently building fully formal proofs of real arithmetic formulas within dL (e.g. through exhibiting appropriate witnesses [39,59]) is an important task for the future.

Conclusion
Among verification practitioners, the amount of manual effort required for formal verification of hybrid systems is one of the chief criticisms leveled against the use of deductive verification tools.Manually crafting continuous invariants often requires expertise and ingenuity, just like manually selecting support function templates for reachability tools [22], and presents the major practical hurdle in the way of wider industrial adoption of this technology.In this article, we describe our development of a system designed to help overcome this hurdle by automating the discovery of continuous invariants.To our knowledge, this work represents the first large-scale effort at combining continuous invariant generation methods into a single invariant generation framework and making it possible to create more powerful invariant generation strategies.The approach we pursue is unique in its integration with a theorem prover, which provides formal guarantees that the generated invariants are indeed correct (in the form of dL proofs, automatically).The results we observe in our evaluation are highly encouraging and suggest that invariant discovery can be improved considerably, opening many exciting avenues for applications and extensions.

1 2 x2Fig. 9 :
Fig. 9: Rational first integral r λ constructed from three Darboux polynomials.Zero sets of the three Darboux polynomials shown in solid green, blue and red.Invariant level sets of the rational first integral shown in dashed black for values r λ = 1 10 , 1, −2, respectively.

Fig. 10 :
Fig.10:(Left) A candidate invariant generated using numerical barrier certificates (in blue) for the safety verification problem of showing that solutions from the green initial state never reach the red unsafe states.(Right) A zoomed out view of the safety verification problem, showing that the candidate invariant is, in fact, not an invariant of the ODE because some states can exit the invariant (highlighted with a dashed red circle).

Fig. 13 :
Fig. 13: Benchmark suite classification unsafe states.In the rightmost (domain constraint x 1 > 0) plot, Pegasus finds a barrier certificate (in blue) that solves the subproblem.
b r a i c P r e c o n

Fig. 15 :
Fig.15: Comparison of invariant generation methods.Each column represents one benchmark problem and the colour encodes duration (lighter is faster).Empty columns are unsolved.Legend: the combined Differential Saturation (DS) strategy against Qualitative Analysis (QA), First Integrals (FI), Darboux Polynomials (DP), and Barrier Certificates (BC), on total proof duration (T), generation duration (G), and checking duration (C).Earlier reported results[78] are also shown for comparison (DS'19).ODE classification is annotated at the top: homogeneous polynomial (H), polynomial (P), linear (L), affine (A), multi-affine (M), dashes indicate same class as the enclosing labels.

Fig. 16 :
Fig. 16: Cumulative logarithmic time (in seconds) taken to solve the fastest n problems (more problems solved and flatter is better) Invariant generation duration in multiples of fastest method DS (C)

Fig. 18 :
Fig. 18: Configuration options: cumulative logarithmic time (in seconds) taken to solve the fastest n problems (more problems solved and flatter is better) No Proof Hints (PH): Several examples vastly increase in checking duration or fail to check entirely.Conclusion: KeYmaera X's checking procedures spend time to rediscover efficient proofs that were already known by construction during the generation.Proof hints should be kept wherever possible, especially since they are inexpensive to produce in Pegasus.-NoAuto-Reduction (C1ÄR): slight increase in proof duration on some examples, but compensated well by decrease in generation duration.Conclu-sion: Auto-Reduction is not an essential technique in invariant generation or proof checking.-No Heuristic Search (C2ḦS): slight decrease in generation and checking duration on most examples at the expense of considerably increasing or entirely failing to generate invariants for some of the examples.Conclusion: option C2 should be configurable by users.-No Budget Redistribution (C3BR): slight decrease in generation duration on some examples at the expense of just not solving others, mostly without effect on proof duration.Conclusion: C3 is a useful technique in invariant generation, and should typically be enabled.-No Subsystem Splitting (C4 & & SS): little effect on both generation and checking duration for solved problems, but solves considerably fewer problems including several that are not product systems.Conclusion: C4 is a useful technique in invariant generation, and should typically be enabled.