Qi Lu

Michael Madsen

Martin Milata
email: xmilata@fi.muni.cz

Søren Ravn

Uli Fahrenberg
email: ulrich.fahrenberg@irisa.fr

Kim G Larsen

Reachability Analysis for Timed Automata using Max-Plus Algebra

Keywords: timed automaton, real-time model checking, data structure, max-plus algebra, max-plus polyhedron

We show that max-plus polyhedra are usable as a data structure in reachability analysis of timed automata. Drawing inspiration from the extensive work that has been done on difference bound matrices, as well as previous work on max-plus polyhedra in other areas, we develop the algorithms needed to perform forward and backward reachability analysis using maxplus polyhedra. To show that the approach works in practice and theory alike, we have created a proof-of-concept implementation on top of the model checker opaal.

Introduction

A real-time system is a system where total completion of a task depends not only on the logical ordering of events, but also on the timing at which these events are performed. Examples of real-time systems include airbags, pacemakers, live video streaming, video game systems, and production lines.

A key problem when developing a real-time system is to ensure correctness of the system. For that purpose, it is useful to construct a model of the system and verify certain properties directly on the model. This is known as model checking.

For real-time systems, the model of timed automata introduced by Alur and Dill [START_REF] Alur | Automata for modeling real-time systems[END_REF] has emerged as a preferred formalism, allowing most relevant real-time aspects to be expressed while leaving model checking decidable. Being essentially finite-state systems extended with a finite number of realvalued clocks used for conditioning when transitions are available, the formalism of timed automata describes systems with infinitely many states. However, the development of data structures for symbolic representation and manipulation of state spaces of timed automata has been subject to significant research effort, leading not only to decidability but also efficient model checking tools, e.g. Uppaal [32] and Kronos [START_REF] Yovine | KRONOS: A verification tool for real-time systems[END_REF].

In particular, the notion of zones (state sets that can be described by constraints on individual clocks and clock-differences) is a main datastructure in these tools. A commonly used representation of zones are Difference Bound Matrices, or DBMs for short [START_REF] Alur | An implementation of three algorithms for timing verification based on automata emptiness[END_REF][START_REF] David | Timing assumptions and verification of finite-state concurrent systems[END_REF], due to their efficient time and space properties. A comprehensive account on DBM-based algorithms for zone-based reachability analysis can be found in [START_REF] Bengtsson | and States in Timed Systems[END_REF]; the algorithms used in tools like Uppaal are based on DBM representations of zones.

There are however a number of open problems regarding DBM-based reachability analysis, and other data structures to complement or replace DBMs have been proposed [START_REF] Behrmann | Efficient timed reachability analysis using clock difference diagrams[END_REF][START_REF] Asarin | Data-structures for the verification of timed automata[END_REF][START_REF] Ehlers | Fully symbolic timed model checking using constraint matrix diagrams[END_REF][START_REF] Jesper | Difference decision diagrams[END_REF][START_REF] Bozga | Some progress in the symbolic verification of timed automata[END_REF]. One particular such problem is that zones are not closed under unions, that is, the union of two zones may not again be a zone. This contributes to the so-called state-space explosion during zone-based reachability analysis, unless over-or underapproximations are used.

It is the purpose of this paper to propose replacing zones and DBMs by max-plus polyhedra for real-time reachability analysis. We develop the algorithms necessary for forward and backward reachability analysis using this new data structure and show that their complexities are comparable to the one of the standard DBM-based algorithms. Additionally, overapproximations using max-plus polyhedra are more precise than the ones one achieves using DBMs, which is of great help in combating state-space explosion.

Related work. The study of max-plus analogues of convex sets in max-plus algebra begins with [START_REF] Zimmermann | A general separation theorem in extremal algebras[END_REF]. Originally motivated by some problems in abstract interpretation, [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF] introduce max-plus polyhedra as a new numerical abstraction which can express some non-convex properties without any disjunctive representations. The theory of max-plus polyhedra is further developed in [START_REF] Xavier Allamigeon | Computing the extreme points of tropical polyhedra[END_REF][START_REF] Xavier Allamigeon | The tropical double description method[END_REF][START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF], where also fundamental algorithms to compute with these structures are given.

The idea of using max-plus polyhedra as a data structure for real-time reachability analysis was communicated to us by Eric Goubault. In a previous work [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF] we provide some basic algorithms for forward reachability. Other related work in abstract interpretation include zones [START_REF] Miné | A new numerical abstract domain based on differencebound matrices[END_REF], classical polyhedra [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], octagons [START_REF] Miné | The octagon abstract domain[END_REF] and disjunctive representations [START_REF] Sriram Sankaranarayanan | Static analysis in disjunctive numerical domains[END_REF][START_REF] Giacobazzi | Compositional optimization of disjunctive abstract interpretations[END_REF].

Structure. We start this paper with three sections reviewing timed automata, zones, max-plus polyhedra, and DBMs, finishing with an example which compares the expressivity of DBMs and max-plus polyhedra. Section 5 then introduces the algorithms for reachability analysis based on max-plus polyhedra, analyses their complexities, and compares them with the standard DBM-based algorithms. In the conclusive Section 6 we sum up on our work and lay out some open problems which are still waiting to be solved.

This paper is based on a Bachelor's thesis [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF] and on a project report [START_REF] Lu | Computations on zones using max-plus algebra[END_REF] by the first four authors and Jesper Dyhrberg. On both projects, the last two authors of this paper have acted as supervisors.

Timed automata

Timed automata are finite automata that are enriched with a number of real-valued clocks. They have proven to be useful in modeling and verification of real-time systems. Timed automata were introduced by [START_REF] Alur | A theory of timed automata[END_REF]; a more thorough overview of relevant notions can be found for example in [START_REF] Aceto | Reactive Systems[END_REF] or [START_REF] Fahrenberg | Verification, performance analysis and controller synthesis for real-time systems[END_REF].

Syntax

First, let X be a finite set of real-valued non-negative variables referred to as clocks. Define B(X) to be set of all clock constraints g generated by the following grammar:

g ::= x 1 ∼ n | x 1 -x 2 ∼ n | g 1 ∧ g 2 ,
where x 1 , x 2 ∈ X, n ∈ ◆ is a natural number and ∼ ∈ {<, ≤, =, ≥, >}.

