The semiring I(R±∞) is a partial order (more specifically, it is a bounded distributive lattice) and therefore does not fit existing monitoring frameworks that consider only linear orders (e.g., the max-min semiring R±∞ of the extended reals and the associated sliding-max/min algorithms).

### **3 Symbolic Quantitative Traces and Languages**

In this section we start with our investigation of how to generalize the "robustness degree" of [22] to our abstract algebraic setting. The result of [22] that relates the robustness degree with the robust semantics is an inequality. For this reason, we focus on idempotent semirings, for which there is a natural partial order ≤ that is induced by semiring addition (x ≤ y iff x + y = y). Since our approach is abstract algebraic (i.e., axiomatic), we have no notion of real-valued distance between elements of D. Moreover, V does not need to be a semiring of real numbers. Instead, we rely on the intuition that for an atomic predicate p : D → V and a data item d ∈ D, the value p(d) gives a degree of truth or falsity. We propose using symbolic traces **x** = p0p1 ...pn−1, which are sequences of atomic predicates, in order to compactly represent sets of concrete traces, which are sequences of data items. If each pi represents a subset Si ⊆ D, then **x** represents the set L = S0 × S1 ×···× Sn−1 = {v0v1 ...vn−1 | vi ∈ Si for each i} of concrete traces. Moreover, given a concrete trace u = u0u1 ...un−1 ∈ Dn, we can use the value p0(u0) · p1(u1)··· pn−1(un−1) ∈ V as a quantitative measure of how close the trace u is to the set of traces L. We propose the interpretation of a formula ϕ as a language of symbolic traces. This will allow us to define the "closeness" of a trace u ∈ Dn to the specification ϕ as a (semiring) sum of all the closeness values w.r.t. each symbolic trace in the symbolic language of ϕ. We will also see that this interpretation of a formula ϕ as a symbolic language is compatible with the standard interpretation of ϕ as a set of concrete traces. Using these definitions we obtain a generalization of the theorem of [22] that relates the robustness degree with the robust semantics. Additionally, we characterize precisely the class of semirings for which this generalization is possible.

Let V be an idempotent semiring. For predicates p, q : D → V we define p ≤ q if p(d) ≤ q(d) for every d ∈ D. The intuition for p ≤ q is that p is a stronger predicate than q. We write F(D, V ) to denote the set of atomic quantitative predicates, which always includes the predicates 1 and 0. For symbolic traces **x**, **y** ∈ F(D, V )∞ with λ = |**x**| = |**y**| we define **x** ≤ **y** if **x**(i) ≤ **y**(i) for every i<λ. These relations ≤ on predicates and traces are partial orders. We define the symbolic satisfaction relation |=, where **x**, i |= ϕ says that the formula ϕ : MTL(D, V ) is satisfied by the symbolic trace **x** ∈ F(D, V )∞ at position i < |**x**|. For atomic formulas, we put **x**, i |= p iff **x**(i) ≤ p. The definition is given by induction on ϕ in the usual way. For a formula ϕ : MTL(D, V ), length λ ∈ ω ∪ {ω} and a position i<λ, we define the symbolic language SL(ϕ, λ, i) = {**x** ∈ F(D, V )λ | **x**, i |= ϕ}. For nonempty finite traces **x** ∈ F(D, V )n and u ∈ Dn of the same length, we define **x**[u] = n i=1**x**(i)(u(i)), where n = |**x**| = |u|. Since the semiring multiplication is monotone w.r.t. ≤, we see that **x** ≤ **y** implies **x**[u] ≤ **y**[u] for every u ∈ Dn. Informally, the value **x**[u] quantifies how close the concrete trace u is to the symbolic trace **x**.

**Example 6.** Let D = R and V = R±∞. For c ∈ R, the predicate p = "x ≥ c" is defined by p(d) = d−c for every d ∈ D. The predicate q = "x ≤ c" is given by q(d) = c−d for every d ∈ D. For the symbolic trace **x** = "x ≥ 1" "x ≤ 5" "x ≥ 2" and the concrete trace u = 3 6 8 we get that **x**[u] = min(2, −1, 6) = −1.

Let c, d ∈ R. For the predicates p = "x ≥ c" and q = "x ≥ d" we have that p ≤ q iff d ≤ c. Similarly, for the predicates p = "x ≤ c" and q = "x ≤ d" it holds that p ≤ q iff c ≤ d. Finally, notice that the predicates "x ≥ c" and "x ≤ d" are incomparable. Consider **y** = "x ≥ 0" "x ≤ 7" "x ≥ 1" and observe that **x** ≤ **y**.

For the formula ϕ = p ∧ F1q, where p and q are atomic predicates, we have that SL(ϕ, 2, 0) = {p q ∈ F(D, V )2 | p ≤ p and q ≤ q}.

The definition of the robustness degree in [22] involves the value −dist(u, L) = − infv∈L dist(u, v) = supv∈L −dist(u, v), where u is a trace, L is a set of traces, and dist is a metric. Notice that this is a supremum over a potentially infinite set. The semirings that we have considered so far have an addition operation that can model a finitary supremum. In order to model an infinitary supremum, we need to consider semirings that have an infinitary addition operation. A complete semiring is an algebraic structure (V, +, -, ·, 0, 1), where - i∈Ixi is the sum of the I-indexed tuple of elements (xi)i∈I , that satisfies: (1) - i∈∅ - xi = 0, i∈{j}xi = xj , - i∈{j,k}xi = xj + xk for j = k, and - k∈K - i∈Ik xi = - i∈Ixi where I = k∈KIk and the index sets (Ik)k∈K are pairwise disjoint, (2) (V, ·, 1) is a monoid, (3) the infinite distributivity properties (- i∈Ixi) · y = - i∈I (xiy) and x · ( - i∈I yi) = - i∈I (xyi) hold for every index set I and all xi, y ∈ V , and (4) 0 is an annihilator for multiplication. A complete semiring V is idempotent if - i∈Ixi = x for every non-empty index set I with xi = x for every i ∈ I. For example, (R±∞, max,sup, min, −∞, +∞) is an idempotent complete semiring. For a formula ϕ : MTL(D, V ), a trace u ∈ D+ and i<n = |u|, we define

$$\mathsf{val}(\varphi, u, i) = \sum\_{\mathbf{x} \in \mathsf{SL}(\varphi, n, i)} \mathbf{x}[u]. \tag{2}$$

Informally, val(ϕ, u, i) is a measure of how close the trace u is to satisfying ϕ at position i. It is an abstract algebraic variant of the robustness degree [22].

**Theorem 7 (Approximation).** Let D be a set of data items and V be an idempotent complete semiring. Then, the following are equivalent:

(1) The multiplication of V is idempotent and 1 is the top element of V .

(2) For every ϕ : MTL(D, V ), u ∈ D+ and i < |u|, val(ϕ, u, i) ≤ ρ(ϕ, u, i).

Proof. Assume that (1) holds.  For a symbolic language L ⊆ F(D, V )n and for u ∈ Dn, we define val(L, u) = - **x**∈L**x**(u). Let {Li}i∈I be a collection of languages with Li ⊆ F(D, V )n. Then,

$$\text{val}(\bigcup\_{i \in I} \mathcal{L}\_i, u) = \sum\_{\mathbf{x} \in \bigcup\_{i \in I} \mathcal{L}\_i} \mathbf{x}[u] \le \sum\_{i \in I} \sum\_{\mathbf{x} \in \mathcal{L}\_i} \mathbf{x}[u] = \sum\_{i \in I} \text{val}(\mathcal{L}\_i, u). \tag{3}$$

For symbolic languages L1,L2 ⊆ F(D, V )n, define L = L1 ∩ L2, L 1 = L1 \ L2 and L 2 = L2 \ L1. Then, L1 = L 1 ∪ L and L2 = L 2 ∪ L. The languages L 1,L 2,L are pairwise disjoint. So, we have that val(L1, u) = x + z and val(L2, u) = y + z, where x = val(L 1, u), y = val(L 2, u) and z = val(L, u). It follows that

$$\mathsf{val}(\mathcal{L}\_1 \cap \mathcal{L}\_2, u) = z = zz \le (x + z)(y + z) = \mathsf{val}(\mathcal{L}\_1, u) \cdot \mathsf{val}(\mathcal{L}\_2, u) \tag{4}$$

by the idempotence of multiplication. This property extends to val(L1 ∩···∩ Lk, u) ≤ val(L1, u)··· val(Lk, u). Now, we will prove (2) by induction on ϕ.

The rest of the cases S, S¯, U, U¯ can be dealt with similarly using (3) and (4). The proof that (2) implies (1) is not too difficult, and we therefore omit it.

Theorem 7 could be considered an abstract algebraic counterpart of the result of [22] (page 4268, Theorem 13) for discrete finite traces. We will discuss later how it can be used to obtain the original result (for the max-min semiring R±∞) as a corollary. Additionally, Theorem 7 gives a precise equational characterization of the class of semirings for which the relationship between the two semantics holds.

Let D be a set of data items, V be a semiring and h : V → B. For a formula ϕ : MTL(D, V ), length λ ∈ ω ∪ {ω} and i<λ, we define the concrete trace language CLh(ϕ, λ, i) = {u ∈ Dλ | u, i |=h ϕ}. For a symbolic trace **x** ∈ F(D, V )λ, we define its (concrete) trace language by CLh(**x**) = {u ∈ Dλ | u |=h **x**}, where u |=h **x** means that u(i) |=h **x**(i) for every i<n. Lemma 8 below establishes a correspondence between the symbolic and concrete language of a formula ϕ, which we need to connect Theorem 7 to the concrete setting of [22].

**Lemma 8 (Concrete and Symbolic Languages).** Let D be a set of data items, V be an idempotent semiring with top element 1, and h : V → B be a semiring homomorphism. For every formula ϕ : MTL(D, V ), length λ ∈ ω ∪ {ω}, and position i<λ, it holds that CLh(ϕ, λ, i) = **x**∈SL(ϕ,λ,i)CLh(**x**).

### **4 Relationship with robust semantics**

In this section, we consider the concrete quantitative setting where V is the maxmin semiring R±∞. We will obtain the result of [22] that relates the robustness degree with the robust semantics as a consequence of Theorem 7.

A metric space is a set M together with a function dist : M×M → R≥0, called metric, satisfying: (1) dist(x, y) = 0 iff x = y for all x, y ∈ M, (2) dist(x, y) = dist(y, x) for all x, y ∈ M, and (3) dist(x, z) ≤ dist(x, y)+dist(y, z) for all x, y, z ∈ M. Given a metric dist on M we define the distance function Dist as follows:

$$\begin{aligned} \mathsf{dist}: &M \times \mathcal{P}(M) \to \mathbb{R}\_{\geq 0}^{\infty} \\ \mathsf{dist}(x, S) = \inf\_{y \in S} \mathsf{dist}(x, y) \\ \mathsf{dist}(x, \emptyset) = \infty \end{aligned} \qquad \begin{aligned} \mathsf{Dist}: &M \times \mathcal{P}(M) \to \mathbb{R}^{\pm \infty} \\ \mathsf{Dist}(d, S) = \begin{cases} -\mathsf{dist}(d, S), & \text{if } d \notin S \\ \mathsf{dist}(d, \sim S), & \text{if } d \in S \end{cases} \end{aligned}$$

where ∼S = M \ S is the complement of S. Notice that Dist(x, ∅) = −∞.

