Patricia Bouyer
email: bouyer@lsv.ens-cachan.fr

Uli Fahrenberg

Kim G Larsen

Nicolas Markey
email: markey@lsv.ens-cachan.fr

Quantitative analysis of real-time systems using priced timed automata

The problems of time-dependent behavior in general, and dynamic resource allocation in particular, pervade many aspects of modern life. Prominent examples range from reliability and efficient use of communication resources in a telecommunication network to the allocation of tracks in a continental railway network, from scheduling the usage of computational resources on a chip for durations of nano-seconds to the weekly, monthly or longer-range reactive planning in a factory or a supply chain.

These problems have been subject to substantial research for decades by different communities such as operational research, computer systems performance evaluation as well as planning and scheduling, witnessed by large ACM communities such as SIGMET-RICS and PERFORMANCE. In this paper we argue that the formalism of timed automata together with recent extensions provides an alternative framework with complementary, yet competitive, results in terms of modeling capabilities and efficiency of analysis. Timing: Twenty years ago, R. Alur and D. Dill introduced the notion of timed automata. As a witness for the importance of the formalism one may consider the 2008 Computer-Aided Verification Award given to Alur and Dill for their seminal 1990 article Automata for mod-eling real-time systems, [START_REF] Alur | Automata for modeling real-time systems[END_REF] which provided the theoretical foundation for the computer-aided verification of real-time systems.

Real-time systems and resource allocation problems have manifested themselves under different names in application domains such as manufacturing, transport, communication networks, embedded systems, and digital circuits, and have been treated using theories and methods in several disciplines. Most of these applications involve distributed, reactive systems of considerable complexity, and with a number of real-time constraints in the sense that correctness not only depends on the logical ordering of events of the systems, but also on the relative timing between these.

State-based models have been the basis of a wide range of successful computer-supported verification methodologies allowing the efficient prediction of functional properties, e.g. absence of deadlock or memory overflow. However, many of the models used in this methodology are purely discrete and their treatment of time is purely qualitative, that is, behaviors are just sequences of events appearing one after the other but without any quantitative timing information about the duration of actions and the time between events. Timed automata allow such timing constraints to be expressed, while being amenable to computer-aided analysis methods such as simulation, verification, optimization and controller synthesis. Performance: In all of the above applications, an explicit constraint on timing is only one of a number of quantitative aspects of importance. Within embedded systems additional key quantities include energy and memory consumption, in communication networks required band-width is a key quantity, and within the factory and supply chain applications need for storage and overall cost for a given production are crucial quantities. The extended notion of priced or weighted timed automata has been put forward as a formalism allowing for such additional and timedependent quantities to be modeled, without hampering efficient analysis and even permitting optimization. Uncertainty: Classical models for scheduling in manufacturing, such as job-shop problems, are somewhat detached from industrial practice and reality. They assume that the duration of every step as well as the arrival times are fixed and known with certainty; in practice however, it is rarely the case that a schedule is executed as planned.

For solving problems related to expected time and performance properties, stochastic process models have been very successful. Several refinements of a model (a) of the working mathematician according to Erdős: After insertion of a coin into the coffee dispenser, coffee can be collected, and the scientist can go back to work. In the timed-automaton model (b), precisely five time units pass between coin insertion and coffee collection, and the time which passes between coin insertion and going back to work is less than ten time units. In the priced timed automaton (c), cost rates (modeling e.g. energy consumption) are associated with the three states. In the timed game (d), uncertainty as to precisely when coffee is delivered is modeled as an uncontrolled edge.

When aiming at guaranteed time and performance properties under uncertainty, so-called timed games may be used instead. They provide efficient off-line algorithms for synthesizing reactive schedulers with performance guarantees. Such algorithms can plan for the best or worst case, but the scheduling strategies they produce are adaptive and can take advantage, for example, of the fact that a task has terminated before it was expected to.

In this paper we present the formalism of timed automata and its priced and game extension as a unifying mathematical framework for the modeling, analysis, optimization and synthesis of real-time related phenomena. Figure 1 shows some simple examples of these formalisms; later we provide more elaborate and realistic examples and case studies.

TIMED AUTOMATA 1.1 A model for time

Timed automata 5 are a powerful model for representing and reasoning about systems where the notion of time is essential. They are an extension of classical finite-state automata with real-valued variables called clocks. These clocks all increase at the same rate, and their values can be used to restrict availability of transitions and how long one can stay in a location (or state). Also, clocks can be reset to zero when a transition is taken. To this end, each transition has associated with it a guard (which must be satisfied for the transition to be enabled) and a set of clocks to be reset, and each location carries an invariant which must be continuously satisfied when the system is in the location. Below we show an example of a timed automaton with two clocks x and y, and label set {a, b, c, d, e}. Note that no time can elapse in location ℓ1 due to the invariant (y = 0); locations with this property are called urgent.

ℓ 0 ℓ 1 (y=0) ℓ 2 ℓ 3 x≤2,a,y:=0 b c x = 2 ,d x = 2 ,e
Guards and invariants are given as comparisons x ≤ α or x < α, or the reverse relations, of a clock value with an integer constant, or as conjunctions of these. Sometimes also so-called diagonal constraints xy ≤ α (or < or other) are allowed, but other extensions quickly lead to undecidability issues, see below.

A configuration of the system is made of a location and a clock valuation (in our case, values for both clocks x and y). A possible execution in our example is:

  ℓ0 0 0   delay -→ 1.3   ℓ0 1.3 1.3   action -→ a   ℓ1 1.3 0   action -→ c   ℓ3 1.3 0   delay -→ 0.7   ℓ3 2 0.7   action -→ e   2 0.7  
where the first component of a configuration is the location and the second and third components give the values of clocks x and y, respectively. This execution corresponds to a delay of 1.3 time units in ℓ0, the firing of transition a (which is enabled because the value of clock x is less than 2; clock y is then set to 0), the firing of transition c (which occurs without delay as ℓ1 is urgent), etc.

In the context of verification, several problems are of interest, like the model-checking of safety properties ("Can a distinguished set of states be avoided?"), reachability/liveness properties ("Can/will a distinguished goal state be reached?"), or more involved properties such as response properties ("Is any request eventually granted?"). As a model for real-time systems, these properties can include quantitative constraints, for instance time-bounded reachability, or time-bounded response properties ("Is any request granted within two minutes?"). It is also relevant to compute optimal time bounds for these properties, e.g. optimal-time reachability ("What is the minimum time required for reaching a distinguished set of states?").

The region abstraction

A timed automaton is a syntactical representation of an infinite transition system, since clocks take (nonnegative) real values. However, there is a way to deal with this infinity of configurations by reasoning symbolically: the main theoretical ingredient for solving problems on timed automata is the notion of regions, [START_REF] Alur | Automata for modeling real-time systems[END_REF] which provide a finite partitioning of the state space such that states within a given region are bisimilar, i.e. behaviorally indistinguishable.

The precise definition of regions is such that inside a region, integral parts of clock values do not change, and also the ordering of clocks according to their values' fractional parts stays the same. Special consideration has to be given to the cases where one or more clock values are integers, and finiteness of the region partitioning is ensured by considering as equivalent all clock values which exceed the maximal constant appearing in guards and invariants of the timed automaton in question. In the left part of Figure 2 we show the fortyfour regions for two clocks x and y with maximal constant equal to 2. In this two-clock case, regions can be points From two equivalent configurations (same location, region equivalent valuations), by delaying or by taking a transition, similar regions will be visited and similar behaviors will be possible. Regions are thus a way to finitely abstract the behaviors of a timed automaton. There are finitely many regions, and by considering as abstract configurations pairs of locations and regions, we get a finite automaton, called the region automaton, which preserves many properties including reachability, liveness and safety. Hence verification of those properties on the original timed automaton can be transferred to the finite region automaton and then checked using standard algorithms.

The limits of the region abstraction

Not all properties can be decided on timed automata using the region abstraction, and problems such as checking inclusion ("Are all real-time behaviors of a timed automaton also behaviors of another timed automaton?") and universality ("Can all real-time behaviors be realized in a given timed automaton?") are undecidable.

Also, the set of real-time behaviors exhibited by timed automata is not closed under complement, and not all timed automata are determinizable. As a counterexample for these properties, one can use the following timed automaton:

ℓ 0 ℓ 1 a a a a,x:=0 x=1,a
It accepts all behaviors with at least two a's separated by one time unit. It can be shown that no deterministic timed automaton exists with exactly the same behaviors, and also that no timed automaton can implement precisely all complementary behaviors. x:=0