A timed automaton is a quintuple (L, X, l 0 , E, I), where

• L is a finite set of locations, • X is a finite set of clocks, • l 0 ∈ L is the initial location, • E ⊆ L × B(X) × 2 X × L is the set of edges and
• I : L → B(X) a function which assigns to every location an invariant.

To denote edges, we will write l 1 g,r -→ l 2 instead of (l 1 , g, r, l 2 ). Here, l 1 is the source location of the edge, l 2 is the destination, g is the guard of the edge and r is the set of clocks to be reset after transition.

Strict constraints. In this paper (excluding the preliminary section) we will only work with closed timed automata, i.e. timed automata on clock constraints which only use ∼ ∈ {≤, =, ≥}. It is well-known that enforcing closedness restricts expressivity, but currently our algorithms (to be presented in Section 5) only work with these non-strict constraints. Unlike DBMs, the representation of max-plus polyhedra we use here does not allow a simple extension to strict inequalities; see Section 6.2 on page 25 for further discussion.

Semantics

We will first try to give a very informal description of the semantics of timed automata. The state of a timed automaton is composed of its control location and values of each clock. In the initial state, the location is set to l 0 and the value of all clocks is zero. Whenever the automaton is in some control location l, it has two choices.

1. Similarly to finite state automata, it can do a transition over one of the edges (l, g, r, l ′ ) that lead from it. However, this is only possible when the current clock values satisfy the guard g of the edge. After the transition, all the clocks in the set r of the edge are reset to zero. Additionally, the invariant of the destination state must be satisfied by the clock valuation after resetting. 2. It can "stay" in location l for some period of time. This means that all the clocks increase by the same amount of time. The values of the clocks must satisfy the location invariant I(l) during this whole period.

To formally define the semantics of a timed automaton, we must first define clock valuations. A clock valuation v is a function v : X → | ≥0 , that assigns a non-negative value to each clock. Let δ ∈ | ≥0 and r ⊆ X. We define v + δ to be the valuation such that

(v + δ)(c) = v(c) + δ and v[r] to be the valuation such that v[r](c) = 0 whenever c ∈ r and v[r](c) = v(c) otherwise.
Formally, the semantics of a timed automaton (L, X, l 0 , E, I) is given in terms of a transition system (S, s 0 , →):

• S = {(l, v) | l ∈ L, v : X → | ≥0 , v I(l)} is the set of states,
• s 0 = (l 0 , v 0 ), where v 0 is the clock valuation that assigns 0 to all clocks, is the initial state,

• → ⊆ S × S are transitions such that:

-(l, v) → (l ′ , v ′ ) if l g,r -→ l ′ , v g, v ′ = v[r], -(l, v) → (l, v + δ) for all δ ∈ | ≥0 such that v + t I(l) for any 0 ≤ t ≤ δ.
Such a transition system is in most cases infinite and even uncountable, which means it cannot be directly used in an algorithm. Fortunately, we can also construct transition systems which are finite -this is achieved by replacing individual states with symbolic states, where each such state consists of a control location and a set of clock valuations. We naturally require those sets to have a finite description and the clock valuations contained in them to be in some way equivalent, which usually means that they have to be untimed bisimilar, see [START_REF] Aceto | Reactive Systems[END_REF].

An example of such symbolic semantics is the region graph, where the sets of clock valuations are divided into equivalence classes based on integer parts, orderings of fractional parts of clock values and whether or not the value is greater than some fixed constant. The number of such classes grows very quickly with the number of clocks, hence the region graph is not really suitable for algorithmic use. Nevertheless, region graphs have been shown to be useful when proving decidability of some properties of timed automata.

There is, however, another variant of symbolic semantics, which is suitable for practical use.

Zones and zone graphs

Zones are sets of clock valuations that satisfy a conjunction of clock constraints. Formally, Z ⊆ | X ≥0 is a zone if there is a g ∈ B(X) such that Z = {v | v g}. Zones can also be thought of as convex subsets of |X|dimensional Euclidean space.

In order to define transitions on symbolic states of the form (l, Z), where l is a location and Z is a zone, we again need to define some operations first. Let Z be a zone, g ∈ B(X) and r ⊆ X.

• Z ∧ g = {v ∈ Z | v g}, • Z ↑ = {v + δ | v ∈ Z, δ ∈ | ≥0 }, • Z[r] = {v[r] | v ∈ Z}.
Lemma 2.1. Let Z be a zone, g ∈ B(X) and r ⊆ X. Then Z ∧ g, Z ↑ and Z[r] are also zones [START_REF] Henzinger | Symbolic model checking for real-time systems[END_REF].

We can now define the symbolic successor relation as follows:

• (l, Z) (l, Z ↑ ∧ I(l)) -delay successor, • (l, Z) (l ′ , (Z ∧ g)[r] ∧ I(l ′ )) if l g,r -→ l ′ -discrete successor.
Let Z 0 be a zone containing just the one valuation which assigns zero to all clocks, and (l 0 , Z 0 ) our initial symbolic state. All this together gives us a transition system on symbolic states, called the zone graph. We are usually only interested in the part that is reachable from (l 0 , Z 0 ), but it still may be the case that even this part is infinite.

The zone graph can be made finite by the process of extrapolation (sometimes also referred to as normalization), which exploits the fact that once a clock value exceeds the maximal constant the clock is compared to in the constraints of the automaton, its precise value becomes irrelevant. There exist several such operations [START_REF] Behrmann | Lower and upper bounds in zone-based abstractions of timed automata[END_REF][START_REF] Bouyer | Diagonal constraints in timed automata: Forward analysis of timed systems[END_REF], which will make the state space finite while preserving its properties (i.e. reachability of a state in our case).

From a practical point of view, the representation of zones is very important. We obviously cannot represent a zone as a list of the clock valuations it contains. A logical approach is to represent zones by the constraints that define them, which is the basic underlying principle of the DBM data structure, which is nowadays most commonly used in tools for timed automata analysis. Section 4 on page 11 is devoted to describing this data structure.

Deciding reachability in zone graphs

A zone graph can be directly used to decide whether a particular state is reachable in a timed automaton, i.e. whether there is a run of the automaton that reaches the state. Although the algorithms is basically a depth-first search on the zone graph, we show the forward reachability algorithm, Algorithm 1 on the following page, as to give the motivation for the operations needed.

