S Akshay
email: akshayss@cse.iitb.ac.in

Paul Gastin
email: paul.gastin@lsv.fr

Karthik R Prakash
email: karthikrprakash@cse.iitb.ac.in

Fast zone-based algorithms for reachability in pushdown timed automata

Keywords: Timed automata, Zone-based abstractions, Pushdown automata, Simulations, Reachability

Given the versatility of timed automata a huge body of work has evolved that considers extensions of timed automata. One extension that has received a lot of interest is timed automata with a, possibly unbounded, stack, also called pushdown timed automata (PDTA). While different algorithms have been given for reachability in different variants of this model, most of these results are purely theoretical and do not give rise to efficient implementations. One main reason for this is that none of these algorithms (and the implementations that exist) use the so-called zone-based abstraction, but rely either on the region-abstraction or other approaches, which are significantly harder to implement. In this paper, we show that a naive extension, using simulations, of the zone based reachability algorithm for the control state reachability problem of timed automata is not sound in the presence of a stack. To understand this better we give an inductive rule based view of the zone reachability algorithm for timed automata. This alternate view allows us to analyze and adapt the rules to also work for pushdown timed automata. We obtain the first zone-based algorithm for PDTA which is terminating, sound and complete. We implement our algorithm in the tool TChecker and perform experiments to show its efficacy, thus leading the way for more practical approaches to the verification of timed pushdown systems.

Introduction

Timed automata [START_REF] Alur | A theory of timed automata[END_REF] are a popular formalism for capturing real-time systems, and of use for instance, in model checking of cyber-physical systems. They extend finite automata with real variables called clocks whose values increase over time; transitions are guarded by constraints over these variables. The main problem of interest is the reachability problem, which asks whether a given state can be reached while satisfying the constraints imposed by the guards. This problem is known to be PSPACE-complete (already shown in [START_REF] Alur | A theory of timed automata[END_REF]). The PSPACE algorithm, uses the so-called region-automaton construction, which essentially abstracts the timed automaton into an exponentially larger finite automaton of regions (collections of clock valuations), which is sound and complete for reachability.

Despite this complexity-theoretic hardness, the model of timed automata has proved to be extremely influential and versatile, resulting in an enormous body of work on its theory, variants and extensions over the past 25 years. Almost since its inception, researchers also began to develop tools to extend from theoretical algorithms to solve practical problems. Such tools range from the classical and richly featured tool UPPAAL [START_REF] Bengtsson | Uppaal-a tool suite for automatic verification of real-time systems[END_REF][START_REF] Kim G Larsen | Uppaal in a nutshell[END_REF] to the more recent open-source tool TChecker [START_REF] Herbreteau | Tchecker[END_REF], which have been used on industry strength benchmarks and perform rather well on many of them. These tools use a different algorithm for reachability, where reachable sets of valuations are represented as zones and explored in a graph. While a naive exploration of zones does not terminate, the algorithms used identify different strategies [START_REF] Behrmann | Lower and upper bounds in zone based abstractions of timed automata[END_REF][START_REF] Frédéric Herbreteau | Better abstractions for timed automata[END_REF][START_REF] Herbreteau | Using nonconvex approximations for efficient analysis of timed automata[END_REF], e.g., subsumption or simulations, extrapolations, for pruning the zone-based exploration graphs, while preserving soundness and completeness of reachability. While this does not change the worst case complexity, in practice, the zone exploration results in much better practical performance as it allows on-the-fly computation of reachable zones. One could even argue that the wider adoption of timed automata paradigm in the verification community has been a result of scalable implementations and tools built on this zone-based approach.

In light of this, zone-based algorithms are often looked for to improve practical performance of extensions of timed automata as well. For instance, for timed automata with diagonal constraints, classical zone-based approaches were shown to be unsound [START_REF] Bouyer | Forward analysis of updatable timed automata[END_REF]12], but recently, an approach has been developed which adapts the existing construction and obtains fast zone-based algorithms [START_REF] Gastin | Fast algorithms for handling diagonal constraints in timed automata[END_REF]. In the present paper, we are interesting in adding a different feature to timed automata, namely an unbounded lifo-stack. This results in a powerful model of pushdown timed automata (PDTA for short), in which the source of "infinity" is both from real-time and the unbounded stack. Unsurprisingly, this model and its variants have been widely studied over the last 20 years with several old and recent results on decidability of reachability, related problems and their complexity, including [START_REF] Parosh | Dense-timed pushdown automata[END_REF][START_REF] Akshay | Timed systems through the lens of logic[END_REF]3,[START_REF] Akshay | Revisiting underapproximate reachability for multipushdown systems[END_REF][START_REF] Akshay | Towards an efficient tree automata based technique for timed systems[END_REF][START_REF] Bouajjani | On the automatic verification of systems with continuous variables and unbounded discrete data structures[END_REF][START_REF] Clemente | Timed pushdown automata revisited[END_REF]14,[START_REF] Clemente | Timed pushdown automata and branching vector addition systems[END_REF][START_REF] Dang | Pushdown timed automata: a binary reachability characterization and safety verification[END_REF]. A wide variety of techniques have been employed to solve these problems, from region-based abstractions, to using atoms and systems of constraints, to encoding into different logics etc. However, except for [START_REF] Akshay | Towards an efficient tree automata based technique for timed systems[END_REF][START_REF] Akshay | Revisiting underapproximate reachability for multipushdown systems[END_REF], to the best of our knowledge, none of the others carry an implementation. In [START_REF] Akshay | Towards an efficient tree automata based technique for timed systems[END_REF], the implementation uses a tree-automaton implicitly based on regions and the focus in [START_REF] Akshay | Revisiting underapproximate reachability for multipushdown systems[END_REF] is towards multi-pushdown systems. A common factor of all these works is that none of them consider zone-based abstractions.

In this paper, we ask whether zone-based abstractions can be used to decide efficiently reachability questions in PDTA. We focus on the problem of wellnested control-state reachability of PDTA, i.e., given a PDTA, an initial and a target state, does there exist a run of the PDTA that starts at the initial state with empty stack and reaches the target state with an empty stack (in between, i.e., during the run, the stack can indeed be non-empty). As with timed automata, our goal here is towards its applicability to build powerful tools which could lead to wider adoption of the PDTA model and showcase its utility to model-checking timed recursive systems. As the first step, we examine the difficulties involved in mixing zones with stacks and point out that a naive adaptation of the zone-based algorithm would not be sound. Then we propose a new algorithm that modifies the zone-based algorithm to work for pushdown timed automata. This is done in three steps.

-First we view the zone-graph exploration at the heart of the zone-based reachability algorithm for timed automata as a least fixed point computation of two inductive rules. When applied till saturation, they compute a sound and complete finite abstraction of the set of all reachable zones. -Next, this view allows us to generalize the approach in the presence of a stack by adding new inductive rules that correspond to push and pop transitions, and hence are specific to the stack operation. There are two main technical difficulties in this. First, we need to ensure termination of the fixed point computation, using a strong enough pruning condition of the (a priori infinite) zone graph to ensure finiteness, while being sound and not adding spurious runs. Second, we want to aggressively prune the graph as much as possible to obtain an efficient zone-exploration algorithm. We show how we can minimally change the condition of pruning in the zone exploration graph to achieve this delicate balance. Indeed, in doing so we use a judicious combination of the subsumption (or simulation) relation and an equivalence relation for obtaining a fixed point computation for PDTA that is terminating, while being sound and complete. -Finally, we build new data structures that allow us to write an efficient algorithm that implements this fixed point computation. While getting a correct algorithm is relatively simple, to obtain an efficient one, we must again encounter and overcome several technical difficulties.

We implement our approach to build the first zone-based tool that efficiently solves well-nested control state reachability for PDTA. Our tool is built on top of existing infrastructure of TChecker [START_REF] Herbreteau | Tchecker[END_REF], an open source tool and benefits from many existing optimizations. We perform experiments to show the practical performance of multiple variants of our algorithm and show how our most optimized version is vastly better in performance than other variants and of course the earlier region-based approach on a suite of example benchmarks.

We note that our PDTA model differs slightly from the model considered in [START_REF] Parosh | Dense-timed pushdown automata[END_REF]3], as there is no age on stack and time spent on stack cannot be compared with clocks. Hence our model is closer to [START_REF] Bouajjani | On the automatic verification of systems with continuous variables and unbounded discrete data structures[END_REF][START_REF] Dang | Pushdown timed automata: a binary reachability characterization and safety verification[END_REF]. However, in [START_REF] Clemente | Timed pushdown automata revisited[END_REF], it was shown that these two models are equivalent, more specifically, the stack can be untimed without loss of expressivity (albeit with an exponential blowup). Thus our approach can be applied to the other model as well by just untiming the stack. There are other more powerful extensions [START_REF] Clemente | Timed pushdown automata and branching vector addition systems[END_REF]14] studied especially in the context of binary reachability, where only theoretical results are known. We also remark that the idea of combining the subsumption relation between zones with an equivalence relation also occurs while tackling liveness, or Buchi acceptance, in timed automata. This has been studied in depth [START_REF] Tripakis | Checking timed büchi automata emptiness on simulation graphs[END_REF][START_REF] Laarman | Multi-core emptiness checking of timed büchi automata using inclusion abstraction[END_REF][START_REF] Frédéric Herbreteau | Why liveness for timed automata is hard, and what we can do about it[END_REF], where the naive zone-based algorithm does not work, forcing the authors to strengthen the simulation relation in different ways. Though these problems are quite different, there are surprising similarities in the issues faced, as explained in Section 3.

The structure of the paper is as follows: we start with preliminaries and move on to the difficulty in using zones and simulation relations in solving reachability in PDTA. Then, we introduce in Section 4 our inductive rules for timed automata and PDTA and show their correctness. In Section 5, we present our algorithm and helpful data-structural advancements. We show the experimental performance in Section 6 and end with a brief conclusion. Proofs that are missing and more experimental results can be found in the long version of the paper available at [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF].

Preliminaries

Timed automata

Timed automata extend finite-state automata with a set X of (non-negative) real-valued variables called clocks. We let Φ(X) denote the set of constraints ϕ that can be formed using the grammar:

ϕ ::= x ∼ c | x -y ∼ c | ϕ ∧ ϕ,
where x, y ∈ X, c ∈ N, ∼ ∈ {≤, ≥, <, >}, where each x ∼ c is called an atomic constraint. A clock valuation is a map v : X → R ≥0 and is said to satisfy ϕ, denoted v |= ϕ, if ϕ evaluates to true when each clock x ∈ X is replaced with v(x). For δ ∈ R ≥0 , we write v + δ to denote the valuation defined as (v + δ)(x) = v(x) + δ for all clocks x. For R ⊆ X, we write [R]v to denote the valuation obtained by resetting clocks in R, i.e., ([R]v)(x) = 0 if x ∈ R, and ([R]v)(x) = v(x) otherwise. Finally, v 0 is the valuation that sets all clocks to 0.