Let D be a metric space of points (data items). Let p be a propositional letter (symbol), and O(p) ⊆ D be its interpretation, that is, the set of points for which p is true. The corresponding quantitative predicate is p : D → R±∞ given by p(d) = Dist(d, O(p)) for every d ∈ D. Given the metric dist on D, we obtain a metric dist : Dλ × Dλ → R∞ ≥0 (on the set of traces of length λ, where λ ∈ ω ∪ {ω}) as follows: dist(u, v) = supi<λ dist(u(i), v(i)). Let CLO(ϕ, n, i) = {u ∈ Dn | u, i |=O ϕ} be the set of traces (of length n) that satisfy ϕ at i (defined using the interpretation function O). Corollary 9 below was proved in [22]. We will give a proof that relies on the algebraic variant that we presented earlier.

**Corollary 9.** Let D be a set of data items, and V = R±∞. Let ϕ : MTL(D, V ), u ∈ Dn and i<n (where n ≥ 1). Then, −dist(u, CLO(ϕ, n, i)) ≤ ρ(ϕ, u, i).

Proof. We will use the semiring R±∞ ±0 ∼= B × R∞ ≥0 instead of R±∞, so that the value 0 is not ambiguous (it can be either true or false when we use R±∞). That is, we will have a positive zero +0 (true) and a negative zero −0 (false). The semiring homomorphism h : R±∞ ±0 → B sends the positive (resp., negative) elements to (resp., ⊥). We will interpret a predicate symbol p as the quantitative predicate p : D → R±∞ ±0 given by p(d) = −dist(d, O(p)) if d /∈ O(p) and p(d)=+dist(d, ∼O(p)) if d ∈ O(p). Using these definitions, the satisfaction relations |=O and |=h are the same, hence CLO and CLh are the same. Now,

$$\begin{split} \text{dist}(u, \mathbb{C}\mathcal{L}\_{h}(\varphi, n, i)) &= \text{dist}(u, \bigcup\_{\mathbf{x} \in \mathbb{S}\mathcal{L}(\varphi, n, i)} \mathbb{C}\mathcal{L}\_{h}(\mathbf{x})) \qquad \text{[Lemma 8]} \\ &= \inf\_{\mathbf{x} \in \mathbb{S}\mathcal{L}(\varphi, n, i)} \inf\_{v \in \mathbb{C}\mathcal{L}\_{h}(\mathbf{x})} \sup\_{i < n} \text{dist}(u(i), v(i)) \qquad \text{[def. of dist]} \\ &\geq \inf\_{\mathbf{x} \in \mathbb{S}\mathcal{L}(\varphi, n, i)} \sup\_{i < n} \inf\_{v \in \mathbb{C}\mathcal{L}\_{h}(\mathbf{x})} \text{dist}(u(i), v(i)) \qquad \text{[sup inf} \leq \inf \text{sup]} \\ &= \inf\_{\mathbf{x} \in \mathbb{S}\mathcal{L}(\varphi, n, i)} \sup\_{i < n} \inf\_{v(i) \in \mathcal{O}(\mathbf{x}(i))} \text{dist}(u(i), v(i)) \quad \text{[def. of } \mathbf{C}\mathcal{L}] \\ &= \inf\_{\mathbf{x} \in \mathbb{S}\mathcal{L}(\varphi, n, i)} \sup\_{i < n} \text{dist}(u(i), \mathcal{O}(\mathbf{x}(i))). \qquad \text{[def. of } \mathbf{dist}] \end{split}$$

By negating the above inequality we get that