The input of the algorithm is a timed automaton together with a description of a state to check for reachability, i.e. a location s and a constraint ϕ ∈ B(X).

The algorithm keeps two sets of symbolic states. The Passed set contains already processed states, and the Waiting set contains the initial state at the beginning, for states yet to be processed. The body of the main loop picks a state from the Waiting set, checks whether it satisfies ϕ, and terminates with positive answer if it does. Otherwise it checks if the state is already covered by the passed states, i.e. whether it is included in an already visited zone with the same control location. If this is not the case, the state is added to the Passed set and its successors to the Waiting set. This is repeated as long as the Waiting set is nonempty. Choose and remove (l, Z) from Waiting if Z Z ′ for all (l, Z ′ ) ∈ Passed then 9:

Passed := Passed ∪ {(l, Z)} 10:

Waiting := Waiting ∪ {(l ′ , Z ′ ) | (l, Z) (l ′ , Z ′ ) ∧ Z ′ = ∅} 11:
end if 12: end while 13: return False We can see from the algorithm that in order to develop a data structure that supports forward reachability analysis, we need it to support following operations:

• decide whether some parts of the zone satisfy a constraint ϕ,

• decide whether it is a subset of another zone,

• compute successors of a state, which can be achieved by the three operations listed in Section 2.3 on page 5 and an additional check whether the zone is empty.

Analogously to the forward reachability algorithm above, also backward reachability analysis, which computes predecessors instead of successors, is used. In order to support this analysis, we additionally need operations for backward delay ↓ Z and freeing clocks [r]Z which are the inverses of delay Z ↑ and reset Z[r].

Max-plus algebra

Let | max denote the set | ∪ {-∞}, and let a ⊕ b = max(a, b) and a ⊗ b = a + b.
The max-plus algebra is the semiring (| max , ⊕, ⊗), that is, the set of real numbers equipped with zero element -∞ with the maximum operation as addition and ordinary addition as multiplication. To further conform to the usual semiring notation, we denote -∞ as ✵ and 0, the neutral element with respect to max-plus multiplication, as ✶. As with ordinary multiplication, we will also use the convention that ab = a ⊗ b.

The definitions of addition and multiplication in max-plus algebra can be extended to vectors and matrices in the usual way -addition:

(v 1 , . . . , v n )⊕ (w 1 , . . . , w n ) = (v 1 ⊕ w 1 , . . . , v n ⊕ w n ), multiplication by scalar: α ⊗ (v 1 , . . . , v n ) = (α⊗v 1 , . . . , α⊗v n ), matrix multiplication: (AB) ij = n k=1 A ik ⊗ B kj .
Note that the max-plus semiring is idempotent, because a ⊕ a = a for any element a, and that it is not a ring, because for every a except a = ✵, there is no b such that a ⊕ b = ✵.

The dual notion to max-plus algebra, used in some literature, is the minplus algebra, also called tropical semiring, 3 in which the maximum operation is replaced by the minimum operation and positive infinity is used as zero element. This has also lead some authors to call max-plus algebra the arctic semiring.

Notation. In the rest of the paper we will mainly use the letters a, b, c . . . to denote elements of | max . Greek letters α, β, . . . will be used for elements of | max in context of scalar multiplication. Whenever speaking of dimension makes sense, we will denote it n. Vectors, i.e. elements of | n max , will be denoted v, w, . . .; we will use v i to mean i-th component of vector v, as is usual. Whenever we work with an indexed family of vectors, the indices of individual elements will be written in superscript, i.e. v 1 , . . . , v m to avoid confusion with selecting the element of a vector.

Polyhedra over max-plus algebra

Convex max-plus polyhedra are the max-plus analogues of classical convex polyhedra. Definition 3.1. A convex max-plus polyhedron is a subset of | n max that satisfies a finite set of (max-plus) linear inequalities.

The word convex refers to the property that any max-plus line segment between two points of the set is contained in the set. As shown in Figure 1 on the following page, which displays convex max-plus polyhedra in 2D and 3D, we can see that this is different from the classical understanding of convexity.

Whenever we mention polyhedra in the paper, we refer to closed convex max-plus polyhedra, unless stated otherwise. Max-plus convex sets were introduced by [START_REF] Zimmermann | A general separation theorem in extremal algebras[END_REF], a general introduction can be found for example in [START_REF] Gaubert | The Minkowski theorem for max-plus convex sets[END_REF].

Representation of max-plus polyhedra

Because polyhedra are usually infinite sets of points, we need to represent them in some finite way, which we furthermore would like to be able to manipulate efficiently. In the following, we consider three representations of max-plus polyhedra: systems of constraints, sets of generators and max-plus cones of higher dimension; the third being a slight variation of the second. The external representation is quite similar to a DBM; indeed a DBM can easily be converted into a polyhedron in external representation using the fact that, e.g., a clock constraint x ix j ≥ n can be written as the max-plus inequality ax ⊕ b ≥ cx ⊕ d with a k = ✵ for k = i and a i = ✶, c k = ✵ for k = j and c j = n, and b = d = ✵; see [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF] for further details.

An interesting difference from ordinary linear algebra is that a system of inequalities can be represented as a system of equalities (and vice versa). This follows from the fact that a ≥ b ⇔ a = a ⊕ b. We can therefore also choose to represent max-plus polyhedra as systems of equalities.

Sets of generators. Let

A, B ⊆ | n max . The Minkowski sum is defined as A⊕B = {a ⊕b | a ∈ A, b ∈ B}. Vector a is a convex combination of vectors b 1 , . . . , b m if a = m i=1 α i b i for some scalars α i such that m i=1 α i = ✶. Let co(A)
denote the set of all convex combinations of points from A, i.e. its convex hull.

Vector a is a linear combination of vectors b 1 , . . . , b m if a = m i=1 α i b i for some scalars α i . A max-plus cone is a set that is closed under linear combination. The max-plus cone generated by A, denoted cone(A), is the set of all linear combinations of members of A. Max-plus cones are the analogues of vector (sub)spaces of vectors over a field, and they are studied in different contexts under the name semimodules [START_REF] Gaubert | The Minkowski theorem for max-plus convex sets[END_REF].