ℓ 1 ℓ 2 y<1, x:=0 Figure 2:
The region abstraction is a finite representation of all possible behaviors of the timed automaton. Consider the timed automaton on top of the picture, and assume we enter location ℓ1 with clock values (x0, y0) for which 0 < y0 < x0 < 1 (a point in the red triangle, see the picture on the left); as clock y has value strictly less than 1, we have the option to switch to location ℓ2, which would reset clock x and end up in the purple region. We also have the option to delay in ℓ1; in that case, we exit the red triangle and reach the orange line. Here again we have two options: switching to ℓ2, or delaying to the yellow clock region. In case we still decide to wait in that region, we reach the green line. From that region on, the transition to ℓ2 is not enabled anymore. This description of the possible behaviors starting from the red region, which has been represented on the picture to the right, does not depend on the precise values of the clocks: region equivalence preserves enough information to encode exactly the behaviors of the underlying timed automaton.

Timed automata in practice

The region abstraction is a powerful tool for showing decidability of a number of interesting properties, but unfortunately, the region-based verification algorithms introduced above are infeasible in practice. The number of regions grows exponentially with the number of clocks, hence these algorithms have exponential time complexity. Using onthe-fly techniques, one can reduce the complexity to polynomial-space, and indeed it can be shown 4 that e.g. the reachability problem is PSPACE-complete.

Algorithms which have shown to be feasible, even efficient, in practice are based on the so-called zone graph abstraction: [START_REF] Th | Symbolic modelchecking for real-time systems[END_REF] A zone is a set of clock valuations defined by a clock constraint and can hence be represented by such; the zone graph has as vertices pairs of locations and zones which satisfy the location's invariant, and its edges are derived from the transitions of the given timed automaton. The number of zones is unbounded, so unlike the region graph, the zone graph is infinite. Finiteness can be enforced using a technique known as normalization; [START_REF] Bouyer | Forward analysis of updatable timed automata[END_REF] however the number of zones is still much larger than the number of regions, and moreover the same zone can be represented using many different clock constraints.

The reason for zone-based algorithms to be efficient in practice is two-fold: First, the algorithms used have no need to explore all of the zone graph (they work on-the-fly), and zones are commonly bigger than regions, hence the part of the zone graph to be explored is smaller. Second, operations on zones can be implemented very efficiently (in time cubic in the number of clocks): Zones are usually rep-resented using difference-bound matrices, or DBMs. The DBM representation of a zone on a set of k clocks has (k + 1) × (k + 1) entries, where an entry ci,j represents a clock constraint xixj ≤ ci,j and an extra clock x0 is added to represent absolute clock constraints xi ≤ ci0. DBMs in turn can be represented as directed weighted graphs; see below for an example of a zone and its DBM (graph) representation. Canonical representations of zones can be obtained using shortest-path closure or shortest-path reduction of their DBM graphs, and delay and reset operations on zones can be efficiently implemented on the DBM representations.

Z =                    x 1 ≤ 3 x 1 -x 2 ≤ 10 x 1 -x 2 ≥ 4 x 1 -x 3 ≤ 2 x 3 -x 2 ≤ 2 x 3 ≥ -5 x 0 x 3 x 1 x 2 3 5 2 2 10 -4

Task graph scheduling: time optimality

A task graph problem involves a number of tasks T1, ..., Tm, a number of machines or processors P1, ..., Pn, and a (partial) mapping d giving, for each task Ti and processor Pj, the time d(i, j) for computing Ti on Pj. In addition there is a partial order on the tasks used for describing dependencies. Figure 3 is an example of a task graph problem.

We want to determine a schedule of when to start the execution of tasks, and on which processors, that minimizes the total execution time while being feasible in respecting the following conditions: (a) a task can be executed only if all its predecessors have completed; (b) each machine can process at most one task at a time; (c) tasks cannot be pre-empted.

Task graph scheduling problems may be easily modeled as networks of timed automata so that every run corresponds to a feasible schedule and the fastest run gives the timeoptimal schedule: For each processor we construct a small timed automaton ablewhen idle-to handle within the appropriate amount of time the requests from the tasks. For the processors of Figure 3, these are as follows:

P1: idle + (x≤2) × (x≤3) x:=0 add 1 x:=0 mult 1 x=2 done 1 x=3 done 1 P2: idle + (y≤5) × (y≤7) x:=0 add 2 x:=0 mult 2 y=5 done 2 y=7 done 2
Each task is modeled as a timed automaton waiting to be served by either of the processors, conditioned by the completion of its predecessors (indicated by Boolean variables t1 through t5). Tasks T4 and T5 of our example can be represented as follows: T4:

t 1 ∧t 2 add i t 4 :=1 done i T5: t 3 add i t 5 :=1 done i
Extensive experiments on benchmarks have demonstrated that the above timed automata approach to task graph scheduling is competitive compared with more Compute D×(C×(A+B))+(A+B)+(C×D) using two processors: Given the execution platform with two processors, P1 and P2, and corresponding computation times for addition and multiplication, as well as their energy consumption, Sch1 to Sch3 provide three feasible schedules, where Sch2 is in fact time-optimal, and Sch3 is energy-optimal.

traditional approaches from operations research (e.g. mixed-integer linear programming) as well as specialized, heuristic algorithms from planning and scheduling. [START_REF] Abdeddaïm | Scheduling with timed automata[END_REF] Furthermore the generic approach of timed automata admits easy incorporation of more specialized features (e.g. release times, deadlines) to the models and scheduling.

Extensions of timed automata

As we shall see in Section 4, timed automata are a rich extension of classical automata with efficient tool support and several successful industrial applications. As such, they are often cited as the model of choice for representing and reasoning about embedded and real-time systems.

This success has led to several extensions of the model, for instance with more general guards or resets being allowed (e.g. additive guards [START_REF] Bérard | Timed automata and additive clock constraints[END_REF] or non-deterministic updates of clocks [START_REF] Bouyer | Forward analysis of updatable timed automata[END_REF] ), or with more involved dynamics measuring other quantities than time. Unfortunately, these extensions quickly lead to undecidability; e.g. for timed automata in which clocks can be stopped (so-called stopwatch automata), even basic properties such as safety or liveness are undecidable. [START_REF] Th | What's decidable about hybrid automata[END_REF] On the other hand, the model of hybrid automata [START_REF] Th | What's decidable about hybrid automata[END_REF] though suffering from the same undecidability problems as mentioned for other classes above, has emerged as a popular formalism for which semi-decision and approximation procedures have been developed. The model of priced timed automata, which we shall discuss in the next section, form an intermediate class between timed and hybrid automata for which some of the good decidability properties of timed automata are retained. Other intermediate classes of models have been investigated, including linear hybrid automata [START_REF] Th | What's decidable about hybrid automata[END_REF] and integration graphs, [START_REF] Kesten | Decidable integration graphs[END_REF] providing semi-decidability in general and decidability under certain restrictions.

PRICED TIMED AU-TOMATA 2.1 A model for resources

Time is not the only quantitative notion of interest when designing embedded systems; other quantities such as energy or memory consumption, required bandwidth, or accumulated cost can be important to measure in such systems.

These notions are intimately connected to time, because the longer the device is operating, the more resources it consumes. This makes timed automata the model of choice to reason about those quantities, and has led to the definition of priced timed automata, [START_REF] Alur | Optimal paths in weighted timed automata[END_REF][START_REF] Behrmann | Minimum-cost reachability for priced timed automata[END_REF] extending timed automata with cost (which is the general name we will use in the sequel to refer to the various quantities which can be modeled within this formalism; in some other literature, this is referred to as reward).

A priced timed automaton is hence a timed automaton with extra information indicating how the cost is evolving in locations and during transitions. To avoid the undecidability problems of hybrid automata, cost information cannot be used to guard transitions; the cost is only an observer variable, and whether a transition is enabled only depends on timing information, not cost value. An example of a priced timed automaton, extending the timed automaton of the previous section, is depicted below (labels omitted):

