Set-Based Invariants over Polynomial Systems Alberto Casagrande1,* , Alessandro Cimatti2 , Luca Dorigo2,3 , Carla Piazza3 and Stefano Tonetta2 1 Dip. di Matematica e Geoscienze, UniversitΓ  di Trieste, via Alfonso Valerio 12/1, 34127, Trieste, Italy 3 Dip. di Matematica, Informatica e Fisica, UniversitΓ  di Udine, via delle Scienze 205, 33100, Udine, Italy 2 Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo, Trento, Italy Abstract Dynamical systems model the time evolution of both natural and engineered processes. The automatic analysis of such models relies on different 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. Keywords Dynamical Systems, Infinite State Systems, Invariants, K-Induction 1. Introduction Invariants are subsets of the state space of dynamical systems that contain all the states reach- able from themself. This notion is crucial in may fields such as control theory [1], software verification [2, 3], and analysis of both continuous and hybrid systems [4, 5]. 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 semi- algebraic invariants [4, 6, 7]. Some others focus on invariant sets representable by using polytopes [8, 5] to avoid the double exponential decision procedure of the semi-algebraic theory. The vast majorities of the works in the literature take advantage of the system continuity to either validate or identify an invariant set. This work focuses on discrete-time transition systems. The applications of this kind of systems go from biology [9, 10] to engineering [11, 12]. We rely on a strengthening of a special case of induction, π‘˜-inductiveness [13], to validate a possible invariant 𝑃 . We start from an initial set of states and we enlarge it to include all the states reachable in π‘˜βˆ’1 discrete transitions. CILC’23: 38th Italian Conference on Computational Logic, June 21–23, 2023, Udine, Italy * Corresponding author. $ acasagrande@units.it (A. Casagrande); cimatti@fbk.eu (A. Cimatti); dorigo.luca001@spes.uniud.it (L. Dorigo); carla.piazza@uniud.it (C. Piazza); tonetta@fbk.eu (S. Tonetta)  0000-0002-8681-1482 (A. Casagrande); 0000-0002-1315-6990 (A. Cimatti); 0000-0002-1396-2167 (L. Dorigo); 0000-0002-2072-1628 (C. Piazza); 0000-0001-9091-7899 (S. Tonetta) Β© 2023 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) This is the candidate invariant and it has to be included in the invariant 𝑃 . If this is the case, we have to prove that if 𝑝 is a path of length π‘˜ starting from 𝑃 and having the first π‘˜ βˆ’ 1 steps in 𝑃 itself, then 𝑝 does not exit from 𝑃 . Other approaches involving π‘˜-inductiveness have been proposed in the literature (see, e.g., [14, 15]). However, they assume the possibility of computing the reverse transition, which is not always possible when dealing with discrete polynomial systems. Indeed our goal was to develop a general method for (possibly) non-reversible discrete transition systems. 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. 2. Preliminaries 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 𝐼 βŠ† 𝑆. The transition relation 𝑅 induces a transition function 𝑇 that maps any set 𝐴 βŠ† 𝑆 into the 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 𝐹𝐡 (𝐴) = 𝐹 (𝐴) ∩ 𝐡. We write 𝐹 π‘˜ (𝐴) to def denote π‘˜ consecutive applications of 𝐹 to 𝐴, i.e., {οΈƒ 𝐴 if π‘˜ = 0 𝐹 π‘˜ (𝐴) = def (οΈ€ π‘˜βˆ’1 )οΈ€ . 𝐹 𝐹 (𝐴) 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 𝐼. 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 𝐼 is a set 𝑃 βŠ† 𝑆 such that ⋃︁ 𝑇 π‘˜ (𝐼) βŠ† 𝑃. (1) π‘˜βˆˆN 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., [7, 8, 5]). 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 difficult 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 efficiency in the case of non-linear systems [16]. K-Inductiveness A standard approach for verifying invariants is the use of k-induction [13] in combination with Bounded Model Checking (BMC) [17, 18]. If TS reaches from 𝐼 a state outside 𝑃 in π‘˜ steps, i.e., 𝑇 π‘˜ (𝐼) ΜΈβŠ† 𝑃 , 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 (𝐼) βŠ† 𝑃 (basis) 𝑇 𝑖 (𝐼) βŠ† 𝑃 β†’ 𝑇 𝑖+1 (𝐼) βŠ† 𝑃 (step) Induction and π‘˜-induction are equivalent over N. Thus, 𝑃 is an invariant if and only if β‹€οΈ€π‘˜βˆ’1 𝑗 (︁⋀︀ )︁ βŠ† 𝑃 𝑗=0 𝑇 (𝐼) (π‘˜-basis) π‘˜βˆ’1 𝑖+𝑗 𝑖+π‘˜ 𝑗=0 𝑇 (𝐼) βŠ† 𝑃 β†’ 𝑇 (𝐼) βŠ† 𝑃 (π‘˜-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 [19]. 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 𝑇 (𝐼) βŠ† 𝑃 (π‘˜-init) 𝑇 (π‘‡π‘ƒπ‘˜βˆ’1 (𝑃 )) βŠ† 𝑃 (π‘˜-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 Set-Based Reachability 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 [20]. In such a context, set-based analysis aims at defining techniques for dealing with possible infinite sets of states, processing each set as a single object. 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 sufficient 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 difference with respect to algorithms such as IC3/PDR [14, 15] which use transition reversibility to compute counter- examples and refine over-approximations of the reachability sets. 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 π‘˜. 3. Set-Based K-Inductiveness 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 (οΈ€ π‘˜βˆ’1 False )οΈ€ 9: if 𝑇 𝑇 CI (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 𝑗 𝑗=0 𝑇 (𝐼). This observation allows us to prove some properties about the results of Algorithm 2. Theorem 3.1 (Correctness of Algorithm 2). If 𝑇 under-approximates 𝑇 and Algorithm 2 re- turns 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 de- pend 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 reachable from⋃︀𝐼 does not converge after a finite⋃︀number of while-loop iterations. of points ⋃︀ 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 finite number of forward transitions. Corollary 3.1.1 (Termination). If 𝑇 (𝐴) exactly represents 𝑇 (𝐴) for any 𝐴 ∈ π’Ÿ, then Al- gorithm 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 𝑃 (b) Algorithm 2 returns False even though 𝑃 is not an invariant for TS and 𝐼 when 𝑇 is an invariant for TS and 𝐼 when 𝑇 over- under-approximates 𝑇 . However, 𝑇 (𝐼) is approximates 𝑇 . However, 𝑇 (𝐼) is a subset 3 not a subset of 𝑃 and, as a consequence, 𝑃 of 𝑃 and 𝑇 (𝐼) = 𝑇 (𝐼). Thus, 𝑃 is an is not an invariant for TS and 𝐼. invariant for TS and 𝐼. Figure 1: Two scenarios in which, when the forward image function is approximated, Algorithm 2 terminates, but does not return the aimed answer. The dashed blue line rectangle depicts 𝑇 (𝐼). 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 π‘˜βˆ’1 π‘˜βˆ’1 𝑇 (𝑇 CI (CI)) βŠ† CI ∨ 𝑇 (𝑇 𝑃 (𝑃 )) βŠ† 𝑃, 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 0 π‘˜ π‘˜βˆ’1 consists in the list [𝑇 (𝐼), . . . , 𝑇 (𝐼)]. Thus, the evaluation of 𝑇 CI (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. Over-Approximation of the Union 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 different 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 [21]. 𝑗 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 1: function join-candidate-inductive(TS = βŸ¨π‘†, π‘…βŸ©, 𝐼, 𝑃 ) 2: if 𝐼 ΜΈβŠ† 𝑃 then 3: return False 4: 𝐽 ←𝐼 5: π‘˜β†1 6: CI ← 𝐽 7: while True do π‘˜ 8: if 𝑇 (𝐽) ΜΈβŠ† 𝑃 then 9: return (οΈ€ π‘˜βˆ’1 False )οΈ€ 10: if 𝑇 𝑇 CI (CI) βŠ† CI then 11: return True π‘˜ 12: CI ← Join(CI, 𝑇 (𝐽)) 13: if CI ΜΈβŠ† 𝑃 then π‘˜ 14: 𝐽 ← 𝑇 (𝐽) 15: π‘˜β†0 16: CI ← 𝐽 17: π‘˜ β†π‘˜+1 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 over- approximates 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 𝑇 (𝐽) 2 𝑇 (𝐽) 𝑃 𝐽 CI β‹ƒοΈ€π‘˜ 𝑖 Figure 2: A scenario in which CI ΜΈβŠ† 𝑃 despite 𝑖=0 𝑇 (𝐽) βŠ† 𝑃 when the Join operator of Algorithm 3 approximates the union as the bounding box. 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. Approximating Intersections The analysis of Algorithm 3 in the previous section relies on the hypothesis that intersections, π‘˜βˆ’1 needed to compute 𝑇 CI (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 affect 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 𝐼. 4. Implementation and Tests over Polynomial Systems 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 [10, 22, 16] implements bounded reachability, verification of Signal Temporal Logic specifications, and parameter synthesis for such systems. The framework supports state sets specified as polytopes, i.e., closed and bounded subsets of R𝑛 enclosed by a finite number of hyperplanes. A direction of the polytope is a unit vector normal to one of the enclosing hyper- plane. Forward images of polytopes are over-approximated exploiting Bernstein’s coefficients. The forward image of a polytope is a new polytope specified by the same directions of the initial one, but with new thresholds. A detailed description of both algorithms and data structures implemented by Sapo is given in [16]. 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 Al- gorithm 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; 2. implement the test 𝐴 βŠ† 𝑃 where 𝐴 is a generic polytope and 𝑃 is the candidate invariant; 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 Β¬(𝐴 ΜΈβŠ† 𝑃 ). The implementation of the test 𝐴 ΜΈβŠ† 𝑃 relies on the solution of up to π‘š linear programming problems where π‘š is the number of linear inequalities in ℒ𝑃 . If the problems are solved by using Karmarkar’s algorithm [23], the asymptotic complexity of the algorithm is polynomial with respect to both the space dimension and the number of constraints in ℒ𝐴 and ℒ𝑃 . 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 effective 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 {οΈ‚ π‘₯π‘˜+1 = π‘¦π‘˜ . π‘¦π‘˜+1 = π‘₯π‘˜ 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. 𝑦 𝑦 π‘₯ + 𝑦 = 5 𝐼 𝐼 π‘₯ + 𝑦 = 5 CI CI 𝑇 (𝐼) 𝑇 (𝐼) 𝑃 :π‘₯+𝑦 ≀5 𝑃 :π‘₯+𝑦 ≀5 π‘₯ π‘₯ (a) Algorithm 3 with the Packaging se- (b) Algorithm 3 with the enhanced mantics cannot validate the invariant Packaging semantics that adds the π‘₯ + 𝑦 ≀ 5. The bounding box of 𝐼 and candidate invariant directions to the 𝑇 (𝐼) violates π‘₯ + 𝑦 ≀ 5 and the algo- packaging can validate the invariant rithm loops forever. π‘₯ + 𝑦 ≀ 5 in 2 while-loop iterations. Figure 3: A graphical representation of Example 3. The system evolution starts from set labelled 𝐼. One step brings the system in the state set 𝑇 (𝐼) and one further step let it go back to its initial condition. 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 different 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 [16]. We now briefly illustrate their effectiveness on the following transition system over R2 {οΈ‚ π‘₯π‘˜+1 = 2π‘₯π‘˜ π‘¦π‘˜ + 𝑏 π‘¦π‘˜+1 = π‘¦π‘˜2 βˆ’ π‘₯2π‘˜ + π‘Ž where π‘Ž, 𝑏 ∈ R. This system represents the complex quadratic generator function 𝑔𝑐 (𝑧) = 𝑧 2 +𝑐 with 𝑧, 𝑐 ∈ C which is used to define Mandelbrot set [24, 25]. The variables π‘₯ and 𝑦 are the real and imaginary components of 𝑧, respectively, and similarly π‘Ž and 𝑏 are the components of 𝑐. The dynamical system iterates the computation of 𝑔𝑐 . Let [0.09, 0.11] Γ— [0.09, 0.11] be the system initial set 𝐼. Different choices of π‘Ž and 𝑏 can produce very different 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). 𝑦 𝑦 0.2 0.2 𝑃2 0.15 0.15 CI 0.1 𝐼 0.1 𝐼 CI 0.05 0.05 𝑃2 0.1 0.2 0.3 0.4 π‘₯ 0.1 0.2 0.3 0.4 π‘₯ (a) Packaging: CI is not included in 𝑃2 at (b) enhanced Packaging: the candidate in- the beginning the second iteration of Algo- variant directions are added to 𝐼 specifica- π‘˜ rithm 3 while-loop, thus, CI is assigned to tion. The set CI and all the 𝑇 (𝐼) with π‘˜ > 0 2 𝑇 (𝐼). After 3 further while-loop iterations, are no more rectangles and 𝑃2 always in- CI becomes the bounding box of the forward cludes them. Thus, CI is never re-initialized images which is 3-inductive. and it eventually becomes 4-inductive. Figure 4: The Mandelbrot system with π‘Ž ∈ [0.19, 0.2] and 𝑏 ∈ [0.29, 0.3], initial set 𝐼 = [0.09, 0.11] Γ— [0.09, 0.11], and 𝑃2 : π‘₯+𝑦 ≀ 0.6 as candidate invariant. The red area represents the candidate invariant, the green region is CI during the 𝑛-th and final iteration of the while-loop, the blue regions are the over-approximation of the set reachable in up to 𝑛 steps. 5. Conclusions and Future Work 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 [24, 25]. As future work we plan to investigate other examples and eventually include heuristics for allowing convergence in limit cases. Acknowledgments 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” . References [1] F. Blanchini, Set invariance in control, Automatica 35 (1999) 1747–1767. [2] X. Deng, M. B. Dwyer, J. Hatcliff, M. Mizuno, Invariant-based specification, synthesis, and verification of synchronization in concurrent programs, in: Proceedings of the 24th international conference on software engineering, 2002, pp. 442–452. [3] B. Becker, D. Beyer, H. Giese, F. Klein, D. Schilling, Symbolic invariant verification for systems with dynamic structural adaptation, in: Proceedings of the 28th international conference on Software engineering, 2006, pp. 72–81. [4] A. Platzer, E. M. Clarke, Computing differential invariants of hybrid systems as fixedpoints, Formal Methods in System Design 35 (2009) 98–120. [5] M. A. B. Sassi, A. Girard, Computation of polytopic invariants for polynomial dynamical systems using linear programming, Automatica 48 (2012) 3114–3121. [6] S. Sankaranarayanan, Automatic invariant generation for hybrid systems using ideal fixed points, in: Proceedings of the 13th ACM international conference on Hybrid systems: computation and control, 2010, pp. 221–230. [7] A. Sogokon, K. Ghorbal, P. B. Jackson, A. Platzer, A method for invariant generation for polynomial continuous systems, in: International Conference on Verification, Model Checking, and Abstract Interpretation, Springer, 2016, pp. 268–288. [8] S. Sankaranarayanan, T. Dang, F. IvančiΔ‡, Symbolic model checking of hybrid systems using template polyhedra, in: TACAS ’08, Springer, 2008, pp. 188–202. [9] A. Veliz-Cuba, A. S. Jarrah, R. Laubenbacher, Polynomial algebra of discrete models in systems biology, Bioinformatics 26 (2010) 1637–1643. [10] T. Dang, T. Dreossi, C. Piazza, Parameter synthesis using parallelotopic enclosure and applications to epidemic models, in: Hybrid Systems and Biology, HSB, 2014, pp. 67–82. [11] H. Ichihara, S. Tanabe, Y. Ebihara, D. Peaucelle, Analysis and synthesis of discrete-time interconnected positive systems, SICE Journal of Control, Measurement, and System Integration 11 (2018) 91–99. [12] T. Dreossi, T. Dang, C. Piazza, Reachability computation for polynomial dynamical systems, Formal Methods in System Design 50 (2017) 1–38. [13] M. Sheeran, S. Singh, G. StΓ₯lmarck, Checking safety properties using induction and a sat-solver, in: W. A. Hunt, S. D. Johnson (Eds.), Formal Methods in Computer-Aided Design, Springer Berlin Heidelberg, Berlin, Heidelberg, 2000, pp. 127–144. [14] A. R. Bradley, Sat-based model checking without unrolling, in: International Workshop on Verification, Model Checking, and Abstract Interpretation (VMCAI 2011), volume 6538 of Lecture Notes in Computer Science, Springer, 2011, pp. 70–87. [15] N. EΓ©n, A. Mishchenko, R. Brayton, Efficient implementation of property directed reachabil- ity, in: International Conference on Formal Methods in Computer-Aided Design (FMCAD 2011), IEEE, 2011, pp. 125–134. [16] A. Casagrande, T. Dang, L. Dorigo, T. Dreossi, C. Piazza, E. Pippia, Parameter synthesis of polynomial dynamical systems, Information and Computation 289 (2022) 104941. [17] A. Biere, A. Cimatti, E. Clarke, Y. Zhu, Symbolic model checking without BDDs, in: TACAS ’99, volume 1579 of Lecture Notes in Computer Science, Springer, 1999, pp. 193–207. [18] A. Biere, A. Cimatti, E. M. Clarke, O. Strichman, Y. Zhu, Bounded model checking, Handbook of satisfiability 185 (2009) 457–481. [19] A. F. Donaldson, L. Haller, D. Kroening, P. RΓΌmmer, Software Verification Using k-Induction, in: E. Yahav (Ed.), Static Analysis, Springer, Berlin, Heidelberg, 2011, pp. 351–368. [20] A. Casagrande, T. Dreossi, C. Piazza, Hybrid automata and πœ–-analysis on a neural oscillator, Electronic Proceedings in Theoretical Computer Science, EPTCS 92 (2012) 58 – 72. [21] P. Cousot, N. Halbwachs, Automatic discovery of linear restraints among variables of a program, in: Proceedings of the 5th ACM SIGACT-SIGPLAN Symposium on Principles of Programming Languages, POPL ’78, Association for Computing Machinery, New York, NY, USA, 1978, pp. 84–96. [22] T. Dang, T. Dreossi, C. Piazza, Parameter synthesis through temporal logic specifications, in: Formal Methods, FM, 2015, pp. 213–230. [23] N. Karmarkar, A new polynomial-time algorithm for linear programming, in: Symposium on Theory of computing, STOC, ACM, 1984, pp. 302–311. [24] R. Brooks, J. P. Matelski, The dynamics of 2-generator subgroups of PSL(2, 𝑐), in: Riemann surfaces and related topics: Proceedings of the 1978 Stony Brook Conference, Ann. of Math. Stud, volume 97, 1981, pp. 65–71. [25] B. B. Mandelbrot, Fractal aspects of the iteration of 𝑧 ← πœ† 𝑧(1 βˆ’ 𝑧) for complex πœ† and 𝑧, Annals of the New York Academy of Sciences 357 (1980) 249–259.