$$-\mathsf{dist}(u,\mathsf{CL}\_h(\varphi,n,i)) \le \sup\_{\mathbf{x}\in\mathsf{SL}(\varphi,n,i)} \inf\_{i$$

which is ≤ - **x**∈SL(ϕ,n,i)**x**[u] = val(ϕ, u, i). From Theorem 7 we get val(ϕ, u, i) ≤ ρ(ϕ, u, i) and therefore −dist(u, CLO(ϕ, n, i)) ≤ ρ(ϕ, u, i).

From Corollary 9 we can also obtain ρ(ϕ, u, i) ≤ dist(u, ∼CLO(ϕ, n, i)). This inequality is equivalent to −dist(u, ∼CLO(ϕ, n, i)) ≤ −ρ(ϕ, u, i), which in turn is equivalent to −dist(u, CLO(∼ϕ, n, i)) ≤ ρ(∼ϕ, u, i). The operation ∼ on formulas is a pseudo-negation, that is, ∼ϕ is the formula that results by "dualizing" all connectives and negating the atomic predicates. This operation is meaningful for the semiring R±∞. The final inequality is an instance of Corollary 9 for ∼ϕ.

Theorem 7 and Corollary 9 are not used later for the monitoring algorithm. The significance of our theorem is that it can be instantiated to give the existing result from [22]. This serves as a sanity check for our algebraic framework and it supports the semiring-based semantics of Sect. 2.

### **5 Online Monitoring**

For an infinite input trace u ∈ Dω, the output of the monitor for the time instant t should be ρ(ϕ, u, t), but the monitor has to compute it by observing only a finite prefix of u. In order for the output value of the monitor to agree with the standard temporal semantics over infinite traces we may need to delay an output item until some part of the future input is seen. For example, in the case of F1p we need to wait for one time unit: the output at time t is given after the input item at time t + 1 is seen. In other words, the monitor for F1p has a delay (the output is falling behind the input) of one time unit. Symmetrically, we can allow monitors to emit output early when the correct value is known. For example, the output value for P1p is 0 in the beginning and the value at time t is already known from time t − 1. So, we also allow monitors to have negative delay (the output is running ahead of the input). The function dl : MTL → Z gives the amount of delay required to monitor a formula. It is defined by dl(p) = 0 and

$$\begin{aligned} \mathsf{dl}(\varphi \wedge \psi) &= \max(\mathsf{dl}(\varphi), \mathsf{dl}(\psi)) & \mathsf{dl}(\varphi \ \mathsf{S}\_{[a,b]} \ \psi) &= \max(\mathsf{dl}(\varphi), \mathsf{dl}(\psi)) - a \\ \mathsf{dl}(\varphi \ \mathsf{S}\_{[a,\infty)} \ \psi) &= \max(\mathsf{dl}(\varphi), \mathsf{dl}(\psi)) - a & \mathsf{dl}(\varphi \ \mathsf{U}\_{[a,b]} \ \psi) &= \max(\mathsf{dl}(\varphi), \mathsf{dl}(\psi)) + b. \end{aligned}$$

The monitor TL(ϕ) for a formula ϕ is a variant of a Mealy machine. If dl(ϕ) = 0, the TL(ϕ) is precisely a Mealy machine (one output item per input item) with inputs D and outputs V . If = dl(ϕ) > 0, then TL(ϕ) emits no output for the first steps and then behaves like a Mealy machine. If = dl(ϕ) < 0, then TL(ϕ) emits items upon initialization and continues to behave like a Mealy machine.

Let A and B be sets. A monitor of type M(A, B) is a state machine G = (St, init, o, next, out), where St is a set of states, init ∈ St is the initial state, o ∈ B∗ is the initial output, next : St × A → St is the transition function, and out : St × A → Opt(B) is the output function, where Opt(B) = B ∪ {nil}.

```
map(op) : M(A, B)
                    St = Unit
                   init = u
                     o = ε
            next(s, a) = s
             out(s, a) = op(a)
                                         aggr(b, op) : M(A, B)
                                                   St = B
                                                 init = b
                                                   o = ε
                                           next(s, a) = op(s, a)
                                            out(s, a) = op(s, a)
                                                                        emit(n, v) : M(A, A)
                                                                                St = Unit
                                                                               init = u
                                                                                o = vn
                                                                         next(s, a) = s
                                                                          out(s, a) = a
ignore(n) : M(A, A)
        St = [0, n]
       init = 0
         o = ε
next(s, a) = s + 1, if s<n
next(s, a) = s, if s = n
 out(s, a) = nil, if s<n
 out(s, a) = a, if s = n
                                wnd(n, v, op) : M(A, A)
                                           St = Buf(A)
                                          init = Buf(n, v)
                                            o = ε
                                   next(s, a) = s.ins(a)
                                    out(s, a) = s.ins(a).agg(op)
                                                                    wndV(n, op) : M(A, A)
                                                                              St = Buf(A)
                                                                             init = Buf()
                                                                               o = ε
                                                                      next(s, a) = s.ins(a)
                                                                       out(s, a) = ε, if size(s) < n − 1
                                                                       out(s, a) = s.ins(a).agg(op), o/w
```
Fig. 3: Basic building blocks for constructing temporal quantitative monitors.

In Fig. 3 we give several examples of simple monitors that can be used as building blocks. The monitor map(op) applies the function op : A → B elementwise. The monitor aggr(b, op) applies a running aggregation to the input trace that is specified by the initial aggregate b : B and the aggregation function op : B×A → B (similar to the fold combinator used in functional programming). The monitor emit(n, v) emits n copies of the value v ∈ A upon initialization and then echoes the input trace. The monitor ignore(n) discards the first n items of the trace and proceeds to echo the rest of the trace. The monitor wnd(n, v, op) performs an aggregation, given by the associative function op : A×A → A, over a sliding window of size n. It initializes the window using the value v : A and emits output at the arrival of every item. The monitor wndV(n, op) is different in that it starts with an empty window and it only starts emitting output when the window fills up with n items. We will combine monitors using the operations serial composition >> and parallel composition par. In the serial composition G >> H the output trace of G is propagated as input trace to H. In the parallel composition par(G, H) the input trace to copied to two concurrently executing monitors G and H and their output traces are combined. Both combinators >> and par are given by variants of the product construction on state machines. In the case of par the output traces of G and H may not be synchronized (one may be ahead of the other), which requires some bounded buffering in order to properly align them. The construction for par is described in [37]. Some variants of the combinators of Figure 3 are part of the StreamQL language [29], which has been proposed for the processing of streaming time series.

The identities of Fig. 2 suggest that MTL monitoring can be reduced to a small set of computational primitives. In fact, the primitives described earlier are sufficient to specify the monitors, as shown in Fig. 4. We write π1 : A × B → A for the left projection and π2 : A × B → B for the right projection.

Let u ∈ D+ and n = |u|. If n>a then ρ(ϕ S[0,a] ψ, u, n − 1) = ρ(ϕ S ψ, v, a), where v is the suffix of u with a + 1 items. If n ≤ a then ρ(ϕ S[0,a] ψ, u, n − 1) =

```
TL(p) = map(p)
    TL(ϕ ∨ ψ) = par(TL(ϕ), TL(ψ)) >> map(+)
       TL(Pϕ) = TL(ϕ) >> aggr(0, +)
      TL(Paϕ) = TL(ϕ) >> emit(a, 0)
  TL(P[a,∞)ϕ) = TL(PaPϕ)
     TL(ϕ S ψ) = par(TL(ϕ), TL(ψ)) >> aggr(0, opS)
          opS : V × (V × V ) → V , where
 opS(s, x, y)=(s · x) + y
TL(ϕ S[a,∞) ψ) = TL(Pa(ϕ S ψ) ∧ H[0,a−1]ϕ)
 TL(ϕ S[0,b] ψ) = par(TL(ϕ), TL(ψ)) >>
                 wnd(b + 1, 0, ⊗S) >> map(π2)
 TL(ϕ S[a,b] ψ) = TL(Pa(ϕ S[0,b−a] ψ) ∧ H[0,a−1]ϕ)
      TL(Faϕ) = TL(ϕ) >> ignore(a)
   TL(F[a,b]ϕ) = TL(FbP[0,b−a]ϕ)
 TL(ϕ U[0,b] ψ) = par(TL(ϕ), TL(ψ)) >>
                 wndV(b + 1, ⊗U) >> map(π2)
TL(ϕ U[a,b] ψ) = TL(G[0,a−1]ϕ ∧ Fa(ϕ U[0,b−a] ψ))
                                                       // fill buffer with v (initial values)
                                                       T[n] buf ← [n; v]
                                                       // calculate partial aggregates
                                                       for i ← n − 2 to 0 do
                                                           buf [i] ← buf [i] ⊗ buf [i + 1]
                                                       // initial total aggregate
                                                       T agg ← buf [0]
                                                       Nat m ← 0 // size of new block
                                                       T z ← nil // aggregate of new block
                                                       Function Add(T d):
                                                           if m = n then // full new block
                                                                // convert new block to old
                                                                for i ← n − 2 to 1 do
                                                                    buf [i] ← buf [i] ⊗ buf [i + 1]
                                                                m ← 0 // empty new block
                                                                z ← nil
                                                           // evict oldest item, replace with d
                                                           buf [m] ← d
                                                           m ← m + 1 // new block enlarged
                                                           z ← z ⊗ d // where nil ⊗ d = d
                                                           if m<n then
                                                                agg ← buf [m] ⊗ z
                                                           else // m = n agg ← z
```
Fig. 4: Online monitors for bounded-future MTL formulas & sliding aggregation.

ρ(ϕ S ψ, 0a+1−nu, a). So, we can implement a monitor for the connective S[0,a] by computing S over a window of exactly a + 1 data items.

**Proposition 10 (Aggregation for** S**,** U**).** Let V be a semiring. For every trace u = u0u1 ...un−1 ∈ (V ×V )+ of length n = |u|, the values ρ(π1 Sπ2, u, n−1) and ρ(π1 Uπ2, u, 0) can be written as aggregates of the form π2(u0 ⊗u1 ⊗···⊗un−1).

Proposition 10 justifies the translation of S[0,b]/U[0,b] into monitors (Fig. 4). Now, we will describe the data structure that performs the sliding aggregation. It is used in Fig. 3 in the monitors wnd and wndV. The implementation is shown in Fig. 4. Suppose that the current window (of size n) is [x0, x1,...,xn−1]. We maintain a buffer of the form [xn−m,...,xn−1, y0,...,yn−1−m], where the part [xn−m,...,xn−1] is the block of newer elements ("new block") and the part [y0,...,yn−1−m] contains aggregates of the older elements ("old block"). They satisfy the invariant yi = xi⊗···⊗xn−1−m for every i = 0,...,n−1−m. We also maintain the aggregate z = xn−m ⊗···⊗ xn−1 of the new block. So, the overall aggregate of the window is agg = y0⊗z. When a new item d arrives, we evict the aggregate y0 corresponding to the oldest item x0 and replace it by d. Thus, the new block is expanded with the additional item d and therefore we also update the aggregates z and agg. When the new block becomes full (i.e., m = n) then we convert it to an old block by performing all partial aggregations from right to left. This conversion requires n−1 applications of ⊗, but it is performed once every n items. So, the algorithm needs O(1) amortized time-per-item.

**Theorem 11.** Let D be a set of data items, V be a semiring, and ϕ : MTL(D, V ) be a bounded-future formula. The monitor TL(ϕ) : M(D, V ) is a streaming algorithm that needs O(2|ϕ| ) space and O(|ϕ|) amortized time-per-item.

Proof. The algorithm needs space that is exponential in the size of ϕ because of the connectives of the form X[a,∞) and X[a,b]. The monitor uses buffers of size a or b − a. Since the constants a, b are written in binary notation, we need space that is exponential in the size. The O(|ϕ|) amortized time per element hinges on the algorithm of Fig. 4, which is used for S[0,b] and U[0,b]. As discussed earlier, this algorithm needs O(1) amortized time-per-item.

### **6 Experimental Evaluation**

We have implemented our semiring-based monitoring framework in Rust. We compare our implementation with the verified lattice-based monitors of [13] and the monitoring tool Reelay [40]. We perform our experiments using the (R±∞, max, min) semiring for truth values, which are approximately represented using 64-bit floating-point numbers.

We have observed that all three tools process items at a roughly constant rate. We summarize the performance of a monitor with the average time it takes to process one data item (i.e., amortized time-per-item). In Fig. 5, we consider formulas X[0,n], Xn, X[n,2n], X[n,∞) where X ∈ {S, P}. We show the time-per-item for the monitors for n = 1, 10, 102, 103, 104, 105, 106. We have also evaluated how the monitors for future temporal connectives scale with respect to the constants in the intervals. In Fig. 6, we benchmark all tools using formulas from the Timescales benchmark [39]. Our monitors are generally more than 100 (resp., 10) times faster than Reelay (resp., the lattice-based tool of [13]).

The profiling tools Valgrind [38] and Heaptrack [41] are used to analyze the memory consumption of the monitors. Our Rust implementation, given a formula, begins by allocating a fixed amount of memory and does not allocate any more memory during the rest of the computation. Reelay allocates and de-allocates memory throughout its execution. The lattice-based monitor is implemented in OCaml (which is a garbage-collected language) and consumes a larger amount of memory. In Fig. 5, we plot the peak memory usage of the monitors. We note that our tool does not seem to be allocating an increasing amount of memory for Pn and similar formulas. This is because the corresponding monitor for Pn emits output as early as possible and therefore does not need to use a buffer. In the case of the lattice-based monitor and our tool, we observe that the memory consumption does not depend on the input trace (it only depends on the formula). In the case of Reelay, it appears that the memory consumption depends on the input trace. We have plotted the behavior for two different input traces: one that consists of an increasing sequence of values ("reelay-ascending"), and another one that is decreasing ("reelay-descending"). We have only measured the memory usage of Reelay for up to n = 213, as the execution becomes very slow beyond this value.

We use *case studies* from the automotive domain, which have been suggested as benchmarks for hybrid system verification [25]. The **A**utomatic Transmission System has two input signals (a throttle and a break) and three output signals: the gear sequence (gi for each gear i), the engine rotation speed (in

Fig. 5: Microbenchmark

rpm, denoted ω) and the vehicle speed (denoted v). Based on the suggestions in [25], we consider five properties: A1 = ω < ω, A2 = (ω < ω) ∧ (v < v), A3 = g1 ∧ Y(g2) → YH[0,2.5]g2 (where Y is notation for P1), A4 = H[0,5](ω < ω) → H[0,2.5](v < v) and A5 = (v > v) S¯[0,1] ((ω > ω) S¯[0,2] ((¬g4) S¯[0,10] ((¬g3) S¯ ((¬g2) S¯ (¬g1))))). All constants in the temporal connectives are in seconds, and we choose the constants v = 120 and ω = 4500. Formula A3 says that before changing from the second to the first gear, at least 2.5 seconds must first pass. Formula A4 says that keeping the engine speed low enough should ensure that the vehicle does not exceed a certain speed. Formula A5 says that changing the gear from the first to the fourth within 10 seconds, and then having the engine speed exceed ω will cause the vehicle speed to exceed v. The other case study is a **F**ault-Tolerant Fuel Control System. We monitor two properties. The first is

Fig. 6: Macrobenchmarks

that the fuel flow rate should frequently become and remain non-zero for a sufficient amount of time. We encode this as F1 = H[0,10]P[0,1](F uelFlowRate > 0). The other property is to ensure that whenever the air-to-fuel ratio goes out of bounds, then within 1 second it should settle back and stay there for a second. This is written as F2 = (H[0,1]airF uelRatio < 1) S¯[0,2] airF uelRatio < 1. The experimental results are shown in Fig. 6.

All of our experiments were executed on a laptop with an Intel Core i7 10610U CPU clocked at 2.30GHz and 16GB of memory. Each value reported is the mean of 20 executions of the experiment. The whiskers in the plots indicate the standard deviation across all executions.

### **7 Related Work**

Fainekos and Pappas [22] define the robustness degree of satisfaction in terms of the distance of the signal from the set of desirable ones (or its complement). They also suggest an under-approximation of the robustness degree which can be effectively monitored. This is called the robust semantics and is defined by induction on STL formulas, by interpreting conjunction (resp., disjunction) as min (resp., max) of R±∞. Our paper explores this robust semantics (and the related approximation guarantee) in the general algebraic setting of semirings.

In [27], the authors study a generalization of the robustness degree by considering idempotent semirings of real numbers. They also propose an online monitoring algorithm that uses symbolic weighted automata. While this approach computes the precise robustness degree in the sense of [22], the construction of the relevant automata incurs a doubly exponential blowup if one considers STL specifications. In [13], it is observed that an extension of the robust semantics to bounded distributive lattices can be effectively monitored. In this paper, we generalize this semantics by considering semirings (bounded distributive lattices are semirings). Semirings are also used in [9], where the authors consider a spatio-temporal logic. They consider the class of constraint semirings, which require the semiring order to induce a complete lattice. Efforts have been made to define notions of robustness that take temporal discrepancies into account. In [20], we see a definition of temporal robustness by considering the effect of shifting the signal in time. The "edit distance" between discretized signals is proposed as a measure of robustness in [26]. Abbas et al. [3] define a notion of (τ,ε) closeness between signals, which considers temporal and value-based guarantees separately. In [2], a metric based on conformance is put forward for applications in cardiac electrophysiology. Averaging temporal operators are used in [5], which assign a higher value to temporal obligations that are satisfied earlier.

A key ingredient for the efficient monitoring of STL is a streaming algorithm for sliding-window maximum [19,15]. The tool Breach [17,18], which is used for the falsification of temporal specifications over hybrid systems, uses the sliding-maximum algorithm of [32]. In contrast, we use a more general sliding aggregation which applies to any associative operation (not only max/min) and does not require the truth values to be totally ordered.

Different approaches for interpreting future temporal connectives in the context of online monitoring have been studied. While [16] assumes the availability of a predictor to interpret future connectives, [21] considers robustness intervals: the tightest intervals which cover the robustness for all possible extensions of the available trace prefix. Reelay [40] exclusively uses past-time connectives. The transducer-based framework of [37] can be used to monitor rich temporal properties which depend on bounded future input by allowing some bounded delay in the output.

There is a large amount of work on formalisms, domain-specific languages and associated tools for quantitative online monitoring and, more generally, for data stream processing. The synchronous language LOLA [14] has served as the basis for the StreamLAB tool [23], which is used for monitoring cyber-physical systems. Quantitative Regular Expressions [36] and associated automata-theoretic models with registers [7,8,6] have been used to express complex online detection algorithms for medical monitoring [1,4]. There are many synchronous languages and models of computation based on Kahn's dataflow model [28] that have been used for signal processing [31] and embedded controller design [12,11,10]. The construction of online monitors described in Sect. 5 relies on a set of combinators that constitute a simple domain-specific language for stream processing. Our focus here, however, is on providing efficient monitors for MTL formulas with a quantitative semantics, rather than designing a general-purpose language for monitor specification. The compositional construction of automata-based monitors from temporal specifications has also been considered in [34,35,24].

### **8 Conclusion**

We have presented a new efficient algorithm for the online monitoring of MTL properties over discrete traces. We have used an abstract algebraic semantics based on semirings, which can be instantiated to the widely-used Boolean (qualitative) and robustness (quantitative) semantics, as well as to other partially ordered semirings. We also provide a theorem that relates our quantitative semantics with an algebraic generalization of the robustness degree of [22]. We have provided an implementation of our algebraic monitoring framework, and we have shown experimentally that our monitors scale reasonably well and are competitive against the state-of-the-art tool Reelay [40].

### **References**

ing of synchronous systems. In: TIME 2005. pp. 166–174. IEEE (2005). https://doi.org/10.1109/TIME.2005.26

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/ 4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Neural Networks**

# **Synthesizing Context-free Grammars from Recurrent Neural Networks**

Daniel M. Yellin1 and Gail Weiss2

1 IBM, Givatayim, Israel<EMAIL_ADDRESS>2 Technion, Haifa, Israel<EMAIL_ADDRESS>
**Abstract.** We present an algorithm for extracting a subclass of the context free grammars (CFGs) from a trained recurrent neural network (RNN). We develop a new framework, *pattern rule sets* (PRSs), which describe sequences of deterministic finite automata (DFAs) that approximate a non-regular language. We present an algorithm for recovering the PRS behind a sequence of such automata, and apply it to the sequences of automata extracted from trained RNNs using the L∗ algorithm. We then show how the PRS may converted into a CFG, enabling a familiar and useful presentation of the learned language.

Extracting the learned language of an RNN is important to facilitate understanding of the RNN and to verify its correctness. Furthermore, the extracted CFG can augment the RNN in classifying correct sentences, as the RNN's predictive accuracy decreases when the recursion depth and distance between matching delimiters of its input sequences increases.

**Keywords:** Model Extraction · Learning Context Free Grammars · Finite State Machines · Recurrent Neural Networks

### **1 Introduction**

Recurrent Neural Networks (RNNs) are a class of neural networks adapted to sequential input, enjoying wide use in a variety of sequence processing tasks. Their internal process is opaque, prompting several works into extracting interpretable rules from them. Existing works focus on the extraction of deterministic or weighted finite automata (DFAs and WFAs) from trained RNNs [18,6,26,3].

However, DFAs are insufficient to fully capture the behavior of RNNs, which are known to be theoretically Turing-complete [20], and for which there exist architecture variants such as LSTMs [14] and features such as stacks [9,23] or attention [4] increasing their practical power. Several recent investigations explore the ability of different RNN architectures to learn Dyck, counter, and other non-regular languages [19,5,28,21], with mixed results.

While the data indicates that RNNs can generalize and achieve high accuracy, they do not learn hierarchical rules, and generalization deteriorates as the length and 'depth' of the input grows [19,5,28]. Sennhauser and Berwick conjecture that © The Author(s) 2021

J. F. Groote and K. G. Larsen (Eds.): TACAS 2021, LNCS 12651, pp. 351–369, 2021. https://doi.org/10.1007/978-3-030-72016-2 19

**Fig. 1.** Overview of steps in algorithm to synthesize the hidden language L

"what the LSTM has in fact acquired is sequential statistical approximation to this solution" instead of "the 'perfect' rule-based solution" [19]. Similarly, Yu et. al. conclude that "the RNNs can not truly model CFGs, even when powered by the attention mechanism" [28]. This is line with Hewitt et. al., who note that a fixed precision RNN can only learn a language of fixed depth strings [13].

Goal of this paper We wish to extract a CFG from a trained RNN. In particular, we wish to find the CFG that not only explains the finite language learnt by the RNN, but generalizes it to strings of unbounded depth and distance.

Our approach Our method builds on the DFA extraction work of Weiss et al. [26], which uses the L∗ algorithm [2] to learn the DFA of a given RNN. As part of the learning process, L∗ creates a sequence of hypothesis DFAs approximating the target language. Our main insight is in treating these hypothesis DFAs as coming from a set of underlying rules, that recursively improve each DFA's approximation of the target CFG by increasing the distance and embedded depth of the sequences it can recognize. In this light, synthesizing the target CFG becomes the problem of recovering these rules.

We propose the framework of pattern rule sets (PRSs) for describing such rule applications, and present an algorithm for recovering a PRS from a sequence of DFAs. We also provide a method for converting a PRS to a CFG, and test our method on RNNs trained on several PRS languages. Pattern rule sets are expressive enough to cover several variants of the Dyck languages, which are prototypical context-free languages (CFLs): the Chomsky–Sch¨utzenberger representation theorem shows that any CFL can be expressed as a homomorphic image of a Dyck language intersected with a regular language[16].

A significant issue we address is that the extracted DFAs are often inexact, either through inaccuracies in the RNN, or as an artifact of the L∗ algorithm.

To the best of our knowledge, this is the first work on synthesizing a CFG from a general RNN (though some works extract push-down automata [23,9] from RNNs with an external stack, they do not apply to plain RNNs). The overall steps in our technique are given in Figure 1.

Contributions The main contributions of this paper are:

**–** An implementation of our technique1, and an evaluation of its success on recovering various CFLs from trained RNNs.

### **2 Definitions and Notations**

#### **2.1 Deterministic Finite Automata**

**Definition 1 (Deterministic Finite Automata).** A deterministic finite automaton (DFA) over an alphabet Σ is a 5-tuple Σ,q0, Q, F, δ such that Q is a finite set of states, q0 ∈ Q is the initial state, F ⊆ Q is a set of final (accepting) states and δ : Q × Σ → Q is a (possibly partial) transition function.

Unless stated otherwise, we assume each DFA's states are unique to itself, i.e., for any two DFAs A, B – including two instances of the same DFA – QA ∩QB = ∅. A DFA A is said to be complete if δ is complete, i.e., the value δ(q, σ) is defined for every q, σ ∈ Q × Σ. Otherwise, it is incomplete.

We define the extended transition function ˆδ : Q × Σ∗ → Q and the language L(A) accepted by A in the typical fashion. We also associate a language with intermediate states of A: L(A, q1, q2) {w ∈ Σ∗ | ˆδ(q1, w) = q2}. The states from which no sequence w ∈ Σ∗ is accepted are known as the sink reject states.

**Definition 2.** The sink reject states of a DFA A = Σ,q0, Q, F, δ are the maximal set QR ⊆ Q satisfying: QR ∩ F = ∅, and for every q ∈ QR and σ ∈ Σ, either δ(q, σ) ∈ QR or δ(q, σ) is not defined.

**Definition 3 (Defined Tokens).** Let A = Σ,q0, Q, F, δ be a complete DFA with sink reject states QR. For every q ∈ Q, its defined tokens are def(A, q) {σ ∈ Σ | δ(q, σ) ∈/ QR}. When the DFA A is clear from context, we write def(q).

All definitions for complete DFAs are extended to incomplete DFAs A by considering their completion - an extension of A in which all missing transitions are connected to a (possibly new) sink reject state.

**Definition 4 (Set Representation of** δ**).** A (possibly partial) transition function δ : Q×Σ → Q may be equivalently defined as the set Sδ = {(q, σ, q ) | δ(q, σ) = q }. We use δ and Sδ interchangeably.

**Definition 5 (Replacing a State).** For a transition function δ : Q × Σ → Q, state q ∈ Q, and new state qn ∈/ Q, we denote by δ[q←qn] : Q × Σ → Q the transition function over Q = (Q \ {q}) ∪ {qn} and Σ that is identical to δ except that it redirects all transitions into or out of q to be into or out of qn.

1 The implementation for this paper, and a link to all trained RNNs, is available at https://github.com/tech-srl/RNN to PRS CFG.

#### **2.2 Dyck Languages**

A Dyck language of order N is expressed by the grammar D ::= ε | L1 D R1 | ... | LN D RN |DD with unique symbols L1,...,LN ,D1,...,DN . A common measure of complexity for a Dyck word is its maximum distance (number of characters) between matching delimiters and embedded depth (number of unclosed delimiters) [19]. We generalize and refer to Regular Expression Dyck (RE-Dyck) languages as languages expressed by the same CFG, except that each Li and each Ri derive some regular expression.

We present regular expressions as is standard, for example: L({a|b}·c) {ac,bc}.

### **3 Patterns**

Patterns are DFAs with a single exit state qX in place of a set of final states, and with no cycles on their initial or exit states unless q0 = qX.

**Definition 6 (Patterns).** A pattern p = Σ,q0, Q, qX, δ is a DFA Ap = Σ,q0, Q, {qX}, δ, satisfying: 1. L(Ap) = ∅, and 2. either q0 = qX, or def(qX) = ∅ and L(A, q0, q0) = {ε}. If q0 = qX then p is called circular, otherwise, it is noncircular. Patterns are always given in minimal incomplete presentation.

We refer to a pattern's initial and exit states as its edge states. All the definitions for DFAs apply to patterns through Ap. We denote each pattern p's language Lp L(p), and if it is marked by some superscript i, we refer to all of its components with superscript i: pi = Σ,qi 0, Qi , qi X, δi .

### **3.1 Pattern Composition**

We can compose two non-circular patterns p1, p2 by merging the exit state of p1 with the initial state of p2, creating a new pattern p3 satisfying Lp3 = Lp1 ·Lp2 .

**Definition 7 (Serial Composition).** Let p1, p2 be two non-circular patterns. Their serial composite is the pattern p1 ◦ p2 = Σ,q1 0, Q, q2 X, δ in which Q = Q1 ∪ Q2 \ {q1 X} and δ = δ1 [q1 X←q2 0] ∪ δ2. We call q2 0 the join state of this operation.

If we additionally merge the exit state of p2 with the initial state of p1, we obtain a circular pattern p which we call the circular composition of p1 and p2. This composition satisfies Lp = {Lp1 ·Lp2 }∗.

**Definition 8 (Circular Composition).** Let p1, p2 be two non-circular patterns. Their circular composite is the circular pattern p1 ◦c p2 = Σ,q1 0, Q, q1 0, δ in which Q = Q1 ∪ Q2 \ {q1 X, q2 X} and δ = δ1 [q1 X←q2 0] ∪ δ2 [q2 X←q1 0] . We call q2 0 the join state of this operation.

Figure 2 shows 3 examples of serial and circular compositions of patterns.

Patterns do not carry information about whether or not they have been composed from other patterns. We maintain such information using pattern pairs.

**Fig. 2.** Examples of the composition operator

**Definition 9 (Pattern Pair).** A pattern pair is a pair P, Pc of pattern sets, such that Pc ⊂ P and for every p ∈ Pc there exists exactly one pair p1, p2 ∈ P satisfying p = p1 0 p2 for some 0 ∈ {◦, ◦c}. We refer to the patterns p ∈ Pc as the composite patterns of P, Pc, and to the rest as its base patterns.

We will often discuss patterns that have been composed into larger DFAs.

**Definition 10 (Pattern Instances).**  We say that pˆ is an instance of p in A if pˆ is isomorphic to p.

A pattern instance in a DFA A is uniquely determined by its structure and initial state: (p, q). If p is a composite pattern with respect to some pattern pair P, Pc, the join state of its composition within A is also uniquely defined.

**Definition 11.** For every pattern pair P, Pc, for each composite pattern p ∈ Pc, DFA A, and initial state q of an instance pˆ of p in A, join(p, q, A) returns the join state of pˆ with respect to its composition in P, Pc.

### **4 Pattern Rule Sets**

For any infinite sequence S = A1, A2, ... of DFAs satisfying L(Ai) ⊂ L(Ai+1), for all i, we define the language of S as the union of the languages of all these DFAs: L(S) = ∪iL(Ai). Such sequences may be used to express CFLs.

In this work we take a finite sequence A1, A2, ..., An of DFAs, and assume it is a (possibly noisy) finite prefix of an infinite sequence of approximations for a language, as above. We attempt to reconstruct the language by guessing how the

sequence may continue. To allow such generalization, we must make assumptions about how the sequence is generated. For this we introduce pattern rule sets.

Pattern rule sets (PRSs) create sequences of DFAs with a single accepting state. Each PRS is built around a pattern pair P, Pc, and each rule application connects a new pattern instance to the current DFA Ai, at the join state of a composite-pattern inserted into Ai at some earlier point. To define where a pattern can be connected to Ai, we introduce an enabled instance set I.

**Definition 12.** An enabled DFA over a pattern pair P, Pc is a tuple A, I such that A = Σ,q0, Q, F, δ is a DFA and I ⊆ Pc × Q marks enabled instances of composite patterns in A.

Intuitively, for every enabled DFA A, I and (p, q) ∈ I, we know: (i) there is an instance of pattern p in A starting at state q, and (ii) this instance is enabled; i.e., we may connect new pattern instances to its join state join(p, q, A).

**Definition 13.** A PRS **P** is a tuple Σ, P, Pc, R where P, Pc is a pattern pair over the alphabet Σ and R is a set of rules. Each rule has one of the following forms, for some p, p1, p2, p3, pI ∈ P, with p1 and p2 non-circular:

(1) ⊥ pI (2) p c (p1 0 p2)◦= p3, where p = p1 0 p2 for 0 ∈ {◦, ◦c}, and p3 is circular (3) p s (p1 ◦ p2)◦= p3, where p = p1 ◦ p2 and p3 is non-circular

A PRS derives sequences of enabled DFAs as follows: first, a rule of type (1) creates A1, I1 according to pI . Then, for every Ai, Ii, each rule may connect a new pattern instance to Ai, specifically at a state determined by Ii.

**Definition 14 (Initial Composition).** D1 = A1, I1 is generated from a rule ⊥ pI as follows: A1 = ApI , and Ii = {(pI , qI 0)} if pI ∈ Pc and otherwise I1 = ∅.

Let Di = Ai, Ii be the enabled DFAat step i and denote Ai = Σ,q0, Q, F, δ. Note that for A1, |F| = 1, and for all Ai+1, F is unchanged (by future definitions).

Rules of type (1) extend Ai by grafting a circular pattern to q0, and then enabling that pattern if it is composite.

**Definition 15 (Rules of type** (1)**).** A rule ⊥ pI with circular pI may extend Ai, Ii at the initial state q0 of Ai iff def(q0)∩def(qI 0) = ∅. This creates the DFA Ai+1 = Σ,q0, Q∪QI \{qI 0}, F, δ∪δI [qI 0←q0] . If pI ∈ Pc then Ii+1 = Ii∪{(pI , q0)}, else Ii+1 = Ii.

Rules of type (2) graft a circular pattern p3 = Σ,q3 0, q3 x, F, δ3 onto the join state qj of an enabled pattern instance pˆ in Ai, by merging q3 0 with qj . In doing so, they also enable the patterns composing ˆp, if they are composite.

**Definition 16 (Rules of type** (2)**).** A rule p c (p1 0 p2)◦= p3 may extend Ai, Ii at the join state qj = join(p, q, Ai) of any instance (p, q) ∈ Ii, provided def(qj ) ∩ def(q3 0) = ∅. This creates Ai+1, Ii+1 as follows: Ai+1 = Σ,q0, Q ∪ Q3 \ q3 0, F, δ ∪ δ3 [q3 0←qj ] , and Ii+1 = Ii ∪ {(pk, qk) | pk ∈ Pc, k ∈ {1, 2, 3}}, where q1 = q and q2 = q3 = qj .

**Fig. 3.** Structure of DFA after applying rule of type 2 or type 3

Example applications of rule (2) are shown in Figures 3(i) and 3(ii).

We also wish to graft a non-circular pattern p3 between p1 and p2, but this time we must avoid connecting the exit state q3 X to qj lest we loop over p3 multiple times. We therefore replicate the outgoing transitions of qj in p1 ◦ p2 to the inserted state q3 X so that they may act as the connections back into the DFA.

**Definition 17 (Rules of type** (3)**).** A rule p s (p1 ◦ p2)◦= p3 may extend Ai, Ii at the join state qj = join(p, q, Ai) of any instance (p, q) ∈ Ii, provided def(qj ) ∩ def(q3 0) = ∅. This creates Ai+1, Ii+1 as follows: Ai+1 = Σ,q0, Q ∪ Q3 \ q3 0, F, δ ∪ δ3 [q3 0←qj ] ∪ C where C = { (q3 X, σ, δ(qj , σ))| σ ∈ def(p2, q2 0)}, and Ii+1 = Ii ∪ {(pk, qk) | pk ∈ Pc, k ∈ {1, 2, 3}} where q1 = q and q2 = q3 = qj .

We call C the connecting transitions. We depict this rule application in example in Fig. 3 (iii), in which a member of C is labeled 'c'.

Multiple applications of rules of type (3) to the same instance pˆ will create several equivalent states in the resulting DFAs, as all of their exit states will have the same connecting transitions. These states are merged in a minimized representation, as depicted in Diagram (iv) of Figure 3.

We write A ∈ G(**P**) if there exists a sequence of enabled DFAs derived from **P** s.t. A = Ai for some Ai in this sequence.

**Definition 18 (Language of a PRS).** The language of a PRS **P** is the union of the languages of the DFAs it can generate: L(**P**) = ∪A∈G(**P**)L(A).

#### **4.1 Examples**

Example 1: Let p1 and p2 be the patterns accepting 'a' and 'b' respectively. Consider the PRS Rab with rules, ⊥ p1 ◦ p2 and p1 ◦ p2 s (p1 ◦ p2)◦= (p1 ◦ p2). This PRS creates only one sequence of DFAs. Once the first rule creates the initial DFA, by continuously applying the second rule we obtain the infinite sequence of DFAs each satisfying L(Ai) = {aj bj : 1 ≤ j ≤ i}, and so L(Rab) = {ai bi : i > 0}. Figure 2(i) presents A1, while A2 and A3 appear in Figure 4(i). We can substitute any non-circular patterns for p1 and p2, creating the language {xi yi : i > 0} for any non-circular pattern regular expressions x and y.

**Fig. 4.** DFA sequences for Rab and RDyck2

Example 2: Let p1,p2,p4, and p5 be the non-circular patterns accepting '(', ')', '[', and ']' respectively. Let p3 = p1 ◦c p2 and p6 = p4 ◦c p5. Let RDyck2 be the PRS containing rules ⊥ p3, ⊥ p6, p3 c (p1 ◦c p2)◦= p3, p3 c (p1 ◦c p2)◦= p6, p6 c (p4 ◦c p5)◦= p3, and p6 c (p4 ◦c p5)◦= p6. RDyck2 defines the Dyck language of order 2. Figure 4 (ii) shows one of its possible DFA-sequences.

### **5 PRS Inference Algorithm**

A PRS can generate a sequence of DFAs defining, in the limit, a context-free language. We are now interested in inverting this process: given a sequence of DFAs generated by a PRS **P**, can we reconstruct **P**? Coupled with an L∗ extraction of DFAs from a trained RNN, solving this problem will enable us to extract a PRS from an RNN – provided the extraction follows a PRS (as we often find it does).

We present an algorithm for this problem, and show its correctness. In practice the DFAs we are given are not "perfect"; they contain noise that deviates from the PRS. We therefore augment this algorithm, allowing it to operate smoothly even on imperfect DFA sequences created from RNN extraction.

In the following, for each pattern instance pˆ in Ai, we denote by p the pattern that it is an instance of. We use similar notation pˆ1, pˆ2, and pˆI to refer to specific instances of patterns p1, p2 and pI . Additionally, for each consecutive DFA pair Ai and Ai+1, we refer by ˆp3 to the new pattern instance in Ai+1.

Main steps of inference algorithm. Given a sequence of DFAs S = A1 ··· An, the algorithm infers **P** = Σ, P, Pc, R in the following stages:

1. Discover the initial pattern instance pˆI in A1. Insert pI into P and mark pˆI as enabled. Insert the rule ⊥ → pI into R.

	- (a) Discover the new pattern instance ˆp3 in Ai+1 that extends Ai.
	- (b) If pˆ3 starts at the state q0 of Ai+1, then it is an application of a rule of type (1). Insert p3 into P, mark ˆp3 as enabled, and add ⊥ p3 to R.
	- (c) Otherwise (pˆ3 does not start at q0), find the unique enabled pattern pˆ = pˆ1 0 pˆ2 in Ai s.t. pˆ3's initial state q is the join state of pˆ. Add p1, p2, and p3 to P, p to Pc, and mark pˆ1,pˆ2, and pˆ3 as enabled. If pˆ3 is noncircular, add p s (p1 ◦ p2)◦= p3 to R; otherwise add p c (p1 0 p2)◦= p3.

We now elaborate on how we determine the patterns ˆpI , ˆp3, and ˆp.

**Discovering new patterns** *p***ˆ***I* **and** *p***ˆ3** A1 provides an initial pattern pI . For subsequent DFAs, we need to identify which states in Ai+1 = Σ,q 0, Q , F , δ are 'new' relative to Ai = Σ,q0, Q, F, δ. From the PRS definitions, we know that there is a subset of states and transitions in Ai+1 that is isomorphic to Ai:

**Definition 19.** (Existing states and transitions) For every q ∈ Q , we say that q exists in Ai with parallel state q ∈ Q iff there exists a sequence w ∈ Σ∗ such that q = ˆδ(q0, w), q = ˆδ (q0, w), and neither is a sink reject state. Additionally, for every q 1, q 2 ∈ Q with parallel states q1, q2 ∈ Q, we say that (q 1, σ, q 2) ∈ δ exists in Ai iff (q1, σ, q2) ∈ δ. We denote Ai+1's existing states and transitions by QE ⊆ Q and δE ⊆ δ , and the new ones as QN = Q \ QE and δN = δ \ δE.

By construction of PRSs, each state in Ai+1 has at most one parallel state in Ai, which can be found in one simultaneous traversal of the two DFAs.

The new states and transitions form a new pattern instance pˆ in Ai+1, excluding its initial and possibly its exit state. The initial state of pˆ is the existing state q s ∈ QE that has outgoing new transitions. The exit state q X of pˆ is identified by the Exit State Discovery algorithm:

Finally, the new pattern instance is p = Σ,q s, Qp, q X, δp, where Qp = QN ∪ {q s, q X} and δp is the restriction of δN to the states of Qp.

**Discovering the pattern** *p***ˆ (step 2c)** In [27] we show that no two enabled pattern instances in a DFA can share a join state, that if they share any non-edge states, then one is contained in the other, and finally that a pattern's join states is never one of its edge states. This makes finding pˆ straightforward: denoting qj

as the parallel of pˆ3's initial state in Ai, we seek the enabled composite pattern instance (p, q) ∈ Ii for which join(p, q, Ai) = qj . If none is present, we seek the only enabled instance (p, q) ∈ Ii that contains qj as a non-edge state, but is not yet marked as a composite. (Note that if two enabled instances share a non-edge state, then the containing one is already marked as a composite: otherwise we would not have found and enabled the other).

In [27] we define the concept of a minimal generator and prove the following:

**Theorem 1.** Let A1, A2, ...An be a finite sequence of DFAs that has a minimal generator **P**. Then the PRS Inference Algorithm will discover **P**.

#### **5.1 Deviations from the PRS framework**

Given a sequence of DFAs generated by the rules of PRS **P**, the inference algorithm given above will faithfully infer **P**. In practice however, we want to apply the algorithm to a sequence of DFAs extracted from a trained RNN using the L∗ algorithm (as in [26]). Such a sequence may contain noise: artifacts from an imperfectly trained RNN, or from the behavior of L∗ . The major deviations are incorrect pattern creation, simultaneous rule applications, and slow initiation.

Incorrect pattern creation Whether due to inaccuracies in the RNN classification, or as artifacts of the L∗ process, incorrect patterns are often inserted into the DFA sequence. Fortunately, these patterns rarely repeat, and so we can discern between them and 'legitimate' patterns using a voting and threshold scheme.

The vote for each discovered pattern p ∈ P is the number of times it has been inserted as the new pattern between a pair of DFAs Ai, Ai+1 in S. We set a threshold for the minimum vote a pattern needs to be considered valid, and only build rules around the connection of valid patterns onto the join states of other valid patterns. To do this, we modify the flow of the algorithm: before discovering rules, we first filter invalid patterns by splitting step 2 into two phases. Phase 1: Mark all the inserted patterns between each pair of DFAs, and compute their votes. Add to P those whose vote is above the threshold. Phase 2: Consider each DFA pair Ai, Ai+1 in order. If the new pattern in Ai+1 is valid, and its initial state's parallel state in Ai also lies in a valid pattern, then synthesize the rule according to the original algorithm. If a pattern is discovered to be composite, add its composing patterns to P.

As almost every DFA sequence produced by our method has some noise, the voting scheme greatly extended the reach of our algorithm.

Simultaneous rule applications In the theoretical framework, Ai+1 differs from Ai by applying a single PRS rule, and therefore q s and q X are uniquely defined. L∗ however does not guarantee such minimal increments between DFAs. In particular, it may apply multiple PRS rules between two subsequent DFAs, extending Ai with several patterns. To handle this, we expand the initial and exit state discovery methods given above.

1. Mark the new states and transitions QN and δN as before.

If Ai+1's new patterns have no overlap and do not create an ambiguity around join states, then they may be handled independently and in arbitrary order. They are used to discover rules and then enabled, as in the original algorithm.

Simultaneous but dependent rule applications – such as inserting a pattern and then grafting another onto its join state – are more difficult to handle, as it is not always possible to determine which pattern was grafted onto which. However, there is a special case which appeared in several of our experiments (examples L13 ad L14 of Section 7) for which we developed a technique as follows.

Suppose we discover a rule r1 : p0 s (pl ◦ pr)◦= p and p contains a cycle c around some internal state qj . If later another rule inserts a pattern pn at the state qj , we understand that p is in fact a composite pattern, with p = p1 ◦ p2 and join state qj . However, as patterns do not contain cycles at their edge states, c cannot be a part of either p1 or p2. We conclude that the addition of p was in fact a simultaneous application of two rules: r 1 : p0 s (pl ◦ pr)◦= p and r2 : p c (p1 ◦ p2)◦= c, where p is p without the cycle c, and update our PRS and our DFAs' enabled pattern instances accordingly. The case when p is circular and r1 is of rule type (2) is handled similarly.

Slow initiation Ideally, A1 directly supplies an initial rule ⊥ pI to our PRS. In practice, the first few DFAs generated by L∗ have almost random structure. We solve this by leaving discovery of the initial rules to the end of the algorithm, at which point we have a set of 'valid' patterns that we are sure are part of the PRS. From there we examine the last DFA An generated in the sequence, note all the enabled instances (pI , q0) at its initial state, and generate a rule ⊥ pI for each of them. This technique has the weakness that it will not recognise patterns pI that do not also appear as extending patterns p3 elsewhere in the sequence, unless the threshold for patterns is minimal.

### **6 Converting a PRS to a CFG**

We present an algorithm to convert a given PRS to a context free grammar (CFG), making the rules extracted by our algorithm more accessible.

 For simplicity, we restrict the PRS so that every pattern p can only appear on the LHS of rules of type (2) or only on the LHS of rules of type (3) but cannot only appear on the LHS of both types of rules. Similarly, we assume that for each rule ⊥→ pI , the RHS patterns pI are all circular or non-circular. This restriction is natural: all of the examples

in Sections 4.1 and 7.3 conform to it. Still, in [27] we show how to remove this restriction.

We create a CFG G = Σ, N, S, P rod. Σ is the same alphabet of **P** and we take S as a special start symbol. For every pattern p ∈ P, let Gp = Σp, Np, Zp, P rodp be a CFG describing L(p). Let PY ⊆ PC be those composite patterns that appear on the LHS of a rule of type (2). Create the nonterminal CS and for each p ∈ PY , create an additional non-terminal Cp. We set N = {S, CS} {Np} {Cp}.

p∈P p∈PY Let ⊥ pI be a rule in **P**. If pI is non-circular, create a production S ::= ZpI . If pI is circular, create the productions S ::= SC , SC ::= SC SC and SC ::= ZpI . For each rule p s (p1 ◦ p2)◦= p3 create a production Zp ::= Zp1Zp3Zp2 . For each rule p c (p1 ◦ p2)◦= p3 create productions Zp ::= Zp1CpZp2 , Cp ::= CpCp, and Cp ::= Zp3 . Let P rod be the all the productions defined by the above process. We set P rod = { p∈P P rodp} ∪ P rod .

**Theorem 2.** Let G and **P** be as above. Then L(**P**) = L(G).

The proof is given in the extended version of this paper [27].

Expressibility Every RE-Dyck language (Section 2.2) can be expressed by a PRS, but the converse is not true; RE-Dyck languages nest delimiters arbitrarily, while PRS grammars may not. For instance, language L12 of Section 7.3 is not a Dyck language. Meanwhile, not every CFL can be expressed by a PRS [27].

Succinctness The construction above does not necessarily yield a minimal CFG G. For a PRS defining the Dyck language of order 2 – which can be expressed by a CFG with 4 productions and 1 non-terminal – our construction yields a CFG with 10 non-terminals and 12 productions. In this case, and often in others, we can recognise and remove the spurious productions from the generated grammar.

### **7 Experimental results**

#### **7.1 Methodology**

We test the algorithm on several PRS-expressible context free languages, attempting to extract them from trained RNNs using the process outlined in Figure 1. For each language, we create a probabilistic CFG generating it, train an RNN on samples from this grammar, extract a sequence of DFAs from the RNN, and apply our PRS inference algorithm. Finally, we convert the extracted PRS back to a CFG, and compare it to our target CFG.

In all of our experiments, we use a vote-threshold s.t. patterns with less than 2 votes are not used to form any PRS rules (Section 5.1). Using no threshold significantly degraded the results by including too much noise, while higher thresholds often caused us to overlook correct patterns and rules.

#### **7.2 Generating a sequence of DFAs**

We obtain a sequence of DFAs for a given CFG using only positive samples[11,1] by training a language-model RNN (LM-RNN) on these samples and then extracting DFAs from it with the aid of the L∗ algorithm [2], as described in [26]. To apply L∗ we must treat the LM-RNN as a binary classifier. We set an 'acceptance threshold' t and define the RNN's language as the set of sequences s satisfying: 1. the RNN's probability for an end-of-sequence token after s is greater than t, and 2. at no point during s does the RNN pass through a token with probability < t. This is identical to the concept of locally t-truncated support defined in [13].

To create the samples for the RNNs, we write a weighted version of the CFG, in which each non-terminal is given a probability over its rules. We then take N samples from the weighted CFG according to its distribution, split them into train and validation sets, and train an RNN on the train set until the validation loss stops improving. In our experiments, we used N = 10, 000. For our languages, we used very small 2-layer LSTMs: hidden dimension 10 and input dimension 4.

In some cases, especially when all of the patterns in the rules are several tokens long, the extraction of [26] terminates too soon: neither L∗ nor the RNN abstraction consider long sequences, and equivalence is reached between the L∗ hypothesis and the RNN abstraction despite neither being equivalent to the 'true' language of the RNN. In these cases we push the extraction a little further using two methods: first, if the RNN abstraction contains only a single state, we make an arbitrary initial refinement by splitting 10 hidden dimensions, and restart the extraction. If this is also not enough, we sample the RNN according to its distribution, in the hope of finding a counterexample to return to L∗ . The latter approach is not ideal: sampling the RNN may return very long sequences, effectively increasing the next DFA by many rule applications. We place a time limit of 1, 000 seconds (∼ 17 minutes) on the extraction.

#### **7.3 Languages**

We experiment on 15 PRS-expressible languages L1 − L15, grouped into 3 classes:

**Table 1.** Results of experiments on DFAs extracted from RNNs

#### **7.4 Results**

Table 1 shows the results. The 2nd column shows the number of DFAs extracted from the RNN. The 3rd and 4th columns present the number of patterns found by the algorithm before and after applying vote-thresholding to remove noise. The 5th column gives the minimum and maximum votes received by the final patterns (we count only patterns introduced as a new pattern p3 in some Ai+1). The 6th column notes whether the algorithm found a correct CFG, according to our manual inspection. For languages where our algorithm only missed or included 1 or 2 valid/invalid productions, we label it as partially correct.

Alternating Patterns Our algorithm struggled on the languages L3, L6, and L11, which contained patterns whose regular expressions had alternations (such as ab|cd in L3, and ab|c in L6 and L11). Investigating their DFA sequences uncovered the that the L∗ extraction had 'split' the alternating expressions, adding their parts to the DFAs over multiple iterations. For example, in the sequence generated for L3, ef appeared in A7 without gh alongside it. The next DFA corrected this mistake but the inference algorithm could not piece together these two separate steps into a single rule. It will be valuable to expand the algorithm to these cases.

Simultaneous Applications Originally our algorithm failed to accurately generate L13 and L14 due to simultaneous rule applications. However, using the technique described in Section 5.1 we were able to correctly infer these grammars. However, more work is needed to handle simultaneous rule applications in general.

Additionally, sometimes a very large counterexample was returned to L∗ , creating a large increase in the DFAs: the 9thiteration of the extraction on L3 introduced almost 30 new states. The algorithm does not manage to infer anything meaningful from these nested, simultaneous applications.

Missing Rules For the Dyck languages L7−L9, the inference algorithm was mostly successful. However, due to the large number of possible delimiter combinations, some patterns and nesting relations did not appear often enough in the DFA

sequences. As a result, for L8, some productions were missing in the generated grammar. L8 also created one incorrect production due to noise in the sequence (one erroneous pattern was generated two times,passing the threshold).

RNN Noise In L15, the extracted DFAs for some reason always forced that a single character d be included between every pair of delimiters. Our inference algorithm of course maintained this peculiarity. It correctly allowed the allowed optional embedding of "abc" strings. But due to noisy (incorrect) generated DFAs, the patterns generated did not maintain balanced parenthesis.

### **8 Related work**

Training RNNs to recognize Dyck Grammars. Recently there has been a surge of interest in whether RNNs can learn Dyck languages [5,19,21,28]. While these works report very good results on learning the language for sentences of similar distance and depth as the training set, with the exception of [21], they report significantly lower accuracy for out-of-sample sentences.

Among these, Sennhauser and Berwick [19] use LSTMs, and show that in order to keep the error rate within a 5 percent tolerance, the number of hidden units must grow exponentially with the distance or depth of the sequences (though Hewitt et. al. [13] find much lower theoretical bounds). They conclude that LSTMs do not learn rules, but rather statistical approximations. Bernardy [5] experimented with various RNN architectures, finding in particular that the LSTM has more difficulty in predicting closing delimiters in the middle of a sentence than at the end. Based on this, he conjectures that the RNN is using a counting mechanism, but has not truly learnt the Dyck language (its CFG). For the simplified task of predicting only the final closing delimiter of a legal sequence, Skachkova, Trost and Klakow [21] find that LSTMs have nearly perfect accuracy across words with large distances and embedded depth.

Yu, Vu and Kuhn [28] compare the three works above, and note that the task of predicting only the closing bracket of a balanced Dyck word is not sufficient for checking if an RNN has learnt the language, as it can be computed by only a counter. In their experiments, they present a prefix of a Dyck word and train the RNN to predict the next valid closing bracket. They experiment with an LSTM using 4 different models, and show that the generator-attention model [17] performs the best, and is able to generalize quite well at the tagging task . However, they find that it degrades rapidly with out-of-domain tests. They also conclude that RNNs do not really learn the Dyck language. These experimental results are reinforced by the theoretical work in [13], who remark that no finite precision RNN can learn a Dyck language of unbounded depth, and give precise bounds on the memory required to learn a Dyck language of bounded depth.

Despite these findings, our algorithm nevertheless extracts a CFG from a trained RNN, discovering rules based on DFAs synthesized from the RNN using the algorithm in [26]. Because we can use a short sequence of DFAs to extract the rules, and because the first DFAs in the sequence describe Dyck words with

increasing but limited distance and depth, we are often able to extract the CFG perfectly even when the RNN does not generalize well. Moreover, we show that our approach works with more complex types of delimiters, and on Dyck languages with expressions between delimiters.

Extracting DFAs from RNNs. There have been many approaches to extract higher level representations from a neural network (NN), both to facilitate comprehension and to verify correctness. One of the oldest approaches is to extract rules from a NN [24,12]. In particular, several works attempt to extract FSAs from RNNs [18,15,25]. We base our work on [26]. Its ability to generate sequences of DFAs providing increasingly better approximations of the CFL is critical to our method.

There has been less research on extracting a CFG from an RNN. One exception is [23], where they develop a Neural Network Pushdown Automata (NNPDA) framework, a hybrid system augmenting an RNN with external stack memory. They show how to extract a push-down automaton from an NNPDA, however, their technique relies on the PDA-like structure of the inspected architecture. In contrast, we extract CFGs from RNNs without stack augmentation.

Learning CFGs from samples. There is a wide body of work on learning CFGs from samples. An overview is given in [10] and a survey of work for grammatical inference applied to software engineering tasks can be found in [22].

Clark et. al. studies algorithms for learning CFLs given only positive examples [11]. In [7], Clark and Eyraud show how one can learn a subclass of CFLs called CF substitutable languages. There are many languages that can be expressed by a PRS but are not substitutable, such as xnbn. However, there are also substitutable languages that cannot be expressed by a PRS (wxwR - see [27]). In [8], Clark, Eyraud and Habrard present Contextual Binary Feature Grammars. However, it does not include Dyck languages of arbitrary order. None of these techniques deal with noise in the data, essential to learning a language from an RNN.

### **9 Future Directions**

Currently, for each experiment, we train the RNN on that language and then apply the PRS inference algorithm on a single DFA sequence generated from that RNN. Perhaps the most substantial improvement we can make is to extend our technique to learn from multiple DFA sequences. We can train multiple RNNs and generate DFA sequences for each one. We can then run the PRS inference algorithm on each of these sequences, and generate a CFG based upon rules that are found in a significant number of the runs. This would require care to guarantee that the final rules form a cohesive CFG. It would also address the issue that not all rules are expressed in a single DFA sequence, and that some grammars may have rules that are executed only once per word of the language.

Our work generates CFGs for generalized Dyck languages, but it is possible to generalize PRSs to express a greater range of languages. Work will then be needed to extend the PRS inference algorithm.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# Automated and Formal Synthesis of Neural Barrier Certificates for Dynamical Models

Andrea Peruffo1(-) , Daniele Ahmed2 , Alessandro Abate1

1 Department of Computer Science, University of Oxford, Oxford, UK<EMAIL_ADDRESS>2 Amazon Inc, London, UK

Abstract. 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 structured as a neural network, and a sound verifier, which either certifies 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 three orders less), whilst needing a far smaller data set (up to three orders less) for the learning part. Beyond improvements over the state of the art, we further challenge the new approach on a hybrid dynamical model and on larger-dimensional models, and showcase the numerical robustness of our algorithms and codebase.

