Fast algorithms for handling diagonal constraints in timed automata

A popular method for solving reachability in timed automata proceeds by enumerating reachable sets of valuations represented as zones. A na\"ive enumeration of zones does not terminate. Various termination mechanisms have been studied over the years. Coming up with efficient termination mechanisms has been remarkably more challenging when the automaton has diagonal constraints in guards. In this paper, we propose a new termination mechanism for timed automata with diagonal constraints based on a new simulation relation between zones. Experiments with an implementation of this simulation show significant gains over existing methods.


Introduction
Timed automata have emerged as a popular model for systems with real-time constraints [2]. Timed automata are finite automata extended with real-valued variables called clocks. All clocks are assumed to start at 0, and increase at the same rate. Transitions of the automaton can make use of these clocks to disallow behaviours which violate timing constraints. This is achieved by making use of guards which are constraints of the form x ≤ 5, x − y ≥ 3, y > 7, etc where x, y are clocks. A transition guarded by x ≤ 5 says that it can be fired only when the value of clock x is ≤ 5. Another important feature is the reset of clocks in transitions. Each transition can specify a subset of clocks whose values become 0 once the transition is fired. The combination of guards and resets allows to track timing distance between events. A basic question that forms the core of timed automata technology is reachability: given a timed automaton, does there exist an execution from its initial state to a final state. This question is known to be decidable [2]. Various algorithms for this problem have been studied over the years and have been implemented in tools [24,5,23,27,26].
Since the clocks are real valued variables, the space of configurations of a timed automaton (consisting of a state and a valuation of the clocks) is infinite and an explicit enumeration is not possible. The earliest solution to reachability was to partition this space into a finite number of regions and build a region graph that provides a finite abstraction of the behaviour of the timed automaton [2]. However, this solution was not practical. Subsequent works introduced the use of zones [13]. Zones are special sets of clock valuations with efficient data structures and manipulation algorithms [5]. Within zone based algorithms, there is a division: forward analysis versus backward analysis. The current industry strength tool UPPAAL [24] implements a forward analysis approach, as this works better in the presence of other discrete data structures used in UPPAAL models [8]. We focus on this forward analysis approach using zones in this paper.
The forward analysis of a timed automaton essentially enumerates sets of reachable configurations stored as zones. Some extra care needs to be taken for this enumeration to terminate. Traditional development of timed automata made use of extrapolation operators over zones to ensure termination. These are functions which map a zone to a bigger zone. Importantly, the range of these functions is finite. The goal was to come up with extrapolation operators which are sound: adding these extra valuations should not lead to new behaviours. This is where the role of simulations between configurations was studied and extrapolation operators based on such simulations were devised [13]. A certain extrapolation operation, which is now known as Extra M [4] was proposed and reachability using Extra M was implemented in tools [13].
A seminal paper by Bouyer [8] revealed that Extra M is not correct in the presence of diagonal constraints in guards. These are constraints of the form x − y ⊳ c where ⊳ is either < or ≤, and c is an integer. Moreover, it was proved that no such extrapolation operation would be correct when there are diagonal constraints present. It was shown that for automata without diagonal constraints (henceforth referred to as diagonal-free automata), the extrapolation works. After this result, developments in timed automata reachability focussed on the class of diagonal-free automata [4,3,20,21], and diagonal constraints were mostly sidelined. All these developments have led to quite efficient algorithms for diagonalfree timed automata.
Diagonal constraints are a useful modeling feature and occur naturally in certain problems, especially scheduling [16] and logic-automata translations [15,22]. It is however known that they do not add any expressive power: every timed automaton can be converted into a diagonal-free timed automaton [6]. This conversion suffers from an exponential blowup, which was later shown to be unavoidable: diagonal constraints could potentially give exponentially more succinct models [9]. Therefore, a good forward analysis algorithm that works directly on a timed automaton with diagonal constraints would be handy. This is the subject of this paper.
Related work. The first attempt at such an algorithm was to split the (extrapolated) zones with respect to the diagonal constraints present in the automaton [5]. This gave a correct procedure, but since zones are split, an enumeration starts from each small zone leading to an exponential blow-up in the number of visited zones. A second attempt was to do a more refined conversion into a diagonal free automaton by detecting "relevant" diagonals [12,25] in an iterative manner. In order to do this, special data structures storing sets of sets of diagonal constraints were utilized. A more recent attempt [17] extends the works [4] and [20] on diagonal-free automata to the case of diagonal constraints. All the approaches suffer from either a space or time bottleneck and are incomparable to the efficiency and scalability of tools for diagonal-free automata.
Our contributions. The goal of this paper is to come up with fast algorithms for handling diagonal constraints. Since the extrapolation based approach is a dead end, we work with simulation between zones directly, as in [20] and [17]. We propose a new simulation relation between zones that is correct in the presence of diagonal constraints (Section 3). We give an algorithm to test this simulation between zones (Section 4). We have incorporated this simulation test in the tool TChecker [18] that checks reachability for timed automata, and compared our results with the state-of-the-art tool UPPAAL. Experiments show an encouraging gain, both in the number of zones enumerated and in the time taken by the algorithm, sometimes upto four orders of magnitude (Section 6). The main advantage of our approach is that it does not split zones, and furthermore it leverages the optimizations studied for diagonal-free automata.
From a technical point of view, our presentation does not make use of regions and instead works with valuations, zones and simulation relations. We think that this presentation provides a clearer perspective -as a justification of this claim, we extend our simulation to timed automata with general updates of the form x := c and x := y + d in transitions (where x, y are clocks and c, d are constants) in a rather natural manner (Section 5). In general, reachability for timed automata with updates is undecidable [11]. Some decidable cases have been proposed for which the algorithms are based on regions. For decidable subclasses containing diagonal constraints, no zone based approach has been studied. Our proposed method includes these classes, and also benefits from zones and standard optimizations studied for diagonal-free automata.