A timed automaton A is a tuple (Q, X, q 0 , ∆, F ), where Q is a finite set of states, X is a finite set of clocks, q 0 ∈ Q is an initial state, F ⊆ Q is the set of final states and ∆ ⊆ Q × Φ(X) × 2 X × Q is a set of transitions. A transition t ∈ ∆ is of the form (q, g, R, q ), where q, q are states, g ∈ Φ(X) is the guard of the transition and R ⊆ X is the set of clocks that are reset at the transition. The semantics of a timed automaton A is given as a transition system TS (A) over configurations. A configuration is a pair (q, v) where q ∈ Q is a state and v is a valuation, with the initial configuration being (q 0 , v 0 ). The transitions are of two types. First, for a configuration (q, v) and δ ∈ R ≥0 , (q, v) δ -→ (q, v + δ) is a delay transition. Second, for t = (q, g, R, q ) ∈ ∆, (q, v)

t - → (q , v ) is a discrete transition if v |= g and v = [R](v).
A run is an alternating sequence of delays and discrete transitions starting from the initial configuration, and is said to be accepting if the last state in the sequence is a final state.

Reachability, Zones and simulations

The problem of control-state reachability asks whether a given timed automaton has an accepting run. This problem is known to be PSPACE-complete [START_REF] Alur | A theory of timed automata[END_REF], originally shown via the so-called region abstraction. Note that, since TS (A) is infinite, some abstraction is needed to get an algorithm. In practice however, the abstraction used to solve reachability, e.g., in tools such as UPPAAL [START_REF] Kim G Larsen | Uppaal in a nutshell[END_REF] or TChecker [START_REF] Herbreteau | Tchecker[END_REF] is the zone abstraction. A zone Z is defined as a set of valuations defined by a conjunction of atomic clock constraints. Given a guard g and reset R, we define the following operations on zones: time elapse

- → Z = {v + δ | v ∈ Z, δ ∈ R ≥0 }, guard intersection g ∩ Z = {v ∈ Z | v |= g} and reset [R]Z = {[R]v | v ∈ Z}.
The resulting sets are also zones. With this, we can define the zone graph ZG(A) as a transition system obtained as follows: the nodes are (state, zone) pairs and (q, Z)

t - → (q , Z ), if t = (q, g, R, q ) is a transition of A and Z = -------→ [R](g ∩ Z)
. The initial node is (q 0 , Z 0 = --→ {v 0 }) and a path in the zone graph is said to be accepting if it ends at an accepting state. The zone graph is known to be sound and complete for reachability, but as the graph may still be infinite, this does not give an algorithm for solving reachability yet.

To obtain an algorithm, one resorts to different techniques such as extrapolation or simulation. Here we focus on simulation relations which will lead to finite abstractions. Given a timed automaton A, a binary relation on configurations is called a simulation if whenever (q, v) (q , v ), we have q = q and for each delay δ ∈ R ≥0 , (q, v + δ) (q, v + δ) and for each t = (q, g, R, q

1 ) ∈ ∆, if v |= g then v |= g and (q 1 , [R]v) (q 1 , [R]v ).
We often simply write v q v instead of (q, v) (q, v ). We can now lift this to sets Z, Z of valuations as Z q Z if for all v ∈ Z there exists v ∈ Z such that v q v . We say that node (q, Z) is subsumbed by node (q, Z ) when Z q Z . As a consequence we obtain the following lemma.

Lemma 1. If (q, Z) t - → (q 1 , Z 1 ) in ZG(A) and Z q Z , then (q, Z ) t - → (q 1 , Z 1 ) and Z 1 q1 Z 1 . Proof. Indeed, let v 1 ∈ Z 1 = -------→ [R](g ∩ Z). We find v ∈ Z and δ ≥ 0 such that v |= g and v 1 = [R]v + δ. Since Z q Z , we find v ∈ Z with v q v . We deduce that v |= g and [R]v q1 [R]v , which implies v 1 q1 v 1 with v 1 = [R]v + δ ∈ Z 1 = -------→ [R](g ∩ Z ).
A simulation is said to be finite if for every sequence of nodes (q 1 , Z 1 ), (q 2 , Z 2 ), . . . there exist two nodes (q i , Z i ) and (q j , Z j ) with i < j such that q i = q j and Z j qi Z i . The importance of the finiteness is that it allows us to stop exploration of zones along a branch of the zone graph: when a node (q j , Z j ) is reached which is subsumed by an earlier node (q i , Z i ), we may cut the exploration since all control states reachable from the latter are already reachable from the former. For a timed automaton A, we call this pruned graph as ZG (A). Thus, if the simulation relation is finite, then ZG (A) is finite, sound and complete for control state reachability. We formalize this algorithm in Section 4, using inductive rules.

Various finite simulations have been shown to exist in the literature, including the famous LU-abstractions [START_REF] Behrmann | Lower and upper bounds in zone based abstractions of timed automata[END_REF], and more recent G-abstractions based on sets of guards [START_REF] Gastin | Fast algorithms for handling diagonal constraints in timed automata[END_REF]. Hence this theory indeed has resulted in better implementations and is used in standard tools in this domain.

We will see that using simulation in the context of pushdown timed automata is not always sound, in some cases we need a stronger condition to stop the exploration. Towards this, we consider the equivalence relation on nodes induced by the simulation relation: Z ∼ q Z if Z q Z and Z q Z. We say that the simulation is strongly finite if the induced equivalence relation ∼ has finite index. Notice that strongly finite implies finite but the converse does not necessarily hold. Fortunately, the usual simulations for timed automata, in particular the LU-simulation and the G-simulation, are strongly finite.

Pushdown timed automata (PDTA)

A Pushdown Timed Automaton A is a tuple (Q, X, q 0 , Γ, ∆, F ), where Q is a finite set of states, X is a finite set of clocks, q 0 ∈ Q is an initial state, Γ is the stack alphabet, F ⊆ Q is the set of final states and ∆ is a set of transitions. A transition t ∈ ∆ is of the form (q, g, op, R, q ), where q, q are states, g ∈ Φ(X) is the guard of the transition and R ⊆ X is the set of clocks that are reset at the transition, op is one of three stack operations: nop or push a or pop a for some a ∈ Γ .

The semantics of a PDTA A is given as a transition system TS (A) over configurations. A configuration here is a tuple (q, v, χ) where q ∈ Q is a state, v is a valuation, χ ∈ Γ * is the stack content, with the initial configuration being (q 0 , v 0 , ε). The transitions are of two types. First, for a configuration (q, v, χ)

and δ ∈ R ≥0 , (q, v, χ) δ -→ (q, v + δ, χ) is a delay transition. Second, for t = (q, g, op, R, q ) ∈ ∆, (q, v, χ) t - → (q , v , χ ) is a discrete transition if v |= g, v = [R](v) and -if op = nop, then χ = χ, -if op = push a then χ = χ • a, -if op = pop a , then χ = χ • a.
A run is an alternating sequence of delays and discrete actions starting from the initial configuration. It is accepting if the last state in the sequence is final.

Our main focus is the well-nested control state reachability problem for PDTA, which asks whether a configuration (q, v, ε) with q ∈ F is reachable, where the stack is empty. Later, in Section 7, we remark how our solution can be extended to solve general control state reachability, i.e., asking whether a configuration (q, v, χ) with q ∈ F is reachable, possibly with a nonempty stack χ.

q0 q1 q2 q3 q4 push a push a push a y ≤ 3, pop a , {x} x ≥ 1, pop a , {x}
Fig. 1: A simple PDTA with 2 clocks {x, y}. Note that if we ignore the push/pop actions we get a TA, say A.

q0, Z0 q1, Z0 q2, Z0 q3, Z0 q4, (0 ≤ y -x ≤ 3) q4, (1 ≤ y -x ≤ 4) push a
push a push a pop a pop a Fig. 2: Zone graph with simulation edges for finiteness. Again ignoring push/pop actions gives us a zone graph for the TA. Z 0 is the initial zone.

3 Zones in PDTA and the problem with simulations

As mentioned earlier, zones are collections of clock valuations defined by conjunctions of timing constraints, and exploring zones reached by a timed automaton gives a sound and complete abstraction for state reachability. To make sure that the exploration is finite we need to prune the graph and one way this is done by simulation, i.e., not exploring paths from some nodes if they are "subsumed" by earlier nodes visited in the graph. Consider Figure 1, in which we ignore the push a and pop a or we can think of them as internal actions. Then the usual zone-graph construction with simulation would give the graph depicted in Figure 2. In this section, just for illustration we instantiate the simulation relation to be the well-known LU-simulation (we do not give the definition here as it is not relevant to what comes later, instead we refer to earlier work [START_REF] Behrmann | Lower and upper bounds in zone based abstractions of timed automata[END_REF]). Using this, we obtain that the rightmost node is subsumed by the previous one, and hence the dotted simulation edge. If we did not do this we immediately observe that we get an infinite graph with increasing sets of zones. Now, our first question is whether this zone exploration with simulation can be lifted to PDTA. In this example, if we were to add back the push/pop edges, we get exactly the same Zone graph with annotations, and further, the final state is indeed reachable. Hence, for this particular example we do obtain a finite, sound and complete graph exploration. However, in general it turns out that the procedure is not sound.

Consider the example in Figure 3. In this example, again considering it as a TA (ignoring the push/pops), we would get the zone graph below, which would be finite, sound and complete for reachability in that TA. But if we consider it as a PDTA, now doing the same does not preserve soundness. In other words, in the PDTA, q 3 is no longer reachable. However, in the zone graph we would conclude that it is reachable due to the simulation edge. If, to fix this, we remove the dotted simulation edge, then we will lose finiteness.

Thus, it seems that we have a difficult situation where zones with the simulation relation, needed for termination, do not preserve soundness. This situation resembles the situation studied in [START_REF] Tripakis | Checking timed büchi automata emptiness on simulation graphs[END_REF][START_REF] Laarman | Multi-core emptiness checking of timed büchi automata using inclusion abstraction[END_REF][START_REF] Frédéric Herbreteau | Why liveness for timed automata is hard, and what we can do about it[END_REF], where the authors study liveness or Buchi-acceptance conditions in timed automata. Again in that situation, the

q0 q2 q3 q1 x ≥ 1, {x} y ≤ 1, push a pop a pop a q0, Z0 q2, Z0 q1, y -x ≥ 1 q3, Z0 q0, y -x = 1 push a pop a pop a
Fig. 3: A PDTA and its zone graph with simulation. With the simulation (dotted) edges, q 3 is reachable in the zone graph, but its not reachable in the PDTA.

naive algorithm with zone simulation does not work and the authors are forced to strengthen the simulation relation in different ways. Surprisingly, it turns out, that even in our very different problem setting of reachability in PDTA, a similar solution works. That is, we replace simulation by equivalence (defined in the previous section) as the pruning criterion. However, there are two issues (i) it is not easy to prove its correctness and (ii) this is far from efficient as shown in the experimental section. Our goal to use zones in the first place was efficiency and hence we would like to prune the zone graph as much as possible, i.e., we would like to use simulation edges as much as possible. In the next two sections, we describe our fix. We first show a different view of the exploration algorithm as a fixed point rule based approach. This allows us to then describe our fix in the same language, which is much easier to understand conceptually. Also as a corollary we will be able to show that using equivalence everywhere also gives a correct algorithm. After proving the correctness of our rule-based algorithm, we then tackle the challenges in implementing it.

