Automated Formal Synthesis of Neural Barrier Certificates for Dynamical Models

We introduce an automated, formal, counterexample-based approach to synthesise Barrier Certificates (BC) for the safety verification of continuous and hybrid dynamical models. The approach is underpinned by an inductive framework: this is structured as a sequential loop between a learner, which manipulates a candidate BC as a neural network, and a sound verifier, which either certifies through algorithmic proofs the candidate's validity or generates counter-examples to further guide the learner. We compare the approach against state-of-the-art techniques, over polynomial and non-polynomial dynamical models: the outcomes show that we can synthesise sound BCs up to two orders of magnitude faster, with in particular a stark speedup on the verification engine (up to five orders less), whilst needing a far smaller data set (up to three orders less) for the learning part. Beyond the state of the art, we further challenge the (verification side of the) approach on a hybrid dynamical model.


Introduction
Barrier Certificates (BC) are an effective and powerful technique to prove safety properties on models of continuous and hybrid dynamical systems [19,20]. Whenever existing, a BC partitions the state space of the model into two parts, ensuring that all trajectories starting from a given initial set, located within one side of the BC, cannot reach a given set of states (deemed to be unsafe), located on the other side. Thus a successful synthesis of a BC (usually not a unique object) represents a formal proof of safety for the dynamical model. This work addresses the safety of systems modelled by non-linear differential equations (ODE), and presents a novel method for the automated and formal synthesis of BC. Our new approach leverages Satisfiability Modulo Theory (SMT) and inductive reasoning (CEGIS, Figure 1, introduced later), to guarantee the formal correctness of the synthesis procedure: this soundness prevents algorithmic or numerical errors related to BC synthesis [9].
Background A few techniques have been developed to synthesise BC. For polynomial models, sum-of-squares (SOS) and semi-definite programming relaxations [9,13,15,26] convert the BC synthesis problem into constraints expressed as linear or bilinear matrix inequalities, which are however numerically solved unsoundly. To increase scalability and to enhance expressiveness, numerous barrier formats have been considered: BC based on exponential conditions are presented in [13]; BC based on Darboux polynomials are outlined in [30]; [27] newly introduces a multi-dimensional generalisation of BC, thus broadening their scope and applicability. BC can also be used to verify safety of uncertain (e.g. parametric) models [18]. BC find various applications spanning robotics, multi-agent systems, and biology [6,29].
Model invariants (namely, regions that provably contain model trajectories) can be employed as BC, though their synthesis is less general, as it does not comprise the unsafe set: [17] introduces a fixpoint algorithm to find algebraicdifferential invariants for hybrid models; invariants can be characterised analytically [3] or synthesised computationally [7]. Invariants can be alternatively studied by Lyapunov theory [4], which provides stability guarantees for dynamical models, and thus can characterise invariants (and barriers) as side products. Whilst Lyapunov theory is classically approached either analytically (explicit synthesis) or numerically (with unsound techniques), relevant for the automated and sound results of this work Lyapunov function synthesis has been recently studied via SMT. In [25] Lyapunov functions are soundly found within a parametric framework, by constructing a system of linear inequality constraints over unknown coefficients. Similarly [21,22,23]   Contributions We introduce a method that efficiently exploits machine learning, whilst guaranteeing formal proofs of correctness via SMT. We follow a CounterExample-Guided Inductive Synthesis (CEGIS) procedure [28], which is structured as an inductive loop between a learner and a verifier (cf. Fig. 1). A learner numerically trains a neural network (NN) to fit the conditions for a BC expressed as a loss function; then a verifier through an SMT solver either formally proves the validity or provides a counter-example, where the barrier conditions are violated, that is passed back to the learner. Our synthesis method for neural BC is formally sound and fully automatic, and is shown to be much faster and to evidently require less data than state-of-the-art cognate results.
Beyond the State of the Art Cognate work [31] presents a method to compute BC using neural networks and to verify their correctness a-posteriori: as such, it does not generate counter-examples within an inductive loop, as in this work. [31] considers large sample sets that are randomly divided into batches and fed to a feed-forward NN; the verification at the end of the (rather long) training either validates the candidate, or invalidates it and the training starts anew on the same dataset. In Section 4 the method in [31] is shown to be slower (both in the training and in the verification), and to require more data than the CEGIS-based approach underpinning this work, which furthermore introduces numerous bespoke optimisations, as outlined in Section 3: our CEGIS-based technique exploits fast learning, fast verification, and an enhanced communication between these two components. Related to the work on BC is the synthesis of Lyapunov functions, mentioned above. The construction of Lyapunov Neural Networks (LNNs) has been studied with approaches based on simulations and numerical optimisation, which are formally unsound [24]. Formal methods for Lyapunov synthesis are introduced in [4], together with a counterexample-based approach using polynomial candidates. The work is later extended in [1], which employs NN as candidates. The generation of control Lyapunov functions using counterexample-based NN is similarly considered in [8], however by means of differing architectural details and with a different SMT solver.