ℓ 0 +5 ℓ 1 (y=0) ℓ 2 +10 ℓ 3 +1 x≤2,y:=0 x = 2 + 1 x = 2 + 7
A decoration +10 on a location indicates that cost increases by 10 units per time unit in the location; a decoration +7 on a transition indicates that taking the transition increases The corner-point abstraction refines the region abstraction by also keeping track of the corner point close to which an execution runs. This is needed to measure costs: for instance, if we are in location ℓ1 and in the red region where 0 < y < x < 1, the price of delaying depends on the value of the clocks. From (a), where both x and y are arbitrarily close to 0, we can let almost one time-unit elapse and reach (b). The resulting cost is arbitrarily close to +3. On the other hand, from (c), where x and y are arbitrarily close to 1 and 0, respectively, letting time elapse takes us to the subsequent region, so that the cost is arbitrarily close to 0. (Notice that for readability, some resetting transitions have been omitted.) overall cost by 7 units (locations and transitions without cost indication have cost 0). The executions of such an automaton are those of the underlying timed automaton. The total cost of the example execution given in Section 1.1 (delaying 1.3 time units in ℓ0, 0.7 time units in ℓ3, and ending in the rightmost location) can be computed as

1.3 × (+5) + 0.7 × (+1) + 7 = 14.2

Optimizing the resources

Natural optimization questions can be posed on that model, e.g. the optimal reachability problem (minimum cost for reaching a given goal), the mean-cost optimization problem (mean cost used in the long run), or the discounted-cost optimization problem (where costs are discounted exponentially as time elapses).

As an example, we compute the minimum cost that is required for reaching location in the previous example. There are two families of executions: those that go through ℓ2 and those that go through ℓ3. Furthermore, in each family, there is a single parameter t: the time elapsed in location ℓ0; everything else is determined by the guards in the automaton. Hence the minimum cost is:

inf 0≤t≤2 min 5t+10(2-t)+1 5t+(2-t)+7 = 9 (1)
where the expressions 5t + 10(2t) + 1 and 5t + (2t) + 7 give the cost of executions going through ℓ2 respectively ℓ3 after delaying t time units in location ℓ0.

The standard region construction is not accurate enough to properly keep track of cost information, and a refinement of the region abstraction, the corner-point abstraction, [START_REF] Bouyer | Optimal infinite scheduling for multi-priced timed automata[END_REF] has to be used to solve the optimiza-tion problems mentioned above. For this abstraction, regions are refined by distinguishing their corner points. As an example, the two-dimensional region depicted below is refined into three region-corner pairs; the meaning of a region-corner pair is that the current clock valuation is arbitrarily close to the distinguished corner:

(a) (b) (c)
Similar to the refinement of regions, the transitions in the region automaton have to be refined to keep track of the corners. In the example above, there is a (delay) transition from region-corner pair (a) to (b), whereas (c) cannot be reached neither from (a) nor from (b). Figure 4 illustrates the corner-point abstraction of an example priced timed automaton. This graph has two types of delay edges: either within a region, from one corner to another one, or from a corner of a region to the corresponding corner in the subsequent region. The first case corresponds to a delay of "almost" one time unit, while the second case corresponds to a delay of "almost" zero time units. In addition, there are edges representing transitions of the timed automaton (which reset clock x in our example of Figure 4). In that case as well, there is a natural mapping between corners.

The edges of the corner-point abstraction are labeled with discrete cost information: if the cost rate in the current location is +3, all one time-unit edges have label +3, and all zero time-unit edges get label 0. Edges coming from discrete transitions are labeled with the cost of the transition (+5 in the example).

The corner-point abstraction can be used to solve many optimization problems, as it can be shown that in these cases, optimal total cost is obtained for runs which always take transitions close to integer clock values. Hence the optimization problem reduces to a problem on a finite graph which can be solved using different standard techniques. This is the case for the mean-cost optimization problem [START_REF] Bouyer | Optimal infinite scheduling for multi-priced timed automata[END_REF] and the discounted-cost problem. [START_REF] Fahrenberg | Discountoptimal infinite runs in priced timed automata[END_REF] For optimal reachability, another technique (priced regions) has been used [START_REF] Behrmann | Minimum-cost reachability for priced timed automata[END_REF] which also extends to a setting of more than one cost variable. [START_REF] Larsen | Optimal reachability for multi-priced timed automata[END_REF] As for algorithm and tool support, the zonebased approach has been successfully extended to solve the optimal reachability problem, [START_REF] Larsen | As cheap as possible: Efficient cost-optimal reachability for priced timed automata[END_REF] by introducing priced zones, and tool support is available in UPPAAL CORA. For mean-cost and discounted-cost optimization, active research is being conducted in developing efficient zone-based algorithms, or alternatively showing that no such algorithms exist.

Task graph scheduling: energy optimality

Reconsidering our running task graph scheduling problem of Section 1.5, costoptimal reachability for priced timed automata may be used to provide energyoptimal schedules. The energy needed for performing computation steps is modeled as cost in the priced timed automaton model, and optimal reachability techniques can be used for finding an energy-optimal schedule.

For the task graph scheduling instance of Figure 3, energy consumption of the two processors is reflected in the respective timed automata by suitable cost-rates in the locations corresponding to the processor being idle or in use. The processors can then be represented by the following two priced timed automata: Figure 5: The resource management problem asks whether it is possible to maintain the cost level within fixed bounds. There can be a lower bound only (a), a lower and an upper bound (b, c), or a lower bound and a soft upper bound above which cost level cannot increase. Figures (a), (b) and (d) represent solutions to the respective problems for the priced timed automaton depicted on the left: there is an infinite run that satisfies the global constraint. In case (a) for instance, we have depicted a possible schedule for the first cycle, and this run can be repeated because at the start of the second cycle, the cost level is larger than at the start of the first cycle. In figure (c), the proposed schedule violates the lower bound, and it can be shown that there exists no infinite run which maintains cost level within the specified bounds.

P1:

idle +10 + (x≤2) +90 × (x≤3) +90 x:=0 add 1 x:=0 mult 1 x=2 done 1 x=3 done 1 P2: idle +20 + (y≤5) +30 × (y≤7) +30 x:=0 add 2 x:=0 mult 2 y=5 done 2 y=7 done 2

Managing the resources

Up to this point we have only employed priced timed automata as a formalism for modeling time-dependent consumption of resources. However, in several situations resources may not just be consumed but also occasionally regained, e.g. in autonomous robots with rechargeable batteries, or in tanks which may not only be emptied but also filled. Extending priced timed automata to allow for both positive (regaining) and negative (consumption) rates provides a natural modeling formalism. [START_REF] Chakrabarti | Resource interfaces[END_REF] However, a new question now emerges related to the appropriate management of resources: "Is it possible to maintain the level of resources within fixed bounds?" Such resource-bound problems are highly relevant to the analysis of several embedded systems, e.g. it is natural to plan the usage of a device with rechargeable batteries so that one never runs out of energy, nor exceeds the maximum capacity for energy storage. Figure 5 shows a priced timed automaton together with some resource management problems.

Few results have been obtained on this problem so far: only the case of oneclock priced timed automata has been investigated. [START_REF] Bouyer | Infinite runs in weighted timed automata with energy constraints[END_REF] This restriction has two important consequences: Cycle detection can be done statically, as each resetting transition leads to a configuration with clock value 0, and the region automaton can be coarsened so that the partition consists of intervals with endpoints given by the constants in the automaton's guards. As a consequence, there are only polynomially many regions.

Under the additional assumption that the cost cannot be updated during transitions (hence cost evolves only in locations), it can be shown [START_REF] Bouyer | Infinite runs in weighted timed automata with energy constraints[END_REF] that for finding runs which satisfy a global lower-bound constraint, with or without soft upper bound, one can restrict oneself to look for runs with integral delays. Hence the corner-point abstraction can be used for this, and the problems are solvable in polynomial time.

For priced timed automata with more than one clock, no results are known, but even for one-clock automata with cost updates during transitions, there are some difficulties which mean that abstractions like the above are insufficient. As an example, consider the following priced timed automaton:

ℓ 0 +2 ℓ 1 +2 ℓ 2 +4 x=1,x:=0 -3 x=1,x:=0
Assuming that we start with initial cost 0, this automaton has exactly one feasible execution in which the cost level remains nonnegative: after spending 1 time unit in location ℓ0, we alternately spend half a time unit in ℓ1 and half a time unit in ℓ2. Any other execution eventually violates the lower bound. Hence in this case, runs satisfying the lower bound cannot be found using the corner-point abstraction.

PRICED TIMED GAMES

A model for uncertainties