Viewing reachability algorithms using rewrite rules

In this section, our goal is to compute a set S of nodes of the zone graph of a PDTA, as a least fixed point of a small set of inductive rules, such that a control state q occurs in S, i.e., (q, Z) ∈ S for some Z iff q is reachable in the PDTA from its initial state. To understand the rules and their correctness it is easier to first visualize this on plain timed automata without any push-pop edges.

Rewrite rules for Timed automata.

Given a TA A = (Q, X, q 0 , ∆, F ), the set S containing all reachable nodes of the zone graph, can be obtained as the least fixed point of the following inductive rules, with a natural deduction style of presentation.

S := {(q 0 , Z 0 )} (q, Z) ∈ S q g,R --→ q Z = ------→ R(g ∩ Z) = ∅ Trans S := S ∪ {(q , Z )}
Let S * denote the set at the fixed point by starting with the start rule and repeatedly applying the trans rule. It is easy to see that this computes the set of all reachable nodes of the zone graph: the start rule starts with the initial node and each application of trans rule takes a reachable node and applies a transition of the automaton and includes the resulting node reached. However, this set S * is a priori infinite since number of zones is infinite.

To make it finite we add a condition under which we will apply the transition rule based on a finite simulation relation (let us denote it ) for A.

(q, Z) ∈ S q g,R --→ q Z = ------→ R(g ∩ Z) = ∅ Trans- S := S ∪ {(q , Z )}, unless ∃(q , Z ) ∈ S, Z q Z
Thus to obtain an algorithm, we would explore all nodes in the Zone graph using a search algorithm (say DFS/BFS) and we would add a node only if it is not subsumed by an already visited node, according to the simulation relation. We explained in Section 2.2 that doing this preserves soundness and completeness and gives a finite exploration.

Lemma 2. Let S * denote any set obtained from the start rule and by repeatedly applying Trans-till a fixed point is reached. Note that depending on the order of applications we may have different sets. Then we have:

1. (finiteness) S * is finite. 2. (soundness and completeness) For all q ∈ Q, a configuration (q, v) is reachable from (q 0 , v 0 ) in the TA A iff (q, Z) ∈ S * for some zone Z.

We do not give the proof here as (i) it is only a reformulation of known results and (ii) it will be subsumed by the much stronger theorem we prove next.

Rewrite rules for PDTA.

Let A = (Q, X, q 0 , Γ, ∆, F ) be a PDTA, we will need not just a set but a tuple of sets. More precisely, we maintain a set of nodes S called root nodes. For each root node (q, Z) ∈ S, we also maintain a set of nodes, denoted S (q,Z) . The intuition is that root nodes are those that can be reached after pushing a symbol to the stack, whereas S (q,Z) will be the set of nodes that can be reached from (q, Z) with a well-nested run, i.e., starting with an empty stack and ending in an empty stack. This is to avoid storing the stack contents in our algorithm, which would be another source of infinity. Again, we use simulations to make the computation finite. So we fix a strongly finite simulation relation for A.

Our inductive rules for the control state reachability of pushdown timed automata are given in Table 1. Note that the internal rule is just the same as for timed automata above. The start rule not only starts the set of nodes computation but also the set of roots computation as described above. So the only interesting rules are the Push and Pop rules. The push rule says that when a push is encountered, then we must start exploring from a new root (i.e., context). So the only complicated rule is the Pop rule. Here the intuition is that if we see

Start S := {(q0, Z0)}, S (q 0 ,Z 0 ) := {(q0, Z0)} (q, Z) ∈ S (q , Z ) ∈ S (q,Z) q g,nop,R -----→ q Z = ------→ R(g ∩ Z ) = ∅ Internal S (q,Z) := S (q,Z) ∪ {(q , Z )}, unless ∃(q , Z ) ∈ S (q,Z) , Z q Z (q, Z) ∈ S (q , Z ) ∈ S (q,Z) q g,push a ,R ------→ q Z = ------→ R(g ∩ Z ) = ∅ Push S := S ∪ {(q , Z )}, S (q ,Z ) = {(q , Z )}, unless ∃(q , Z ) ∈ S, Z ∼ q Z (q, Z) ∈ S (q , Z1) ∈ S (q , Z ) ∈ S (q,Z) (q 1 , Z 1 ) ∈ S (q ,Z 1 ) q g,push a ,R ------→ q q 1 g 1 ,pop a ,R 1 -------→ q2 Z = ------→ R(g ∩ Z ) ∼ q Z1 Z2 = --------→ R1(g1 ∩ Z 1 ) = ∅ Pop S (q,Z) := S (q,Z) ∪ {(q2, Z2)}, unless ∃(q2, Z 2 ) ∈ S (q,Z) , Z2 q 2 Z 2
Table 1: Inductive rules for control state reachability of PDTA a push at a node and from a root equivalent to the root created from it, (i.e., its context) we see a matching pop reaching a new node, then this push-pop context is complete, and we can add this new node to the set of reachable nodes. This is precisely the point where we need equivalence rather than simulation and this will be made clear in the proof of the theorem below.

Theorem 1. Let S * and (S (q,Z) ) (q,Z)∈S * denote any tuple of sets obtained from the start rule and by repeatedly applying the rules in Table 1 till a fixed point is reached 3 . Note that we always have (q 0 , Z 0 ) ∈ S * . The following statements hold:

1. (finiteness) S * is finite and for each (q, Z) ∈ S * , S (q,Z) is finite.

2. (completeness) For each (q, Z) ∈ S * , if there exists a run (q, v, ε) * -→ (q , v , ε) of A with {v} q Z, then there exists (q , Z ) ∈ S (q,Z) such that {v } q Z . 3. (soundness) For each (q, Z) ∈ S * , (q , Z ) ∈ S (q,Z) and v ∈ Z , there exists a run in PDTA from (q, v, ε) to (q , v , ε) with v ∈ Z and v q v . Proof. 1. Note that only the Push rule creates new root nodes and the red condition states that a new root node is added only if there isn't already an equivalent node in S * . Since the simulation relation is strongly finite, the set of roots S * must be finite. Also, before adding a node to some S (q,Z) with the internal rule or the pop rule, we check that the node is not subsumed by an existing one. Since the simulation relation is finite, this ensures that each set S (q,Z) is finite.

2. Let (q, Z) ∈ S * and assume that (q , v , ε) is reachable from some (q, v, ε) with v q Z, i.e., there exists a run (q, v, ε)

= (q 1 , v 1 , χ 1 ) → • • • → (q n , v n , χ n ) = (q, v, ε) * -→ (qi, vi, ε) t -→ (qi+1, vi+1, a) * -→ (qn-1, vn-1, a) t1 -→ (qn, vn, ) v q Z vi qi Zi vi+1 qi+1 Zi+1 vn-1 qn-1 Zn-1 vn qn Zn (q, Z) ∈ S (qi, Zi) ∈ S (q,Z) Zi+1 = ------→ R(g ∩ Zi) (qn-1, Zn-1) ∈ S (qi+1,Z i+1 ) Zn = ----------→ R1(g1 ∩ Zn-1) Zi+1 ∼q i+1 Z i+1 Zn qn Z n (qi+1, Z i+1 ) ∈ S (qn, Z n ) ∈ S (q,Z)
Fig. 4: Construction for the completeness-push-pop last sub-case.

(q , v , ε). We will then show that v n qn Z for some (q n , Z ) ∈ S (q,Z) . The proof is by induction on n. Base case: For n = 1 we have q = q and v = v. The result is obtained by taking Z = Z. Notice that (q, Z) ∈ S (q,Z) follows immediately from the start rule if q = q 0 , Z = Z 0 or from the push-create rule.

Let us then assume that the statement holds for runs of length at most n -1. Consider any run of the form (q, v, ε) = (q 1 , v 1 , χ 1 ) → • • • → (q n , v n , χ n = ε) with v q Z. Notice that its last transition (q n-1 , v n-1 , χ n-1 ) → (q n , v n , χ n ) cannot be a push transition (in the PDTA) since χ n = ε. Hence, we have three subcases, depending on the last transition.

-Time elapse.

χ n-1 = χ n = ε, q n-1 = q n = q , v n = v n-1 + δ for some δ ∈ R ≥0 .
Applying induction hypothesis, we have v n-1 q Z for some (q , Z ) ∈ S (q,Z) . Since zones are closed under time elapse, we get Z = -→ Z and by definition of the simulation relation v n = v n-1 + δ q -→ Z = Z . This completes the case.

-Discrete internal transition. In this case χ n-1 = χ n = ε, t = q n-1 g,nop,R -----→ q n , v n-1 |= g and v n = [R]v n-1 . Then applying induction hypothesis, there exists (q n-1 , Z ) ∈ S (q,Z) such that v n-1 qn-1 Z . Now let Z = ------→ R(g ∩ Z ). From the definition of the simulation relation we get v n qn Z . Then, applying the Internal rule, there exists (q n , Z ) ∈ S (q,Z) such that Z qn Z , with possibly Z = Z . Hence, v n qn Z qn Z , which completes the case.

-Pop transition. Then there exists 1 ≤ i < n -1 such that the run has the form:

(q 1 , v 1 , ε) → . . . → (q i , v i , χ i = ε) push a ----→ (q i+1 , v i+1 , χ i+1 = a) → . . . → (q n-1 , v n-1 , χ n-1 = a) pop a ---→ (q n , v n , χ n = ε)
, where the push and pop are matching transitions, i.e., |χ j | ≥ 1 for all i < j < n -1 (see Figure 4). Then by induction hypothesis at i, we have v i qi Z i for some (q i , Z i ) ∈ S (q,Z) .

(

) 1
From the push transition we have

∃t = q i g,push a ,R ------→ q i+1 ∈ ∆ with v i |= g and v i+1 = [R]v i . (2)
Let

Z i+1 = ------→ R(g ∩ Z i )
. By definition of the simulation relation, we deduce from v i qi Z i that v i+1 qi+1 Z i+1 . We can apply the Push rule to obtain

(q i+1 , Z i+1 ) ∈ S * for some Z i+1 ∼ qi+1 Z i+1 (3)
possibly with Z i+1 = Z i+1 as a special case. Further the segment of run (q i+1 , v i+1 , a) → . . . (q n-1 , v n-1 , a) in the PDTA never pops the symbol a (by choice, since otherwise the push and pop would not be matching). Hence we will also have the same sequence of transitions forming a run (q i+1 , v i+1 , ε) → . . . (q n-1 , v n-1 , ε). Using v i+1 qi+1 Z i+1 ∼ qi+1 Z i+1 , we deduce that v i+1 qi+1 Z i+1 . By induction hypothesis,

v n-1 qn-1 Z n-1 for some (q n-1 , Z n-1 ) ∈ S (qi+1,Z i+1 ) . (4)
Finally, we have the pop transition

t 1 = q n-1 g1,pop a ,R1 -------→ q n ∈ ∆ with v n-1 |= g 1 and v n = [R 1 ]v n-1 . (5)
We let