### 1 Introduction

Barrier Certificates (BC) are an effective and powerful technique to prove safety properties on models of continuous dynamical systems, as well as hybrid models (featuring both continuous and discrete states) [21,22]. Whenever found, 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 (which is in general not a unique object) represents a formal proof of safety for the dynamical model. BC find various applications spanning robotics, multi-agent systems, and biology [7,32].

This work addresses the safety of dynamical systems modelled in general by non-linear ordinary differential equations (ODE), and presents a novel method for the automated and formal synthesis of BC. The approach leverages Satisfiability Modulo Theory (SMT) and inductive reasoning (CEGIS, Figure 1, introduced later), to guarantee the correctness of the automated synthesis procedure: this rules out both algorithmic and numerical errors related to BC synthesis [10].

c The Author(s) 2021

J. F. Groote and K. G. Larsen (Eds.): TACAS 2021, LNCS 12651, pp. 370–388, 2021. https://doi.org/10.1007/978-3-030-72016-2\_20

*Background and Related Work* A few techniques have been developed to synthesise BC. For polynomial models, sum-of-squares (SOS) and semi-definite programming relaxations [14,16,29] convert the BC synthesis problem into constraints expressed as linear or bilinear matrix inequalities: these are numerically solved as a convex optimisation problem, however unsoundly. To increase scalability and to enhance expressiveness, numerous barrier formats have been considered: BC based on exponential conditions are presented in [14]; BC based on Darboux polynomials are outlined in [33]; [30] 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 [20]. Let us remark that SOS approaches are typically *unsound*, namely they rely on iterative and numerical methods to synthesise the BC. [10] a-posteriori verifies SOS candidates via computer-aided design (CAD) techniques [15].