Preliminaries
Let N be the set of natural numbers, R ≥0 the set of non-negative reals and Z the set of integers. Let X be a finite set of variables ranging over R ≥0 , called clocks. Let Φ(X) denote the set of constraints ϕ formed using the following grammar: Constraints of the form x ⊳ c and c ⊳ x are called non-diagonal constraints and those of the form x − y ⊳ c are called diagonal constraints. We have adopted a convention that in non-diagonal constraints x ⊳ c and c ⊳ x, the constant c is restricted to N. A clock valuation v is a function which maps every clock x ∈ X to a real number v(x) ∈ R ≥0 . A valuation is said to satisfy a guard g, written as v |= g if replacing every x in g with v(x) makes the constraint g true. For δ ∈ R ≥0 we write v + δ for the valuation which maps every x to v(x) + δ. Given a subset of clocks R ⊆ X, we write [R]v for the valuation which maps each x ∈ R to 0 and each x ∈ R to v(x). A timed automaton A is a tuple (Q, X, q 0 , T, F ) where Q is a finite set of states, X is a finite set of clocks, q 0 ∈ Q is the initial state, F ⊆ Q is a set of accepting states and T ∈ Q × Φ(X) × 2 X × Q is a set of transitions. Each transition t ∈ T is of the form (q, g, R, q ′ ) where q and q ′ are respectively the source and target states, g is a constraint called the guard, and R is a set of clocks which are reset in t. We call a timed automaton diagonal-free if guards in transitions do not use diagonal constraints.
A configuration of A is a pair (q, v) where q ∈ Q and v is a valuation. The semantics of a timed automaton is given by a transition system S A whose states are the configurations of A. Transitions in S A are of two kinds: delay transitions are given by (q, v) δ − → (q, v + δ) for all δ ≥ 0, and action transitions are given by (q, v) −→ for a sequence of delay δ followed by action t. A run of A is an alternating sequence of delay-action transitions starting from the initial state q 0 and the initial valuation 0 which maps every clock to 0: (q 0 , 0) A run of the above form is said to be accepting if the last state q n ∈ F . The reachability problem for timed automata is the following: given an automaton A, decide if there exists an accepting run. This problem is known to be PSPACE-complete [2]. Since the semantics S A is infinite, solutions to the reachability problem work with a finite abstraction of S A that is sound and complete. Before we explain one of the popular solutions to reachability, we state a result which allows to convert every timed automaton into a diagonal-free timed automaton.
◮ Theorem 1. [6] For every timed automaton A, there exists a diagonal-free timed automaton A df s.t. there is a bijection between runs of A and A df . The number of states in A df is 2 d · n where d is the number of diagonal constraints and n is the number of states of A.
Proof. (Sketch) We give the construction to remove one diagonal. This can be applied iteratively to remove all diagonal constraints. Let x−y ⊳ c be a diagonal constraint appearing in A. Initially, the value of x − y is 0 and hence x − y ⊳ c is true iff 0 ⊳ c. The value x − y does not change by time elapse. It can change only due to resets. Therefore, the idea is to use the states of the automaton to maintain if the constraint is true or false. Let A ′ be the new automaton which does not contain x − y ⊳ c.
The states of A ′ are of the form Q × {0, 1}. For every transition (q, g, R, q ′ ) of A we have the following transitions in and ≥ otherwise. The guards −y ⊳ c and −y ⊲ 1 c can be suitably removed if they are either trivially true or false, and if not, rewritten to the form y ∼ d where d ≥ 0. if x ∈ R and y ∈ R, the construction is similar to the above case.
The initial state is (q 0 , 0 ⊳ c) and the accepting states are F × {0, 1}. Note that there is a one-one correspondence between runs of A and A ′ and moreover A ′ does not have the diagonal constraint x − y ⊳ c. ◭ The above theorem allows to solve the reachability of a timed automaton A by first converting it into the diagonal free automaton A df and then checking reachability on A df . However, this conversion comes with a systematic exponential blowup (in terms of the number of diagonal constraints present in A). It was shown in [9] that such a blowup is unavoidable in general. We will now recall the general algorithm for analyzing timed automata, and then move into specific details which depend on whether the automaton has diagonal constraints or not.