The systems we have considered so far are closed in the sense that we have a complete description of the system. This is not sufficient to model embedded systems where interaction with the environment is crucial, or systems with some imprecision. These can be modeled using (two-player) timed games, [START_REF] Asarin | Controller synthesis for timed automata[END_REF] in which some actions are triggered by the environment (we can think of signals received by sensors, or of unexpected events). The aim is to control, or guide, the system so that it will be safe or correct regardless of the way the environment interferes. An example of a timed game is depicted below:

ℓ 0 ℓ 1 (y=0) ℓ 2 ℓ 3 x≤2,y:=0 x = 2 x = 2
Dashed edges belong to the environment (they are uncontrollable): when they are fireable, the system cannot prevent (nor force) them to be fired. Here, the system cannot decide whether it goes through ℓ2 or through ℓ3.

For simple correctness criteria, e.g. reachability or safety, the set of winning states (i.e., states from which the system can be controlled under the safety constraint) and also winning strategies (i.e., policies for how to control the system) can be computed using the region abstraction. [START_REF] Asarin | Controller synthesis for timed automata[END_REF] Also computability of time-optimal strategies, 7 as well as strategies under partial observability, has been demonstrated. For the latter, decisions are based on discrete observations giving only partial information of the system state, depending on the availability and precision of sensors. [START_REF] Cassez | Timed control with observation based and stuttering invariant strategies[END_REF] For efficient algorithms, a zone-based approach for solving timed games with reachability and safety objectives has been developed, [START_REF] Cassez | Efficient on-the-fly algorithms for the analysis of timed games[END_REF][START_REF] Tripakis | Controller synthesis for discrete and dense-time systems[END_REF] and tool support is available in UPPAAL-TIGA.

Task graph scheduling: timing uncertainty

Returning to our running task graph scheduling example, we can use the formalism of timed games to model uncertainty in precisely how much time a certain computation on a given processor takes. In Section 1.5 we modeled computation times by precise numbers, whereas we now can make the model more realistic by only providing interval bounds within which computation times are prescribed to lie. The timed game models below provide version of the processors P1 and P2 from Figure 3 in which computation times are prescribed to lie in the intervals [1, 2] for addition and [1, 3] for multi-plication on P1, and similarly for P2.

P1:

idle

+ (x≤2) × (x≤3) x:=0 add 1 x:=0 mult 1 x≥1 done 1 x≥1 done 1 P2: idle + (x≤2) × (x≤3) x:=0 add 2 x:=0 mult 2 y≥3 done 2 y≥2 done 2
Using these models, a computed timeoptimal schedule will no longer be a simple fixed assignment of tasks and time slots to processors, but rather a flexible dynamic assignment, where task scheduling can be adapted on-line according to actual completion times of previous tasks. (Hence we cannot display the solution here.)

Cost-optimal strategies

It is natural to extend the timed game framework with cost information, hence making it possible to model uncertainty as well as resource use, and to ask for controllability under resource constraints, or for optimal controllability. The model of priced timed games is a synthesis of priced timed automata and timed games; we show an example below:

ℓ 0 +5 ℓ 1 (y=0) ℓ 2 +10 ℓ 3 +1 x≤2,y:=0 x = 2 + 1 x = 2 + 7
In this example we may e.g. want to compute the minimum cost for reaching location regardless of the moves of the environment (which is in charge of the edges out of ℓ1 as before). As the system cannot control whether execution goes through ℓ2 or ℓ3, the minimum cost is given by the term

inf 0≤t≤2 max 5t+10(2-t)+1 5t+(2-t)+7 = 14 + 1 3
where t is the delay spent in location ℓ0. Solving this, one arrives at a minimum cost of 14 1 3 which is attained for t = 4 3 . As this is not an integer, one sees that techniques based on the corner-point abstraction are not sufficient for computing optimal-reachability strategies, even in case of one-clock priced timed games.

Generally, priced timed games are much more difficult to analyze than priced timed automata. Using reductions from the Halting problem for two-counter machines, one can show that cost-optimal strategies are undecidable, [START_REF] Th | On optimal timed strategies[END_REF] even when restricted to priced timed games with only three clocks.

Decidability has been shown for classes of priced timed games with strong conditions on the cost evolution 3 and for one-clock priced timed games. [START_REF] Bouyer | Almost optimal strategies in one-clock priced timed automata[END_REF] The reason for the latter is the same as for one-clock priced timed automata above: resetting the clock leads to a configuration with a known clock valuation.

APPLICATIONS & TOOLS

Timed automata and their extensions have been applied to the modeling, analysis and optimization of numerous real-time applications. In this section we give a few examples, not aiming at being exhaustive but rather to illustrate the wide range of application domains.

A variety of mature tools are available which provide important computer-aided support for applications. Well-known tools include UPPAAL, KRONOS, and HYTECH, but there is a large number of other tools available. The electronic version of this article contains an extra section which aims to give an overview together with references to the individual tools.

The timed automata formalism is now routinely applied to the modeling and analysis of real-time control programs, including a wide class of Programmable Logic Controller (PLC) control programs [START_REF] Dierks | PLC-automata: a new class of implementable real-time automata[END_REF][START_REF] Mader | Timed automaton models for simple programmable logic controllers[END_REF] and timing analysis and code generation of vehicle control software, [START_REF] Tripakis | Timing analysis and code generation of vehicle control software using TAXYS[END_REF] and the timed automaton approach has also demonstrated its viability to the timing analysis of certain classes of asynchronous circuits. [START_REF] Bozga | Verification of asynchronous circuits using timed automata[END_REF] Similarly, numerous real-time communication protocols have been analyzed using timed automata technology, often with inconsistencies being revealed: e.g. using real-time model checking, the cause of a 10 year old bug in the IR-link protocol used by Bang&Olufsen was identified and corrected. [START_REF] Havelund | Formal modeling and analysis of an audio/video protocol: An industrial case study using UPPAAL[END_REF] Most recently, real-time model checking has been applied to the clock synchronization algorithm currently used in a wireless sensor network that has been developed by the Dutch company CHESS. [START_REF] Schuts | Modelling clock synchronization in the Chess gMAC WSN protocol[END_REF] Here it is shown that in certain cases a static, fully synchronized network may eventually become unsynchronized if the current algorithm is used, even in a setting with infinitesimal clock drifts.

During the last years, timed automata modeling of multitasking applications running under real-time operating systems has received substantial research effort. Here the goals are multiple: to obtain less pessimistic worst-case response time analysis compared with classical methods for single-processor systems; [START_REF] Waszniowski | Formal verification of multitasking applications based on timed automata model[END_REF] to relax the constraints of period task arrival times of classical scheduling theory to task arrival patterns that can be described using timed automata; [START_REF] Fersman | Schedulability analysis of fixed-priority systems using timed automata[END_REF] to allow for schedulability analysis of tasks in terms of concurrent objects executing on multiprocessor or distributed platforms (e.g. MPSoC). [START_REF] David | Model-based framework for schedulability analysis using UPPAAL 4[END_REF] Just as symbolic reachability checking of finite-state models has lead to very efficient planning and scheduling algorithms, reachability checking for (priced) timed automata has demonstrated competitive and complementary performance with respect to classical approaches such as MIPL on optimal scheduling problems involving real-time constraints, e.g. job-shop and task-graph scheduling [START_REF] Abdeddaïm | Scheduling with timed automata[END_REF][START_REF] Behrmann | Efficient guiding towards cost-optimality in UPPAAL[END_REF] and aircraft landing problems. [START_REF] Larsen | As cheap as possible: Efficient cost-optimal reachability for priced timed automata[END_REF] In fact a translation of the variant PDDL3 of PDDL (Planning Domain Definition Language) into priced timed automata has been made [START_REF] Dierks | Finding optimal plans for domains with continuous effects with UPPAAL CORA[END_REF] allowing optimal planning questions to be answered by cost-optimal reachability checking. Industrial applications include planning a wafer scanner from semiconductor industry [START_REF] Hendriks | Model checker aided design of a controller for a wafer scanner[END_REF] and computation of optimal paper paths for printers. [START_REF] Igna | Formal modeling and scheduling of datapaths of digital document printers[END_REF] Most recently, computation of winning strategies for timed games has been applied to controller synthesis for embedded systems, including synthesis of most general nonpreemptive online schedulers for real-time systems with sporadic tasks, 2 synthesis of climate control for pig stables provided by the company Skov A/S, [START_REF] Jessen | Guided controller synthesis for climate controller using UPPAAL-TIGA[END_REF] and automatic synthesis of robust and near-optimal controllers for industrial hydraulic pumps. [START_REF] Cassez | Automatic synthesis of robust and optimal controllers -an industrial case study[END_REF]