Safety Analysis with Barrier Certificates
We address the safety verification of continuous-time dynamical models by designing barrier certificates (BC) over the continuous state space X of the model. We consider n-dimensional dynamical models described bẏ where f : X → R n is a continuous vector field, X ⊆ R n is an open set defining the state space of the system, and X 0 represents the set of initial states. Given model (1) and an unsafe set X u ⊂ X, the safety verification problem concerns checking whether or not all trajectories of the model originating from X 0 reach the unsafe region X u . BC offer a sufficient condition asserting the safety of the model.

Definition 1.
The Lie derivative of a continuously differentiable scalar function B : X → R, with respect to a vector field f , is defined as followṡ Intuitively, this derivative denotes the rate of change of function B along the model trajectories.
Proposition 1 (Barrier Certificate for Safety Verification, [19]). Let the model in (1) and the sets X, X 0 and X u be given. Suppose there exists a function B : X → R that is differentiable with respect to its argument and satisfies the following conditions: then the safety of the model is guaranteed. That is, there exists no trajectory of the model contained in X, starting from an initial state in X 0 , that ever enters set X u .
Consider a trajectory x(t) starting in x 0 ∈ X 0 and the evolution of B(x(t)) along this trajectory. Whilst the first of the three conditions guarantees that B(x 0 ) < 0, the last condition asserts that the value of B(x(t)) along a trajectory x(t) cannot become positive. Hence such a trajectory x(t) cannot enter the set X u , where B(x) > 0 (second condition), thus ensuring the safety of the model.

Synthesis of Neural Barrier Certificates via Learning and Verification
We introduce an automated and formal approach for the construction of barrier certificates (BC) that are expressed as feed-forward neural networks (NN). The procedure leverages CEGIS (see Fig. 1) [28], an automated and sound procedure for solving second-order logic synthesis problems, which comprises two interacting parts. The first component is a learner, which provides candidate BC functions by training a NN over a finite set of sample inputs. The candidate is passed to the second component, a verifier, which acts as an oracle: either it proves that the solution is valid, or it finds one (or more) instance (called a counter-example) where the candidate BC does not comply with required conditions. The verifier consists of an SMT solver [14], namely an algorithmic decision procedure that extends Boolean SAT problems to richer, more expressive theories, such as linear and non-linear arithmetics. More precisely, the learner trains a NN composed of n of input neurons (this matches the dimension of the model f ), k hidden layers, and one output neuron (recall that B(x) is a scalar function): this NN candidate B is required to closely match the conditions in Eq. (3) over a discrete set of samples S, which is initialised randomly. The verifier checks whether the candidate B violates any of the conditions in Eq. (3) over the entire set X and, if so, produces one (or more) counter-examples c. We add c to the samples set S as the loop restarts, hence forcing the NN to be trained also over the generated counter-examples c. This loop repeats until the SMT verifier proves that no counter-examples exist or until a timeout is reached. CEGIS offers a scalable and flexible alternative for BC synthesis: on the one hand, the learner does not require soundness, and ensures a rapid synthesis exploiting the training of NN architectures; on the other, the algorithm is sound, i.e. a valid output from the SMT-based verifier is provably correct; of course we cannot claim any completeness, since CEGIS might in general not terminate with a solution because it operates over a continuous model and its domain X.
The performance of the CEGIS algorithm in practice hinges on the effective exchange of information between the learner and the verifier [2]. A core contribution of this work is to tailor the CEGIS architecture to the problem of BC synthesis: we devise several improvements to NN training, such as a bespoke loss function and a multi-layer NN architecture, together with an informative counter-example generation by the SMT verifier that is adapted to the candidate BC and the underlying dynamical model. These tailored architectural details generate in practice a rapid and efficient CEGIS loop, which is shown in this work to clearly outperform state-of-the-art methods.