Zones and simulations
Fix a timed automaton A with clock set X for the rest of the discussion in this section. As the space of valuations of A is infinite, algorithms work with sets of valuations called zones. A zone is set of clock valuations given by a conjunction of constraints of the form x − y ⊳ c, x ⊳ c and c ⊳ x where c ∈ Z and ⊳ ∈ {<, ≤}, for example the solutions of x − y < 5 ∧ y ≤ 10 is a zone. The transition relation over configurations (q, v) is extended to (q, Z) where Z is a zone. We define the following operations on zones given a guard g and a set of clocks R: time It can be shown that all these operations result in zones. Zones can be efficiently represented and manipulated using Difference Bound Matrices (DBMs) [14].
The zone graph ZG(A) of timed automaton A is a transition system whose nodes are of the form (q, Z) where q is a state of A and Z is a zone. For each transition t := (q, g, R, q ′ ) of A, and each zone (q, Z) there is a transition (q, The initial node is (q 0 , Z 0 ) where q 0 is the initial state of A and Z 0 = {0 + δ | δ ≥ 0} is the zone obtained by elapsing an arbitrary delay from the initial valuation. A path in the zone graph is a sequence (q 0 , Z 0 ) ⇒ t0 (q 1 , Z 1 ) ⇒ t1 · · · ⇒ tn−1 (q n , Z n ) starting from the initial node. The path is said to be accepting if q n is an accepting state. The zone graph is known to be sound and complete for reachability.
◮ Theorem 2. [13] A has an accepting run iff ZG(A) has an accepting path.
This does not yet give an algorithm as the zone graph ZG(A) is still not finite. Moreover, there are examples of automata for which the reachable part of ZG(A) is also infinite: starting from the initial node, applying the successor computation leads to infinitely many zones. Two different approaches have been studied to get finiteness, both of them based on the usage of simulation relations. A (time-abstract) simulation relation between configurations of A is a reflexive and transitive relation such that for all states q. The simulation relation can be extended to zones: The simulation relation is said to be finite if the function mapping zones Z to the down sets ↓Z has finite range. We now recall a specific simulation relation which was shown to be finite and correct for diagonal-free timed automata. Current algorithms and tools for diagonal-free automata are based on this simulation. The conditions required for v LU v ′ ensure that when all lower bound constraints c ⊳ x satisfy c ≤ L(x) and all upper bound constraints x ⊳ c satisfy c ≤ U (x), whenever v satisfies a constraint, v ′ will also satisfy it.
◮ Definition 3 (LU-bounds and the relation LU [4,20]). An LU -bounds function is a pair of functions L : X → N ∪ {−∞} and U : X → N ∪ {−∞} that map each clock to either a non-negative constant or −∞. Given an LU -bounds function, The LU simulation is extended to zones as mentioned above. Let us write

Reachability in diagonal-free timed automata
A natural method to get finiteness of the zone graph is to prune the zone graph computation through simulations Z LU Z ′ : do not explore a node (q, Z) if there is an already visited node (q, Z ′ ) such that Z LU Z ′ . Since these simulation tests need to be done often during the zone graph computation, an efficient algorithm for performing this test is crucial. Note that Z LU Z ′ iff Z ⊆ ↓ LU Z ′ . However, it is known that the set ↓ LU Z ′ is not necessarily a zone [4], and hence no simple zone inclusions are applicable. The first algorithms for timed automata followed a different approach, which we call the extrapolation approach. A function Extra + LU [4] on zones was defined, which had two nice properties -(1) Extra + LU (Z) ⊆ ↓ LU Z and (2) Extra + LU (Z) is a zone for all Z. These properties give rise to an algorithm that performs only efficient zone operations: successor computations and zone inclusions.
Reachability algorithm using zone extrapolation. The input to the algorithm is a timed automaton A. The algorithm maintains two lists, Passed and Waiting. Initially, the node (q 0 , Extra + LU (Z 0 )) is added to the Waiting list (recall that (q 0 , Z 0 ) is the initial node of the zone graph ZG(A)). Wlog. we assume that q 0 is not accepting. The algorithm repeatedly performs the following steps: Step 1. If Waiting is empty, then return "A has no accepting run"; else pick (and remove) a node (q, Z) from Waiting. Add (q, Z) to Passed.
The correctness of this algorithm follows from Theorem 4 and the fact that Extra + LU (Z) ⊆ ↓ LU Z [4]. Intuitively, the additional valuations on top of Z in Extra + LU (Z) belong to ↓ LU Z, and hence are simulated by some valuation in Z. This therefore entails that there are no new states visited due to these extra valuations. More recently, an O(|X| 2 ) algorithm for Z LU Z ′ was proposed [20]. This algorithm has the same complexity as the test for inclusion between zones. Additionally, the LU test was a simple extension of the zone inclusion algorithm, to account for the LU bounds. This gives an algorithm that can avoid explicit extrapolations, and also use a better set ↓ LU Z ′ to prune the search.
Reachability algorithm using simulations. The initial node (q 0 , Z 0 ) is added to the Waiting list. Wlog. we assume that q 0 is not accepting. The algorithm repeatedly performs the following steps: Step 1. If Waiting is empty, then return "A has no accepting run"; else pick (and remove) a node (q, Z) from Waiting. Add (q, Z) to Passed.
Step 2. For each transition t := (q, g, R, q 1 ), compute the successor (q, Z) ⇒ t (q 1 , Z 1 ): if Z 1 = ∅ perform the following operations -if q 1 is accepting, return "A has an accepting run"; else check if there exists a node (q 1 , Z ′ 1 ) in Passed or Waiting such that Z 1 LU Z ′ 1 : if yes, ignore the node (q 1 , Z 1 ), otherwise add (q 1 , Z 1 ) to Waiting.

Reachability in the presence of diagonal constraints
The LU relation is no longer a simulation when diagonal constraints are present. Moreover, it was shown in [8] that no extrapolation operator (along the lines of Extra + LU ) can work in the presence of diagonal constraints. The first option to deal with diagonals is to use Theorem 1 to get a diagonal free automaton and then apply the methods discussed previously. One problem with this is the systematic exponential blowup introduced in the number of states of the resulting automaton. Another problem is to get diagnostic information: counterexamples need to be translated back to the original automaton [5]. Various methods have been studied to circumvent the diagonal free conversion and instead work on the automaton with diagonal constraints directly. We recall the approach used in the state-of-the-art tool UPPAAL below.
Zone splitting [5]. The paper introducing timed automata gave a notion of equivalence between valuations v ≃ M v ′ parameterized by a function M mapping each clock x to the maximum constant M among the guards of the automaton that involve x. This equivalence is a finite simulation for diagonal-free automata. Equivalence classes of ≃ M are called regions. This was extended to the diagonal case by [5] M relation splits the regions further, such that each region is either entirely included inside g, or entirely outside g for each g. The next step is to use this notion of equivalence in zones. The paper [5] follows the extrapolation approach: to each zone Z, an extrapolation operation Extra M (Z) is applied; this adds some valuations which are ≃ M equivalent to valuations in Z; then it is further split into multiple zones, so that each small zone is either inside g or outside g for each diagonal constraint g. If d is the number of diagonal constraints present in the automaton, this splitting process can give rise to 2 d zones for each zone Z. From each small zone, the zone graph computation is started. Essentially, the exponential blow-up at the state level which appeared in the diagonal-free conversion now appears in the zone level.
Finding relevant diagonal constraints [12]. This approach wants to perform the diagonalfree conversion, however restricted to a subset of diagonal constraints. Initially, none of the diagonal constraints are removed and the extrapolated zone graph (using Extra M ) is computed. Since this graph contains more valuations, this is an over-approximation: if there is a path over these extrapolated zones to an accepting state, it might not be an actual path over concrete zones. Hence, each time an accepting state is reached, a concrete zone computation is done to check if it is spurious. If it is indeed spurious, the path was disabled in between by some guard. The goal is to now find all diagonal constraints along the path that were responsible for this disability: the authors propose to do this by maintaining for each zone Z, and each pair of clocks x, y, a set of diagonal constraints that were used to get the current value of x − y in Z. In order to get smaller sets of diagonal constraints, the authors in fact maintain a set-of-sets so that any set in this set-of-sets can be picked as a candidate. Once this candidate set S is picked, the diagonal free conversion restricted to S is done. This is the refinement. Again, an extrapolated zone graph is computed, and the process is continued till there are no false positives. Since each refinement removes a diagonal constraint, in the end the algorithm will compute the full diagonal-free conversion. This method works well when not many diagonal constraints result in false positives.
LU extended to diagonals. Recently, [17] gave an extension of LU simulation, called d LU , to the case of diagonals. This was parameterized by two constants LU associated to each x − y for x, y ∈ X. An algorithm Z d LU Z ′ was developed. However, the complexity of this test is unsatisfactory:it was shown that deciding Z d LU Z ′ is NP-hard. Although this results in a tremendous decrease in the number of nodes explored, the algorithm performs bad in terms of timing due to the implementation of the simulation test.
In this paper, we propose a new simulation to handle diagonal constraints. This has two advantages -using this avoids the blow-up in the number of nodes arising due to zone splitting, and the simulation test between zones has an efficient implementation and is significantly quicker than the simulation of [17].

A new simulation relation
In this section, we start with a definition of a relation between timed automata configurations, which in some sense "declares" upfront what we need out of a simulation relation that can be utilized in a reachability algorithm. As we proceed, we will make its description more concrete and give an effective simulation algorithm between zones, that can be implemented. Fix a set of clocks X. This generates constraints Φ(X).
Our goal is to utilize the above relation in a simulation (as defined in Page 5) for a timed automaton. Directly from the definition, we get the following lemma which shows that the ⊑ G relation is preserved under time elapse.
The other kind of transformation over valuations is resets. Given sets of guards G 1 , G and a set of clocks R, we want to find conditions on To do this, we need to answer this question: what guarantees should we This motivates the next definition. ◮ Definition 7 (weakest pre-condition of ⊑ G over resets). For a constraint ϕ and a set of clocks R, we define a set of constraints wp(⊑ ϕ , R) as follows: Note that the relation ⊑ G is parameterized by a set of constraints. Additionally, we desire this set to be finite, so that the relation can be used in an algorithm. We need to first link an automaton A with such a set of constraints. One way to do it is to take the set of all guards present in the automaton and to close it under weakest pre-conditions with respect to all possible subsets of clocks. A better approach is to consider a set of constraints for each state, as in [3] where the parameters for extrapolation (the maximum constants appearing in guards) are calculated at each state.
◮ Definition 8 (State based guards). Let A = (Q, X, q 0 , T, F ) be a timed automaton. We associate a set of guards G(q) for each state q ∈ Q, which is the least set of guards (for the coordinate-wise subset inclusion order) such that for every transition (q, g, R, q 1 ): the guard g and the set wp(⊑ G(q 1 ) , R) are present in G(q). More precisely, {G(q)} q∈Q is the least solution to the following set of equations written for each q ∈ Q: All constraints present in the set wp(⊑ G(q 1 ) , R) contain constants which are already present in ⊑ G(q 1 ) . The least solution to the above set of equations can therefore be obtained by a fixed point computation which starts with G(q) set to (q,g,R,q1)∈T {g} and then repeatedly updates the weakest-preconditions. Since no new constants are generated in this process, the fixed point computation terminates. We now have the ingredients to define a simulation relation over configurations of a timed automaton with diagonal constraints.
◮ Definition 9 (A-simulation). Let A = (Q, X, q 0 , T, F ) be a timed automaton and let the set of guards G(q) of Definition 8 be associated to every state q ∈ Q. We define a relation ◮ Lemma 10. The relation A is a simulation on the configurations of timed automaton A.
Proof. Pick two configurations such that (q, v) A (q, v ′ ). We need to show the two conditions required for a simulation as given in Page 5. The first condition follows directly from Lemma 6. For the second condition, pick a transition t := (q, g, R, q 1 ) such that Moreover, the guard g ∈ G(q), by Definitions 9 and 8. This implies that v ′ satisfies the guard g (taking δ = 0 in Definition 5) thereby giving us a transition (q, v ′ ) This follows directly from the definition of the weakest precondition wp(⊑ G(q 1 ) , R). ◭ As pointed before, Definition 5 gives a declarative description of the simulation and it is unclear how to work with it algorithmically, even when the set of constraints G is finite. The main issue is with the ∀δ quantification, which is not finite. We will first provide a characterization that brings out the fact that this ∀δ quantification is irrelevant for diagonal constraints (essentially because value of v(x)−v(y) does not change with time elapse). Given a set of constraints G, let G − ⊆ G be the set of non-diagonal constraints in G.
Proof. The left-to-right implication is direct from Definition 5. For the right-to-left implication, it is sufficient to note that v(x) − v(y) is invariant under time elapse, for every pair of clocks x, y. ◭ It now amounts to solving the ∀δ problem for non-diagonals. It turns out that the LU simulation achieves this, almost. Since the LU bounds cannot distinguish between guards with < and ≤ comparisons, the LU simulation does not characterize v ⊑ G − v ′ completely (converse of Lemma 13 does not hold). Although we are aware of the (rather technical) modifications to LU simulation that are needed for this characterization, we choose to use the existing LU directly as this has already been implemented in tools. This gives us a finer simulation than v ⊑ G − v ′ . Let us first define the LU bounds corresponding to a set of general constraints G.
◮ Definition 12 (LU -bounds from G). Let G be a finite set of constraints. We define LU (G) to denote the pair of functions L G and U G defined as follows: Hence v, v ′ (and thereby v + δ, v ′ + δ) satisfy all constraints of the form c ⊳ x from G, giving the ⊑ G − condition for c ⊳ x guards. Analogous reasoning works for the case when v(x) < v ′ (x), giving the ⊑ G − condition for the remaining guards. ◭ The above observations call for the next definition and subsequent lemmas.
◮ Definition 14 (approximating ⊑ G ). Let G be a finite set of constraints. We define a relation

◮ Lemma 15. The relation LU
A is a finite simulation on the configurations of A.

Proof. The fact that LU
A is a simulation follows from definitions of LU , the state based guards, and the way we have constructed the weakest pre-condition over resets. We now focus on showing that it is finite. To do this, we will show that for each q, the simulation ⊑ LU G is finite. Recall that a simulation is finite, if the function mapping zones Z to down sets ↓Z has finite range. From Theorem 4, simulation LU (G) is finite. Now consider the We will now show that this simulation is also finite. Let G + ⊆ G be the set of diagonal constraints in G; now for each subset S ⊆ G + consider the set R |X| ≥0 ∩ ϕ∈S ϕ ∩ ϕ∈G + \S ¬ϕ. These sets partition the space of all valuations into atmost 2 |G + | classes. The downset ↓Z for a zone with respect to will be the union of these classes that intersect with Z. This shows that is also finite. Since LU (G) and are finite, it follows that ⊑ LU G is also finite. ◭ The ⊑ LU G relation can be extended to an algorithm for simulation between zones, which we do in the next section. With this, we get the following theorem.

4
Algorithm for Z ⊑ G Z ′ Fix a finite set of guards G. Restating the definition of ⊑ G extended to zones: In this section, we will view the characterization of ⊑ G as in Proposition 11 and give an algorithm to check Z ⊑ G Z ′ that uses as an oracle a test Z ⊑ G − Z ′ . Using the discussion in the previous section, we later approximate Z ⊑ G − Z ′ with Z LU (G) Z ′ in the implementation. We start with an observation following from Proposition 11.
Firstly, note that the set of non-diagonals in G and G ′ are the same, since G and G ′ differ only by the diagonal constraint ϕ. Secondly, note that both v |= ϕ and v ′ |= ϕ. These two observations coupled with The last part of the Lemma follows from Proposition 11. ◭ This leads to the following algorithm consisting of two mutually recursive procedures. This algorithm is essentially an implementation of the above lemma, with two optimizations: we start with the non-diagonal check in Line 6 of Algorithm 1 -if this is already violated, then the algorithm returns false; suppose Z ⊑ G − Z ′ , the next task is to perform the checks in the first statement of Lemma 17 -this is done by Algorithm 2; note however that when Algorithm 2 is called, we already have Z ⊑ G − Z ′ , hence Z ∩ ¬ϕ ⊑ G − Z ′ . Therefore we use an optimization in Line 7 by calling Algorithm 2 directly (as the check in Line 6 of Algorithm 1 will be redundant).

Algorithm 2
As mentioned before, we approximate Z ⊑ G − Z ′ with Z LU (G) Z ′ . In this case, each zone operation inside a call can be done in O(|X| 2 ) [28,20].

◮ Theorem 18. When using
Proof. The algorithm makes two recursive calls per diagonal constraint. In each call the operations are either Z LU (G) Z ′ or an intersection of Z with a single constraint. From [20], the check Z LU (G) Z ′ can be done in O(|X| 2 ). Intersection with a single constraint can also be done in O(|X| 2 ) (similar to approach in [28]). This gives the complexity. ◭ From a complexity viewpoint, this algorithm is not efficient since it makes an exponential number of calls in the number of diagonal constraints (in fact this may not be avoidable due to Lemma 22, which follows from the NP-hardness result in [17]). Although the above algorithm does involve many calls, the internal operations involved in each call are simple zone manipulations. Moreover, the preliminary checks (for instance line 6 of Algorithm 1) cut short the number of calls. This is visible in our experiments which are very good, especially with respect to running time, as compared to other methods. A similar hardness was shown for a different simulation in [17], but the implementation there indeed witnessed the hardness, as the time taken by their algorithm was unsatisfactory.

Hardness of checking
In this section we prove that the problem of checking Z ⊑ LU G Z ′ is NP-complete. In order to do this, we will give a polynomial time reduction from the problem of checking Z d LU Z ′ , which was shown to be NP-hard [17], to our problem. Let us first recall the L d U d bounds function (superscripting with d to distinguish from the diagonal free LU bounds) and the simulation relation d LU of [17]. if The following two results hold when Z, Z ′ and the L d U d bounds satisfy certain properties.
◮ Lemma 19. [17] Under certain assumptions on Z, Z ′ and L d U d : We will use this to show ⊑ G is NP-hard. Given L d U d , we construct G to be the set containing: for every x, y ∈ X ∪ {0} the constraint x − y < L d (x − y) and the constraints ◮ Lemma 20. Given L d U d , let G be constructed as above. Then for all Z, Z ′ satisfying the conditions needed for Lemma 19 Then from Lemma 19 we get that there exists an integral valuation that serves as a witness. We choose that integral valuation and call it v. Then there exists some pair of clocks x, y ∈ X ∪ {0} such that for every v ′ ∈ Z ′ , one of the following two cases holds: In the second case, from the construction of G we get that c On the other hand, if Z LU (G) Z ′ there is v ∈ Z such that for every v ′ ∈ Z ′ one of the following two cases happen: . Use similar argument as previous case.
This completes the proof. ◭ Before proceeding to prove NP-completeness, we present a new characterization of the simulation relation ⊑ LU G that will help in proving the upper bound.
Suppose there exists S such that (Z ∩ ϕ∈S ϕ)

Simulations for updateable timed automata
In the timed automata considered so far, clocks are allowed to be reset to 0 along transitions. We consider in this section more sophisticated transformations to clocks in transitions. These are called updates. An update up : R |X| ≥0 → R |X| is a function mapping non-negative |X|dimensional reals (valuations) v to general |X|-dimensional reals (which may apriori not be valuations as the coordinates may be negative). The syntax of the update function up is given by a set of atomic updates up x to each x ∈ X, which are of the form x := c or x := y + d where c ∈ N, d ∈ Z and y ∈ X (possibly equal to x). Note that we want d to be an integer, since we allow for decrementing clocks, and on the other hand c ∈ N since we have non-negative clocks. Given a valuation v and an update up, the valuation up(v) is: Note that in general, due to the presence of updates x := y + d, the update up(v) may not yield a clock valuation. However, when it does give a valuation, it can be used as a transformation in timed automata transitions. We say An updateable timed automaton A = (Q, X, q 0 , T, F ) is an extension of a classic timed automaton with transitions of the form (q, g, up, q ′ ) where up is an update. Semantics extend in the natural way: delay transitions remain the same, and for action transitions We allow the transition only if the update results in a valuation. The reachability problem for these automata is known to be undecidable in general [11]. Various subclasses with decidable reachability have been discussed in the same paper. Decidability proofs in [11] take the following flavour, for a given automaton A: (1) divide the space of all valuations into a finite number of equivalence classes called regions (2) to build the parameters for the equivalence, derive a set of diophantine equations from the guards of A; if they have a solution then construct the quotient graph of the equivalence (called region graph) parameterized by the obtained solution and check reachability on it; if the equations have no solution, output that reachability for A cannot be answered. Sufficient conditions on the nature of the updates that give a solution to the diophantine equations have been tabulated in [11]. When the automaton is diagonal-free, the "region-equivalence" can be used to build an extrapolation operation which in turn can be used in a reachability algorithm with zones. When the automaton contains diagonals, the region-equivalence is used to only build a region graphno effective zone based approach has been studied.
We use a similar idea, but we have two fundamental differences: (1) we want to obtain reachability through the use of simulations on zones, and (2) we build equations over sets of guards as in Definition 8. The advantage of this approach is that this allows the use of coarser simulations over zones. Even for automata with diagonal constraints and updates, we get a zone based algorithm, instead of resorting to regions which are not efficient in practice.
The notion of simulations as in Page 5 remains the same, now using the semantics of transitions with updates. We will re-use the simulation relation ⊑ G . We need to extend Definition 7 to incorporate updates. We do this below. Here is a notation: for an update function up, we write up(x) to be c if up x is x := c, and up(x) to be y + c if up x is x := y + c. Some examples: wp(x ≤ 5, x := x + 10) is empty, since up(x) is x + 10, and the guard x + 10 ≤ 5 is not satisfiable; wp(x ≤ 5, x := x − 10) is x ≤ 15, wp(x ≤ 5, x := c) is empty, wp(x − y ≤ 5, x := z 1 , y := z 2 + 10 ) will be z 1 − (z 2 + 10) ≤ 5, giving the constraint

and is empty otherwise.
◮ Definition 24 (State based guards). Let A = (Q, X, q 0 , T, F ) be an updateable timed automaton. We associate a set of constraints G(q) for each state q ∈ Q, which is the least set of constraints (for the coordinate-wise subset inclusion order) such that for every transition (q, g, up, q 1 ): the guard g and the set wp(⊑ G(q 1 ) , up) are present in G(q), and in addition constraints that allow the update to happen are also present in G. The last condition is given by the weakest precondition of the set of constraints {x ≥ 0 | x ∈ X}. Overall, {G(q)} q∈Q is the least solution to the following set of equations, for each q ∈ Q: The least solution {G(q)} q∈Q is said to be finite if each G(q) is a finite set of constraints.
In contrast to the simple reset case, the above set of equations may not have a finite solution. Consider a self-looping transition: (q, x ⊳ c, x := x − 1, q). We require x ⊳ c ∈ G(q). Now, wp(x ⊳ c, x := x − 1) is x ⊳ c + 1 which should be in G(q) according to the above equation. Continuing this process, we need to add x ⊳ d for every natural number d ≥ c. Indeed this is consistent with the undecidability of reachability when subtraction updates are allowed. We deal with the subject of finite solutions to the above equations later in this section. On the other hand, when the above system does have a solution with finite G(q) at every q, we can use the A simulation of Definition 9 and its approximation LU A to get an algorithm. Proof. Pick two configurations such that (q, v) A (q, v ′ ). We need to show the two conditions required for a simulation as given in Page 5. The proof follows along the lines of Lemma 10, except for one change. Consider a transition t : (q, g, up, q 1 ) and let . Pick a constraint ϕ ∈ G(q 1 ) such that v |= ϕ. Suppose ϕ is of the form x − y ⊳ c. If both up(x) and up(y) contain the same clock z, then in the difference x − y the role of z gets cancelled, and we get that , and furthermore this difference is invariant over time elapse. Hence the ⊑ G(q 1 ) property is true for ϕ. The next case is when up(x) and up(y) do not contain the same clock. Then, up(x) − up(y) ⊳ c is a guard (diagonal or non-diagonal) which according to Definition 23 is present in G. Since v ⊑ G v ′ , we get the ⊑ G(q 1 ) condition for ϕ. Similar reasoning follows when ϕ is a non-diagonal. One however needs to use the fact that since v |= ϕ, if ϕ is of the form x ⊳ c then up(x) ⊳ c is a valid guard.

Finite solution to the state-based guards equations
The least solution to the equations of Definition 24 can be obtained by a standard Kleene iteration for fixed points computation. For each i ≥ 0 and each state q, define: The iteration stabilizes when there exists a k satisfying G k+1 (q) = G k (q) for all q. At stabilization, the values G k (q) satisfy the equations of Definition 24, and give the required G(q). However, as we mentioned earlier, this iteration might not stabilize at any k. We will now develop some observations that will help detect after finitely many steps if the iteration will stabilize or not.
Suppose we colour the set G i+1 (q) to red if either there exists a diagonal constraint x−y ⊳ c ∈ G i+1 (q)\G i (q) (a new diagonal is added) or there exists a non-diagonal constraint x ⊳ c or c ⊳ x in G i+1 (q) \ G i (q) such that the constant c is strictly bigger than c ′ for respectively every non-diagonal x ⊳ c ′ or c ′ ⊳ x in G i (q) (a non-diagonal with a bigger constant is added). If this condition is not applicable, we colour the set G i+1 (q) green. The next observations say that the iteration terminates iff we reach a stage where all sets are green. Intuitively, once we reach green, the only constraints that can be added are non-diagonals having smaller (non-negative) constants and hence the procedure terminates.
is green for all q, then G i+1 (q) is green for all q.
Proof. Pick an i > 0 such that G i (q) is green for all q. Suppose there is a state p such that G i+1 (p) is red. By definition, there is a new constraint ϕ ∈ G i+1 (p) \ G i (p). This means that there is a transition (p, g, up, p 1 ) and a constraint ψ ∈ G i (p 1 ) \ G i−1 (p 1 ) (which is new at the i th stage) such that ϕ is wp(⊑ ψ , up).
If ϕ is a diagonal constraint, then ψ is also a diagonal constraint. This contradicts G i (p 1 ) being green and hence ϕ cannot be a diagonal constraint. Suppose ϕ is a non-diagonal x ⊳ c. Then, ψ is of the form z ⊳ d and the update up is of the form z := x + c 1 . Therefore ϕ is would have been added to G i (p). We know that is not the case as G i+1 (p) has been coloured red due to ϕ, and hence d − c 1 > d ′ − c 1 . This gives d > d ′ and hence G i (p 1 ) should be coloured red due to the added constraint z ⊳ d. Yet again, this contradicts that G i (p 1 ) is green. The case when ϕ is of the form c ⊳ x is similar. ◭ ◮ Lemma 28. Let K = 1 + |Q| · |X| · (|X| + 1). If there is a state p such that G K (p) is red, then there is no i such that G i (q) is green for all q.
Proof. Suppose there is a state p K such that G K (p K ) is red. This is due to some constraint with the fact that G K (p K ) is red due to ϕ K . The proof for the case ϕ K−1 = c ⊳ x is similar. Iterating this construction we obtain a sequence (t j = (p j , g j , up j , p j−1 )) 2≤j≤K of transitions and a sequence of constraints (ϕ j ) 1≤j≤K such that ϕ j ∈ wp(⊑ ϕ j−1 , up j ) for 2 ≤ j ≤ K and G j (p j ) is red due to the addition of ϕ j for 1 ≤ i ≤ K. Since K = 1 + |Q| · |X| · (|X| + 1) there exist two indices a and b with b > a such that p a = p b = p and both ϕ a and ϕ b have the same "support", that is, either ϕ a := x − y ⊳ c and ϕ b := x − y ⊳ c ′ , or ϕ a := x ⊳ c and ϕ b := x ⊳ c ′ , or ϕ a := c ⊳ x and ϕ b := c ′ ⊳ x.
Suppose ϕ a and ϕ b are non-diagonals. Since G b (p) is red due to ϕ b and ϕ a ∈ G a (p) ⊆ G b (p), we deduce that c ′ > c. Since p a = p b = p, the weakest pre-condition computations for the sequence of updates up a+1 , . . . , up b can be applied starting from ϕ b ∈ G b (p b ) and results in a constraint x ⊳ c ′′ or c ′′ ⊳ x with c ′′ − c ′ > c ′ − c > 0. This can be iterated infinitely often to give red sets at each stage.
Assume now that ϕ a := x − y ⊳ c and ϕ b := x − y ⊳ c ′ . Since ϕ b is a new diagonal constraint in G b (p) ⊇ G a (p) we get c ′ = c. As above, the weakest pre-condition computations for the sequence of updates up a+1 , . . . , up b can be applied starting from ϕ b ∈ G b (p b ) and results in a constraint x − y ⊳ c ′′ with c ′′ − c ′ = c ′ − c = 0. This can be iterated infinitely often giving infinitely many new diagonal constraints.
Both cases imply that there can be no i where all G i (q) are green. ◭ ◮ Proposition 29. The least solution of the local constraint equations for an updateable timed automaton is finite iff G K (q) is green for all q and where K = 1 + |Q| · |X| · (|X| + 1). All decidable classes of [11] can be shown decidable with our approach, by showing stabilization of the G(q) computation. Proof. With the above combination of guards and updates, each constraint added in the G(q) computation apart from the guards of the automaton is either a diagonal constraint with a constant appearing in the guards of the automaton, or is a non-diagonal constraint with a constant c ≤ 2M where M = max{|d| | d occurs in a guard or an update of the automaton}. Hence the G(q) computation terminates. ◭

Experiments
We have implemented the reachability algorithm for timed automata with diagonal constraints (and only resets as updates) based on the simulation approach (Page 6) using the LU A simulation (Definition 14) for pruning zones. The algorithm for Z ⊑ LU G Z ′ comes from Section 4. Our implementation is over a prototype tool TChecker [18]. The core reachability algorithm based on the simulation approach is already implemented in TChecker for timed automata without diagonal constraints. We have added the state based guards computation, and the ⊑ LU G simulation algorithm. In the rest of this section, we explain the timed automata models and the experimental results (Table 1). Experiments are reported in Table 1. We take model Cex from [7,25] and Fischer from [25]. We are not aware of any other "standard" benchmarks containing diagonal constraints. In addition to these two models, we introduce the following new benchmark. Figure 1 UPPAAL model capturing the timing requirements in a scheduling problem. The state marked C is called a committed location, which is a convenient modeling feature in UPPAAL. Time is not allowed to elapse in such a state.

Job-shop scheduling with dependent tasks.
This is an extension of the job-shop scheduling using (diagonal-free) timed automata [1]. The benchmark models the following situation: there are n jobs J 1 , J 2 , . . . , J n and k machines m 1 , m 2 , . . . , m k . Each job J i consists of two tasks given by triples where m, a, b (with appropriate indices) denote respectively the machine on which the task needs to be executed and [a, b] is an interval such that the time needed to execute the task is in this interval. In addition to this, each job has an overall deadline D i . The question now is if the tasks can be scheduled so that all jobs finish within their deadline. This problem can be modeled using diagonal-free automata. We now add a further constraint (excluding the indices for clarity): for a job J with tasks T 1 and T 2 , let t 1 and t 2 be the execution times. We know that t 1 ∈ [a 1 , b 1 ] and t 2 ∈ [a 2 , b 2 ]. Now, we have a dependency: if t 1 ≥ c 1 , then t 2 ≤ c 2 and if t 1 ≤ d 1 then t 2 ≥ d 2 for suitable constants c 1 , c 2 , d 1 , d 2 . This says that if the first task is executed for a longer time, the second task needs a shorter time to finish and vice versa. This kind of a dependency has a natural modeling using diagonal constraints. We illustrate in Figure 1 a part of the UPPAAL-style automaton for a job J capturing the timing requirements -we make use of the feature of committed locations in UPPAAL [5]. Time is not allowed to elapse in this state, and in a product the components with committed locations execute first. This modeling template of the job scheduling problem with dependent tasks can easily be extended when jobs have more tasks and there are similar dependencies between them. In addition to this automaton, we model the mutual exclusion between machines (by use of extra boolean variables): each machine can have atmost one task at a time. We vary the number of jobs and get different timed automata with diagonal constraints. Note that although this model is acyclic, there are "cross simulations" which help in pruning the search: the same state q can be reached by multiple paths and simulations help in cutting out new searches.
Each model considered above is a product of a number of k timed automata. In the table we write the name of the model and the number k of automata involved in the product. Along with this, we report the number of diagonal constraints in each of them.
Experimental results. We considered four algorithms as mentioned in the caption of Table 1. Under each algorithm, we report on the number of zones enumerated and the time taken. The first algorithm gives a huge gain over the second algorithm (upto four orders of magnitude in the number of nodes, and even better for time) and gives a less marked, but still significant, gain over the third and fourth algorithms. We provide a brief explanation of this phenomenon. The performance of the reachability algorithm is dependent on three factors: the parameters on which the extrapolation or simulation is based on: the maximum constant based M -simulations which use the maximum constant appearing in the guards, versus the LU -simulations which make a distinction between lower bound guards c ⊳ x and upper bound guards x ⊳ c (refer to [4] for the exact definitions of extrapolations based on these parameters, and [20] for simulations based on these parameters); LUsimulations are superior to M -simulations. the way the parameters are computed: global parameters which associate a bound to each clock versus the more local state based parameters as in Definition 8 which associate a set of bounds functions to each state [3]; local bounds are superior to global bounds. when diagonal constraints are present, whether zones get split or not: each time a zone gets split, new enumerations start from each of the new nodes; clearly, a no-splitting-ofzones approach is superior to zone splitting. Algorithm of column 1 uses the superior heuristic in all the three optimizations above. The no-splitting-of-zones was possible thanks to our simulation approach, which temporarily splits zones for checking Z ⊑ LU G Z ′ , but never starts a new exploration from any of the split nodes. The algorithm of column 2, which is implemented in the current version UPPAAL 4.1 uses the inferior heuristic in all the three above. In particular, it is not clear how the extrapolation approach can avoid the zone splitting in an efficient manner. The superiority of our approach gets amplified (by multiplicative factors) when we consider bigger products with many more diagonals. In the third algorithm, we give a diagonal free equivalent of the original model (c.f. Theorem 1) and use the UPPAAL engine for diagonal free timed automata. The UPPAAL diagonal free engine is highly optimized, and makes use of the superior heuristics in the first two optimizations mentioned above (the third is not applicable now as it is a diagonal free automaton). The third algorithm can be considered as a good approximation of the zone splitting approach to diagonal constraints using LU -abstractions and local guards.
The second and the third methods are the only possibilities of verifying timed models coming with diagonal constraints in UPPAAL. Both these approaches are in principle prone to a 2 #D blowup compared to the first approach, where D gives the number of diagonal constraints. The table shows that a good extent of this blowup indeed happens. The UPPAAL diagonal free engine uses "minimal constraint systems" [5] for representing zones, whereas in our implementation we use Difference Bound Matrices (DBMs) [14]. This explains why even with more nodes visited, UPPAAL performs better in terms of time in some cases. We have not included in the table the comparison with two other works dealing with the same problem: the refined diagonal free conversion [25] and the LU simulation extended with LU d simulation for diagonals [17]. However, our results are better than the tables reported in these papers.

Conclusion
We have proposed a new algorithm for handling diagonal constraints in timed automata, and extended it to automata with general updates. Our approach is based on a simulation relation between zones. From our preliminary experiments, we can infer that the use of simulations is indispensable in the presence of diagonal constraints as zone-splitting can be avoided. Moreover, the fact that the simulation approach stores the actual zones (as opposed to abstracted zones in the extrapolation approach) has enabled optimizations for diagonalfree automata that work with dynamically changing simulation parameters (LU -bounds), which are learnt as and when the the zones are expanded [19]. Working with actual zones is also convenient for finding cost-optimal paths in priced timed automata [10]. Investigating these in the presence of diagonal constraints is part of future work. Currently, we have not implemented our approach for updateable timed automata. This will also be part of our future work.
Working directly with a model containing diagonal constraints could be convenient (both during modeling, and during extraction of diagnostic traces) and can also potentially give a smaller automaton to begin with. We believe that our experiments provide hope that diagonal constraints can indeed be used.