CONCLUSIONS

Timed automata and their priced and game extensions provide a uniform and expressive formalism for dynamic resource allocation problems with hard real-time constraints, i.e. timing constraints that must be satisfied under all circumstances. This is in contrast to soft real-time constraints, which only need to be met with a certain probability, .999 say, and which require stochastic modeling formalisms such as discrete-time or continuoustime Markov chains, queueing models, etc. While hard real-time focuses on worst-case analysis, soft real-time addresses more refined properties such as average-case performance.

However within the setting of hard realtime, timed automata and their extensions allow for analysis of a wide collection of performance and optimization problems, with results competitive with respect to more traditional approaches such as mixed-integer linear programming or others.

Particularly challenging problems remaining to be settled include decidability of synthesis for priced timed games under partial observability, as well as a range of resource management problems in the setting of priced timed automata and games with both consumption and regaining of resources.

Figure 1 :

 1 Figure1: Several refinements of a model (a) of the working mathematician according to Erdős: After insertion of a coin into the coffee dispenser, coffee can be collected, and the scientist can go back to work. In the timed-automaton model (b), precisely five time units pass between coin insertion and coffee collection, and the time which passes between coin insertion and going back to work is less than ten time units. In the priced timed automaton (c), cost rates (modeling e.g. energy consumption) are associated with the three states. In the timed game (d), uncertainty as to precisely when coffee is delivered is modeled as an uncontrolled edge.

  (both clocks have integer values), open line segments (one clock has integer value, or their fractional parts are equal), open triangles, or open unbounded rectangles.