Training of the Barrier Neural Network
The learner instantiates the candidate BC using the hyper-parameters k and h (depth and width of the NN), trains it over the N samples in the set S, and later refines its training whenever the verifier adds counter-examples to the set S. The class of candidate BC comprises multi-layered, feed-forward NN with polynomial activation functions. Unlike most learning applications, the choice of polynomial activations comes from the need for interpretable outputs from the NN, whose analytical expression must be readily processed by the verifier. Whilst this work emphasises the novel use of polynomial activation functions, we remark that alternatives are possible: in particular, as displayed in one case study, tanh are well-suited to our objective, being universal function approximators. The order γ of the polynomial activations is a hyper-parameter fed at the start of the CEGIS procedure. Specifically, we split the i-th hidden layer into γ portions and apply polynomial activations of order j to the neurons of the j-th portion, as shown next.
Example 1 (Polynomial Activations). Assume a NN composed of an input x, 5 hidden neurons and 1 activation-free output, with γ-th order polynomial activation, γ = 5. We split the hidden layer in γ sub-vectors, each containing one neuron. The hidden layer after the activation results in are the i-th row of the first-layer weight matrix, and the b i form the bias vector.
The learning process updates the NN parameters to improve the satisfaction of the BC conditions in (3): B(x) ≤ 0 for x ∈ X 0 , B(x) > 0 for x ∈ X u , and a negative Lie derivativeḂ (eq. (2)) over the set implicitly defined by B(x) = 0.
The training minimises a loss comprising thee terms, namely where s i , i = 1, . . . , N are the samples taken from the set S. The constants τ 0 , τ u , τ d are offsets, added to improve the numerical stability of the training. Notably, B(x) = 0 is a measure-zero set, thus it is highly unlikely that a single sample s will satisfy B(s) = 0. We then relax this last condition and consider a belt B around B(s) = 0, namely B = |B(x)| ≤ β, which depends on the hyper-parameter β. Note that we must use continuously differentiable activations throughout, as we require the existence of Lie derivatives (cf. Eq. (2)), and thus cannot leverage simple ReLUs.