Model *invariants* (namely, regions that provably contain model trajectories, such as *basins of attractions* [28]) can be employed as BC, though their synthesis is less general, as it does not comprise an unsafe set and tacitly presupposes the initial set to be "well placed" within the state space (that is, within the aforementioned basin): [19] introduces a fixpoint algorithm to find algebraic-differential invariants for hybrid models; invariants can be characterised analytically [4] or synthesised computationally [8]. Invariants can be alternatively studied by *Lyapunov theory* [5], which provides *stability* guarantees for dynamical models, and thus can characterise invariants (and barriers) as side products: however this again requires that initial conditions are positioned around stable equilibria, and does not explicitly encompass unsafe sets in the synthesis. Whilst Lyapunov theory is classically approached either analytically (explicit synthesis) or numerically (with unsound techniques), an approach that is relevant for the results of this work looks at automated and sound Lyapunov function synthesis: in [27] Lyapunov functions are soundly found within parametric templates, by constructing a system of linear inequality constraints over unknown coefficients. [23,24,25] employ a counterexample-based approach to synthesise control Lyapunov functions, which inspires this work, using a combination of SMT solvers and convex optimisation engines: however unlike this work, SMT solvers are never used for verification, which is instead handled by solving optimisation problems that are numerically unsound. As argued above, let us emphasise again that the BC synthesis problem, as studied in this work, cannot in general be reduced to a problem of Lyapunov stability analysis, and is indeed more general.