P 1 (Figure 3 :

 13 Figure 3: Task graph problem with 6 tasks, where each task corresponds to the computation of a given sub-expression of the term(D × (C × (A + B)) + ((A + B) + (C × D)).Given the execution platform with two processors, P1 and P2, and corresponding computation times for addition and multiplication, as well as their energy consumption, Sch1 to Sch3 provide three feasible schedules, where Sch2 is in fact time-optimal, and Sch3 is energy-optimal.

Figure 4 :

 4 Figure4: The corner-point abstraction refines the region abstraction by also keeping track of the corner point close to which an execution runs. This is needed to measure costs: for instance, if we are in location ℓ1 and in the red region where 0 < y < x < 1, the price of delaying depends on the value of the clocks. From (a), where both x and y are arbitrarily close to 0, we can let almost one time-unit elapse and reach (b). The resulting cost is arbitrarily close to +3. On the other hand, from (c), where x and y are arbitrarily close to 1 and 0, respectively, letting time elapse takes us to the subsequent region, so that the cost is arbitrarily close to 0. (Notice that for readability, some resetting transitions have been omitted.)

Acknowledgments

The authors are partly supported by the European project QUASIMODO (FP7-ICT-STREP-214755). The French authors are supported by project DOTS (ANR-06-SETI-003). The Danish authors are supported by the Danish Center of Excellence MT-LAB.


Alessio Colombo

Daniele Fontanelli

Axel Legay

Luigi Palopoli

Sean Sedwards

Motion Planning in Crowds using Statistical Model Checking to Enhance the Social Force Model

Crowded environments pose a challenge to the comfort and safety of those with impaired ability. To address this challenge we have developed an efficient algorithm that may be embedded in a portable device. The algorithm anticipates undesirable circumstances in real time, by verifying simulation traces of local crowd dynamics against temporal logical formulae. The model incorporates the objectives of the user, pre-existing knowledge of the environment and real time sensor data. The algorithm is thus able to suggest a course of action to achieve the user's changing goals, while minimising the probability of problems for the user and others in the environment.

To demonstrate our algorithm we have implemented it in an autonomous computing device that we show is able to negotiate complex virtual environments. The performance of our implementation demonstrates that our technology can be successfully applied in a portable device or robot.

I. INTRODUCTION

With unimpaired ability, pedestrians are able to negotiate crowded areas with few problems. With reduced ability or under panic conditions [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], finding a good strategy to proceed can be challenging. As a result, people afflicted by a decline in physical or cognitive abilities can be discouraged from attending crowded places, with a consequent negative impact on their physical condition (reduced exercise), on the quality of their nutrition (reduced fresh food) and on their psychological wellbeing (reduced social contact). Motivated by these considerations, the DALi project [START_REF]DALi project web site[END_REF] aims to devise an intelligent 'walker' (an assistive wheeled device) that detects the presence of other pedestrians in the environment, anticipates their intent and plans an appropriate path that is suggested to the user via a combination of audio, visual and haptic interfaces. In this work we present an efficient algorithm that employs advanced modelling and verification techniques to address the path planning problem in a crowded and unfamiliar environment.

Succinctly, the problem is one of devising an online motion planning algorithm for an autonomous agent (the user) in a dynamic environment. The position of most fixed objects (e.g., buildings and rooms) are known a priori, but the algorithm must account for the possibility of changes, such as temporary obstructions. The environment contains moving objects (i.e., other pedestrians), whose positions and velocities cannot be known before they are encountered. The overall goal is to allow the user to visit pre-defined locations in the environment, while avoiding collisions, crowding and delays. The output of the algorithm is a suggested trajectory, so the algorithm must be reactive to the potentially uncooperative response of the user. Practically, the algorithm will be implemented in a low power embedded computing device and must be sufficiently efficient to make course corrections in a time of the order of seconds. This time scale is dictated by the typical velocities of pedestrians and by the fact that frequent readings help to reduce the random errors produced by sensors.

Our solution is a two-tiered algorithm, comprising a low level predictive mathematical model of pedestrian dynamics, managed by a statistical model checking (SMC) engine that checks temporal logical properties expressing the high level goals and constraints of the user. The algorithm uses dynamic input from sensors to reconstruct the user's position from fixed objects and to account for non-fixed objects, such as other pedestrians and temporary obstructions.

A. Related work

Our work is related to sampling methods (e.g., [START_REF] Barraquand | A monte-carlo algorithm for path planning with many degrees of freedom[END_REF]- [START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF]) and to recent methods using temporal logic (.e.g, [START_REF] Loizou | Automatic synthesis of multi-agent motion tasks based on ltl specifications[END_REF]- [START_REF] Karaman | Sampling-based motion planning with deterministic µ-calculus specifications[END_REF]). It is also related to methods that predict behaviour based on models parametrised with data from sensors (e.g., [START_REF] Kuwata | Motion planning in complex environments using closed-loop prediction[END_REF]).

In common with existing sampling methods, our algorithm uses randomisation to cover an intractably large configuration space. In contrast to many existing uses of sampling, however, we do not assume a fixed environment. In our application the environment contains both fixed and dynamic elements, such that a single optimal path cannot be defined a priori. Hence, the problem we solve by sampling is not one of creating an optimal global plan, but one of finding an optimal local plan given a changing environment.

Model checking is an automatic process to verify that a system satisfies a property specified in temporal logic. In the present context, temporal logic can express complex dynamical properties such as "the user will visit all the desired locations in a specified sequence, within the specified time" and "the user will never get too close to any other pedestrian". If the notion of an optimal path can be so defined, the principles of model checking can be used to directly synthesise a 'correct' motion planner or to prove that an existing motion planner is correct [START_REF] Loizou | Automatic synthesis of multi-agent motion tasks based on ltl specifications[END_REF]- [START_REF] Karaman | Sampling-based motion planning with deterministic µ-calculus specifications[END_REF]. The use of probabilistic model checking in combination with the theory of stochastic hybrid automata [START_REF] Hu | Towars a theory of stochastic hybrid systems[END_REF] is particularly appealing for control and robotic applications where a nonzero probability of failing the mission can be tolerated. For example, in [START_REF] Hu | Aircraft conflict prediction in the presence of a spatially correlated wind field[END_REF], [START_REF] Prandini | A probabilistic framework for aircraft conflict detection[END_REF] for air traffic control or in [START_REF] Asaula | Safety provisions for human/robot interactions using stochastic discrete abstractions[END_REF], [START_REF]A probabilistic methodology for predicting injuries to human operators in automated production lines[END_REF] for industrial robotics. Combining model checking with sampling, algorithms can be constructed which provably converge to optimal schedulers [START_REF] Karaman | Sampling-based motion planning with deterministic µ-calculus specifications[END_REF]. Standard model checking algorithms are computationally intensive, hence existing applications have used model checking offline. By using statistical model checking, we are able to perform online verification. We do not prove correctness, but find a local plan that maximises the probability of success.

We believe our work is novel in using formal verification in real time to solve an adaptive motion planning problem in a dynamic environment.

B. Organisation of the paper

The rest of the paper is organised as follows. Section II-A gives an overview of our approach. Section II-B introduces our mathematical model in detail and Section II-C introduces the basic notions of statistical model checking. Section III gives a detailed description of our algorithm and implementation, while Section IV describes the results of a number of experiments that demonstrate the utility of our approach. Section V discusses our choices and highlights areas of ongoing development.

II. PRELIMINARIES

A. Overview of the approach

Fig. 1 gives a high level overview of the algorithm. At each iterative step the algorithm acquires the state of the system, comprising the position of static objects and the position and velocity of the user and of other people in the environment.

Given the current state, the algorithm hypothesises alternative courses of action using the social force model. Each hypothesised trajectory is formally verified (modelchecked) against properties that express goals and constraints required for the user's trajectory (i.e., where the user wants to go, obeying the appropriate social rules). This leads to a statistical distribution of potentially successful trajectories. The algorithm uses this distribution to choose an immediate action that maximises the probability of achieving the user's objectives and minimises the probability of problems. In this probabilistic context, the measurement noise is considered as an additional source of stochasticity.

The social force model may be programmed with the user's objectives and is an efficient way to describe the continuous interactions that allow pedestrians to avoid collisions. The model also includes stochasticity to model the typical unpredictability of human behaviour. We use the stochasticity to generate a random sample of possible futures and choose the course of action that maximises the probability of success. Such trajectories respect the basic social and physical laws of pedestrian interactions and include the possibility of unpredicted behaviour. Their distribution allows the algorithm to choose a course of action that maximises the probability of success.

The stochasticity, while realistic, places an upper bound on the predictive accuracy of the model. Moreover, the model alone cannot account for the overall "mission" of the user. The predictive model needs to be managed reactively. Fortunately, the field of statistical model checking (SMC) encapsulates the technologies that we require to do this. SMC provides efficient algorithms to verify hypothesised trajectories against the user's constraints and objectives expressed in temporal logic. SMC can estimate the probability of success and bound the error of the estimation.

The key elements of our approach are (i) the social force model to hypothesise trajectories that respect low level social and physical "forces"; (ii) temporal logic to express the high level goals of the user and (iii) a statistical model checker to verify the traces with respect to the goals.

B. The social force model

The social force model [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF]- [START_REF] Helbing | Simulation of Pedestrian Crowds in Normal and Evacuation Situations[END_REF] combines real and psychological forces to predict the behaviour of pedestrians in crowds, under normal and panic situations. The model recognises that pedestrians are constrained by the physical laws of motion and also by social 'laws' that can be modelled by external forces. The model considers an environment comprising fixed objects (walls) and moving agents (pedestrians) that respond to attractive and repulsive forces, originated by social and physical interactions.

The model is constructed in two dimensions [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF]- [START_REF] Helbing | Simulation of Pedestrian Crowds in Normal and Evacuation Situations[END_REF], with agents represented as circular discs. In what follows we adopt the convention of denoting vectors in bold type. Thus, agent i has mass m i centred at position x i ∈ R 2 in the environment, radius r i and velocity v i ∈ R 2 . The linear model for the i-th agent is given by

   ẋi = v i vi = v 0 i -v i τ i + f i + ξ i m i (1)
v 0 i is the driving (desired) velocity of agent i, represented by a product of speed amplitude v 0 i and normalised direction e 0 i , which is given by the direction of the line joining the initial and desired configurations. τ i is the time taken to react to the difference between desired and actual velocity, while ξ i is a noise term (a random variable) that models random fluctuations not accounted for by the deterministic part of the model. The inclusion of the noise term makes the model stochastic, such that a different trajectory is generated each time (1) is solved. This allows the application of SMC and serves to avoid deadlocks that might arise if, by chance, some of the deterministic forces are equal and opposite.

f i is the force acting on agent i resulting from other objects in the environment and, hence, given by

f i = ∑ j̸ =i [f soc ij + f att ij + f ph ij ] + ∑ b [f soc ib + f ph ib ] + ∑ c f att ic . ( 2
)
The first term on the right-hand side of ( 2) includes all the forces on agent i resulting from interactions with other agents: f soc ij is the repulsive social force that inhibits agents getting too close, f att ij is the attractive social force that brings friends together, f ph ij is the physical force that exists when two agents touch. The second summation includes the forces acting on agent i as a result of the boundaries (walls): f soc ib is the social force that inhibits agent i from getting too close to boundaries, f ph ib is the physical force that exists when agent i touches boundary b. Finally, f att ic is the attractive social force that draws agent i towards fixed objects of incidental interest (shops, cafs, toilets, etc.).

In general, the force acting on any agent is calculated with respect to the distance between its centre of mass and all other visible objects. Since the model mixes both notional (social) and real forces, the mass m i is notionally the real mass of agent i. Other parameters can be used to model the unique characteristics of individual agents. For example, the latency factor τ i can be used to model the possibly reduced mobility of agent i. Full details of these and other parameters can be found in [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF]. In [START_REF] Colombo | Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments[END_REF] we show how the model may be parametrised from captured motion.

C. Statistical and probabilistic model checking

Model checking is an automatic technique to verify that a system satisfies a property [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]. Typically, the system has discrete states and the property is specified in temporal logic. The output of standard model checking algorithms is either true or false, with the possibility to give corresponding examples or counter-examples. The logics used are desired to be expressive, but must be decidable and tractable. Typical logics for standard model checking include LTL and CTL [START_REF] Baier | Principles of Model Checking[END_REF]. To give a result with certainty, the algorithms effectively perform an exhaustive exploration of the state space of the system. The number of states scales exponentially with the number of interacting components in the system, leading to a 'state explosion problem' [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF] that can make standard model checking slow or intractable.

Probabilistic model checking extends the standard notion to include probabilistic or stochastic transitions. These can express the uncertainties of modelling and reality. The output of probabilistic model checking algorithms is the probability that an arbitrary execution of the system will satisfy a given property. Such properties are specified in probabilistic or stochastic logics, such as PCTL and CSL [START_REF] Baier | Principles of Model Checking[END_REF]. The probabilistic model checking problem is solved with numerical model checking algorithms. These calculate the notionally exact probability by considering all the states. As such, they suffer the same state explosion problem as standard model checking algorithms.

Statistical model checking (SMC) is a type of probabilistic model checking that avoids the state explosion problem by estimating the probability of a property ϕ from executions (simulations) of the system. Given N independent simulation traces ω i , i ∈ {1 . . . N }, and a model checking function 1(ω i |= ϕ) ∈ {0, 1} that indicates whether ω i |= ϕ (read "ω i satisfies ϕ"), the probability γ that an arbitrary execution satisfies ϕ is estimated using γ ≈ 1/N ∑ N i=1 1(ω i |= ϕ). This reduces the probabilistic model checking problem to estimating the parameter of a Bernoulli random variable, hence the confidence of the estimate can be guaranteed by standard statistical bounds (e.g., the Chernoff bound [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF]). In general, the confidence of the estimate increases with increasing N . In comparison to standard and numerical model checking, SMC does not require decidable logics nor a finite state space, making it particularly suitable for the present application that considers continuous time and space.

Bounded Linear Temporal Logic: Our model checking engine is based on the PLASMA-lab library [START_REF]PLASMA-lab web site[END_REF]. PLASMAlab implements the function 1(ω i |= ϕ) using bounded linear temporal logic (BLTL [START_REF] Jha | A bayesian approach to model checking biological systems[END_REF]) to express the property ϕ:

ϕ = ϕ ∨ ϕ | ϕ ∧ ϕ | ¬ϕ | F ≤t ϕ | G ≤t ϕ | ϕ U ≤t ϕ | Xϕ | α
∨, ∧ and ¬ are the standard logical connectives and α is a Boolean constant or an atomic proposition constructed from numerical constants, state variables and relational operators. X is the next temporal operator: Xϕ means that ϕ will be true on the next step. F, G and U are temporal operators bounded by time interval [0, t], relative to the time interval of any enclosing formula. We refer to this as a relative interval. F is the finally or eventually operator: F ≤t ϕ means that ϕ will be true at least once in the relative interval [0, t]. G is the globally or always operator: G ≤t ϕ means that ϕ will be true at all times in the relative interval [0, t]. U is the until operator: ψU ≤t ϕ means that in the relative interval [0, t], either ϕ is initially true or ψ will be true until ϕ is true. Combining these temporal operators creates complex properties with interleaved notions of eventually (F), always (G) and one thing after another (U). A detailed description of the semantics of BLTL is given in [START_REF] Jha | A bayesian approach to model checking biological systems[END_REF].

III. SMC-BASED MOTION PLANNER

Our motion planner is based on the scheme depicted in Fig. 1 and Algorithm 1. The planner assumes the existence of a pre-calculated global plan (GlobalPlan) that visits the user's objectives in an a priori optimal way, that is, considering all things known in advance. Typically, the global plan Algorithm 1 The planning algorithm is computed with respect to a map of the static objects in the environment, the user's objectives and predicted anomalies (e.g., known crowded areas). Any contradiction of the a priori assumptions (e.g., an unforeseen blockage) triggers a recalculation of the global plan. We do not describe this recalculation in the present work.

The sensor board provides the current state of the local environment, located with respect to the global plan: the position and velocity of the user (state user ); the positions and velocities of other pedestrians (state ped1 , state ped2 , . . . ); the position of static objects (Map). The algorithm calculates a local way point w, which is the user's point of greatest straight line progress along the global plan within the sensor range. w is used to calculate the user's driving velocity v 0 , assuming a constant desired speed. The driving velocities of the other pedestrians are estimated from their current velocities.

The algorithm uses the above information to construct social force models (1) of the local environment. Specific characteristics (e.g., τ i ) of other pedestrians are unknown to the algorithm, so it assumes the default values given in [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF]. In the current implementation we construct the noise term ξ i from two normal distributions; one for the magnitude and one for the direction. The results presented here are based on heuristically estimated parametrisations of these distributions, which appear to be adequate. In [START_REF] Colombo | Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments[END_REF] we present a way of obtaining better parameters from captured motion.

The motion planner assumes the user will follow the global plan, but need to temporarily deviate to avoid collisions. The output of the algorithm is a suggested deviation, α best , in the range ±90 degrees relative to the user's direct path to w. To find α best , the algorithm constructs models for each hypothesised deviation in the set {0, ±25, ±50, ±75, ±90}. These values are chosen to span ±75 degrees using a tractable number of different values, with ±90 included in case the user needs to sidestep an obstacle (see [START_REF] Colombo | Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments[END_REF]). Each model is then investigates using statistical model checking.

The algorithm sets α curr ∈ {0, ±25, ±50, ±75, ±90} and calls function SMC with arguments N and Formula. SMC estimates the probability of success P curr for a particlar deviation α curr by the proportion of N simulation traces that satisfy the BLTL property Formula. The value of α curr is used as the intial deviation: the user's driving velocity is initially rotated by α curr , but at each successive step of the simulation the deviation from a direct path to w is reduced to zero. This ensures that the user will eventually be close to the global plan.

BLTL is expressive enough to define complex sequences of high and low level requirements. For the results presented here, Formula merely expresses the basic constraints of the user:

(G [0,T horizon ] ∧ i̸ =u ∥x u -x i ∥ > 0.5)∧ (F [0,T horizon ] ∥x u -w∥ < 0.2) (3)
x u denotes the position of the user and ∥•∥ denotes Euclidean distance. Intuitively, (3) means that "in the next T horizon time units the user will get no closer than 0.5m to any other pedestrian and will eventually be less than 0.2m from the global plan". T horizon is chosen to be the expected time for the user to walk a distance equivalent to the range of the sensors. Using a higher value might produce impossible trajectories that pass through unseen fixed objects; using a lower value might exclude possible collisions. In our implementation we use T horizon = 4s.

For each hypothesised deviation α curr , the SMC function returns the probabillity of success P curr and the expected distance from the global plan, d curr . These are used by function is better to decide α best . is better chooses the smallest | α curr | which maximises P curr . Ties are resolved by choosing the alpha curr with smallest d curr or randomly if d curr also ties. If P best == 0 the user is required to stop (the global plan will be recalculated).

T decision is the actaul time the algorithm takes to make its predictions and must be less than the time period it is predicting, i.e., T horizon . In practice T decision is bounded below by the performance of the hardware, the complexity of the environment (fixed and moving objects) and the confidence required (controlled by the number of simulations, N ). In our implementation, T decision ≈ 1s.

At each decision point α best is suggested to the user. The user may ignore this suggestion and move in a different direction, but the operation of the algorithm in the next decision period remains the same: α best is calculated according to the global plan and the actual positions and velocities of the user and other pedestrians. Since the user specifies the global plan, when generating hypothesised traces the algorithm assumes that the user is compliant, however ξ user may be used to model a lack of compliance. The present work does not consider how the user's non-compliance might affect predictions.

Given an accurate stochastic model of the behaviour of pedestrians, the Chernoff bound [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF] predicts that with N = 10 simulation runs the estimate of the probability of success has a maximum error of ±0.3 with probability 0.7. With N = 50 the probability of success has a maximum error of ±0.2 with probability 0.90. In general, the statistical confidence of the estimate increases with increasing N , but this only increases the probability of choosing the correct α best . The predictive power of the model is bounded by its stochasticty. Thus, given finite computational power, we choose a value of N that balances the reactive and predictive aspects of the algorithm. That is, we choose a value of N that allows us to make T decision sufficiently small.

The algorithm solves (1) using a standard ODE solver [START_REF] Ahnert | Odeint -Solving Ordinary Differential Equations in C++[END_REF], which produces traces comprising a sequence of states at discrete time points. Since the model given in Section II-B is based on continuous time and space, to guarantee properties that rely on the distance between objects it is necessary to choose time points that are sufficiently close. This is achieved by the ODE solver using adaptive time steps.

Simulating the traces accounts for most of the computational cost of the algorithm. We have found our chosen ODE solver to be efficient and presume its performance scales in a standard way with respect to the number of visible moving agents M and the complexity of their interactions. Since the forces in the model are dependent on the distances between agents, there is an additional O(M 2 ) cost, however M is bounded by the range of the sensors.

IV. SIMULATIONS

To demonstrate our algorithm we have implemented a prototype on a off-the-shelf, low power embedded system, the Beagleboard xM 1 . It is a portable device that may run from battery power and provides performance comparable to a small computer. We use PLASMA-lab [START_REF]PLASMA-lab web site[END_REF] as the statistical model checking library. To test the algorithm we have created a virtual environment that evolves according to the Social Force Model and contains fixed objects and other pedestrians that react to the user's presence.

The pedestrians are assigned individual global plans to simulate their objectives and individual parameters that reflect the variation seen in reality. The values of the parameters are based on the ones estimated in [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF] and two different but correlated sets, one for the planner and one for the virtual environment, have been defined in order to increase the sense of reality. The noise term ξ has been differentiated as well, the standard deviation of the two normal distributions in the planner has been set as the double of the one in the virtual environment.

In this way we simulate pedestrians that are reactive to the user and each other, with behaviour that is realistically unpredictable. Moreover, the simulated device has limited omnidirectional sensing range, we suppose it is able to detect agents moving within a radius of 4 meters with respect to the current position of the user. In the final application, a sensor board connected to the single board computer will provide the real (estimated) positions and velocities of the user and nearby pedestrians.

We compared three different strategies:

• SMC with the Social Force Model (SM C +SF M ): our novel approach, where Algorithm 1 computes the local 3). We also used the following parameters: T horizon = {1, 2, 4, 6, 8}, T decision = 1 and N = 50. We performed 500 independent runs for SF M and 500 for every combination of T horizon for SM C+SF M and SM C + LIN . Our objective was to demonstrate that 1) the higher complexity of our approach leads to valuable payoff in terms of performance and 2) it can be implemented online on an embedded device with limited computing power.