Enhanced Loss Functions
The loss function in Eq. (4) experimentally yields possible drawbacks, which suggest a few ameliorations. Terms L 0 and L u solely penalise samples with incorrect value of B(x) without further providing a reward for samples with a correct value. The NN thus stops learning when the samples return correct values of B(x) without further increasing the positivity of B over X u or the negativity over X 0 . As such, the training often returns a candidate B(x) with values just below τ 0 in X 0 or above τ u in X u . These candidates are easily falsified, thus potentially leading to a large number of CEGIS iterations. We improve the learning by adopting a (saturated) Leaky ReLU, hence rewarding samples that evaluate to a correct value of B(x). Noting that where α is a small positive constant, we rewrite term L 0 as where satReLU is the saturated ReLU function. Term L u is similarly modified. In view of the composite nature of our training objective, incorrect samples account for the major contribution to the loss function, leading the NN to correct those first. At a second stage, the network finds a direction of improvement by following the leaky portion of the loss function. This is saturated to prevent the training from following only one of these directions, without improving the other loss terms.
Another possible drawback of the loss function in (4) derives from the term L d : it solely accounts for a penalisation of the sample points within B. To quickly and myopically improve the loss function, the training can generate a candidate BC for which no samples are within B -we experimentally find that this behaviour persists, regardless of the value of β. Similarly to L 0 and L u , we reward the points within a belt fulfilling the BC condition: namely, we solely apply the satReLU function to reward samples s with a negativeḂ(s), whilst not penalising valueṡ B(s) ≥ 0. The training is driven to include more samples in B, guiding towards a negativeḂ(s), and finally enhancing learning. The expression of L d results in Finally, we choose an asymmetric belt B = −β 1 ≤ B(s) ≤ β 2 , with β 2 > β 1 > 0 to both ensure a wider sample set and a stronger safety certificate.
Multi-layer Networks Polynomial activation functions generate interpretable barrier certificates with analytical expressions that are readily verifiable by an SMT solver. However, when considering polynomial networks, the use of multilayer architectures quickly increase the order of the barrier function: a k-layer network with γ-th order activations returns a polynomial of kγ degree. We have experienced that deep NN provide numerical robustness to our method, although the verification complexity increases with the order of the polynomial activation functions used and with the depth of the NN. As a consequence, our procedure leverages a deep architecture whilst maintaining a low-order polynomial by interchanging linear and polynomial activations over adjacent layers. Linear activations positively affect the training, providing the robustness needed to the synthesis of BC, without increasing the order by new polynomial terms.

Learning in Separate Batches
The structure of the conditions in (3) and the learning loss in (4) naturally suggests a separate, parallel approach to training. We then split the dataset S into three batches S 0 , S u and S x , each including samples belonging to X 0 , X u and the complement X \ {X 0 , X u }, respectively. For training, we compute the loss function in a parallel fashion. Similarly, for the verifier, generated counter-examples are added to the relevant batch.

Certification of the Barrier Neural Network, or Falsification via Counter-examples
Every candidate BC function B(x) which the learner generates requires to be certified by the verifier. Equivalently, in practice the SMT-based verifier aims at finding states that violate the barrier conditions in (3) over the continuous domain X. To this end, we express the negation of such requirements, and formulate a nonlinear constrained problem over real numbers, as The verifier searches for solutions of the constraints in Eq. (8): this requires manipulating non-convex functions, in view of our choice of activations. In particular, tanh activations result in a more complex task: whilst still decidable, non-polynomial clauses render the computation NP-hard [14]. On the one hand, the soundness of our CEGIS procedure heavily relies on the correctness of SMT solving: an SMT solver never fails to assert the absence of solutions for (8). As a result, when it states that formula (8) is unsatisfiable, i.e. returns unsat, B(x) is formally guaranteed to fulfil the BC conditions in Eq. (3). On the other hand, the CEGIS algorithm offers flexibility in the choice of the verifier, hence we implement and discuss two SMT solvers: dReal [12] and Z3 [10]. dReal is a δ-complete solver, namely the unsat decision is correct [11], whereas when a solution for (8) is found, this comes with a δ-error bound. The value of δ characterises the procedure precision. In our setting, it is then possible to return spurious counter-examples: nevertheless, these are then used as additional samples and do not invalidate the sound outcomes of the procedure, but rather help synthesising a more robust barrier candidate. dReal is capable of handling non-polynomial terms, such as as exponentials or trigonometric vector fields f for some of the models considered in Section 4. Z3 is a powerful, sound and complete SMT-solver, namely its conclusions are provably correct both when it determines the validity of a BC candidate and when it provides counter-examples. The shortcoming of Z3 is that it is unable to fully handle non-polynomial formulae.