Z n = ----------→ R 1 (g 1 ∩ Z n-1
). From v n-1 qn-1 Z n-1 and the definition of the simulation relation we obtain v n qn Z n . Then, combining all the above equations (1-5), and applying the Pop-rule we obtain some (q n , Z n ) ∈ S (q,Z) with Z n qn Z n (possibly Z n = Z n ). Finally we get v n qn Z n qn Z n . This completes the proof. 3. We will show that the following property is invariant by rule applications:

∀(q, Z) ∈ S, ∀(q , Z ) ∈ S (q,Z) , ∀v ∈ Z , there is a run (q, v, ε) * -→ (q , v , ε) with v ∈ Z and v q v (Inv)
The invariant holds initially, i.e., after application of the start rule. Indeed, in this case we have S = {(q 0 , Z 0 )} and S (q0,Z0) = {(q 0 , Z 0 )}. Hence (q , Z ) = (q, Z) = (q 0 , Z 0 ) and for all v ∈ Z 0 we can choose the empty run (q 0 , v, ε) 0 -→ (q 0 , v, ε). We show below that (Inv) is preserved by application of an internal/push/pop rule. Therefore, the invariant still holds when reaching the fixed point, which proves the soundness. Let us write S -and S - (q,Z) for the sets before the application of the rule and S and S (q,Z) for the sets after the application of the rule. Internal rule. Let (q, Z) ∈ S = S -, (q , Z ) ∈ S (q,Z) and v ∈ Z . If (q , Z ) ∈ S - (q,Z) then we get the result since (Inv) holds before applying the internal rule. Otherwise, there is some (q 1 , Z 1 ) ∈ S - (q,Z) and a transition t

= q 1 g,nop,R -----→ q with Z = ------→ R(g ∩ Z 1 ). By definition, there exists v 1 ∈ Z 1 such that v 1 |= g and v = [R]v 1 + δ for some δ ≥ 0. Hence we have a run (q 1 , v 1 , ε) t - → δ -→ (q , v
, ε). Since the invariant holds before the internal rule, there is a run (q, v, ε)

* -→ (q 1 , v 1 , ε) with v ∈ Z and v 1 q1 v 1 . Now since is a simulation we obtain that (q 1 , v 1 , ε) t - → δ -→ (q , v , ε)
with v q v and we are done. Push rule. Let (q, Z) ∈ S, (q , Z ) ∈ S (q,Z) and v ∈ Z . If (q, Z) ∈ S -then we get the result since (Inv) holds before applying the Push rule. Otherwise, we must have (q , Z ) = (q, Z) and we can choose the empty run (q, v , ε) Pop rule. Let (q, Z) ∈ S = S -, (q 2 , Z 2 ) ∈ S (q,Z) and v ∈ Z 2 . Again, if (q 2 , Z 2 ) ∈ S - (q,Z) then we get the result since (Inv) holds before applying the Pop rule. Otherwise, by definition of the pop rule we have:

0 -→ (q, v , ε). (q 1 , v4, a) t 1 -→ δ -→ (q2, v , ) (q , v3, a) * -→ (q 1 , v 4 , a) (q , v2, ε) t -→ δ -→ (q , v 3 , a) (q, v, ε) * -→ (q , v 2 , ε) t -→ δ -→ (q , v 3 , a) * -→ (q 1 , v 4 , a) t 1 -→ δ -→ (q2, v , )
1. some (q , Z ) ∈ S (q,Z) , 2. some push transition t = q g,push a ,R ------→ q , 3. some (q , Z 1 ) ∈ S with

Z 1 ∼ q Z = ------→ R(g ∩ Z ), 4. some (q 1 , Z 1 ) ∈ S (q ,Z1) , 5. some pop transition t 1 = q 1 g1,pop a ,R1 -------→ q 2 , with Z 2 = --------→ R 1 (g 1 ∩ Z 1 ). The construction below is illustrated in Figure 5. Since v ∈ Z 2 , we get some v 4 ∈ Z 1 such that v 4 |= g 1 and v = [R 1 ]v 4 + δ
for some δ ≥ 0. Hence we have a run (q 1 , v 4 , a) t1 -→ δ -→ (q 2 , v , ε). Now, applying the invariant to (q , Z 1 ) ∈ S, (q 1 , Z 1 ) ∈ S (q ,Z1) and v 4 ∈ Z 1 , we get a run (q , v 3 , ε) * -→ (q 1 , v 4 , ε) with v 3 ∈ Z 1 and v 4 q 1 v 4 . Hence, we also have a run (q , v 3 , a) * -→ (q , v 4 , a).

Let v 3 ∈ Z = ------→ R(g ∩ Z ) ∼ q Z 1 with v 3 q v 3 . we get some v 2 ∈ Z such that v 2 |= g and v 3 = [R]v 2 + δ for some δ ≥ 0. Hence we have a run (q , v 2 , ε) t - → δ -→ (q , v 3 , a).
Finally, we apply the invariant to (q, Z) ∈ S, (q , Z ) ∈ S (q,Z) and v 2 ∈ Z , we get a run (q, v, ε) * -→ (q , v 2 , ε) with v ∈ Z and v 2 q v 2 . By repeatedly applying the property of simulation , we may extend the run from (q , v 2 , ε) with (q , v 2 , ε)

t - → δ -→ (q , v 3 , a) * -→ (q 1 , v 4 , a) t1 -→ δ -→ (q 2 , v , ε)
where v 3 q v 3 q v 3 and v 4 q 1 v 4 q 1 v 4 . Finally v q2 v . Therefore, the invariant holds after the pop rule.

Algorithm for PDTA Reachability via Zones

In this section, we describe Algorithm 1 implementing the fixed point computation defined by the inductive rules in Table 1. We describe the structure of the algorithm and its main data-structures.

Notice first that the sets S and S (q,Z) for (q, Z) ∈ S can be alternatively represented as a single set of pairs of nodes: S = {[(q, Z), (q , Z )] | (q, Z) ∈ S and (q , Z ) ∈ S (q,Z) } .

We can recover S as the first projection of S and S (q,Z) as the second projection of S filtered by the first component being (q, Z). We use both notations below depending on which is more convenient. The start rule initializes S to {[(q 0 , Z 0 ), (q 0 , Z 0 )]}.

Let us consider first the rule for internal transitions. For each already discovered pair of nodes [(q, Z), (q , Z )] ∈ S (or (q , Z ) ∈ S (q,Z) with (q, Z) ∈ S), we have to consider each possible internal transition q g,nop,R -----→ q and check whether the node (q , Z ) with Z = ------→ R(g ∩ Z ) should be added to S (q,Z) or is subsumed by an existing node. This is like a graph traversal. The set S stores the already discovered pairs of nodes, and we will use a ToDo (unordered) list to store the newly discovered nodes from which outgoing transitions should be considered. The ToDo list should also consist of pairs [(q, Z), (q , Z )] so that when a new node (q , Z ) is discovered by an internal transition from (q , Z ) we know to which set S (q,Z) it should be added.

As we can see from Theorem 1-soundness, given (q, Z) ∈ S, the set S (q,Z) should consist of nodes reachable from (q, Z) via a well-nested run. Hence, when dealing with a pair [(q, Z), (q , Z )] ∈ S and we see a push transition q g,push a ,R ------→ q with Z = ------→ R(g ∩ Z ), we should not try to add the pair (q , Z ) to S (q,Z) since the corresponding run would not be well-nested. Instead, we should search for a matching pop transition which could be taken after a well-nested run starting from (q , Z ). This is why the push rule adds the new root (q , Z ) to S (unless it is equivalent to an existing root). The pair of nodes [(q , Z ), (q , Z )] is newly discovered and added to the ToDo list for further exploration.

The push transition may be matched with several pop transitions (which could be already discovered or yet to be discovered by the algorithm). To avoid revisiting the push transition many times, it will be stored by the algorithm in an additional set S push . More precisely, we will store in S push the tuple [(q, Z), a, (q , Z )] meaning that the root node (q , Z ) may be reached from the root node (q, Z) via a well-nested run reaching some (q , Z ) followed by a transition pushing a onto the stack.

Finally, assume that, when dealing with a pair [(q 1 , Z 1 ), (q 1 , Z 1 )] ∈ S, we see

a pop transition q 1 g1,pop a ,R1 -------→ q 2 with Z 2 = --------→ R 1 (g 1 ∩ Z 1
). We will check whether it can be matched with an already visited push transition, stored in the set S push as a pair [(q, Z), a, (q , Z )] with (q , Z ) = (q 1 , Z 1 ). If this is the case, the pop rule may be applied and the node (q 2 , Z 2 ) added to S (q,Z) (unless it is subsumed by an existing node). The newly discovered pair of nodes [(q, Z), (q 2 , Z 2 )] is also added to the ToDo list for further exploration. Once again, the pop transition may also be matched with push transitions that will be discovered later by the algorithm. To avoid revisiting the pop transition many times, we store the tuple [(q 1 , Z 1 ), a, (q 2 , Z 2 )] in a new set S pop . Data structures. We use a data structure TLM to store the triple of sets (S, S push , S pop ) and which is accessed with the following methods.

-TLM.create() creates the data structure with the three sets empty. -TLM.add(q, Z, q , Z ) adds [(q, Z), (q , Z )] to S. -TLM.addPush(q, Z, a, q , Z ) adds [(q, Z), a, (q , Z )] to S push . -TLM.addPop(q, Z, a, q , Z ) adds [(q, Z), a, (q , Z )] to S pop . -TLM.isNewRoot(q, Z) returns [false, Z ] if there exists some (q, Z ) ∈ S with Z ∼ q Z, and returns [true, Z] otherwise. -TLM.isNewNode(q, Z, q , Z ) returns false if ∃[(q, Z), (q , Z )] ∈ S with Z q Z , and returns true otherwise. -TLM.isNewPop(q, Z, a, q , Z ) returns false if ∃[(q, Z), a, (q , Z )] ∈ S pop with Z q Z , true otherwise. -TLM.isNewPush(q, Z, a, q , Z ) returns false if [(q, Z), a, (q , Z )] ∈ S push , and returns true otherwise. -TLM.iterPop(q, Z, a) returns the list of (q , Z ) with [(q, Z), a, (q , Z )] ∈ S pop .

-TLM.iterPush(a, q , Z ) returns the list of (q, Z), s.t. [(q, Z), a, (q , Z )] ∈ S push .

Concretely, the data structure should store sets of nodes (q, Z) and be able to search or iterate through such sets. In order to make the algorithm slightly faster, we will segregate our sets of nodes, with the name of the state. We will use a hashmap in order to accomplish this task. See Figure 6 where the concrete data structure is depicted.

We will use a first level hashmap to store the set of roots S. To implement TLM.isNewNode(q, Z, q , Z ), we first search for (q, Z) in the first level map, then a pointer TLM[q][Z][0] will lead to a second level hashmap for the set of nodes S (q,Z) and we search for (q , Z ) in this second level map. See Figure 6(b).

To implement TLM.isNewPop(q, Z, a, q , Z ) and TLM.iterPop(q, Z, a), we first search the root node (q, Z) in the first level map, then a pointer TLM[q][Z] [START_REF] Akshay | Timed systems through the lens of logic[END_REF] will lead to a second level hashmap storing the set of triples (a, q , Z ) such that [(q, Z), a, (q , Z )] ∈ S pop . To speed up the access, this second level pop map is segregated first on the key a, then on the key q to get the list of corresponding zones Z . See Figure 6(c,d).

