Dynamical systems model the time evolution of both natural and engineered processes. The automatic analysis of such models relies on diferent techniques ranging from reachability analysis, model checking, theorem proving, and abstractions. In this context, invariants are subsets of the state space containing all the states reachable from themself. The verification and synthesis of invariants is still a challenging problem over many classes of dynamical systems, since it involves the analysis of an infinite time horizon. In this paper we propose a method for computing invariants through sets of trajectories propagation. The method has been implemented and tested in the tool Sapo which provides reachability methods over discrete time polynomial dynamical systems.
Invariants are subsets of the state space of dynamical systems that contain all the states
reachable from themself. This notion is crucial in may fields such as control theory [
Many studies investigating the verification and synthesis of invariants have been published
in the literature so far. Some of them deals with polynomial systems and investigated
semialgebraic invariants [
This work focuses on discrete-time transition systems. The applications of this kind of
systems go from biology [
Other approaches involving -inductiveness have been proposed in the literature (see, e.g.,
[
When all operations are computed without approximations the proposed method converges if and only if reachability converges within a finite number of steps. Surprisingly, when reachability, unions and intersections are over-approximated, convergence is achieved on a larger class of systems.
The paper is organized as follows. Section 2 introduces the problem, the notation, and the basic tools. Section 3 incrementally presents our proposal for a set-based algorithm deciding -inductiveness. In Section 4 we detail the approach over polynomial systems as implemented in Sapo and discuss an example. Section 5 draws concluding remarks.
In this section, we present some central notions for this work.
Definition 1 (Transition System). A transition system is a pair ⟨, ⟩, where is the set of states and ⊆ × is the transition relation. An initial set for a transition system TS = ⟨, ⟩ is a set of states ⊆ .
forward image of itself, i.e. () =def { ∈ | ∃′ ∈ (′, ) ∈ }.
For any function : 2 → 2 mapping a subset of the states into another subset of the states, we define the image of restricted to as () =def () ∩ . We write () to denote consecutive applications of to , i.e., () =def {︃ (︀ − 1())︀ if = 0 if > 0 .
So, (), () and () stand for the forward image of restricted to , consecutive applications of the forward image, and consecutive applications of the forward image restricted to , respectively.
A set under-approximates a set ′ when ⊆ ′. A set over-approximates a set ′ when ⊇ ′. A function ′ : → under-approximates (over-approximates) a function ′′ when ′() ⊆ ′′() ( ′() ⊇ ′′(), respectively) for all ∈ . It is easy to see that, if over/under-approximates , then , , and do the same for , , and , respectively. Definition 2 (Reachability). Given a transition system TS and an initial set , the set of states reachable from in -steps is the set () and the set of states reachable from is ⋃︀∈N ().
A property for a system is a set of states of the system. An invariant for a system and an initial set of states, , is a property that contains all the states reachable from . is a set ⊆ such that Definition 3 (Property and Invariant).
Let TS = ⟨, ⟩ be a transition system and let ⊆ be an initial set of states for TS. A property of TS is a set of states ⊆ . An invariant for TS and ⋃︁ () ⊆ . ∈N (1)
A well-know problem in the context of verification is that of checking whether a given
property is an invariant for a system (see, e.g., [
One could consider the problem of constructing the smallest invariant for a system. The problem is well defined since if both 1 and 2 are invariants for TS and , then also 1 ∩ 2 is an invariant. In other terms, a given property is an invariant for the system if and only if it includes the smallest invariant. However, as one could imagine finding the smallest invariant could be far more dificult than checking whether a given property is an invariant.
In this paper, we will describe a method for checking whether a property is an invariant,
while trying to construct the smallest invariant. The method has to work also on infinite state
systems assuming that only a restricted number of operations are computable. For instance,
we will be able to over-approximate the forward image computation, but not the backward
one. This is a reasonable compromise in term of computability and eficiency in the case of
non-linear systems [
A standard approach for verifying invariants is the use of k-induction [
, then is not an invariant and we have a counter-example of length . If this is not the case, then we try to prove that is -inductive, which implies that it is an invariant.
We now briefly explain the meaning of -inductiveness. By Equation 1 and by induction principle over the naturals, is an invariant if and only if 0() ⊆
() ⊆ → +1() ⊆ (basis) (step) Induction and -induction are equivalent over N. Thus, is an invariant if and only if ⋀︀− 1 () ⊆
=0 ︁( ⋀︀− 1 + () ⊆
︁)
=0
→ +() ⊆
(-basis)
(-step)
Unfortunately, the inductive -step has to be proven for a generic ∈ N, which makes the
technique not algorithmic for infinite systems. A stronger approach, named -inductiveness,
can be used to automatize invariant verification. Intuitively, -inductiveness strengthens -steps
by requiring that all the paths laying inside for −
1 steps stay inside also after steps,
no matter whether they started from or not. Considering all the possible paths starting from
and not only those from allows to drop the dependence on . The following definition is a
set-based version of the definition given in [
Definition 4 ( -Inductive Property). Given a transition system TS and an initial set , a inductive property for TS and is a set of states ⊆ such that ⋀︀− 1 ( ) ⊆
=0 (− 1( )) ⊆ (-init) (-ind) It is well known that -inductiveness implies invariance.
Lemma 2.1. Let TS be a transition system and be and initial set for TS. If is -inductive for TS and , then is an invariant for TS and .
The reader should be aware that this is a specialized use of induction: a property could be an invariant and not -inductive. Moreover, a set could be + 1-inductive and not -inductive. Example 1. Let TS be the transition system ⟨, ⟩, with = {0, 1, 2, 3} and = {(0, 0), (1, 3), (3, 2), (2, 2)}, and let be the set {0}. The set {0, 1} is an invariant for TS and , but it is not -inductive for any . Moreover, = {0, 1, 2} is 2-inductive because 1 ( ) = ( ) ∩ = {0, 2, 3} ∩ {0, 1, 2} = {0, 2} and ({0, 2}) = {0, 2} ⊆ . However, is not 1-inductive. As a matter of the fact, 0 ( ) = and ( ) = {0, 2, 3} ̸⊆ .
If -init does not hold, not only is not -inductive, but it is not an invariant. If both -init and -ind succeed, is -inductive and it is an invariant. If -init succeeds, but -ind does not, then is not -inductive, but it could be + 1-inductive, i.e., has to be incremented. This forces Algorithm 1 to return True if and only if there exists such that is -inductive. Using two auxiliary variables, we can implement Algorithm 1 to require only 2 forward images, 1 intersection, and 2 subset tests at each iteration. Unfortunately, Algorithm 1 is not complete and, if is not -inductive for any , but it is an invariant, it does not terminate. Algorithm 1 K-Inductiveness 1: function -inductive(TS = ⟨, ⟩, , ) 2: if ̸⊆ then 3: return False 4: ← 1 5: while True do 6: if ( ) ̸⊆ then 7: return False 8: if (− 1( )) ⊆ then 9: return True 10: + 1
←
In formal verification of infinite state systems and more in general in the context of symbolical
representations of transition systems, traces representing punctual evolutions are not always
feasible because of the infinite domain. Processing sets of traces, potentially infinite in number,
even though involving a finite set of states, is certainly computationally more expensive than
computing the evolution of sets consisting in an infinite number of states. Moreover, it is usually
not reasonable to work on single points/traces when we do not have knowledge of the system
at infinite precision level [
In particular, in the case of infinite states systems and reachability computation a generic set-based technique depends on the choice of: • A domain ⊆ 2 providing finite representations for the set of points to be manipulated.
The most used domains are ellipsoids, polytopes, zonotopes. • Methods for computing/approximating basic set operations and tests over , e.g., union, intersection, and inclusion. • A family of admissible transition relations. For instance, linear or polynomial systems could be considered. • For all the transition functions ∈ , a method : → that computes/approximates the image of any set ∈ through .
Such tools are suficient to obtain an implementation of Algorithm 1 instantiated on a specific domain. In particular, as far as basic set operations are concerned, it only requires the computation/approximation of intersections and the tests of inclusion. If both forward images and intersections are over-approximated and inclusion tests never provide false positive answers, then the procedure returns True if and only if is -inductive for some .
It is worth to notice that computing backward images, which could be problematic in the
case of non-linear transition relations, is not required. This is a crucial diference with respect
to algorithms such as IC3/PDR [
Algorithm 1 has two main flaws: 1) if is an invariant, but it is not -inductive for any , Algorithm 1 does not terminate; 2) the algorithm is not able to refine and discover possible -inductive invariant included in . Section 3 tries to overcome both issues by presenting algorithms that, along the computation, search for the “smallest” invariant included in that is also -inductive for some .
We propose a set-based algorithm that takes in input a system TS, a set of initial states and a property and aims at constructing a -inductive property CI, candidate inductive, included in . Intuitively, if TS reaches a state outside in transition steps from , then is not an invariant and the algorithm returns False. If this is not the case, then the algorithm considers the candidate invariant CI built so far. CI has two main properties: it is included in and it includes all the states reachable from in − 1 steps. If CI is -inductive, then is an invariant and the algorithm returns True. Otherwise, CI is augmented with the states reachable in -steps and is incremented. If is an invariant, but it does not includes any -inductive set, then the algorithm does not terminate. These ideas are formalized in Algorithm 2 which describes sets by using a finite representation domain and the method to compute/approximate the images of as discussed in Section 2. In particular, while maps any subsets of in another subset of according to , : → with ⊆ 2 .
Algorithm 2 Set–based K–inductiveness 1: function candidate-inductive(TS = ⟨, ⟩, , ) 2: if ̸⊆ then 3: return False 4: ← 1 5: CI ← 6: while True do 7: if () ̸⊆ then 8: return False 9: if (︀ C−I 1(CI))︀ ⊆ CI then 10: return True 11: CI ← CI ∪ () 12: ← + 1
It is easy to see that the variable CI in Algorithm 2 accumulates the reachability set of TS and as the computation proceeds. In particular, if under/over-approximates , then CI under/over-approximates ⋃︀− 1 (). This observation allows us to prove some properties =0 about the results of Algorithm 2.
Theorem 3.1 (Correctness of Algorithm 2). If under-approximates and Algorithm 2 returns False, then is not an invariant for TS and . If over-approximates and Algorithm 2 returns True after while-loop iterations, then CI is -inductive for TS and and is an invariant for them. Furthermore, if () = () for all ∈ , Algorithm 2 returns False if and only if is not an invariant for TS and and, when it returns True, any ⊂ CI is not an invariant for TS and .
Let us underline that both the termination of Algorithm 2 and the meaning of its outcomes depend on the approximation of done by . Figure 1 depicts two examples in which Algorithm 2 terminates, but does not return the aimed answer.
When instead exactly represents , Algorithm 2 does not terminate only when is an invariant of the system by Theorem 3.1. Moreover, when the algorithm loops forever, the set of points reachable from does not converge after a finite number of while-loop iterations. Indeed, if ⋃︀∈N () = ⋃︀≤ ℎ () for some ℎ, then CI = ⋃︀∈N () after ℎ iterations, CI is ℎ-inductive, and the algorithm returns True. Hence, if we are assuming that all set-operations are implemented without approximations and exactly represents , Algorithm 2 terminates if and only if either is not an invariant or the set of reachable states can be computed with a ifnite number of forward transitions.
Corollary 3.1.1 (Termination). If () exactly represents () for any ∈ , then Algorithm 2 terminates if and only if either is not an invariant or there exists ℎ such that ⋃︀∈N () = ⋃︀≤ ℎ (). 2()
(()) 2() (a) Algorithm 2 returns True even though is not an invariant for TS and when under-approximates . However, () is not a subset of and, as a consequence, is not an invariant for TS and .
(b) Algorithm 2 returns False even though is an invariant for TS and when overapproximates . However, () is a subset of and 3() = (). Thus, is an invariant for TS and .
While the number of while-loop iterations of Algorithm 2 cannot be upper-bounded in the general case, we can provide a cost depending on the number of iterations for it. Each repetition of the while-loop in Algorithm 2 requires Θ() forward image computations, 2 inclusion tests, Θ() intersection and 1 union. Such operations are denoted as set-operations, since they transform sets of states.
Theorem 3.2 (Set Complexity). If Algorithm 2 terminates exactly after iterations of the while-loop, then it requires Θ(2) set operations.
Without increasing the asymptotic complexity of the algorithm we could modify line 9 by considering the test
( C−I 1(CI)) ⊆ CI ∨ ( − 1( )) ⊆ , i.e., at each step we also check whether is -inductive. This would guarantee termination also in the case in which is -inductive. However, we are interested in classes of systems where it is not reasonable to assume that can be represented/approximated in the domain of sets, for instance, and can be a half-space and the set of all the ellipsoids in the same space, respectively. This was one of the reasons for which we moved from Algorithm 1 to Algorithm 2.
Theorem 3.2 provides a complexity bound for Algorithm 2 in terms of set operations. However, the cost of each operation actually depends on the complexity of the involved sets. For instance, if and are the union of 1 and 2 basic sets, respectively, e.g., = ⋃︀=1 1 and = ⋃︀2
=1 , the evaluation of ∩ may consists in up to 1 * 2 basic sets: the non-empty intersections between each and each . In the context of Algorithm 2, if we assume that both and its forward image are basic sets, after iterations of the algorithm while-loop, CI consists in the list [ 0(), . . . , ()]. Thus, the evaluation of C−I 1(CI) at line 9 may produce a set consisting of up to basic sets and, as a consequence, requires () basic set operations.
The complexity of Algorithm 2 in terms of basic set operations pushes us to either avoid the test of line 9 or replace the exact union performed at line 11 with an approximate version of it that reduces the number of basic sets in CI.
In order to overcome the limits pointed out by both Corollary 3.1.1 and other complexity issues discussed at the end of the previous section, we consider a diferent approach that aims to reduce the number of basic sets in CI. This is done by replacing the line 11 in Algorithm 2 with the assignment:
CI ←
Join(CI, ( )),
where the operator Join over-approximates the set union and returns one single basic set and
is a new variable initially set to . Please notice that we are not forcing Join to be a specific
function, but instead we are requiring it to return one single basic set over-approximating the
union of its parameters. For instance, depending on the adopted set representation, Join may
be implemented by the convex-hull function, the bounding box function, or the Halbwachs
widening operator [
Due to the proposed change, CI may be not included in even when ( ) ⊆ for all ∈ N (e.g., see Figure 2). When this happens we can either return Unknown, since we do not have a counter-example, but our -inductive set is growing outside , or re-initialize , CI, and . Algorithm 3 implements this approach.
Algorithm 3 Set-based Join K-inductiveness
Figure 2 depicts a situation in which and CI have to be re-initialized, i.e., CI ̸⊆ at line 13 of Algorithm 3. In this example the domain consists in rectangles in R2 and the Join overapproximates the union of two rectangles with the smallest rectangle including both of them, i.e., their bounding box. The violation of is clearly an artifact of the over-approximation. Hence, the algorithm can safely re-initialize and CI.
As done for Algorithm 2, we aim to investigate the correctness of Algorithm 3 when is an under-approximation, an over-approximation, and the exact representation of the forward image function . In order to achieve this goal, we must notice two main properties of the algorithm. First of all, CI is a subset of at the beginning of each while-loop iteration. Moreover, if under-approximates (over-approximates) the forward image function , then, after ℎ iterations of the while-loop, ( ) under-approximates (over-approximates) ℎ(). These properties allow us to prove the following correctness result.
Theorem 3.3 (Correctness of Algorithm 3). If under-approximates and Algorithm 3 returns False, then is not an invariant for TS and . If over-approximates and Algorithm 3 returns True, then is an invariant for TS and . Let () equal () for any ∈ . If Algorithm 3 returns True, then is an invariant for TS and . Moreover, Algorithm 3 returns False if and only if is not an invariant for TS and .
The set complexity of Algorithm 3 equals that of Algorithm 2 where the cost of exact unions is replaced by that of the Join operations. Nevertheless, the number of basic sets represented in CI and the costs of line 10 depend on Join. If Join always returns one single basic set, then the evaluation of the condition at line 10 requires to compute forward images of a single basic set and − 1 intersections and 1 inclusion test between two basic sets. Thus, in this case, the asymptotic complexity of the test at line 10, whose costs in Algorithm 2 was (), becomes (). In Section 4, we present two possible semantics for the Join operator and one of them achieves this complexity.
As far as termination is concerned, a lot depend once more on the choice of Join. Algorithm 3 can loop forever only when is an invariant as Algorithm 2 did too. Intuitively, if there exists a -inductive set which includes the forward images of ℎ() for a given ℎ ≥ 0 and such set can be obtained through successive Join operations, then the algorithm terminates. However, it is no more necessary that the forward images converge within a finite number of steps.
The analysis of Algorithm 3 in the previous section relies on the hypothesis that intersections, needed to compute C−I 1(CI) at line 10, are computed without approximations. However, this is not always feasible and depends on the choice of the domain . For instance, let us consider the domain of ellipsoids; the intersection of two incomparable ellipsoids is not an ellipsoid. So, we now discuss the properties of Algorithm 3 under the hypothesis that the intersections between sets in are approximated.
As already observed, the intersections are exclusively used to evaluate the condition at line 10 and do not afect the correctness of the result when False is returned. Because of this, we focus on the cases in which Algorithm 3 returns True which, by Theorem 3.3, implies is an invariant for TS and only when is an over-approximation of . It is easy to see that if we over-approximate ∩, the algorithm still has no false positive outcomes. Intuitively, if CI contains the over-approximation of the reachable sets it also contains the exact ones, and if the over-approximation of the states reachable from CI after − 1 steps inside it are still inside CI, then this also holds for the exact reachable sets.
Corollary 3.3.1 (Over-approximation of and ∩). Let () and () over-approximate () and () ∩ , respectively, for any , ∈ . If Algorithm 3 returns True, then is an invariant for TS and .
We now consider transition systems defined through discrete time polynomial dynamical systems, i.e., dynamical systems described by equations of the form x+1 = f (x), where x ∈ R and f : R → R is a vector of -variate polynomial functions. The transition system TS = ⟨, ⟩ associated to a dynamical system of this form has = R and = {(x, x+1) | x, x+1 ∈ R, x+1 = f (x)}.
Sapo [
Since Sapo can exclusively over-approximate forward images, we can only aim to validate candidate invariants as proved in Theorem 3.3. Indeed, nothing can be deduced when Algorithm 3 returns False as highlighted by Figure 1b which depicts a scenario in which the algorithm returns False and is an invariant.
Sapo natively represents exact unions of polytopes as list of polytopes. Moreover, it supports exact intersections and inclusion test of exact union polytopes. Thus, in order to implement Algorithm 3 we need to: 1. choose a specification for candidate invariants to be tested. Sapo already supports polytopes, but we aimed to deal with a broader class of properties;
3. define a function Join to over-approximate the union of polytopes.
As far as the candidate invariant specification is concerned, we decided to admit convex linear sets, i.e., sets corresponding to the solutions of a linear system. This specification includes polytopes, but it also covers an infinite number of non-bounded sets.
Example 2. A possible candidate invariant for the Sapo implementation of Algorithm 3 is the set of unbounded solutions of the linear system : ︂{ + 2 * ≤ 5 .
3 * ≥ + 1
We implemented the test ̸⊆ , where is a generic polytope and is the candidate invariant, by using linear programming. Any polytope is specified as a system of linear inequalities, ℒ. Testing ̸⊆ is equivalent to test whether there exists an inequality (x) ≤ in the specification of such that ℒ enriched with (x) > has a feasible solution. As a matter of the fact, such a solution does exist if and only if there exists a state, i.e., the solution itself, that belongs to , because it satisfies ℒ, and does not belong to , because it satisfies (x) > that is the negation of one of the constraints of . Once ̸⊆ has been implemented, ⊆ can be trivially computed as ¬( ̸⊆ ).
problems where is the number of linear inequalities in ℒ . If the problems are solved by
using Karmarkar’s algorithm [
We consider two alternative semantics for the function Join: Listing and Packaging. The most naïve semantics for Join is the exact union, i.e., Join(, ) = ∪ . Since exact unions are usually represented as lists of basic sets, we refer to this approach as Listing. An alternative semantics for the function Join(, ) consists in finding the smallest polytope, having the same directions of , that includes both and . We call this approach Packaging because it packs many polytopes into a new one. Figures 2 and 3 depicts CI during an evaluation of Algorithm 3 when the Join implements the Packaging semantic.
Packaging usually produces coarser union approximations than Listing or, for instance, convex hull. However, it preserves the number of directions needed to represent sets and, as a consequence, it guarantees that both the memory required to represent a single set and the cost of set functions do not change during the computation. Neither Listing nor convex hull ensure these features: the former increases the number of polytopes per set as the computation proceeds, the latter produces only one polytopes per set, but tends to increase the number of directions required to represent it.
While Packaging is usually more efective than Listing, there are cases in which the former leads to an infinite computation, while the latter produces an answer. Example 3. Consider the following system If the system initial state is (, ), this system keeps flipping and values at each evolution step. For instance, if the initial set is = [2.9, 3.1] × [0.9, 1.1], after one step the system reaches () = [0.9, 1.1] × [2.9, 3.1], and one further step brings back the system to its initial condition since 2() = . So, if we select as candidate invariant + ≤ 5, the packaging of any pair of consecutive forward images violates it and Algorithm 3 resets CI and to the last one of the two. Despite this, the packaging of the new pair of consecutive forward images still does not satisfy + ≤ 5 and the algorithm loops forever (see Figure 3a).
Intriguingly, Algorithm 3 with Listing can easily prove that + ≤ 5 is an actual invariant for the system and = [2.9, 3.1] × [0.9, 1.1]. The reachability computation converges after 1 step and all reachable points satisfies the candidate invariant, hence, in this scenario, Listing is preferable to Packaging.
In Example 3, Listing and Packaging produce diferent outcomes because the packaging of two polytopes and may not satisfy a property even when and individually do it (see Figure 3a). We can mitigate this issue by adding the property directions to the packaging and minimizing and maximizing them on and . This approach is called enhanced Packaging. Since we focus on convex linear properties, and satisfy a property if and only if the enhanced packaging of and satisfies too. Indeed, the packaging does not satisfy if and only if there exists an inequality (x) ≤ in the specification of that is not satisfied. However, the enhanced packaging of and is also defined by a set of inequalities among which there is (x) ≤ max(, ) where and are the maxima of (x) in and , respectively. Thus, if the packaging does not satisfy (x) ≤ either or do not satisfy it too. The enhanced Packaging may increase the complexity of the packaging representation with respect to and . However, this only happens when the property directions are not yet present in the representation of the two sets. Figure 3b shows enhanced Packaging on Example 3.
The algorithms presented in Sections 3 and 4 have been implemented in Sapo [
Let [0.09, 0.11] × [0.09, 0.11] be the system initial set . Diferent choices of and can produce very diferent behaviours. We focus on the behaviour obtained with ∈ [0.19, 0.2] and ∈ [0.29, 0.3] and we consider the candidate invariant 1 : ≤ 0.3. Sapo proved that 1 is an invariant for the system by using Listing in less than 1 second on a MacBook Pro 2020 with 16GB of RAM. The proof ended after 12 iterations of the Algorithm 3 while-loop. Packaging can achieve the same result after only 4 iterations. Packaging can also be used to prove that 2 : + ≤ 0.6 is an invariant for the investigated system in 5 while-loop iterations; enhanced Packaging requires only 4 while-loop iterations for the same goal (see Figure 4b).
We presented a set-based method for checking invariant properties over discrete time dynamical systems. The method exploits the notion of -inductiveness together with over-approximation procedures for reachability, unions, and intersections. When it converges to a positive answer, it computes the “smallest” invariant that is a subset of the given one. Our approach achieves convergence on a larger class of systems when set-based operations are over-approximated rather then when exact computations are performed.
The approach has been implemented and tested in Sapo providing two possible semantics
for the over-approximation of unions. The implementation was able to prove two invariants
over a well-known and deeply investigated fractal example [
As future work we plan to investigate other examples and eventually include heuristics for allowing convergence in limit cases.
This work is partially supported by PRIN project NiRvAna CUP G23C22000400005 and National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 1.4 - Call for tender No. 3138 of 16 December 202, rectified by Decree n.3175 of 18 December 2021 of Italian Ministry of University and Research funded by the European Union - NextGenerationEU; Project code CN_00000033, Concession Decree No. 1034 of 17 June 2022 adopted by the Italian Ministry of University and Research, CUP G23C22001110007, Project title “National Biodiversity Future Center - NBFC” .