Prioritisation and Relaxation of Constraints
The effectiveness of the CEGIS framework is underpinned by rapid exchanges between the learner and the verifier, as well as by quick NN training and SMT verification procedures. We have experienced that the bottleneck resides in the handling of the constraint η d = (B(x) = 0 ∧Ḃ(x) > 0) by the SMT solver, since the formula contains the high-order expressionḂ(x) and because it is defined over the thin region of the state space implicitly characterised by B(x) = 0. As a consequence, we have prioritised constraints η 0 = (x ∈ X 0 ∧ B(x) > 0) and η u = (x ∈ X u ∧ B(x) ≤ 0): that is, if either clauses is satisfied, i.e. a counter-example is found for at least one of them, the verifier omits testing η d whilst the obtained counter-examples are passed to the learner. The constraint η d is thus checked solely if η 0 and η u are both deemed to be unsat. Whenever this occurs, and the verification of η d times out, the solver searches for a solution of a relaxed constraint (|B(x)| < τ v ∧Ḃ(x) > 0), similarly to the improved learning conditions discussed in Eq. (7). Whilst this constraint is arguably easier to solve in general, it may generate spurious counterexamples, namely a samplex that satisfy the relaxed constraint, but such that B(x) = 0. The generation of these samples does not contradict the soundness of the procedure, and indeed are shown to improve the robustness of the next candidate BC -this of course comes with the cost of increasing the number of CEGIS iterations.

Increased Information from Counter-examples
The verification task encompasses an SMT solver attempting to generate a counter-example, namely a (single) instance satisfying Eq. (8). However, a lone sample might not always provide insightful information for the learner to process. Naïvely asking the SMT solver to generate more than one counter-example is generally expensive, as it is done sequentially. Specifically, the verifier solves Eq. (8) to find a first counter-examplex; then, to find any additional sample, we include the statement (x =x) and solve again for the resulting formula. We are interested in finding numerous points invalidating the BC conditions and feed them to the learner as a batch, or perhaps in increasing the information generated by the verifier by finding a sample that maximises the violation of the BC conditions. To this end, we randomly generate a cloud of points around the generated counter-example: in view of the continuity of the candidate function B, samples around a counterexample are also likely to invalidate the BC conditions. Secondly, for all points in the cloud, we compute the gradient of B (or ofḂ) and follow the direction that maximises the violation of the BC constraints. As such, we follow the B (resp.Ḃ) maximisation when considering x ∈ X 0 (x s.t. |B(x)| < τ v ), and viceversa when x ∈ X u . This gradient computation is extremely fast as it exploits the neural architecture, and it provides more informative samples for further use by the learner.