Finally, we also store the set S push to implement TLM.isNewPush(q, Z, a, q , Z ) and TLM.iterPush(a, q , Z ). Notice that S push consists of triples [(q, Z), a, (q , Z )] where both (q, Z) and (q , Z ) are root nodes from S. Notice also that for the iteration we fix the second node (q , Z ). To get an efficient implementation, we first search the root node (q , Z ) in the first level map, then a pointer TLM[q ][Z ] [START_REF] Parosh | Dense-timed pushdown automata[END_REF] will lead to a second level hashmap storing the set of triples (a, q, Z) such that [(q, Z), a, (q , Z )] ∈ S push . To speed up the access, this second level push map is segregated first on the key a, then on the key q to get the list of corresponding zones Z. See Figure 6(c,d (d) For pushes/pops map, this is a map corresponding to root node (q1, Z1), and symbol a2 (say). The (q, Z) stored here is constructed using equivalence (pushes map), or using simulation (pops map).

Fig. 6: Two level map implementing the data structure TLM storing the sets S, S push , S pop .

We now show correctness of Algorithm 1. Note that TLM encodes a triple of sets (S, S push , S pop ) defined by:

S = {[(q, Z), (q , Z )] | (q , Z ) ∈ TLM[q][Z][0]} S push = {[(q, Z), a, (q , Z )] | (a, q, Z) ∈ TLM[q ][Z ][1]} S pop = {[(q, Z), a, (q , Z )] | (a, q , Z ) ∈ TLM[q][Z][2]}
Recall also the correspondence explained at beginning of Section 5 between a set S of pairs of nodes, and the set of roots S together with the sets of nodes S (q,Z) for (q, Z) ∈ S. Theorem 2. The set S encoded by TLM computed by Algorithm 1 is a fixed point obtained starting from the empty set by applying the inductive rules in Table 1. Therefore, Algorithm 1 terminates and is sound and complete for wellnested control state reachability of pushdown timed automata.

Proof (sketch).

1. For termination, if we look at our algorithm, we can clearly see that before adding a pair of nodes to the ToDo list, we add the same pair to S with TLM.add, and before that, we always check whether the pair is already in our TLM or not (isNewNode or isNewRoot). Since the size of the TLM is always bounded because we check either the first level map or the second level map before adding, the outer while loop will be called only a finite number of times. From this we can conclude that the algorithm will terminate. 2. For soundness we have prove that any change to the TLM is equivalent to applying one of the rewrite rules to (S, S push , S pop ), which is already known to be sound from Theorem 1. The changes to the TLM occur in lines 3, 9, 14, 21, 31. Since line 3 simply adds [(q 0 , Z 0 ), (q 0 , Z 0 )] to S, it simulates the start rule. For line 9, we can see that the pre conditions of internal rule Table 1 are met, with (q, Z) ∈ S, (q , Z ) ∈ S (q,Z) , the ifstatement (just above the line) stating that there is an nop transition from q to q , and Z = φ. Using all these we can see that indeed the operation can be performed. Similar arguments can be made for line 14, which simulates the push rule, and line numbers 21, 31 both for the pop-rule. 3. For completeness we have to prove that after termination of the algorithm, using (S, S push , S pop ) to encode TLM, we cannot use any of the rules mentioned in Table 1, to add anything extra to the TLM. Then from Theorem 1completeness we can conclude. For the start rule we can simply say that it was definitely executed (Line 3), so it cannot be executed again. For the internal rule we argue that if it can be applied after termination, then it should have been applied during execution. Since all transitions are considered in the for-loop, and the conditions before line 9 checks all the preconditions of the internal rule, it is certainly the case that a node (q , Z ) could not be added because either it was already added, or (q , Z ) ∈ S (q,Z) , Z q Z . The argument for the push rule is similar. For the pop rule to be applied we argue that there must be a push transition and a pop transition satisfying the pre-conditions in the pop rule. Since both of these are already present for zones in the TLM, we say that they must have been added to S push and S pop . We then concern ourselves with the order, arguing that if the push transition was discovered later the node either must already have been added (Line 21) or another node simulating the node must have been present in the TLM (Line 20). A similar argument is made in case the order of discovery is reversed.

For the full proof details, we refer the reader to the long version [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF].

Experiments and Results

Implementation We build on the existing architecture of an open-source tool for analysis of timed automata, TChecker [START_REF] Herbreteau | Tchecker[END_REF]. Our tool along with the benchmarks we used is available at https://github.com/karthik-314/PDTAreachability and more details can be found [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF]. The input for our implementation are PDTA, rather than TA so we modify TChecker in order to run our experiments. While most of the TChecker file format will remain the same, the only place where we make a change to the syntax of the input, will be the edges. TChecker uses the following format, for its transitions, edge:<Process>:<src>:<tgt>:<label>{ do:<Reset1(x=0)> ; <Reset2(y=0)> : provided: <guard1(x==0)> && <guard2(y>=1)>}

The new format in order to incorporate the pushes and pops will be, edge:<Process>:<src>:<tgt>:<label>{ do:<Reset1(x=0)> ; <Reset2(y=0)> : provided: <guard1(x==0)> && <guard2(y>=1)>} [<push/pop>:<symbol>]

In case the operation is nop, the square brackets are left empty.

We have implemented two variants of Algorithm 1 for PDTA and we will compare these between each other and also with a region-based approach. More precisely, we consider the following 3 algorithms:

-Simulation Based Approach ( LU ): Direct implementation of Algorithm 1.

-Equivalence Based Approach (∼ LU ): This is a variation of Algorithm 1, with two methods changed, • TLM.isNewNode(q, Z, q , Z ): Returns false if ∃[(q, Z), (q , Z )] ∈ S with Z ∼ q Z , and true otherwise. • TLM.isNewPop(q, Z, a, q , Z ): Returns false if ∃[(q, Z), a, (q , Z )] ∈ S pop with Z ∼ q Z , and true otherwise. As mentioned in Section 4, if instead of simulation, we just use equivalence everywhere, we do obtain a correct algorithm for reachability in PDTA. Hence it is interesting to compare it with the above approach.

-Region Based Implementation (RB): A previous implementation [START_REF] Akshay | Towards an efficient tree automata based technique for timed systems[END_REF], uses a region based approach in order to solve the non-emptiness problem in PDTA. We note two features of the algorithm. First, it uses a tree-automaton based approach for efficiency and correctness, but underlying it is the region (rather than zone) construction. Second, it works only with closed guards, while our approach works with closed and open guards. We note the following important points regarding our implementation: 1. The used in our implementation will be LU [START_REF] Behrmann | Lower and upper bounds in zone based abstractions of timed automata[END_REF], without extrapolation and with global clock bounds. 2. The ToDo list used currently uses LIFO (stack) ordering for popping of elements. This corresponds to a DFS exploration of the zone-graph. But we can use other data structures for this purpose as well, e.g., changing it to FIFO would give us a BFS exploration etc. 3. Both the simulation based and equivalence based approach are tested on PDTA with empty and non-empty languages, but we have ensured that both of them return an answer only after the entire exploration has been completed. In other words, we do not stop the exploration when we reach a final state. This is to make fair comparisons, where we do not terminate because of being "lucky" in encountering the final state early. Of course in practice we would not do this. In contrast, we note that the RB approach is an on the fly approach which returns non-empty as soon as the final state turns out to be reachable. All experiments are run on Intel-i5 10th Generation processor, with an 8GB RAM, with a timeout of 120 seconds.

Benchmarks. We used a total of 10 benchmarks in our experiments, but parameterized several of them in order to test the scalability and to give us more insight into performance comparisons. The benchmark and their parameterizations are explained in [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF]. We highlight only some salient points here. The benchmark B 1 is the PDTA from Figure 1. B 2 (k) is directly adapted from Figure 3 with the constant y ≤ 1 parametrized to y ≤ k, and k + 1 pops between q 0 and q 2 . Note that q 3 is unreachable regardless of the value of k. Benchmarks B 3 , B 4 are adapted from [START_REF] Akshay | Towards an efficient tree automata based technique for timed systems[END_REF] with B 3 involving untiming of a stack age into normal clocks. B 5 , B 6 involve significant interplay of push/pop edges and clocks and B 6 , B 7 also have open guards. More details can be found in [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF]. We also note that automata B 1 , B 3 (3, 4), B 5 (k 1 , k 2 ), B 8 , B 9 (k 1 , k 2 ) accept a nonempty language, while the rest are empty. As described earlier this does not change the performance of the simulation and equivalence based approaches, but may significantly change the performance of the Region Based Approach.

Results Table 2 contains a selection of our experimental results; more can be found in [START_REF] Akshay | Fast zone-based algorithms for reachability in pushdown timed automata[END_REF]. From the table, we conclude first that the zone based approach is indeed faster than the Region Based Approach for all examples. Second, the simulation based approach runs faster than the equivalence based approach for all examples if the ToDo priority for removal remains the same. In fact, the performance of the simulation based approach depends mostly on the size of the PDTA, but the equivalence based approach is dependant on the constants used in guards as well, which is even more the case for the region based approach. Finally, our approach can easily handle closed and open guards. Most of the timeouts that occurred during the experiments are due to Out of Memory (OoM) kills, especially for larger sized PDTAs. For smaller sized PDTA such as B 2 (100), the recorded number of nodes before timeout was 154700.

Regarding the performance, we would like to emphasize that B1, B2, B3, B4, B7 were designed to compare the Zone approach to the region (RB) approach. As a consequence these models are very simple and the number of explored nodes remains almost the same regardless of whether we use ∼ or to prune, which reflects in the times/sizes not being too different. However, the other examples B5, B6 are more complex and have nodes that get pruned during exploration (both using ∼ and ). Here we can see the clear improvement of over ∼ both in terms of time taken and also of number of explored nodes.

Discussion and Future work

In this paper, we examined how an unbounded stack can be integrated seamlessly with zone-abstractions in timed automata. We would like to point out that two easy extensions of our work are possible. First, as remarked earlier, our algorithm checks for well-nested reachability, i.e., it requires to reach a final state with empty stack for acceptance. But we can generalize this to general control-state reachability by showing that a control state q is reachable in the PDTA (with possibly a non-empty stack) iff some node (q, Z) is discovered by our algorithm and added to some S (q ,Z ) (and not just to S (q0,Z0) as in the well-nested case). While this idea is simple and requires only minor edits to the existing algorithm, the proof of correctness requires more work and we leave this for future work.

Secondly, we can handle the model with ages in stack as in [3,[START_REF] Parosh | Dense-timed pushdown automata[END_REF] with an exponential blowup (thanks to [START_REF] Clemente | Timed pushdown automata revisited[END_REF]). However, an open question is whether this blowup can be avoided in practice. As noted earlier, there exist extensions [START_REF] Clemente | Timed pushdown automata and branching vector addition systems[END_REF]14] studied especially in the context of binary reachability, which are expressively strictly more powerful, for which decidability results are known. It would be interesting to see how we can extend the zone-based approach to those models.

Finally, it seems interesting to examine further the link to the liveness problem, possibly allowing us to transfer ideas and obtain faster implementations. Another possibility would be to use the extrapolation operator (rather than, or in addition to, simulation), which we have not considered in this work.

  start

Fig. 5 :

 5 Fig. 5: Construction for the soundness.

  ). First level map constructed using equivalence ∼q for controlling size. Keys will be state names, values will be lists of quadruplets, each of which has four pointers to second level maps, second level pushes maps, second level pops maps, and zones.

	q1 q2 q3 q4 q5 q6 (a) q1 Z1 . . . Z2 . . . Z31 . . . Z32 . . . Z33 . . . Z41 . . . Z42 . . . Z43 . . . Z51 . . . Z52 . . . Z53 . . . Z61 . . . Z62 . . . Z63 . . . q2 q3 q4 q5 q6 (b) Second level map corre-Z1 Z 2 Z 3 Z42 Z45 Z57 Z61
				sponding to S (q 1 ,Z 1 ) . Each first
				level map node will have its
				own second level map.
				q1	Z1
	a1	.		q2	Z2
	a2	.		q3	Z32
	a3	.		q4	Z42	Z43
	a4	.	.	q5	Z53
	a5	.		q6	Z61
	(c) Pushes/Pops map cor-
	responding to root node
	(q1, Z1). Each pointer points
	to a different map where
	(q, Z) are stored.