The following theorem is an analogue of a similar theorem for classical polyhedra.

Theorem 3.2. Any max-plus polyhedron can be represented as Minkowski sum of a bounded convex set and a cone, both of which are finitely generated. Additionally, the Minkowski sum of some bounded convex set and some cone, both of which are finitely generated, represents some max-plus polyhedron. [START_REF] Gaubert | The Minkowski theorem for max-plus convex sets[END_REF] In other words, max-plus polyhedra can be represented as co(V ) ⊕ cone(W ), where V and W are finite sets of points. This can be written explicitly as

α 1 v 1 ⊕ • • • ⊕ α p v p ⊕ β 1 w 1 ⊕ • • • ⊕ β q w q ,
where p i=1 α i = ✶ and β 1 , . . . , β q ∈ | max . The representation is also referred to as internal. Furthermore, every polyhedron has a unique minimal V, the elements of which are called extreme points [START_REF] Gaubert | The Minkowski theorem for max-plus convex sets[END_REF][START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF][START_REF] Butkovič | Generators, extremals and bases of max cones[END_REF]. Extreme points have the property that they cannot be written as a convex combination of any other points in V . There is no unique minimal W , because any generator of the minimal set can be replaced by a scalar multiple. The set of all scalar multiples of a vector is called a ray, and every cone has a unique minimal set of extreme rays that generate it. Therefore the minimal W is only unique up to the choice of representative for each ray.

The internal representation is radically different from the external representation, and the conversion between these representations involves solving max-plus matrix equations, an operation which is computationally rather expensive [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF]. Consequently, throughout this paper, we will only work with the internal representation.

Homogeneous coordinates for max-plus polyhedra. Max-plus polyhedra in n dimensions can be also represented as max-plus cones in n+1 dimensions, which can be thought of as using homogeneous coordinates. Let P = co(V )⊕ cone(W ) be the max-plus polyhedron generated by the sets [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF]. The requirement that the last component of the vector must be equal to ✶ enforces that the scalars used to multiply elements that were in V sum to ✶, while no restriction is placed on elements of W . The advantage of representing polyhedra as cones is that we don't have to distinguish between two kinds of generators, which allows the algorithms to be considerably simpler.

V, W ⊆ | n max . Now, let Z ⊆ | n+1 max be defined as Z = {(v, ✶) | v ∈ V } ∪ {(w, ✵) | w ∈ W }. It can be seen that P = {x | (x, ✶) ∈ cone(Z)}, cf.

Choice of representation.

Owing to the high cost of converting between external and internal representation, it is necessary for our work to choose one representation and then stick to it. We have here chosen to use the internal representation using sets of generators, and we convert freely back and forth between the convex-cone representation and the one using homogeneous coordinates.

One advantage of using the internal representation is that overapproximating unions of max-plus polyhedra can be computed easily, cf. Section 5.3.6 on page 22. The external representation on the other hand resembles much more DBMs (see the next section), hence one might be able to re-use many DBM algorithms; also, strict constraints can be handled more easily in the external representation (see Section 6.2 on page 25). By the resemblance of the external representation to DBMs however, we conjecture that one cannot obtain major speedups using external-representation based reachability analysis vs. DBM-based analysis.

Difference Bound Matrices

Difference Bound Matrices [START_REF] David | Timing assumptions and verification of finite-state concurrent systems[END_REF] is currently one of the most efficient data structures for representing zones [START_REF] Bengtsson | and States in Timed Systems[END_REF]. A DBM is, as the name suggests, a matrix with entries representing the difference between clocks. To be able to do this in a uniform way for both regular difference constraints as well as comparing just one clock to a constant, a zero clock, 0, with the constant value 0, is introduced. This approach benefits from the fact that, as mentioned in Section 2. To compute the DBM for a zone, every clock in C 0 is numbered, assigning one column and one row to the clock. Every entry in the matrix, D, now represents the bound x i -x j n where x i , x j are clocks, i is the row index of the matrix, j the column index. This means that rows and columns encode lower and upper bounds respectively.

To be able to handle strictness, an entry in the DBM is not just a value, but rather the tuple (n, ) where n ∈ ) and ∈ {≤, <}, representing the bound x ix j n. As we will only be concerned with non-strict inequalities here however, we can simplify representations and only use values.

When no bound is present for the given clock difference, ∞ is used, since everything is less than or equal to infinity. Additionally, since all clocks are positive, the implicit constraints 0x i ≤ 0 are added, and since the difference between a clock and itself should always be 0, x ix i ≤ 0 is added as well.

Finally, to be able to manipulate DBMs, comparison and addition of bounds must be defined. Both are straight-forward extensions of integer comparison and addition; see [START_REF] Bengtsson | and States in Timed Systems[END_REF] for the details.

Since there can be infinitely many different conjunctions of constraints representing the same solution set and thereby the same actual zone, a canonical representation is necessary. The canonical representation for DBMs is the one representing the tightest constraints, without altering the solution set. This canonical representation can be computed by converting the DBM into a directed graph, where clocks are represented by nodes and difference constraints are labelled edges between the appropriate nodes. Now all there is to do is compute the shortest path between nodes, e.g. by using the Floyd-Warshall algorithm [START_REF] Floyd | Algorithm 97: Shortest path[END_REF], and then converting back to matrix form [START_REF] Bengtsson | and States in Timed Systems[END_REF].

DBM operations

As mentioned above, DBMs provide efficient algorithms for reachabilityanalysis operations. This includes linear-time algorithms for the basic operations of delaying and resetting as well as freeing a clock. For checking inclusion of DBMs and intersecting with one clock constraint, quadratic algorithms are provided. This is also the case for backward delay and extrapolation. For comparison of the DBM algorithms to our proposed algorithms on max-plus polyhedra we refer the reader to Section 5 on page 15.

Example

Since DBMs are the current industry standard for performing real-time model checking, it is useful to provide a more direct comparison between max-plus polyhedra and DBMs. Figure 2 shows an example of a small TA of two clocks which benefits from using max-plus polyhedra rather than DBMs. Imagine that this is just the tiny initial part with start location l 0 of a much bigger TA connected by the transition going out of location l 3 .