Case Studies and Experimental Results
We have implemented the discussed new procedure using the PyTorch library. Corroborating the flexibility and scalability of the approach, all experiments are performed on an unassuming laptop workstation, running Ubuntu 18.04 with 8 GB RAM. We demonstrate that the proposed method finds provably correct BCs on benchmarks from literature comprising both polynomial and non-polynomial dynamics (vector fields f ), and we newly tackle a hybrid model (this is a model with two different operating models and associated dynamics) that is challenging for the verification engine. To confirm the flexibility of our architecture in accepting different SMT solvers, we verify generated candidate BC using dReal in the first four benchmarks, whereas we study the hybrid model using the Z3 tool. In all the examples, we use a learning rate of 0.1 for the NN and the loss function described in Section 3.1 with α = 10 −4 , τ 0 = τ u = τ d = 0.1. The region in Eq. (7) is limited by β 1 = 0.1 whilst β 2 = ∞. We set a verification parameter τ v = 0.05 (cf. Sec. 3.2), and the verification precision for dReal to δ = 10 −6 . Table 1 summarises the outcomes. For the first four benchmarks, we compare our procedure, denoted as CEGIS, with the results from [31], which however does not handle the hybrid model in the last benchmark. We have run the algorithm in [31] and reported the cumulative synthesis time under the 'Learn' column, however the verification (we emphasise this is done only once, a-posteriori, in [31]) is not included in the repeatability package, hence we report the results from [31] (these are generated with much more powerful hardware). The outcomes suggest that we obtain much faster synthesis and verification times (up to five orders of magnitude), whilst requiring up to only 0.02% (see Obstacle Avoidance Problem) of the training data: [31] performs a uniform sampling of the space X, hence suffers especially in the 3-d case, where the learning runs two orders of magnitude faster. It is instructive to notice that all the case studies are solved with a small number of iterations (up to 8) of the CEGIS loop, which is thus promising to tackle synthesis problems over more complex models. The analytical expression of the models under study, together with a detailed analysis of the CEGIS iterations, are reported in the Supplementary Material. Darboux Model This 2-dimensional model is approached using polynomial BC. The initial and unsafe sets are depicted in green and red, respectively, whereas the level set B(x) = 0 is outlined in black. The work [30] reports that LMI-based methods fail to verify this model using polynomial templates of degree 6. Our approach generates the BC shown in Fig. 2 (left) in less than 30 seconds, roughly half as much as in [31], and using only 500 initial samples vs more than 65000. The BC is derived from a 3-layer architecture of 10 nodes each, with linear, polynomial with γ = 3, and linear activations respectively. Exponential Model This model from [16] shows that our approach extends to non-polynomial systems encompassing exponential and trigonometric functions. Our algorithm provides a valid BC in 14 seconds, around 5% of the results in [31], again using solely 500 initial samples. The BC, depicted in Fig.2 (centre), results from a single-layer neural architecture of 10 nodes, with polynomial activation function with γ = 3.
Obstacle Avoidance Problem This 3-dimensional model, originally presented in [5], describes a robotic application: the control of the angular velocity of a two-dimensional airplane, aimed at avoiding a still obstacle. The BC is obtained from a single-layer NN comprising 10 neurons, using polynomial activations with γ = 2. Fig. 2 (right) plots the vector field on the plane z = 0. Our procedure takes 0.6% of the computational time in [31], providing a valid BC with 1 iteration starting from an initial dataset of 500 samples.
Polynomial Model This model describes a polynomial system [20] and presents initial and unsafe sets with complex, non convex shapes [31]. Approaches from literature, such as (unsound) SOS-based procedures [15,26], have required highorder polynomial templates, which has suggested the use of alternative activation functions. The BC, shown in Fig. 3, is generated using a 10-layer NN with tanh activations. Needing just around 1 min and only 1000 initial samples, the overall procedure is 30 times faster than that in [31]. Hybrid Model We finally challenge our procedure with a 2-dimensional hybrid model, which extends beyond the capability of the results in [31]. This hybrid framework partitions the set X into two non-overlapping subsets, X 1 and X 2 .
Each subset is associated to different model dynamics, respectively f 1 and f 2 . In other words, the model trajectories evolve according to the f 1 dynamics when in X 1 , and according to f 2 when in X 2 . The structure of this model represents a non-trivial task for the verification engine, for which we employ the Z3 SMT solver: notice the dimensionally larger computation times. The learning phase has instead been quite fast. The BC (Fig.3) is obtained from a single-layer NN comprising 3 neurons, using polynomial activations with γ = 2, overall in less than 30 seconds, starting with an initial dataset of 500 samples.

Conclusions and Future Work
We have presented a new inductive, formal, automated technique to synthesise neural-based barrier certificates for polynomial and non-polynomial, continuous and hybrid dynamical models. Thanks to a number of architectural choices for the new procedure, our method requires less training data and thus displays faster learning, as well as quicker verification time, than state-of-the-art techniques.
Beyond improving the scalability of the approach, future work includes the development of an automatic selection of activation functions that are tailored to the dynamical models of interest.