Table 2 :

 2 Results on the Benchmarks. Time recorded in ms, and timeout (T.O.) used is 120 seconds. OoM stands for Out of Memory kill. Results rounded up to 1 decimal. # nodes refers to the number of nodes in the zone/region graph explored. In case of timeout ≥ n, refers to recorded number of nodes n before timeout occurred. NA in RB columns represents that the region based approach does not handle open guards in transitions (B 6 , B 7 have open guards.)

As before, there could be several such sets depending on the order in which the rules are applied.

This work was partly supported by ReLaX CNRS IRL 2000, DST/CEFIPRA/INRIA project EQuaVE and SERB Matrices grant MTR/2018/00074.


domain: scco spi

text:
Waldez A Gomes

Pauline Maurice

Serena Ivaldi

Dyadic Motor Behavior During Co-manipulation: A Study on Humans

Introduction Human-Robot Collaboration

• Control a robot to act efficiently with a human during a task that requires co-manipulation;

• Human and robot co-manipulating a shared object • Good collaboration is not only fast and/or accurate, but it is also legible

Results

Trajectory in Co-Manipulation

P1 Leader Natural P2 Leader P1 P2

• Trajectory seems to deviate towards the position of the participant that is the leader (Hypothesis 2)

• When the pipe is closer to the participant, it is harder for the participant to move the hand in different directions (low arm manipulability values are expected)

• Each participant may have its own desired trajectory, that is chosen according to several criteria, including arm manipulability.

0.05 0.10 0.15 0.20 0.25 0.30 0.35

RMS

Contact Sensors

• Less contacts when there is no pre-assigned leader (Hypothesis 3)

• No clear difference between conditions in which there is a leader

• Arm stiffness is modulated within the same participant for different conditions (Hypothesis 1)

• Trajectory deviation may be caused by the arm stiffness modulation.

• The ICC is an indirect measure of the arm stiffness, in here, calculated from a pair of forearm antagonist muscles EMG signals:

Future Work

• Rigorous analysis of the trajectories

• Define manipulability metrics for a dyad, and apply them to all dyads

• Use the stiffness modulation from the human arm to modulate the stiffness of the robot in a human-robot collaboration

•

  • Qualisys optical motion capture: 6 markers on each participant's right arm, 1 marker on the pipe → Arm kinematics → Pipe 3D trajectory • Wireless EMG Delsys Trigno on 3 pairs of antagonistic muscles: forearm, arm, shoulder → Muscle activation → Index of Co-Contraction (ICC) • Contact sensors on each tube (Quantitative measurement of accuracy) Task • Collaborative extraction and insertion of a pipe in a tube • Contact between the pipe and the tube should be execution are a responsibility of all agents•The task requires coordination and synchronization• Is there a Leader, or a Follower?• So we look into natural motor behavior in human-human dyads!Goal• Design a task to Investigate how human-human (HH) dyads motor behavior is influenced by different leadership conditions (leader or follower) during an object co-manipulation:• Hypothesis 1: Arm stiffness is influenced by leadership conditions;• Hypothesis 2: Object's trajectory is influenced by leadership conditions • Hypothesis 3: Accuracy of the task execution is influenced by leadership conditions• Investigate the same task, but with one human executing the task alone.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 731540 (An.Dy)


Norhane Benkahla
email: norhane.benkahla@supcom.tn

Hajer Tounsi
email: hajer.tounsi@supcom.tn

Ye-Qiong Song
email: ye-qiong.song@loria.fr

Mounir Frikha
email: m.frikha@supcom.tn

VHMM-based E-ADR for LoRaWAN networks with unknown mobility patterns

Keywords: IoT, LoRaWAN, ADR, Spreading factor, VHMM(2, 1), mobility

 aiming to maximize both battery life of the end-devices and overall network capacity. ADR performs adaptive tuning of radio configurations of the nodes by adjusting bandwidth, spreading factor, coding rate and transmission power parameters whenever the signal quality changes. The ADR algorithm was established for stable radio channel environments and is not efficient when conditions dramatically change (e.g. mobility). So, we have previously proposed an Enhanced-ADR (E-ADR) [2] that deals mobile nodes in case of predefined mobility patterns. However, several Internet Of Thing (IoT) applications, such as smart cattle ranching in smart farms [3], require sensors travelling with unknown or undefined trajectories. So, this paper extends E-ADR to unknown mobility pattern. This E-ADR extension, called VHMM-based E-ADR, is based on a Variable order Hidden Markov Model (VHMM) to predict the node trajectory. It has been implemented on Waspmote SX1272 hardware platform. Experimental results show its high efficiency in terms of the packet loss rate (PLR) and the energy consumption.

I. INTRODUCTION

Among the most popular communication technologies used for Internet of Things (IoT), Long Range (LoRa) technology [START_REF] Augustin | A study of LoRa: Long range low power networks for the Internet of Things[END_REF] and its Long Range Wide Area Network (LoRaWAN) MAC protocol [START_REF] Alliance | LoRaWAN Specification[END_REF] are widely used for long-range and low power sensor data transfer. IoT is applied in various domains where both static and mobile nodes are used, such as smart farms, habitat monitoring, logistics, and product life-cycle management.The main challenge for LoRaWAN is to ensure a good data transmission QoS (Quality of Service) while minimising nodes energy consumption. For this purpose, the adaptive data rate mechanism (ADR) is used to efficiently set the transmission parameters of each node according to the network state.

A LoRaWAN node can run with different data rates, resulting in different TOA (Time On Air), different power consumption, so, different use of duty cycle. This can be achieved by combining four parameters: Spreading factor (SF), Bandwidth (BW), Coding rate (CR), and Transmission power (TP). There are 6720 possible configurations by combining the different values of these parameters [START_REF] Bor | LoRa Transmission Parameter Selection[END_REF]. However, [START_REF]LoRa Alliance Technical committee[END_REF] was limited to only 8 data rates (DR0 -DR7) among these configurations. In our work we adopted the 10 modes defined in [START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF] and presented in Table I.

Table I shows that for a given TP, each mode defines the range of RSSI (Received Signal Strength Indicator) values [B i , B s ] ensuring the correct reception of the packets. An RSSI lower than B i may result in packet loss, whilst an RSSI higher than B s induces useless high power consumption.

How to dynamically choose the best data rate is not trivial. The ADR mechanism of LoRaWAN has been designed for a stable network with fixed nodes. Its convergence speed is relatively low and does not cope well with high dynamic networks [START_REF] Alliance | LoRaWAN Specification[END_REF]. Several studies and experiments revealed its limitation in mobile node scenarios [START_REF] Alvear | Assessing the Impact of Mobility on LoRa Communications[END_REF]- [START_REF] Kousias | Empirical Analysis of LoRaWAN Adaptive Data rate for Mobile Internet of Things Applications[END_REF] without offering solutions for these cases. So, to take into account mobile nodes, we have previously proposed in [START_REF] Benkahla | Enhanced ADR for LoRaWAN networks with mobility[END_REF] an Enhanced ADR (E-ADR). E-ADR aims at self-adapting the node mobility through quickly re-configuration of the parameters face to dynamic network condition changes. It determines the best configuration corresponding to the RSSI value of the sending node, estimated as a function of the its new position. The position calculation is based on the linear regression prediction method. Our experimental tests have confirmed that E-ADR outperforms Basic ADR thanks to its RSSI-based proactive transmission parameter adjustment, that can also use Right RSSI to prevent packet loss. However, E-ADR being designed for known mobility patterns is not suitable for applications where some nodes follow random/unknown trajectories (e.g., cattle ranching). To cope with applications with unknown trajectories, we propose in this paper to extend E-ADR to VHMM-based E-ADR which consists to predict the next position of the node using Variable order Hidden Markov Models (VHMM).

The main new contribution of this paper is:

• The exploration of HMM (Hiden Markov Model) to predict the next position of a mobile node and the validation of the optimal order (HMM-2) according to several experimentation. • The proposition of a variable VHMM-(2,1) model combining both HMM-1 and HMM-2 models • The integration of VHMM- [START_REF] Benkahla | Enhanced ADR for LoRaWAN networks with mobility[END_REF][START_REF] Alliance | LoRaWAN Specification[END_REF] to the E-ADR to adjust the LoRaWAN node configuration according to its new position. The rest of the paper is organised as follows. Section II describes the background of the work. VHMM-based E-ADR is detailed in Section III. Implementations and results are discussed in Section IV by showing how the proposal is adapted to mobile nodes travelling an unknown trajectory. Finally, we conclude in Section V.

II. BACKGROUND

The use of ADR (also called hereafter Basic ADR) is optional but highly recommended by the LoRaWAN specification [START_REF] Alliance | LoRaWAN Specification[END_REF]. ADR runs in two parts using two independent algorithms on the node-side and on the server-side. The nodeside reacts to the packet loss (no receipt of acknowledgements) by increasing the SF to increase its radio coverage range after reaching a retransmission threshold (m). The server-side aims to increase the DR (therefore decrease the SF) in order to save energy based on the best RSSI of the last 20 received packets. The main objective is to converge to an "optimal" SF that can ensure the packet reception with a high data rate on the one hand, and minimise the energy consumption on the other hand.