Running the forward reachability algorithm on this example will give us the following results in locations l 1 and l 2 , shown as a zone Z, a max-plus polyhedron P and a DBM D. A visual indication of the two zones is shown in Figure 3. In both locations, the status is that we are ready to take the transition out.

l 0 l 1 l 2 l 3 R [ x = 0 ] R [ y = 0 ] y ≥ 2 x ≥ 2
l 1 . Z = x ≥ 0 ∧ y ≥ 0 ∧ x -y ≤ 0 P = co({( 0 0 )}) ⊕ cone({( 0 0 ) , -∞ 0 }) D =    0 0 0 ∞ 0 0 ∞ ∞ 0    l 2 . Z = x ≥ 0 ∧ y ≥ 0 ∧ y -x ≤ 0 P = co({( 0 0 )}) ⊕ cone({( 0 0 ) , 0 -∞ }) D =    0 0 0 ∞ 0 ∞ ∞ 0 0   
Continuing the reachability algorithm we calculate the state space going from l 1 and l 2 to l 3 . At first these are stored as two separate states, containing the zones seen in Figure 4 on the next page.

l 3 via l 1 . Z x ≥ 0 ∧ y ≥ 2 ∧ x -y ≤ 0 P = co({( 0 2 )} ⊕ cone({( 0 0 ) , -∞ 0 }) D =    0 0 -2 ∞ 0 0 ∞ ∞ 0    l 3 via l 2 . Z = x ≥ 2 ∧ y ≥ 0 ∧ y -x ≤ 0 P = co({( 2 0 )}) ⊕ cone({( 0 0 ) , 0 -∞ }) D =    0 -2 0 ∞ 0 ∞ ∞ 0 0   
Notice that P and D are unchanged by a delay operation, hence we are ready to take a transition out of l 3 . This allows us to see the difference between the DBM and max-plus approaches. We have two different states for the location, l 3 . Moving on from here, because neither max-plus polyhedra nor DBMs can do exact union, we must decide whether to do all subsequent operations on both instances of the state space, risking state space explosion, or we make an overapproximating union. Opting for the overapproximation, we get the following state for l 3 . Figure 5 shows the overapproximating union both for DBM and maxplus polyhedra. As we can see the DBM overapproximation includes the entire state space, which is not exact, and therefore introduces the risk of false positives throughout the remainder of the analysis. However, for the max-plus polyhedron, the union in this case is actually exact, thus there is no risk of false positives in this case.

l 3 union. P = co({( 2 0 ) , ( 0 2 )}) ⊕ cone({ -∞ 0 , 0 -∞ }) D =    0 0 0 ∞ 0 ∞ ∞ ∞ 0    (
Note that it is not necessarily always the case that exact union and max-plus union are the same.

Algorithms on Max-plus Polyhedra

With inspiration from [START_REF] Bengtsson | and States in Timed Systems[END_REF], this section elaborates on the different algorithms needed for forward and backward reachability analysis. The algorithms accommodate the checking of properties of max-plus polyhedra and the different transforming operations.

Whereas [START_REF] Bengtsson | and States in Timed Systems[END_REF] consider DBMs, we will work with and suggest equivalent algorithms for max-plus polyhedra represented as the Minkowski sum of generators. Using the internal representation, a polyhedron P is given by a convex set V and a linear set W of generators. Additionally note that n always denotes the number of clocks for a system, and p the number of generators for a max-plus polyhedron.

Conversion algorithms

We will sometimes need to convert forth and back between polyhedra expressed as Minkowski sums of a convex set and a cone, and the homogeneous coordinates described in Section 3.2 on page 10. These conversions are not difficult and can be seen in Algorithms 2 and 3; their complexities are O(p).

Algorithm 2: poly-to-cone(P )

G := ∅ for all v ∈ V do G := G ∪ {(v, 0)} end for for all w ∈ W do G := G ∪ {(w, -∞)} end for return G Algorithm 3: cone-to-poly(G) V, W := ∅ for all g ∈ G do if g n = -∞ then W := W ∪ {g} else for i := 1 to n do g i := g i + -g n end for V := V ∪ {g} end if end for return V, W

Property checking

The algorithms in this section do not alter a polyhedron, but are used in reachability analysis to determine whether a given state is reached, as well as determining if some state has previously been visited.

Emptiness testconsistent(P )

To check for consistency of a polyhedron is to check whether it contains a legal clock valuation. By the definition of max-plus polyhedra, a consistent polyhedron is one that contains at least one generator in the convex set. This is because the set of scalars for the convex set must sum to max-plus one, however if there are no generators, there will be no scalars and the sum of nothing is max-plus zero [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF].

Additionally, clocks will always have a value in positive Euclidean space. Hence, for the polyhedron to be consistent, there must exist at least one combination of scalars such that all dimensions are positive or 0.

However, if we start with a polyhedron in positive space and only use the operations described here (for reset, we additionally require it to be used with positive constants), this will preserve the polyhedron in positive space. This makes the test for emptiness O(1) as we only need to check if we have a non-empty set of convex generators according to the definition.

Membership testcontains-point(G, x)

Testing whether a point x is contained in a polyhedron P is not directly used in reachability analysis, but is used as a subroutine in contains and cleanup algorithms.

Checking whether one point x can be generated by a polyhedron P is done by converting P to a cone G of | n+1 max as shown in Algorithm 2. Similarly, x is converted to x ′ of | n+1 max with the last coordinate set appropriate according to whether it is an actual point or a ray representative.

Then the only thing left to do is to see if Gy = x ′ admits a solution, where G is the matrix containing all generators of G as columns. The following algorithm is provided in [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF]; see also [START_REF] Butkovič | Generators, extremals and bases of max cones[END_REF]: The equation Gy = x ′ may not have a solution, but the inequality Gy ≤ x ′ always does. We can compute the maximal solution y of this inequality according to the formula y i = min 1≤j≤n+1 (x ′ j -G ji ). The equation Gy = x ′ then has a solution if and only if G y = x ′ . If so, then x ′ is in P , otherwise it is not. The algorithm runs in O(pn) time.

Algorithm 4: contains-point(G, x)

for all g i ∈ G do y i := min 1≤j≤n+1 (x jg i j ) end for for all 1 ≤ j ≤ n + 1 do z j := max g i ∈G (y i + g i j ) end for if x = z then return true end if return false

Subset testcontains(P, P ′ )

Inclusion checking is a crucial operation when doing state space exploration. It is needed to determine whether a given state has already been visited, and hence does not need to be traversed again. To do this for maxplus polyhedra is simple: given two polyhedra P and P ′ , to determine if P contains P ′ , we just need to check whether all generators of P ′ can be generated by P .

The complexity of the contains algorithm is O(pp ′ n), where n is the number of clocks, p is the number of generators for P and p ′ the number of generators for P ′ , as every call to contains-point takes O(pn), and there are p ′ of them. For comparison, the DBM algorithm for inclusion checking is quadratic in the number of clocks, O(n 2 ).

Algorithm 5: contains(P, P ′ )

G := poly-to-cone(P ) G ′ := poly-to-cone(P ′ ) for all g ′ ∈ G ′ do if ¬contains-point(g ′ , G) then return false end if end for return true

Constraint satisfactionsatisfied(P, ϕ)

It is important to be able to check, non-destructively, whether a polyhedron partially satisfies a given constraint. This is used to determine whether a state satisfies some timing constraints given in the initial query.

For constraints of the form ϕ = x i -x j ∼ c where ∼ ∈ {≤, =, ≥} and x i ∈ C, x j ∈ C ∪{0}, a simple satisfied algorithm is to intersect the polyhedron with the constraint and then check for emptiness. The correctness of this is trivial and will not need more discussion. Complexity-wise, however, the algorithm is O(p 2 n) due to the use of constraint intersection, described in Section 5.3.1. We believe this is not the fastest approach, but improving this is left for future work.

The DBM algorithm for constraint satisfaction uses the same approach, namely adding the constraint to the zone and checking for consistency, which is O(n 2 ).

Transformations

The transformations are the algorithms which modify a polyhedron in order to determine the reachable state space for a TA. This includes the basic algorithms of delay and reset among others.

Constraint intersectionand(P,

x i -x j ∼ c)
Adding a constraint is done by intersecting the polyhedron P with the difference constraint x ix j ∼ c. An algorithm (Algorithm 7 on the next page) for computing the intersection of a max-plus cone with the half-space satisfying a set of max-plus inequalities is given by [START_REF] Xavier Allamigeon | Computing the extreme points of tropical polyhedra[END_REF]. As any difference constraint can be expressed as max-plus inequality, the only thing that remains is to convert the constraint of the form x ix j ∼ c to two vectors a and b that represent the max-plus half-space {x | a ⊗ x ≤ b ⊗ x}, and use the aforementioned algorithm on the cone representation of the polyhedron.

Without loss of generality, we can assume that the constraint is x i -x j ≤ c, because x ix j ≥ c is equivalent to x jx i ≤ -c and intersection with x ix j = c can be computed by simply computing the intersection with x ix j ≤ c and x ix j ≥ c. This constraint can be rewritten as x i ≤ x j + c, or in the max-plus notation, 0 ⊗x i ≤ c ⊗ x j . This means that a will be a vector the components of which will be -∞, except for the i'th value, which will be 0. Similarly, b will consist of negative infinities except at the j'th place, which will be c.

G ≤ := {g ∈ G | a ⊗ g ≤ b ⊗ g} G > := {g ∈ G | a ⊗ g > b ⊗ g} H := G ≤ for all (g, h) ∈ G ≤ × G > do H := H ∪ {((a ⊗ h) ⊗ g) ⊕ ((b ⊗ g) ⊗ h)} end for return H
The complexity of intersect-halfspace is O(p 2 n) -the evaluation of the expression inside the loop takes O(n) and may be performed O(p 2 ) times, which is also the upper bound on the number of computed generators. As the conversions from and to the cone representation can be naively implemented in O(pn), the overall time complexity of intersection with a constraint is in O(p 2 n). The complexity of this operation for DBMs is O(n 2 ).

Delayup(P )

The delay operation is one of the basic algorithms used for forward exploration. Delay is easily done on a max-plus polyhedron, simply by copying all points in the convex combination to the linear combination [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF]. The complexity is O(pn), as we need to loop though all points and for each point we need to copy its value in every dimension. For DBMs the delay algorithm is linear in the number of clocks, O(n).

Algorithm 8: up(P )

W := W ∪ V

Backward delaydown(P )

Backward delay is the algorithm for determining all the states that could have brought us into a given state by delay. It is used when doing backward state-space exploration rather than forward exploration. To do backward delay on a polyhedron we need to add the generator (-1, . . . , -1) to the convex set and then intersect with the polyhedron for positive Euclidean space. Since and is complexity O(p 2 n), the entire complexity of down is O(p 2 n 2 ). Backward delay for DBMs is O(n 2 ).

Algorithm 9: down(P )

V := V ∪ {(-1, . . . , -1)} for i := 1 to n do and(P,

x i ≥ 0) end for Proof of correctness. Let P = co(V )⊕cone(W ) ⊆ | n
≥0 be a polyhedron. The polyhedron that we would like to get after applying this algorithm is

↓ P = {v | v ∈ | n ≥0 , v + d ∈ P, d ∈ | ≥0 }. Let P = (co(V ∪ {f }) ⊕ cone(W )) ∩ | n
≥0 , where f is the vector which has -1 as all components. We need to prove that ↓ P = P .

• ↓ P ⊆ P : Let p be a point in ↓ P . This means that there is a point p + d ∈ P for some d ∈ | ≥0 , from which we can derive a possible representation of p if we consider the scalars of p + d as max-plus multiples of d:

p + d = p i=1 α i dv i ⊕ q i=1 β i dw i = p i=1 α i v i ⊕ q i=1 β i w i ⊗ d,
where α i , β i ∈ | max , v i ∈ V , w i ∈ W and p i=1 α i d = 0. We can now show that p ∈ P -indeed:

p = p i=1 α i v i ⊕ q i=1 β i w i = p i=1 α i v i ⊕ 0f ⊕ q i=1 β i w i .
In the last expression, we have ( p i=1 α i ) ⊕ 0 = 0 and as we know that p ∈ | n ≥0 , the addition of 0 ⊗f will not change the resulting point. Therefore p ∈ P .

• P ⊆ ↓ P : Let p ∈ P :

p = p i=1 α i v i ⊕ γf ⊕ q i=1 β i w i , where α i , β i , γ ∈ | max , v i ∈ V , w i ∈ W and ( p i=1 α i ) ⊕ γ = 0. Let d = -max p i=1 α i . Now we can show that p + d ∈ P : p + d = p i=1 α i v i ⊕ q i=1 β i w i ⊗ d = p i=1 α i dv i ⊕ q i=1 β i dw i , since p i=1 α i d = 0.
According to the definition, this also means that p ∈ ↓ P.

Resetting clocksreset(P,

x i = c)
Reset is the operation of resetting one clock to a given value. In effect it is an affine projection to the hyperplane equivalent to a given constant for the given dimension. It is done by iterating through all generators, setting the given dimension to the given reset value for the convex generators, and setting the given dimension to -∞ for the linear generators [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF]. This operation can be done in linear time in the number of generators O(p). The reset operation on DBMs is done in linear time in the number of clocks O(n). Freeing a clock on a polyhedron is done by removing all constraints on that particular clock. It is used in combination with constraint intersection to handle resets when exploring the state-space backward. Performing the free operation on a polyhedron can be done by resetting the polyhedron with respect to the clock being freed, and then adding the generator containing -∞ for all clocks except x i , where the value should be 0, to the set W. Complexity-wise, free is O(n + p) since reset is linear in the number of points, and creating a vector is linear in the number of clocks. The free operation on DBMs has complexity O(n).

∈ P. ∀1 ≤ i ≤ n, d = i. v i = v ′ i }.
Let P be the polyhedron obtained by running our algorithm on P , i.e. by resetting the clock x d to zero and adding a generator with zero at the d'th place and negative infinities at all others to the set of linear generators. Let ✵ d denote this generator. We shall now prove that P * d = P .

• P * d ⊆ P : Let p be a point in P * d . By the definition of P * d , there must be a point q ∈ P R[x d =0] such that it has the same values of all components except the d'th. From the definition of the reset operation, we also know that q d = 0. We can see that q ⊕ (p d ⊗ ✵ d ) is equal to p and is still contained in P , because of the added generator. In other words, q ⊕ (p d ⊗ ✵ d ) = p ∈ P .

• This operation returns the smallest polyhedron that contains both P 1 and P 2 . Such an operation is useful when performing reachability analysis with convex hull overapproximation of the state space. The algorithm which is simply taking the union of the generators of P 1 and P 2 returns the smallest overapproximating polyhedron. As max-plus polyhedra can represent any given DBM, the max-plus convex hull of two given zones can never be bigger than the DBM convex hull. On the other hand, as seen in the example in Section 4.2 on page 12 there exists convex-hull polyhedra which are tighter overapproximations than convex-hull DBMs. For more details cf. [START_REF] Dyhrberg | Computations on zones using max-plus algebra[END_REF].

The time complexity of taking the union of two sets is O(pn). For DBMs, the overapproximation algorithm is O(n 2 ) as it simply involves searching through the matrices, picking the higher value for every constraint from them.

Cleaning up

All transformations on max-plus polyhedra may introduce redundant generators, i.e. generators that are not extreme points and can therefore be expressed as a combination of the extreme points. Because the time complexity of most of the algorithms depends on the number of generators, it is desirable to remove such redundant generators to store only the minimal number of generators.

Removing redundant generatorscleanup(P )

Cleaning up results in a combination of a unique and minimal representation for max-plus polyhedra. The uniqueness part is to be understood as a minimal unique convex part, and a minimal but not necessarily unique linear part. Making the linear part unique would require a slight alteration of cleanup which normalizes the linear vectors in some way. However, we do not believe this will give us any significant advantage, so we have decided not to do this. A cleaned up polyhedron can be considered to be canonical in the sense that a cleaned up polyhedron is the optimal input for the different algorithms in a similar way as a canonical DBM is for the DBM algorithms. Therefore the clean algorithm can be compared to the canonicalization algorithm for DBMs, which is the standard Floyd-Warshall algorithm and therefore O(n 3 ).

Summary

Table 1 on the next page shows a comparison of the complexity of the algorithms for max-plus polyhedra and for DBMs. The number of clocks is denoted n, while p denotes the number of generators of the polyhedron. Noting that in general, p will be greater than n (especially when we are using max-plus polyhedra to represent complex shapes), we may conclude from the comparison that the complexity of our polyhedra-based algorithms is comparable to the one of the standard DBM algorithms, with a slight advantage to the DBM side.

The main argument in favor of the max-plus approach is thus the fact that overapproximation using max-plus polyhedra is tighter than the DBM overapproximations; max-plus polyhedra allow more complex shapes to be represented by a single instance than DBMs, which can indeed be a major advantage in real-time model checking.

Table 1 notes.

a. This does require the zone to always be canonical and that every other operation checks whether the upper bound is lower than the corresponding lower bound whenever a bound is changed, setting D 00 to a negative value if this is the case. b. We believe that this can be done in O(p). However, the verification of this is left as future work, see Section 6.1.

Max-plus polyhedra DBM

Emptiness test O(1) O(1) a Inclusion test O(p 2 n) O(n 2 ) Constraint satisfaction O(p 2 n) b O(n 2 ) Constraint intersection O(p 2 n) O(n 2 ) Delay O(pn) O(n) Backward delay O(p 2 n 2 ) O(n 2 ) Resetting clocks O(p) O(n) Removing constraints O(p + n) O(n) Union overapproximation O(pn) O(n 2 ) Removing redundant generators O(p 2 n) O(n 3 ) c
Table 1: Complexity of algorithms for max-plus polyhedra and DBMs.

c. As mentioned in Section 5.4.1, we have decided to compare this to DBM canonicalization. Note that some of the DBM algorithms accommodate the possibility of a slight altering which preserves canonicity, allowing canonicalization to be skipped.

Conclusion and future work

We have shown that convex max-plus polyhedra are a suitable data structure for real-time reachability analysis. We have developed the necessary algorithms for both forward and backward reachability analysis, and we have seen that their complexity is comparable to the one of the standard DBM algorithms.

We have also shown that overapproximations using max-plus polyhedra are more tight than the ones one gets in the DBM-based approach, in that max-plus polyhedra can express some non-convex properties without resorting to disjunction. As any disjunction essentially doubles the state space to be checked, hence leads to state-space explosion, we believe that max-plus polyhedra can provide a useful tool for combating state-space explosion.

Performance

To show that our max-plus approach works in practice, we have created a proof-of-concept implementation on top of the Python-based model checker opaal [START_REF] Chr | opaal[END_REF]. This works as expected, but is very slow compared to stateof-the-art tools as e.g. Uppaal. The main reason for this is that neither our implementation nor the opaal model checker itself are optimized for performance, but also some of our algorithms use a rather naïve approach and can certainly be optimized. One example of such optimization is the satisfied algorithm from Section 5.2.4: instead of using constraint intersection, one could compare the constraint against the generators one by one, which would cut complexity down to O(p).

We are also working on an optimized implementation of a max-plus library which can be used by real-time model checkers such as Uppaal as a drop-in replacement for their DBM library. As we are replacing the basic data structures however, this is not an easy task.

Strict constraints

As we are using them here, max-plus polyhedra can only represent sets of non-strict constraints. To extend our approach to strict constraints poses some difficulties which are not present in the DBM approach, essentially because DBMs represent zones by their codimension-one faces, whereas we represent max-plus polyhedra by their corners (generators).

We currently see no way to represent strict constraints using the internal representation of max-plus polyhedra, other than carrying around extra information as to strictness of their codimension-one faces, which so to speak is a step in the wrong direction. We know of some approaches to solve the same problem for classical polyhedra in [START_REF] Bagnara | Possibly not closed convex polyhedra and the parma polyhedra library[END_REF], but they do not seem to carry over to the max-plus setting.

Federation data structure

The forward reachability algorithm described in Section 2.4 on page 6 needs to check whether newly discovered symbolic states have already been processedpassed. The set of passed states can be implemented as a list of such states, but this is very inefficient. It has been shown that a more efficient data structure, such as clock decision diagrams [START_REF] Kim | Clock difference diagrams[END_REF] in the zonebased case, provides a considerable reduction of the space needed to store passed states.

This data structure (also referred to as federation) needs to represent finite union of symbolic states, i.e. max-plus polyhedra in our case, and support following two operations efficiently:

1. Checking whether a polyhedron P in internal representation is a subset of the federation. 2. Adding a new polyhedron P in internal representation to the federation.

Therefore a new BDD-based data structure to be used for federations of max-plus polyhedra is called for.

Extrapolation operation

In general, the symbolic state space explored during reachability checking may be infinite. One solution is to apply an extrapolation operation to each computed successor, which will make the state space finite while keeping the properties of the original one. Several such operations exist for zones [START_REF] Behrmann | Lower and upper bounds in zone-based abstractions of timed automata[END_REF][START_REF] Bouyer | Diagonal constraints in timed automata: Forward analysis of timed systems[END_REF], but as with strict constraints above, the problem is that they are manipulating the codimension-one faces of zones, whereas we need algorithms which work on the corners (generators). Hence new algorithms for extrapolation are necessary.

Algorithm 1 : 2 :

 12 Forward reachability 1: Waiting := {(l 0 , Z 0 )} Passed := ∅ 3: while Waiting = ∅ do 4:

5 :

 5 if l = s and Z ∩ ϕ = ∅ then

  (a) 2D polyhedron. (b) 3D polyhedron.

Figure 1 :

 1 Figure 1: Examples of polyhedra in 2D and 3D.

  3 on page 5, zones are defined by conjunctions of constraints. With a little rewriting these zone constraints can be converted into difference constraints on the form xy n, where x, y ∈ C ∪ {0}, ∈ {≤, <} and n ∈ ) where C is the set of clocks. This is exactly what a DBM represents, hence one DBM encodes exactly one zone. Since we are only interested in the tightest constraints and every constraint is concerned with two clocks, every zone is defined by at most |C 0 | • (|C 0 | -1) constraints, where C 0 = C ∪ {0}. By defining the upper bound on the difference between two clocks as xy n and the lower bound as yx -n, zones in systems with |C| clocks can be stored as |C 0 | × |C 0 | matrices.