Fig. 1. Schematic representation of the CEGIS loop.

*Core approach* We introduce a method that efficiently exploits machine learning, whilst guaranteeing formal proofs of correctness via SMT. We leverage a CounterExample-Guided Inductive Synthesis (CEGIS) procedure [31], which is structured as an inductive loop between a *Learner* and a *Verifier* (cf. Fig. 1). A learner numerically (and unsoundly) trains a neural network (NN) to fit over a finite set of samples the requirements for a BC, which are expressed through a loss function; then a verifier either formally proves the validity of the BC or provides (a) counter-example(s) through an SMT solver: the counter-examples indicate where the barrier conditions are violated, and are passed back to the learner for further training. This synthesis method for neural BC is formally sound and fully automated, and thanks to its specific new features, is shown to be much faster and to clearly require less data than state-of-the-art results.

*Contributions beyond the State of the Art* Cognate work [34] 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. [34] 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 [34] is shown to be slower (both in the training and in the verification), and to require more data than the CEGIS-based approach of this work, which furthermore introduces numerous bespoke optimisations, as outlined in Section 3: our CEGIS-based technique exploits fast learning, verification simplified by the candidates passed by the Learner, and an enhanced communication between Learner and Verifier. Our approach further showcases numerical robustness and scalability features.

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 in general unsound [26]. Formal methods for Lyapunov synthesis are introduced in [5], together with a counterexample-based approach using polynomial candidates. The work is later extended in [2], which employs NN as candidates over polynomial dynamical models. The generation of control Lyapunov functions using counterexample-based NN is similarly considered in [9], however this is done by means of differing architectural details and does not extend to BC synthesis. Beyond the work in [5], this contribution is not limited to a specific polynomial template, since it supports more general mixtures of polynomial functions obtained through the NN structure, as well as the canonical tanh, sigmoid, ReLU activations (we provide one example of BC using tanh activations). Compared to [5], where we use LP programming to synthesise Lyapunov functions, in this work: a) we use a template-free procedure, thanks to the integration of NNs - these are needed since template-based SOS-programming approaches are not sufficient to provide BCs for several of the presented benchmarks (see Section 4 and [34]); b) we provide an enhanced loss function (naturally absent from [5]), enriched counter-example generation, prioritised check of the verification constraints, and c) we newly synthesise verified barrier certificates for hybrid models, which are

