=Paper=
{{Paper
|id=Vol-3428/paper8
|storemode=property
|title=Set-Based Invariants over Polynomial Systems
|pdfUrl=https://ceur-ws.org/Vol-3428/paper8.pdf
|volume=Vol-3428
|authors=Alberto Casagrande,Alessandro Cimatti,Luca Dorigo,Carla Piazza,Stefano Tonetta
|dblpUrl=https://dblp.org/rec/conf/cilc/CasagrandeCDPT23
}}
==Set-Based Invariants over Polynomial Systems==
<pdf width="1500px">https://ceur-ws.org/Vol-3428/paper8.pdf</pdf>
<pre>
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.

</pre>