a) Algorithm performance analysis: We have devised two scenarios that challenge our algorithm and highlight significant features of its performance. In the first one (namely, scenario 1, depicted in Fig. 2(a)) the user moves on a straight line close to a fixed obstacle, while two agents are moving towards him following a straight line as well. In the second one (namely, scenario 2, showed in Fig. 2(b)) the user attempts to visit a market stall at the end of the market while some pedestrians (Agent 1-6) block the user's progress by entering the scenario and moving from one market stall to another. The user's global plan is a straight line from the left to the right of the market. Figure 3 depicts the distances over time with respect to the user, respectively, for one particular run of scenario 2.

We defined 4 indicators to measure performance: 1) the time needed for the user to reach the right side of the scenario (T exit ), 2) the measured probability of respecting the minimum safety distance to agents (P saf e ), 3) the average deviation in position from the global plan (ϵ x ) and 4) the average deviation of the orientation of the user with respect to the ideal orientation of a user perfectly following the global plan (ϵ θ ). These indicators are formally defined as follows.

Let x(t) represent the cartesian coordinates of the position of the user after the planning for each time t, θ(t) represent its orientation with respect to a fixed frame, x(t) the long term plan and θ(t) the orientation decided according to the long term plan. The integral error of the difference between the corrected plan and the global plan can be defined as