generated using counterexample-based, neural architectures. Finally, beyond [5] the new approach is endowed with numerical robustness features.

SOS programming solutions [14,16,29] are not quite comparable to this work. Foremostly, they are not sound, i.e. do not offer a formal guarantee of numerical and algorithmic correctness. The exception is [10], which verifies SOS candidates a-posteriori by means of CAD [15] techniques that are known not to scale well. Furthermore, they can be hardly embedded within a CEGIS loop - we experimentally show that SOS candidates are handled with difficulty by SMT solvers. Finally, they hardly cope with the experiments we have considered, as already observed in [34]. We instead use SMT solvers (Z3 [11] and dReal [13]) within CEGIS to provide sound outcomes based on NN candidates, proffering a new approach that synthesises and formally verifies candidate BCs altogether, with minimum effort from the user.

*Organisation* The remainder of the paper is organised as follows: Section 2 presents preliminary notions on BCs and outlines the problem. Section 3 describes the approach: training of the NN in Sec. 3.1 and verification in Sec. 3.2. Section 4 presents case studies, Section 5 delineates future work.

### 2 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 by

$$
\dot{x}(t) = \frac{dx}{dt} = f(x), \quad x(0) = x\_0 \in X\_0 \subset X,\tag{1}
$$

where f : X → Rn is a continuous vector field, X ⊆ Rn is an open set defining the state space of the system, and X0 represents the set of initial states. Given model (1) and an unsafe set Xu ⊂ X, the safety verification problem concerns checking whether or not all trajectories of the model originating from X0 reach the unsafe region Xu. BC offer a sufficient condition asserting the safety of the model, namely when no trajectory enters the unsafe region.

Definition 1. *The Lie derivative of a continuously differentiable scalar function* B : X → R*, with respect to a vector field* f*, is defined as follows*

$$\dot{B}(x) = \nabla B(x) \cdot f(x) = \sum\_{i=1}^{n} \frac{\partial B}{\partial x\_i} \frac{dx\_i}{dt} = \sum\_{i=1}^{n} \frac{\partial B}{\partial x\_i} f\_i(x). \tag{2}$$

*Intuitively, this derivative denotes the rate of change of function* B *along the model trajectories.*

Proposition 1 (Barrier Certificate for Safety Verification, [21]). *Let the model in* (1) *and the sets* X*,* X0 *and* Xu *be given. Suppose there exists a function* B : X → R *that is differentiable with respect to its argument and satisfies the following conditions:*

$$B(x) \le 0 \,\,\forall x \in X\_0, \quad B(x) > 0 \,\,\forall x \in X\_u, \quad \dot{B}(x) < 0 \,\,\forall x \in X \,\,\,s.t. \,\,B(x) = 0,\tag{3}$$

*then the safety of the model is guaranteed. That is, there is no trajectory of the model contained in* X*, starting from the initial set* X0*, that ever enters set* Xu*.*

Consider a trajectory x(t) starting in x0 ∈ X0 and the evolution of B(x(t)) along this trajectory. Whilst the first of the three conditions guarantees that B(x0) ≤ 0, the last condition asserts that the value of B(x(t)) along a trajectory x(t) must decrease. Hence such a trajectory x(t) cannot enter the set Xu, where B(x) > 0 (second condition), thus ensuring the safety of the model.