Figure 2 :

 2 Figure 2: Timed automaton used for this example.

Figure 3 :

 3 Figure 3: Location l 1 and l 2 after delay.

  (a) l3 via l1 (b) l3 via l2

Figure 4 :

 4 Figure 4: The two different possible states in location l 3 .

  a) l3 DBM union (b) l3 MPP and exact union

Figure 5 :

 5 Figure 5: Union of zones from the two l 3 states.

Algorithm 6 :

 6 and(P, x ix j ≤ c) a := (-∞, . . . , -∞) b := (-∞, . . . , -∞) a i := 0 b j := c return cone-to-poly(intersect-halfspace(poly-to-cone(P ), a, b)) Algorithm 7: intersect-halfspace(G, a, b)

Algorithm 10 :

 10 reset(P, x i = c) for all v ∈ V do v i := c end for for all w ∈ W do w i := -∞ end for 5.3.5. Removing constraintsfree(P, x i )

Algorithm 11 :

 11 free(P, x i ) reset(P,x i = 0) g := (-∞, . . . , -∞) g i := 0 W := W ∪ {g}Proof of correctness. Let P be a polyhedron. Removing a constraint on clock x d in P should result in the set P * d = {v | ∃v ′

Algorithm 12 :

 12 cleanup(P ) G := poly-to-cone(P ) for all g ∈ G do if contains-point(g, G {g}) then G := G {g} end if end for return cone-to-poly(G) This is done by checking for each point whether it can be generated by the other points. If it does, it is removed. Every such check costs O(pn), which makes the time complexity of this algorithm O(p 2 n).

  P ⊆ P * d : Let p ∈ P . Because the reset and addition of ✵ d only affected the d'th coordinate, we know that there is a point p ′ ∈ P that is equal to p except for the d'th coordinate. Because the d'th coordinate was first set to zero and then possibly incremented by some scalar multiple of ✵ d , it cannot be negative and therefore satisfies the second condition of membership in P * d . Thus, p ∈ P * d . 5.3.6. Union overapproximationconvex-union(P 1 , P 2 )

They are called tropical in honour of Imre Simon, who pioneered the field, apparently because he was from Brazil.