ϵ x = E { √ 1 T ∫ T 0 |x(t) -x(t)| 2 dt } . A similar performance
indicator ϵ θ is defined for the orientation θ(t).

Indicators ϵ θ and P saf e can be used to quantify the "comfort" of the user. Indeed, frequent changes in the direction reduce the user experience, especially if elderly, and so does the probability of accidents. Table I and Table II reports the performance we obtained for scenario 1 and scenario 2, respectively, using different values for T horizon .

Scenario 1 is the most problematic for SF M due to the limitations of this model we discussed in [START_REF] Colombo | Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments[END_REF]. The SM C-based strategies exhibit a higher P saf e and a lower ϵ θ when T horizon ≥ 6. SM C + SF M , in turn, outperform SM C + LIN on all indicators.

In scenario 2, from the safety and comfort point of view of the user, SM C + SF M approach obtains a higher P saf e and a lower ϵ θ with respect to SF M and SM C + LIN , when T horizon ≤ 6. Nonetheless, P saf e decreases and ϵ θ increases when T horizon > 6. This is motivated by the fact that the tested temporal logic formula is less likely to be satisfied over a large horizon in a crowded environment. As a consequence, the planning algorithm suggests the user to stop and/or to change direction in order to avoid the unfeasible path, thus raising ϵ θ .

The SF M strategy exhibits the lower ϵ x because it tends to keep the user closer to the global plan. However, this reflects negatively on the "comfort" of the user, especially on P saf e , because the model doesn't have an explicit notion of minimum safety distance to agents. This behaviour is more evident in scenario 1.

b) Timings on the Beagleboard xM: We measured the time needed by the Beagleboard xM to execute the scenarios presented in the previous section. We ran 500 simulations each and we timed the execution of every single decision step for the SM C + SF M strategy, that is, the time needed for a single run of Algorithm 1 using the Social Force Model as the model for the agents. We also set N = 50. We computed both the average µ 1 = 228.9 ms and µ 2 = 2026.1 ms and standard deviations σ 1 = 392.1 ms and σ 2 = 2432.1 ms of the timings for scenario 1 and scenario 2, respectively. If we allow a maximum 1000 ms latency to compute the decision step, the current implementation is able to satisfy it in 93.4% of the cases for scenario 1 and in 40.9% of the cases for scenario 2. Fig. 3. Scenario 2. Distances of the agents with respect to the user during one particular run of scenario 2. In this case the safety distance has been set to 0.5 m (dashed line) and has been violated once.

V. CONCLUSION AND ONGOING WORK

The social force model is a generative model of pedestrian behaviour, which we have embedded in an efficient online motion planning algorithm. We parametrise the model with data from sensors in real time and thus hypothesise future trajectories. The model's stochastic element limits its predictive ability but allows us to consider a distribution of possible futures. The reactive part of the algorithm is provided by statistical model checking technology. Thus, the algorithm verifies the hypothesised trajectories against the user's goals and constraints expressed in temporal logic. In this way the algorithm finds the immediate course of action that maximises the user's probability of success.

The apparently random behaviour of pedestrians is often the result of deterministic choices on their part. We hope to to improve the performance of our algorithm by recognising these choices and thus replacing some of the stochasticity. To this end, in [START_REF] Colombo | Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments[END_REF] we have identified behavioural templates that may be incorporated into the social force model. Also, equation (2) includes the possibility to explicitly model incidental attractive and repulsive forces that might, for example, arise from interesting shops and areas with high probability of crowding, respectively. Such forces apply to pedestrians in general and could enrich the existing model.

As part of our larger project [START_REF]DALi project web site[END_REF], we also propose to include advanced sensor techniques to recognise known interesting or hostile people (e.g., using facial recognition) and to generally avoid people exhibiting hostile behaviour. Such forces apply asymmetrically and would obviously have to be included in an anisotropic version of the social force model [START_REF] Helbing | Simulation of Pedestrian Crowds in Normal and Evacuation Situations[END_REF].

A significant part of the challenge of our motion planning application is the performance of its implementation. Current hardware performance forces us to accept the necessity of multiple boards to handle the overall computational burden, but there is a clear advantage if a portable device can be made to work on a single board. The embedded computing boards we have chosen for our implementation include high performance graphical processor units (GPUs) that can be used for general purpose parallel computing. Since statistical model checking requires multiple independent simulation runs, we propose to exploit the GPU to gain a significant increase in performance.. T decision is necessarily less than T horizon , hence the algorithm predicts traces in time periods that overlap from one iteration to the next. While the predictions of older simulations are likely to be less accurate with respect to the current reality, data from the previous iterations may be employed, suitably weighted, to build a probabilistic map of the good and bad locations in the local environment. This map can be used to avoid simulations that explore directions that are unlikely to be successful and to provide haptic feedback if the user chooses to diverge from the proposed path.

Fig. 1 .

 1 Fig. 1. Diagrammatic overview of the motion planning framework. The sensor board detects the current state of objects in the environment. This state is used by the social force model to generate plausible future paths of the user and other pedestrians. The distribution of paths is verified against the global objectives of the user in order to suggest an optimal course.

Fig. 2 .

 2 Fig. 2. Scenarios used to test the algorithm, scenario 1 (a) and, scenario 2 (b).

1 :

 1 function FINDLOCALPATH(stateuser, state ped 1 , state ped 2 , . . .

		, Map,
		GlobalPlan, Formula, N )
	2:	Real Pcurr, dcurr, P best , d best ;
	3:	[P best , d best ]=[0,∞];
	4:	for αcurr ∈ {0, ±25, ±50, ±75, ±90} do
	5:	[Pcurr, dcurr] = SMC(N, Formula);
	6:	if is better([Pcurr, dcurr], [P best , d best ]) then
	7:	α best = αcurr;
	8:	[P best , d best ] = [Pcurr, dcurr];
	9:	end if
	10:	end for
	11:	if P best == 0 then
	12:	return STOP;
	13:	else
	14:	return α best ;
	15:	end if
	16: end function

TABLE I SCENARIO 1 :

 I1 PERFORMANCE FOR SM C + SF M , SM C + LIN AND SF M STRATEGIES. 500 SIMULATIONS EACH WERE CONDUCTED. PERFORMANCE FOR SM C + SF M , SM C + LIN AND SF M STRATEGIES. 500 SIMULATIONS EACH WERE CONDUCTED.

		Unit		SM C + SF M			SM C + LIN		SF M
	T horizon	[s]	1	2	4	6	8	1	2	4	6	8	-
	T exit	[s]	23.08	23.38	22.72	21.68	21.11	23.17	24.63	24.55	24.42	24.19	23.56
	P saf e	-	0.7444 0.8923	0.9933	0.9981	0.9985	0.7511	0.8709	0.9565	0.9989	0.9925	0.7386
	ϵx	[m]	0.3504 0.9914	1.4377	1.6131	1.7386	0.3322	0.9832	1.4761	1.9007	2.0384	0.3062
	ϵ θ	[DEG]	53.19	37.22	13.93	10.84	9.36	55.89	48.03	40.11	24.66	22.47	36.86
							TABLE II
	SCENARIO 2: Unit		SM C + SF M			SM C + LIN		SF M
	T horizon	[s]	1	2	4	6	8	1	2	4	6	8	-
	T exit	[s]	26.82	23.89	24.08	21.12	20.27	27.33	31.48	36.03	29.98	23.71	23.16
	P saf e	-	0.9908 0.9998	0.9993	0.9977	0.9316	0.9882	0.9965	0.9977	0.9925	0.9486	0.9665
	ϵx	[m]	0.7677 0.7927	0.6497	0.5701	0.5430	0.7902	1.4279	1.2607	0.8282	0.6760	0.3825
	ϵ θ	[DEG]	9.97	5.50	9.20	10.59	19.97	10.62	14.25	20.33	21.56	21.52	13.67

http://www.beagleboard.org

This work is partially funded by the Devices for Assisted Living (DALi) European project, reference 288917.