### 3 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) [31], 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 network is then translated into a logical formula in an appropriate theory, by evaluating it with symbolic inputs, instead of canonical floating point numbers. The details of this conversion are outlined in [2]. This encoded 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 [15], namely an algorithmic decision procedure that extends Boolean SAT problems to richer, more expressive theories, such as non-linear arithmetics.

More precisely, the learner trains a NN composed of n 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, as in this work) 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 counterexamples c. Note that the NN retains its old weights, and restarts the training from the weights obtained at the end of the previous session. 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.

The performance of the CEGIS algorithm in practice hinges on the effective exchange of information between the learner and the verifier [3]. 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 that ensures robustness and outputs a function that is tailored to the verification engine. Over consecutive loops, the verifier may return similar counter-examples: we thus propose a more informative counter-examples 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, efficient, and robust CEGIS loop, which is shown in this work to clearly outperform state-of-the-art methods.

#### 3.1 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* and non-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. The order γ of the polynomial activations is a hyper-parameter fed at the start of the procedure: we split the i-th hidden layer into γ portions and apply polynomial activations of order j to the neurons of the j-th portion.

*Example 1 (Polynomial Activations).* Assume a NN composed of an input x, 3 hidden neurons and 1 activation-free output, with γ-th order polynomial activation, γ = 3. We split the hidden layer in γ sub-vectors, each containing one neuron. The hidden layer after the activation results in

$$z = \begin{bmatrix} W\_1^{(1)}x + b\_1 & \ & (W\_2^{(1)}x + b\_2)^2 & \ & (W\_3^{(1)}x + b\_3)^3 \end{bmatrix}^T, \ q$$

where the W(1) i are the i-th row of the first-layer weight matrix, and the bi 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 ∈ X0, B(x) > 0 for x ∈ Xu, and a negative Lie derivative B˙ (Eq. (2)) over the set implicitly defined by B(x)=0. The training minimises a loss comprising three terms, namely

$$L = L\_0 + L\_u + L\_d = \frac{1}{N} \sum\_{i=1}^{N} \left( \max\_{s\_i \in X\_0} \{ \tau\_0, B(s\_i) \} + \max\_{s\_i \in X\_u} \{ \tau\_u, -B(s\_i) \} \right)$$

$$+ \max\_{s\_i: B(s\_i) = 0} \{ \tau\_d, \dot{B}(s\_i) \}, \quad (4)$$

where si, 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 can be a set with small volume, thus it is highly unlikely that a single sample s will satisfy B(s)=0. We thus 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 suggests a few ameliorations. Terms L0 and Lu 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 Xu or the negativity over X0. As such, the training often returns a candidate B(x) with values just below τ0 in X0 or above τu in Xu. 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

$$\text{LeakyReLU}(\alpha, x) = \text{ReLU}(x) - \alpha \,\text{ReLU}(-x),\tag{5}$$

where α is a small positive constant, we rewrite term L0 as

$$L\_0 = \frac{1}{N} \sum\_{s\_i \in X\_0} \text{ReLU}(B(s\_i) - \tau\_0) - \alpha \cdot \text{satReLU}(-B(s\_i) + \tau\_0), \tag{6}$$

where satReLU is the saturated ReLU function3. The term Lu is similarly modified. The composite loss function works as follows. Incorrect samples account for the main contribution to the loss function, leading the NN to correct those first via the ReLU term in Eq. (6). 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 Ld: 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 L0 and Lu, 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 B˙(s), whilst not penalising values B˙(s) ≥ 0. The training is driven to include more samples in B, guiding towards a negative B˙(s), and finally enhancing learning. The expression of Ld results in

$$L\_d = -\frac{1}{N} \sum\_{s \in \mathcal{B}} \text{satReLU}(-\dot{B}(s) + \tau\_d). \tag{7}$$

3 Let us define M to be an arbitrary upper bound, then satReLU(x) = min(max(0, x), M).

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 increases 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. We have observed that the use of linear activations, particularly in the output layer, positively affects the training: they provide robustness that is needed to the synthesis of BC (see Experimental results), without increasing the order of the network with new polynomial terms.

*Learning in Separate Batches* The structure of the conditions in (3) and the learning loss in (4) naturally suggests a separate approach to training. We then split the dataset S into three batches S0, Su and Sx, each including samples belonging to X0, Xu and X, respectively. For training, we compute the loss function in a parallel fashion. Similarly, for the verifier, generated counterexamples are added to the relevant batch.

#### 3.2 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

$$(x \in X\_0 \land B(x) > 0) \lor (x \in X\_u \land B(x) \le 0) \lor (B(x) = 0 \land \dot{B}(x) \ge 0). \tag{8}$$

The verifier searches for solutions of the constraints in Eq. (8), which in general requires manipulating non-convex functions. This can be cumbersome and timeconsuming, hence simple expressions of B can enhance the verification procedure. 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 [13] and Z3 [11]. dReal is a δ-complete solver, namely the unsat decision is correct [12],

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 acceptable to return spurious counter-examples: indeed, 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 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 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 ∧ B˙(x) ≥ 0) by the SMT solver, since the formula contains the high-order expression B˙(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 ∈ X0 ∧ B(x) > 0) and ηu = (x ∈ Xu ∧ 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 ∧B˙(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 sample x¯ 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 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 can be in general expensive. Specifically, the verifier solves Eq. (8) to find a first counter-example x¯; 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 in increasing the information generated by the verifier by finding a sample that maximises the violation of the BC conditions. To this end, firstly 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 counter-example are also likely to invalidate the BC conditions. Secondly, for the original counter-example, we compute the gradient of B (or of B˙) and follow the direction that maximises the

violation of the BC constraints. As such, we follow the B (resp. B˙) maximisation when considering x ∈ X0 (x s.t. |B(x)| < τv), and vice versa when x ∈ Xu. This gradient computation is extremely fast as it exploits the neural architecture, and it provides more informative samples for further use by the learner.

### 4 Case Studies and Experimental Results

All experiments are performed on a laptop workstation with 8 GB RAM, running on Ubuntu 18.04. We demonstrate that the proposed method finds provably correct BCs on benchmarks from literature comprising both polynomial and non-polynomial dynamics: we compare our approach against the work [34], as this is the only work on sound synthesis of BCs with NNs to the best of our knowledge, and against the SOS optimisation software SOSTOOLS [18]. Beyond the benchmarks proposed in [34], we newly tackle a hybrid model as well as larger, (up to) 8-dimensional models, which push the boundaries of the verification engine and display a significant extension to the state of the art. To confirm the flexibility of our architecture, we employ SMT-solver dReal in the first four benchmarks, whereas we study the last four using Z3. In all the examples, we use a learning rate of 0.1 for the NN and the loss function 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 = ∞. Accordingly, the training over a large set B results in a candidate B with a negative derivative over this large region, which validity is more likely to be certified by the verifier. We set a verification parameter τv = 0.05 (cf. Sec. 3.2), a timeout (later denoted as OOT) of 60 seconds and the precision for dReal to δ = 10−6. Table 1 summarises the outcomes. We emphasise that our approach supports any network depth and width. The presented results seek a tradeoff between speed (low order, small networks) and expressiveness (high order, larger networks): a different architecture may result in a slower or faster synthesis.

For the first four benchmarks, we compare our procedure, denoted as CEGIS, with the repeated results from [34], which however does not handle the hybrid model in the fifth benchmark. We have run the algorithm in [34] and reported the cumulative synthesis time under the 'Learn' column. However the verification is not included in the repeatability package, hence we report the results from [34], which are generated with much more powerful hardware. Due to this issue of lack of repeatability, we have not run [34] on the larger models. Compared to [34], the outcomes suggest that we obtain *much faster* synthesis and verification times, whilst requiring up to only 0.1% (see Obstacle Avoidance Problem) of the training data: [34] 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. Evidently this gap in performance derives from the different synthesis procedure: it appears to be more advantageous to employ a smaller, randomly sampled initial dataset that is progressively augmented with counter-examples, rather than to uniformly sample the state space to then train the neural network.

Next, we have implemented the SOS optimisation problems in [10] within the software SOSTOOLS [18] to generate barrier candidates, which are polynomials up to order 4 (this is the maximum order of the polynomial candidates generated by our Learner). In a few instances we ought to conservatively approximate the expression of X0 or Xu in order to encode them as SOS program - this makes their applicability less general. SOSTOOLS has successfully found BC candidates for five of the eight benchmarks, and they were generated consistently fast, in view of the convex structure of the underlying optimisation problem. However, recall that these techniques lack soundness (also due to numerical errors), which is instead a core asset of our approach. Consequently, we have passed them to the Z3 SMT solver, which should easily handle polynomial formulae: only one of them ('Hybrid Model') has been successfully verified; instead, the candidate for the 'Polynomial Model' has been invalidated (namely Z3 has found a counterexample for it), whereas the verification of the remaining BC candidates has run out of time. For the latter instances, we have experienced that SOSTOOLS generally returns numerically ill-conditioned expressions, namely candidates with coefficients of rather different magnitude, with many decimal digits: even after rounding, expressions with this structure are known to be hardly handled by SMT solvers [2,5], which results in long time needed to return an answer - this explains the experienced timeouts. These experiments suggest that the use of SOS programs within a CEGIS loop appears hardly attainable.

Notice that all the case studies are solved with a *small number* of iterations (up to 9) of the CEGIS loop: this feature, along with the limited runtimes, is promising towards tackling synthesis problems over larger models.

For the eight case studies, we report below the full expressions of the dynamics of the models, the spatial domain X (as a set of constrains), the set of initial conditions X0 ⊂ X, and the unsafe set Xu ⊂ X. We add a detailed analysis of the CEGIS iterations involved in the synthesis of the corresponding BCs.

Table 1. Outcomes of the case studies: Cumulative time for Learning and Verification steps are given in seconds; 'Samples' indicates the size of input data for the Learner (in thousands); 'Iters' is the number of iterations of the CEGIS loop (which is specific to our work); × indicates a synthesis or verification failure; OOT denotes a verification timeout. The Hybrid and the three ODE Models are newly introduced in this work.

*Darboux Model* This 2-dimensional model is approached using polynomial BCs. Its analytical expression is

$$\begin{cases} \dot{x} = y + 2xy, & X = \{-2 \le x, y \le 2\}, \\ \dot{y} = -x + 2x^2 - y^2, & \text{with domains} \end{cases} \quad \text{with domains} \quad \begin{aligned} X &= \{-2 \le x, y \le 2\}, \\ X\_0 &= \{0 \le x \le 1, 1 \le y \le 2\}, \\ X\_u &= \{x + y^2 \le 0\}. \end{aligned}$$

The work [33] reports that methods based on linear matrix inequalities fail to verify this model using polynomial templates of degree 6. Our approach generates the BC shown in Fig. 2 (left) in approximately 30 seconds, roughly half as much as in [34], and using only 500 initial samples vs more than 65000. The initial and unsafe sets are depicted in green and red, respectively, whereas the level set B(x)=0 is outlined in black. The BC is derived from a single-layer architecture of 10 nodes, with linear activations.

*Exponential Model* This model from [17] shows that our approach extends to non-polynomial systems encompassing exponential and trigonometric functions:

$$\begin{cases} \dot{x} = e^{-x} + y - 1, & X = \{-2 \le x, y \le 2\}, \\ \dot{y} = -\sin^{2}x, & X\_{0} = \{(x + 0.5)^{2} + (y - 0.5)^{2} \le 0.16\}, \\ & X\_{u} = \{(x - 0.7)^{2} + (y + 0.7)^{2} \le 0.09\}. \end{cases}$$