Analysis and implementations of the ADR have previously been undertaken. Studies in [START_REF] Alvear | Assessing the Impact of Mobility on LoRa Communications[END_REF]- [START_REF] Kousias | Empirical Analysis of LoRaWAN Adaptive Data rate for Mobile Internet of Things Applications[END_REF] have shown that the current ADR algorithm (Basic ADR) [START_REF] Alliance | LoRaWAN Specification[END_REF] is effective in stable channel conditions where nodes are static but is not suitable in dynamic networks, especially in the case of mobile nodes or the presence of transient obstacles between a node and its gateway. In addition, it has been identified that the performance of Basic ADR is limited in scale because of the duty cycle restriction (1% or not exceeding 36s TOA every cycle of one hour). Also, in a lossy channel, the nodeside algorithm convergence time to an optimal configuration mode (SF) is very slow (only after ADRACKLIM IT = 64 retransmission attempts). Moreover, tuning the transmission parameters based on the best RSSI of the 20 received packets does not lead to the optimal choice. Therefore, different ADR enhancements have been proposed in [START_REF] Bor | LoRa Transmission Parameter Selection[END_REF], [START_REF] Slabicki | Adaptive configuration of loRa networks for dense IoT deployments[END_REF]- [START_REF] Kim | Adaptive data rate control in low power wide area networks for long range IoT services[END_REF]. However, in those proposals, the server-side algorithm can only decrease the SF even in some cases (due to mobility) it should increase the SF to have a better coverage (Low acceptable RSSI range).

In a previous work [START_REF] Benkahla | Enhanced ADR for LoRaWAN networks with mobility[END_REF], we have proposed an Enhanced ADR (E-ADR) adapted to node mobility, which consists in fine-tuning the configuration by both decreasing or increasing the SF according to the next node position definition based on a predefined mobility model. Unlike Basic ADR and the other proposed enhancements in the literature that are reactive algorithms, E-ADR is a proactive algorithm characterised by a quick re-adaptivity to RSSI changes. In E-ADR, each node specifies its mobility pattern when joining the server by indicating its application ID in the join request message. Then, the node transmits n packets using mode1 in Table I. The server uses the trilateration method to calculate the node position for each received packet. Then, it calculates the average variation Avg variation (displacement) of the node [START_REF] Benkahla | Enhanced ADR for LoRaWAN networks with mobility[END_REF].

Based on this average variation, the server estimates the n next positions based on the Linear regression algorithm. When the position is estimated, the server calculates the distance separating the node and the gateways and deduces the estimated RSSI [START_REF] Goldoni | Experimental data set analysis of RSSI-based indoor and outdoor localization in LoRa networks[END_REF].

To decide which configuration the node will be allocated, the server calculates a transition rate R a→b , where a is the old mode and b is the next estimated mode, according to equation [START_REF] Alliance | LoRaWAN Specification[END_REF]. R a→b defines how much RSSI estim belongs to the interval [B i , B s ] of each configuration in {mode1, ...mode10} different to mode a (See Table I).

R a→b = 1 /2 × |Bs(b) -Bi(a)| |Bs(b) -RSSI (estim) | + |Bi(b) -RSSIestim| (1)
The greatest element of the vector corresponds to the new configuration mode.

Our previous experimental tests [START_REF] Benkahla | Enhanced ADR for LoRaWAN networks with mobility[END_REF] have shown that E-ADR outperforms Basic ADR and the other enhancements ( [START_REF] Bor | LoRa Transmission Parameter Selection[END_REF], [START_REF] Slabicki | Adaptive configuration of loRa networks for dense IoT deployments[END_REF]- [START_REF] Kim | Adaptive data rate control in low power wide area networks for long range IoT services[END_REF]) thanks to its RSSI-based proactive adjustment strategy. However, this E-ADR version does not cope with applications where nodes travelling through unknown trajectories. So, in what follows we propose an E-ADR extension to include unknown random mobility pattern. As a motivating example, we will consider the cattle monitoring application in a smart farm [START_REF] Semtech | Revolutionizing Smart Agriculture Using Semtech's LoRa Technology[END_REF].

III. VHMM(2,1)-BASED E-ADR

Predicting trajectories of moving objects has recently received increasing attention. In [START_REF] Zhang | Mobility Prediction : A Survey on State-of-the-Art Schemes and Future Applications[END_REF], the authors presented the hidden markov model (HMM) for the reservation of resources and for the choice of access points in the case of cellular mobile networks. HMM allows the mobile node to learn the environment and update the information itself, thus improving performance. The authors in [START_REF] Song | Evaluating next-cell predictors with extensive wi-fi mobility data[END_REF] compared four different families of location predictors that were tested on a set of data collected by Dart-mouth College from Wi-Fi users between April 2001 and March 2003. From the obtained results, the authors concede that the more complex predictors (such as LZ (Lempe and Ziv algorithm) [START_REF] Bhattacharya | An information-theoretic approach to track mobile users in PCS networks[END_REF], PPM (Prediction by Partial Matching) [START_REF] Bell | Text Compression[END_REF] and SPM (Sampled Pattern Matching) [START_REF] Jacquet | An universal predictor based on pattern matching,preliminary results[END_REF]) are not necessarily the more accurate than HMM predictors. They also established that HMM predictors beyond the second-order are less accurate.

In our work, we are interested in using HMM [START_REF] Zhang | Mobility Prediction : A Survey on State-of-the-Art Schemes and Future Applications[END_REF] for predicting the position of a mobile node by combining secondorder HMM (HMM-2) with first-order HMM (HMM-1). The The principle idea of our model is that it switches between two HMM orders (γ = 1, 2) to predict the next node position or point of interest (POI), on which the prediction of the next configuration will be based. In fact, as depicted in Fig. 1, our model named VHMM(2,1) acts primarily as an HMM-2. This latter switches to an HMM-1 when its outcome is an empty set (Failed match) in order to avoid a random prediction.

For -π = the probability of the initial state of trajectory, which is known in our case.

-A 1 = {A t,(t+1) }: A t,(t+1) =P(P OI (t+1) | P OI (t) ) is the probability that there is a 1-order sequence: it means that there is P OI (t) followed by P OI (t+1) in a training (observation) sequence.

-A 2 = {A (t-1)t,(t+1) }: A (t-1)t,(t+1) =P(P OI (t+1) | P OI (t) ,P OI (t-1) ) is the probability that there is a 2order sequence: it means that there is a sequence P OI (t-1) →POI (t) followed by a P OI (t+1) in a training (observation) sequence.

-

B 1 = {B t,(t+1) }: B t,(t+1) =P(O (t+1) (l) | P OI (t)
) is the probability of occurrence that defines how many times a P OI (t) is followed by each P OI (t+1) (l) , where l is one of the possible POIs (l ∈ {D, A, R, L, AL, P L, AR, P R}).

-B 2 = {B (t-1)t,(t+1) }: B (t-1)t,(t+1) =P(O (t+1) (l) | P OI (t) , P OI (t-1) ) is the probability of occurrence that defines how many times a sequence P OI (t-1) → P OI (t) is followed by each P OI (t+1) (l) , where l is one of the possible POIs (l ∈ {D, A, R, L, AL, P L, AR, P R}).

The training trajectories of the node's movement in the observation bloc (Fig. 1) are divided into groups where each group is composed of 3 POIs in case of HMM-2 and 2 POIs in case of HMM-1. Then, we define the probability that a P OI (t+1) (l) (l is the POI number) is the succession of the sequence P OI (t-1) → P OI (t) in the case of HMM-2 and the succession of only P OI (t) in the case of HMM-1.

In case the probability P (P OI (t+1) (l) ) presented in Definition 2 corresponds to an empty set (∅) -meaning that there is no sequence corresponding to [P OI (t-1) → P OI (t) → P (P OI (t+1) (l) )] in the training sets-, the server switches to the HMM-1 model and will be based on the probability of P (P OI (t+1) (l) ) presented in Definition 3.

Definition 2: The probability that (P OI (t+1) (l) ) is the succession of P OI (t-1) → P OI (t) (γ=2) is: Definition 3: The probability that (P OI (t+1) (l) ) is the succession of P OI (t) (γ=1) is:

P (P OI (t+1) (l) ) = B (t-
P (P OI (t+1) (l) ) = B t,(t+1) A t,(t+1) = P (O (t+1) (l) |P OI (t) ) P (P OI (t+1) |P OI (t) )
The greatest probability P (P OI (t+1) (l) ) will be considered and the corresponding l (l ∈ {D, A, R, L, AL, P L, AR, P R}) will be the predicted POI. Based on this POI, the server estimates the next configuration mode using [START_REF] Alliance | LoRaWAN Specification[END_REF]. Because the node travels an unknown trajectory, we propose to update the configuration after each n = 1 received packet.

IV. EXPERIMENTAL EVALUATION

To evaluate our proposed VHMM-based E-ADR, we have conducted several experiments using Waspmote SX1272 [START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF], STM32 nodes [START_REF] Corporation | Discovery kit for LoRaWAN, Sigfox, and LPWAN protocols with STM32L0 B-L072[END_REF] and gateways considering a smart farm application. In Waspmote SX1272, there are 10 predefined modes (Table I) that will be used in our experimental tests. It is however worth noting that we can choose any parameter combinations without being constrained to only the 10 predefined modes in Table I.

A. Smart farm: A motivating application

As a motivating domain, in this paper, we emulate a smart farm scenario and we consider a cattle monitoring application. Let's consider a LoRaWAN network presented in Fig. 2. It is composed of 5 sensors carried by 5 animals and dedicated for monitoring and tracking, and 3 gateways (G a , G b and G c ). We suppose that these animals are following random trajectories in an open field S: S=([X min , X max ] = [-120m, 140m], [Y min , Y max ] = [0, 140m]). In order to emulate these animals, nodes are manually handled and moved at an average speed of 9m/min.

We suppose that 3 sensors (D 2 , D 3 and D 4 ) are sending 21 packets, D 1 is sending 36 packets and D 5 28 packets. Each packet contains 25Bytes and all nodes send packets every 2 minutes (P sending (i) ). As presented in Fig. 2, we consider the previously presented L = 8 possible states (POIs) for our VHMM model. During the test phase, the reported results represent the average of 5 tests for each node. We will predict the next POI (P OI (t+1) ) of each moving node based on the VHMM(2,1) model and the trajectory set experienced during the training phase.

In our experiments, we define α i,j in (2) as the maximum number of allowed retransmissions of the i th packet before the next scheduled packet (i + 1) should be sent by D j . In (2), T delay (i, j) corresponds to the transmission delay of packet (i) sent by D j at instant T i and represents the amount of time required to receive successfully the i th packet (including retransmissions and the reception windows). T delay (i, j) can very depending on the used configuration mode and the number of performed retransmissions. α i,j depends also on the sending interval P sending (j) (fixed to 1pkt/2min in our scenario) allowing retransmissions to take place only during [T i , T i+1 ) (which is an interval of 2min, in our case).

α i,j = P sending (j) T delay (i, j) (2)
During these α i,j retransmissions, if the node does not receive an acknowledgement after a threshold m (m < α i,j ), it will increment its SF to increase the coverage area and by thus the probability to reach the gateway. In [?], the threshold m is set to 64 by default. In addition, during our experiments, if a packet i sent at T i is still not acknowledged at instant T i+1 (after α i,j retransmission attempts), it will be considered as "lost" because of bad SF allocation. Moreover, when D j exceeds the duty cycle limitation of 36s per hour due to a high number of retransmissions, the non-transmitted packet will be discarded.

During the experiments, the n = 2 first transmissions are established using M ode1 to guarantee the history of 2 POIs to perform the HMM-2. Moreover, mode1 is chosen to ensure to reach the gateways. Then, each process assigns a parameter combination according to the used strategy. The configuration will be updated after every n = 1.

During the n = 2 first transmissions, we recorded all the RSSI values. For D 1 , for example, we recorded an RSSI values of -112.62dB and -115.7dB From these RSSI values, the server allocates the configuration mode to be used for the next transmission. According to the received RSSI values, using Basic ADR, the server will allocate M ode10 for D 1 , for example. However, using VHMM-based E-ADR, based on the next predicted POI for the configuration allocation, the server will allocate mode8 for D 1 in order to ensure the needed radio coverage range.

Firstly, we compare the VHMM(2,1) model to (HMM-1, HMM-2, and HMM-3) to validate its prediction efficiency by considering the prediction accuracy criterion defined in (3) as the rate of the number of correctly predicted POIs over the total number of travelled POIs.

Accuracy = N P OI correct /N P OI travelled

Hereafter, the VHMM-based E-ADR is compared to Basic ADR in terms of PLR and energy consumption. In addition, as retransmission has a direct impact on PLR and energy consumption, we will also evaluate the impact of the retransmission threshold m.

C. Evaluation of the prediction model

In this section, we evaluate the prediction accuracy of VHMM(2,1). We have firstly tested multiple random trajectories using (HMM-1, HMM-2, HMM-3 and VHMM(2,1)) models.

Fig. 3 presents the obtained predicted trajectories using each model for one of these random trajectories (using D 2 ). According to Fig. 3, the VHMM(2,1) model results are better than those obtained using the HMM-1, HMM-2, and HMM-3 models. We note that among the 3 models, HMM-2 is positioned after VHMM(2,1), which is closer to the real covered pattern compared to HMM-1 and HMM-3.

We have registered the accuracy ratio for each tested trajectory. Table II presents the average accuracy prediction using the different HMM models for these tests. The POI is considered to be correctly predicted after checking the real covered trajectory.

In fact, HMM-1 is based only on the current POI to predict the next one, giving a less accurate prediction compared to VHMM(2,1), or HMM-2 which are based on 2 historical POIs to predict the next position (Table II). Furthermore, the HMM-3 model is the less efficient compared to all other models since it is based on 3 historical POIs requiring more learning in the training phase. Moreover, the HMM predictors beyond the second-order are in general less accurate [START_REF] Song | Evaluating next-cell predictors with extensive wi-fi mobility data[END_REF].

In what follows, we evaluate our VHMM-based E-ADR proposal and compare it to Basic ADR in terms of PLR and energy consumption. To fairly compare them by focusing on the algorithmic efficiency of the strategy without the influence of the algorithm execution frequency, we execute the ADR update every n = 1 received packet instead of every n = 20 packets for Basic ADR.

D. PLR evaluation

We define PLR as the number of the lost packets (due to bad SF or a lot of interference) and non-transmitted packets (due to cycle limitation) over the total planned number of packets that a node D j should send during the experiment. Table III presents the average PLRs of 5 measurements for each node D j using the VHMM-based E-ADR and Basic ADR.

We observe that the nodes experience a high PLR using Basic ADR. In fact, as the server-side in Basic ADR is responsible of only the SF decrease (only the node-side can increase it after a m = 64 retransmissions), even if the RSSI value indicates the need of higher SF, the server can not respond to the node requirement inducing by thus an important PLR.

Because of high number of retransmissions, all nodes suffer from a PLR due to duty cycle exceeding in addition to the bad SF allocation. So, to better understand when and how the packets are lost (bad SF) or non-transmitted (duty cycle exceeding), Fig. 4 traces the PLR evolution in time of D 2 , by defining separately the number of lost packets (bad SF) over the number of packets already transmitted (blue curves) and the number of non-transmitted packets (exceeding the duty cycle limitation) over the number of the already transmitted ones (orange curves).

We notice in Fig. 4 that during the transmission of the first packets (packets from 1 to 14), the duty cycle is still valid and Fig. 4: PLR evolution for D 2 loss is rather due to the bad choice of SF (blue curve). In fact, when nodes are close to the gateways and begin to move away from them, high RSSIs are in the history list, resulting in the attribution of low SFs, inducing an under-optimal allocation, so, a high number of loss. Thus, the retransmission of these packets on the one hand and the sending of the rest of the packets on the other hand will cause the duty cycle to be exceeded and therefore the increase in packet loss (orange curve).

However, using the VHMM-based E-ADR, the packet loss due to bad SF is decreased thanks to the pro-activity of the strategy based on the prediction of the future node locations, inducing less retransmissions and therefore the risk of exceeding the duty cycle.

So, although retransmissions give packets better chance to be successfully received, this does not necessarily mean low PLR, since these retransmissions lead to duty cycle exceeding inducing high PLR because of the not transmitted (discarded) packets.

The VHMM-based E-ADR model allows an optimal configuration allocation allowing to reduce the PLR. In some cases (D 1 , D 2 , and D 5 for example), nodes may loss a few packets because of the inaccurate POI prediction.

E. Energy consumption evaluation

Retransmission in LoRaWAN increases ToA and so the energy consumption which is a crucial factor for IoT devices. In the following, we will evaluate the energy consumption that is calculated according to [START_REF] Bouguera | Energy Consumption Model for SensorNodes Based on LoRa and LoRaWAN[END_REF] for Basic ADR and VHMMbased E-ADR.

Table IV presents the average energy consumption per received packet of each node D j using both Basic ADR and VHMM-based E-ADR. These values are defined as: the total consumed energy over the number of received packets. By allocating a bad SF, nodes using Basic ADR are faced a large number of retransmissions and therefore an increase in energy consumption.

F. Impact of retransmission threshold (m)

As we have noticed, retransmissions may sometimes decrease the PLR at the expense of increasing energy consumption in transitive network condition changes. However, with m = 64 by default, a high number of retransmissions may also lead the node to exceed its authorised duty cycle (increasing PLR) even before the node can reach a suitable SF by linearly incrementing it every m retransmissions. Table V presents the impact of m on respectively the average energy per received packet and the PLR by varying m (m = {8, 64}).

The obtained results in Table V show that a low m value (m = 8) allows more chance to reach quickly the optimal SF, as the node increments its SF every m retransmissions, avoiding a high number of retransmissions, a duty cycle exceeding, an important PLR and high energy consumption. However, even if with a small retransmission threshold (eg., m = 8) the Basic ADR still has higher power consumption compared to VHMM-based E-ADR, since it still has recourse to retransmissions. Furthermore, using VHMM-based E-ADR, nodes did not have to re-transmit packets thanks to the optimal allocation using the prediction of the next POI based on the VHMM(2,1) model. So, VHMM-based E-ADR offers a good trade-off between energy consumption and PLR.

V. CONCLUSION

ADR may be used to re-adapt the node configuration in LoRaWAN networks, but it is facing to low adaptation speed and low performance issues in case of mobility. VHMM-based E-ADR solves those issues in case of mobility with undefined trajectories. We have shown through our experimental tests that our proposed strategy based on the prediction of the device displacement outperforms Basic ADR by offering a fair trade-off between PLR and energy consumption. Moreover, by reducing PLR, the number of retransmissions is minimised, and thus the probability to exceed the duty cycle limitation. Our experiments have shown that our prediction model "VHMM(2,1)" presents 87% accuracy. In Future work, we will be interested to improve that by increasing the number of experimented devices and by exploring other ML techniques. Furthermore, we will study the scalability of our proposed mechanism in a large network.

Fig. 1 :

 1 Fig. 1: VHMM(2,1) prediction model

  our model, we consider l = 8 possible states (POIs): Direct (D), Back (A), Right (R), Left (L), Left Lean back (AL), Left Lean forward (PL), Right Lean back (AR), Right Lean forward (PR) (l = {D, A, R, L, AL, P L, AR, P R}) . The relevant definitions of VHMM(2,1) are as follows. Definition 1: Variable HMM-(2,1) uses λ = (π, A 1 , A 2 , B 1 , B 2 ) to describe a model. Supposing that POI{} (POI ∈ l) represents a sequence state of node positions and O{} (O ∈ l) represents the observation state sequence. VHMM-(2,1) is characterised by the following notations:

Fig. 2 :

 2 Fig. 2: Test scenario

Fig. 3 :

 3 Fig. 3: Real trajectory Vs. HMM-1 Vs. HMM-2 Vs. HMM-3 Vs. VHMM(2,1)

TABLE I :

 I Configuration modes for TP=14dBm[START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF]

	mode	BW (kHz)	CR	SF	RSSI(dB) [Bi, Bs]	DR
						(Kbps)
	1	125	4/5	12	[-134, -131]	0,293
	2	250	4/5	12	[-131, -129]	0,585
	3	125	4/5	10	[-129, -128]	0,976
	4	500	4/5	12	[-128, -126]	1,718
	5	250	4/5	10	[-126, -125.5]	1,953
	6	500	4/5	11	[-125.5, -123]	2,148
	7	250	4/5	9	[-123, -120]	3,515
	8	500	4/5	9	[-120, -117]	7,031
	9	500	4/5	8	[-117, -114]	12,50
	10	500	4/5	7	[≥ -114]	21,875

  |P OI (t-1) , P OI (t) ) P (P OI (t+1) |P OI (t-1) , P OI (t) )

	1)t,(t+1) A (t-1)t,(t+1)	=	P (O (t+1) (l)

TABLE II :

 II Prediction accuracy

	Node	VHMM(2,1)	HMM-1	HMM-2	HMM-3
	D 1	86.08%	54.87%	69.71%	42.25%
	D 2	89.14%	57.12%	71.48%	47.51%
	D 3	88.52%	55.58%	65.89%	44.47%
	D 4	88.09%	52.38%	70.35%	48.02%
	D 5	87.15%	54.24%	67.19%	47.87%

TABLE III :

 III Average PLR of D j

	Node	Basic ADR	VHMM-based E-ADR
	D 1	41.66%	5.55%
	D 2	42.85%	14.28%
	D 3	52.38%	0%
	D 4	38.09%	0%
	D 5	42.85%	7.14%

TABLE IV :

 IV Average energy per received packet (J)

	Node	Basic ADR	VHMM-based E-ADR
	D 1	0.97	0.341
	D 2	1.7	0.73
	D 3	2.04	0.26
	D 4	1.56	0.26
	D 5	1.27	0.34

TABLE V :

 V m impact on energy and PLR

		Impact of m on energy (J)
	Node	Basic ADR	VHMM-based E-ADR
		m = 8	m = 64	m = 8	m = 64
	D 1	0.63	0.97	0.28	0.341
	D 2	0.73	1.7	0.58	0.73
	D 3	1.05	2.04	0.26	0.26
	D 4	0.75	1.56	0.26	0.26
	D 5	0.94	1.27	0.29	0.34
			Impact of m on PLR
	Node	Basic ADR	VHMM-based E-ADR
		m = 8	m = 64	m = 8	m = 64
	D 1	11.11%	41.66%	0%	5.55%
	D 2	14.28%	42.85%	0%	14.28%
	D 3	19.04%	52.38%	0%	0%
	D 4	4.7%	38.09%	0%	0%
	D 5	25%	42.85%	0%	7.14%