Kim G Larsen

Uli Fahrenberg

Claus Thrane

Metrics for Weighted Transition Systems: Axiomatization and Complexity

Simulation distances are essentially an approximation of simulation which provide a measure of the extent by which behaviors in systems are inequivalent. In this paper, we consider the general quantitative model of weighted transition systems, where transitions are labeled with elements of a finite metric space. We study the so-called point-wise and accumulating simulation distances which provide extensions to the well-know Boolean notion of simulation on labeled transition systems.

We introduce weighted process algebras for finite and regular behavior and offer sound and (approximate) complete inference systems for the proposed simulation distances. We also settle the algorithmic complexity of computing the simulation distances.

Introduction

The need for an extension of the state-of-the art modeling and verification techniques to encompass systems with quantitative information has long been recognized; see [START_REF] Henzinger | The embedded systems design challenge[END_REF] for a recent position paper on this subject. Classical modeling formalisms for concurrent and reactive systems have focused on describing qualitative aspects of systems with a range of behavioral equivalences and preorders used for the so-called implementation verification, see e.g. the survey provided in [START_REF] Van Glabbeek | The Linear Time-Branching Time Spectrum I. The Semantics of Concrete, Sequential Processes[END_REF]. This approach requires a model of the systems and specifications, as well as a procedure for checking whether the two are related with respect to the given equivalence or preorder.

During more than a decade, classical modeling formalisms have been extended with quantitative aspects such as real time, probabilistic or continuous (so-called hybrid) information. Despite successful generalization of several behavioral preorders and equivalences they largely remain qualitative, e.g. two (quantitative) system models either are, or are not, equivalent. To properly take account of robustness, it is advocated in [START_REF] Henzinger | The embedded systems design challenge[END_REF] that in the quantitative setting, equivalences and preorders are replaced by real-valued distances: i.e. from deciding on the Boolean truth of equivalence P ∼ Q between two models P and Q, the problem becomes that of computing their distance |P, Q| = ǫ ∈ | ≥0 . It is argued that exact behavioral equivalence for quantitative models is unrealistic -as it typically requires exact matching of all quantitative aspects -whereas in practical application matching up to some error margin given by the distance ǫ suffices.

During the last years, substantial progress has been made towards defining suitable metrics or distances for various types of quantitative models including probabilistic models [START_REF] Desharnais | Metrics for labelled Markov processes[END_REF], real-time systems [START_REF] Henzinger | Quantifying similarities between timed systems[END_REF] and metrics for linear and branching systems in general [START_REF] Černý | Simulation distances[END_REF][START_REF] Davoren | Topologies, continuity and bisimulations[END_REF][START_REF] Luca De Alfaro | Linear and branching system metrics[END_REF][START_REF] Fahrenberg | A quantitative characterization of weighted Kripke structures in temporal logic[END_REF][START_REF] Fahrenberg | Distances for weighted transition systems: Games and properties[END_REF][START_REF] Henrik Johansson | On the regularization of zeno hybrid automata[END_REF][START_REF] Kennaway | Transfinite reductions in orthogonal term rewriting systems[END_REF][START_REF] Thrane | Quantitative simulations of weighted transition systems[END_REF][START_REF] Franck Van Breugel | A theory of metric labelled transition systems[END_REF].

In this paper, which is the third in a series of papers on general system distances [START_REF] Fahrenberg | A quantitative characterization of weighted Kripke structures in temporal logic[END_REF][START_REF] Thrane | Quantitative simulations of weighted transition systems[END_REF], we consider the general quantitative model of weighted transition systems where transitions are labeled with elements from a finite metric space ❑. We consider two different distances on states of such transition systems, point-wise and accumulating simulation distance, and provide sound and complete axiomatizations for these distances on weighted process algebras, akin to the aximoatization of bisimulation for finite and regular process algebra in Milner's seminal paper [START_REF] Milner | A complete inference system for a class of regular behaviours[END_REF]. Note that the maximum-lead distance from [START_REF] Henzinger | Quantifying similarities between timed systems[END_REF][START_REF] Thrane | Quantitative simulations of weighted transition systems[END_REF] is not treated here; we leave this for future work.

We also consider the algorithmic complexity of computing point-wise and accumulating simulation distance for finite-state weighted transition systems. Whereas point-wise simulation distance is shown decidable in polynomial time -similar to that of ordinary simulation and bisimulation as shown by Smolka and Kanellakis [START_REF] Kanellakis | CCS expressions, finite state processes, and three problems of equivalence[END_REF] -we show that the problem of accumulating simulation distance is polynomial-time equivalent to that of computing payoff for discounted games and hence in NP ∩ coNP.

Weighted Transition Systems

We need to fix some terminology and notation which we will use heavily:

A mapping d : X×X → [0, ∞] = | ≥0 ∪{∞}
from a set X to the non-negative reals together with positive infinity is called a hemimetric if d(x, x) = 0 for all x ∈ X and d(x, y) + d(y, z) ≥ d(x, z) (the triangle inequality) for all x, y, z ∈ X; it is called a metric if additionally, d(x, y) > 0 for all x = y ∈ X and d(x, y) = d(y, x) for all x, y ∈ X.

A sequence (x j ) in a metric space X is a Cauchy sequence if it holds that for all ǫ > 0 there exist

N ∈ ◆ such that d(x m , x n ) < ǫ for all n, m ≥ N . X
is said to be complete if every Cauchy sequence in X converges in X.

A continuous function f : X → X is called a contraction if there exists 0 ≤ α < 1 (its Lipschitz constant) such that d(f (x), f (y)) ≤ αd(x, y) for all x, y ∈ X. Finally, we recall the Banach fixed-point theorem: Any contraction on a complete metric space has a unique fixed point.

Throughout this article we fix a finite metric space ❑ of weights with a metric d ❑ : ❑ × ❑ → |. We also fix a discounting factor λ with 0 ≤ λ < 1, which will be used in the definition of accumulating distance below. Definition 1. A weighted transition system is a tuple (S, T ), where S is a finite set of states and T ⊆ S × ❑ × S is a set of (weighted) transitions.

Note that all transition systems in this paper are indeed assumed finite, hence requiring finiteness of the metric space ❑ does not add extra restrictions.

Simulation distances

In this section we fix a weighted transition system (S, T ) and introduce simulation distance between states in (S, T ). We concentrate on two types here, accumulating and point-wise distance, but other kinds may indeed be defined.

Accumulating distance

Definition 2. For states s, t ∈ S, accumulating simulation distance from s to t is defined to be the least fixed point to the set of equations

s, t = max s n -→s ′ min t m -→t ′ d ❑ (n, m) + λ s ′ , t ′ (1)
To justify this definition, we need to show that the equations (1) indeed have a least solution. To this end, write S = {s 1 , . . . , s p } and assume for the moment that the transition system (S, T ) is non-blocking such that every s i ∈ S has an outgoing transition s i n -→ s k for some s k ∈ S. Define a function

F : | p×p ≥0 → | p×p ≥0 by F (x) i,j = max s i n -→s k min s j m -→s ℓ d ❑ (n, m) + λx k,ℓ
Here we are using the standard linear-algebra notation | p×p ≥0 for p × pmatrices with entries in | ≥0 and x k,ℓ for the entry in their k'th row and ℓ'th column.

Lemma 3. With metric on | p×p ≥0 defined by d(x, y) = max p i,j=1 |x i,jy i,j |, F is a contraction with Lipschitz constant λ.

Proof. (Cf. also the proof of Theorem 5.1 in [START_REF] Zwick | The complexity of mean payoff games[END_REF].) We can partition | p×p ≥0 into finitely many (indeed at most 2 p 2 q 2 with q = |❑|) closed polyhedral regions R i,j (some of which may be unbounded) such that for x, y ∈ R i,j in a common region, the p 2 max-min equations get resolved to the same transitions. In more precise terms, there are mappings n, m, k, ℓ : {1, . . . , p} × {1, . . . , p} → {1, . . . , p} such that F (x) i,j = d ❑ (n(i, j), m(i, j)) + λx k(i,j),ℓ(i,j) for all x ∈ R i,j . Now if x, y ∈ R i,j are in a common region, then

d(F (x), F (y)) ≤ λ max i,j |x k(i,j),ℓ(i,j) -y k(i,j),ℓ(i,j) | ≤ λ max i,j |x i,j -y i,j | = λd(x, y) If x ∈ R i 1 ,j 1 , y ∈ R i 2 ,j 2
are in different regions, a bit more work is needed. The straight line segment between x and y admits finitely many intersection points with the regions R i,j ; denote these x = z 0 , . . . , z q = y. We have

d(F (x), F (y)) ≤ d(F (z 0 ), F (z 1 )) + • • • + d(F (z q-1 , z q )) ≤ λ d(z 0 , z 1 ) + • • • + d(z q-1 , z q ) = λd(x, y)
Note that the last equality only holds because all z i are on a straight line.

Using the Banach fixed-point theorem and completeness of | p×p ≥0 we can hence conclude that F has a unique fixed point. In the general case, where (S, T ) may not be non-blocking,

F is a function [0, ∞] → [0, ∞] with (extra) fixed point [∞, . . . , ∞]. Hence as a function [0, ∞] → [0, ∞],
F has at most two fixed points. Now we can write the equation set from the definition as

     s 1 , s 1 s 1 , s 2 • • • s 1 , s p s 2 , s 1 s 2 , s 2 • • • s 2 , s p . . . . . . . . . . . . s p , s 1 s p , s 2 • • • s p , s p      = F      s 1 , s 1 s 1 , s 2 • • • s 1 , s p s 2 , s 1 s 2 , s 2 • • • s 2 , s p . . . . . . . . . . . . s p , s 1 s p , s 2 • • • s p , s p     
hence (1) has indeed a unique least fixed point.

Point-wise distance

For point-wise simulation distance we follow a lattice-theoretic rather than a contraction approach. Definition 4. For states s, t ∈ S, point-wise simulation distance from s to t is defined to be the least fixed point to the set of equations With partial order on [0, ∞] p×p defined by x ≤ y iff x i,j ≤ y i,j for all i, j, G is (weakly) increasing.

s, t • = max s n -→s ′ min t m -→t ′ max d ❑ (n, m), s ′ , t ′ • Let G : [0, ∞] p×p → [0, ∞] p×p be the function defined by G(x) i,j = max s i n -→s k min s j m -→s ℓ max d ❑ (n, m), x k,ℓ s 1 s 2 s 5 s 4 3 5 11 t 1 t 2 t 3 t 4 t 5 3
Proof. Trivial. Now the Tarski fixed-point theorem allows us to conclude that G has a unique least fixed point, hence the above definition is justified.

Example

We show a computation of the two simulation distances between states s 1 and t 1 in the weighted transition system in Figure 1. Here ❑ ⊆ ◆, d(n, m) = |n -m|, edges without specified weight have weight 0, and the discount factor is λ = .90.

Repeated application of the definition yields the following fixed-point equation for s 1 , t 1 (note that there is only one transition from s 1 , t 2 and t 3 , respectively):

s 1 , t 1 = min |3 -3| + .90 s 2 , t 2 |3 -4| + .90 s 2 , t 3 = min            .90 max |11 -5| + .90 s 4 , t 4 |5 -5| + .90 s 5 , t 4 1 + .90 max |11 -15| + .90 s 4 , t 5 |5 -15| + .90 s 5 , t 5 = min            .90 max 6 + .90 2 s 1 , t 1 .90 2 s 1 , t 1 1 + .90 max 4 + .90 2 s 1 , t 1 10 + .90 2 s 1 , t 1 = min .90 6 + .90 2 s 1 , t 1 , 1 + .90 10 + .90 2 s 1 , t 1 = 5.4 + .90 3 s 1 , t 1
Hence s 1 , t 1 ≈ 19.9. For the pointwise distance s 1 , t 1 • we have accordingly, Proof. To show that s, s = s, s • = 0 is trivial. The triangle inequalities can be shown inductively; we prove the one for •, • : For s, t, u ∈ S, we have

s 1 , t 1 • = min max |3 -3|, s 2 , t 2 • max |3 -4|, s 2 , t 3 • = min max 0, |11 -5|, s 4 , t 4 • , |5 -5|, s 5 , t 4 • max 1, |11 -15|, s 4 , t 5 • , |5 -15|, s 5 , t 5 • = min max 6, s 1 , t 1 • max 10, s 1 , t 1 • which has least fixed point s 1 , t 1 • = 6.
s, t + t, u = max s n -→s ′ min t m -→t ′ d ❑ (n, m) + λ s ′ , t ′ + max t m -→t ′ min u z - →u ′ d ❑ (m, z) + λ t ′ , u ′ ≥ max s n -→s ′ min t m -→t ′ min u z - →u ′ d ❑ (n, m) + d ❑ (m, z) + λ s ′ , t ′ + t ′ , u ′ ≥ max s n -→s ′ min u z - →u ′ d ❑ (n, z) + λ s ′ , u ′ = s, u
assuming the triangle inequality has been proven for the triple (s ′ , t ′ , u ′ ).

In the next proposition we take the standard liberty of comparing different (weighted) transition systems by considering their disjoint union.

Proposition 7. The weighted transition systems ❖ = ({s 1 }, ∅) and U = ({s 1 }, {(s 1 , n, s 1 ) | n ∈ ❑}) are respectively minimal and maximal elements with respect to both •, • and

•, • • , that is, ❖,A = ❖,A • = A, U = A, U • = 0 for any WTS A.
Proof. For ❖,A and ❖,A • , the maximum max s 1 n -→s ′ 1 is taken over the empty set and hence is 0. For A, U and A, U • , any transition s n -→ s ′ in A can be matched by s 1 n -→ s 1 in U, hence the distance is again 0.

Our distances are related to the usual notion of simulation in the following way. Recall [START_REF] Thrane | Quantitative simulations of weighted transition systems[END_REF]

that a relation R ⊆ S × S is called a -weighted simulation if whenever (s, t) ∈ R and s n -→ s ′ , then also t n -→ t ′ with (s ′ , t ′ ) ∈ R; -unweighted simulation if whenever (s, t) ∈ R and s n -→ s ′ , then also t m -→ t ′ with (s ′ , t ′ ) ∈ R and d ❑ (n, m) < ∞.
State s is weighted simulated, denoted s ≤ t, respectively unweighted simulated, denoted s ≤ u t, if there exists a weighted simulation relation, respectively an unweighted simulation relation, R with (s, t) ∈ R.

Proposition 8. For s, t states in a WTS (S, T ),

-s ≤ t implies s, t = s, t • = 0, -s, t = 0 or s, t • = 0 imply s ≤ t, -s ≤ u t implies that s, t < ∞ and s, t • < ∞, and -s, t < ∞ or s, t • < ∞ imply s ≤ u t. Proof. For the first claim, if (s, t) ∈ R for some weighted simulation relation R, then any s n -→ s ′ can be matched by a transition t n -→ t ′ with (s ′ , t ′ ) ∈ R, hence s, t = 0+λ s ′ , t ′ and s, t • = max 0, s ′ , t ′
• , and we can proceed by induction. Conversely, if s, t = 0, then max

s n -→s ′ min t m -→t ′ d(n, m) + λ s ′ , t ′ = 0, hence d(n, m) = 0,

Axiomatizations for Finite Weighted Processes

We now turn to a setting where our weighted transition systems are generated by finite or regular (weighted) process expressions. We construct a sound and complete axiomatization of simulation distance in a setting without recursion first and show afterwards how this may be extended to a setting with recursion.

Let P be the set of process expressions generated by the following grammar:

E ::= ❖ | n.E | E + E | n ∈ ❑
Here ❖ is used to denote the empty process, cf. Proposition 7.

The semantics of finite process expressions is a weighted transition system generated by the following standard SOS rules:

n.E n -→ E E 1 n -→ E ′ 1 E 1 + E 2 n -→ E ′ 1 E 2 n -→ E ′ 2 E 1 + E 2 n -→ E ′ 2
We can immediately get the following equalities

E + ❖,E = 0 n.E, m.F = d ❑ (n, m) + λ E, F (2)
E 1 + E 2 , F = max( E 1 , F , E 2 , F ) n.E, F 1 + F 2 = min( n.E, F 1 , n.E, F 2 ) (3)
For the point-wise distance, we again need only exchange [START_REF] Černý | Simulation distances[END_REF] with

n.E, m.F • = max(d ❑ (n, m), E, F • )
In order to see e.g. Equation (3) we simply need to apply the definitions:

n.E, F 1 + F 2 = inf F 1 +F 2 m -→F ′ d ❑ (n, m) + λ E, F ′ = min inf F 1 m -→F ′ d ❑ (n, m) + λ E, F ′ , inf F 2 m -→F ′ d ❑ (n, m) + λ E, F ′ = min n.E, F 1 , n.E, F 2
For Equation ( 2), the sup-inf expression ranges over singleton sets, hence the result is easy; the remaining equalities may shown in a similar way.

The inference system F as given in Figure 2 axiomatizes accumulating simulation distance for finite processes, as we shall prove below. Its sentences are inequalities of the form [E, F ] ⊲⊳ r where ⊲⊳ ∈ {=, ≤, ≥} and 0 ≤ r ≤ ∞. Whenever [E, F ] ⊲⊳ r may be concluded from F , we write ⊢ F [E, F ] ⊲⊳ r.

In addition to reflexivity and transitivity, we will need the following standard properties of ⊲⊳ in latter proofs of soundness and completeness: Whenever a ⊲⊳ b then, for all c: a+c ⊲⊳ b+c, a•c ⊲⊳ b•c, and max{a, c} ⊲⊳ max{b, c}.

We also remark that the left process indeed needs to be guarded in rule (R3) above, i.e. the following proposed rule (R3 ′ ) leads to an unsound inference system:

⊢ [E, F 1 ] ⊲⊳ r 1 ⊢ [E, F 2 ] ⊲⊳ r 2 (R3 ′ ) min(r 1 , r 2 ) ⊲⊳ r ⊢ [E, F 1 + F 2 ] ⊲⊳ r
Indeed, using this rule we can derive the following (incomplete) proof tree with a contradictory conclusion; the reason behind is that with E = E 1 + E 2 non-deterministic as below, both F 1 and F 2 may be needed to answer the challenge posed by E:

(A1) 0 ⊲⊳ r ⊢ [❖, E] ⊲⊳ r (A2) ∞ ⊲⊳ r ⊢ [n.E, ❖] ⊲⊳ r ⊢ [E, F ] ⊲⊳ r 1 (R1) d ❑ (n, m) + λr 1 ⊲⊳ r ⊢ [n.E, m.F ] ⊲⊳ r ⊢ [E 1 , F ] ⊲⊳ r 1 ⊢ [E 2 , F ] ⊲⊳ r 2 (R2) max(r 1 , r 2 ) ⊲⊳ r ⊢ [E 1 + E 2 , F ] ⊲⊳ r ⊢ [n.E, F 1 ] ⊲⊳ r 1 ⊢ [n.E, F 2 ] ⊲⊳ r 2 (R3) min(r 1 , r 2 ) ⊲⊳ r ⊢ [n.E, F 1 + F 2 ] ⊲⊳ r Figure 2: The F proof system. ⊢ [1.❖, 1.❖] ≥ 0 ⊢ [2.❖, 1.❖] ≥ 1 ⊢ [1.❖ + 2.❖, 1.❖] ≥ 1 ⊢ [1.❖, 2.❖] ≥ 1 ⊢ [2.❖, 2.❖] ≥ 0 ⊢ [1.❖ + 2.❖, 2.❖] ≥ 1 ⊢ [1.❖ + 2.❖, 1.❖ + 2.❖] ≥ 1 Theorem 9 (Soundness). If ⊢ F [E, F ] ⊲⊳ r, then E, F ⊲⊳ r.
Proof. By an easy induction in the proof tree for ⊢ F [E, F ] ⊲⊳ r, with a case analysis for the applied proof rule: (R2) By induction hypothesis, E 1 , F ⊲⊳ r 1 and E 2 , F ⊲⊳ r 2 , hence

E 1 + E 2 , F = max( E 1 , F , max{E 2 , F }) ⊲⊳ max(r 1 , r 2 ). (R3) By induction hypothesis, n.E, F 1 ⊲⊳ r 1 and n.E, F 2 ⊲⊳ r 2 , hence n.E, F 1 + F 2 = min( n.E, F 1 , n.E, F 2 ) ⊲⊳ min(r 1 , r 2 ). Theorem 10 (Completeness). If E, F ⊲⊳ r, then ⊢ F [E, F ] ⊲⊳ r.
Proof. By an easy structural induction on E:

(E = ❖) We have ❖,F = 0 ⊲⊳ r. By Axiom (A1), also ⊢ [❖, F ] = 0. (E = n.E ′ )
We use an inner induction on F :

Subcase F = ❖. Here E, F = n.E ′ , ❖ = ∞ ⊲⊳ r. By Axiom (A2), also ⊢ [n.E ′ , ❖] = ∞. Subcase F = m.F ′ . Here E, F = n.E ′ , m.F ′ = d ❑ (n, m)+λ E ′ , F ′ ⊲⊳ r, hence with r ′ = λ -1 (r -d ❑ (, nm)
), E ′ , F ′ ⊲⊳ r ′ . By induction hy- pothesis it follows that ⊢ [E ′ , F ′ ] ⊲⊳ r ′ , and we can use Axiom (R1) to conclude that ⊢ [E, F ] ⊲⊳ r.

Subcase F = F 1 +F 2 . Using Equation (3), we have E, F = n.E ′ , F 1 + F 2 = min n.E ′ , F 1 , n.E ′ , F 2 . Let n.E ′ , F 1 ⊲⊳ r 1 , n.E ′ , F 2 ⊲⊳ r 2 . By the previous case, ⊢ [n.E ′ , F 1 ] ⊲⊳ r 1 . As min{r 1 , r 2 } ⊲⊳ r it follows using (R3) that ⊢ F [n.E, F 1 + F 2 ] ⊲⊳ r. (E = E 1 + E 2 )
By an argument similar to the one in the preceding subcase, we have E, F = max E 1 , F , E 2 , F . If E 1 , F ⊲⊳ r 1 and E 2 , F ⊲⊳ r 2 with max(r 1 , r 2 ) ⊲⊳ r, we can use the induction hypothesis to conclude

⊢ [E 1 , F ] ⊲⊳ r 1 and ⊢ [E 1 , F ] ⊲⊳ r 1 , whence ⊢ [E, F ] ⊲⊳ r by Axiom (R2).

Point-wise distance

We can devise a sound and complete inference system F • for point-wise distance (instead of accumulating) by replacing inference rule (R1) in System F by the rule

⊢ [E, F ] ⊲⊳ r 1 (R1 • ) max(d ❑ (n, m), λr 1 ) ⊲⊳ r ⊢ [n.E, m.F ] ⊲⊳ r
As before, we write

⊢ F • [E, F ] ⊲⊳ r if [E, F ] ⊲⊳ r can be proven by F • . Theorem 11 (Soundness & Completeness). ⊢ F • [E, F ] ⊲⊳ r if and only if E, F • ⊲⊳ r
Proof. The proof is similar to the one for F .

Simulation distance zero

We show here that for distance zero, our inference system F specializes to a sound and complete inference system for simulation. The inference system F 0 is displayed in Figure 3.

Theorem 12 (Soundness & Completeness). ⊢ F 0 E ≤ F if and only if E ≤ F
Proof. Soundness follows immediately from the soundness of Proof system F , and for completeness we note that the arguments one uses in the inductive proof of Theorem 10 all specialize to distance zero. We remark that, contrary to the situation for general distance above, we may indeed replace the guarded process n.E in (R3 ′ 0 ) and (R3 0 ) by a plain E without invalidating the rules. Note also that F 0 may similarly be obtained as a specialization F • 0 of the axiomatization F • of point-wise distance above.

(A1 0 ) ⊢ ❖ ≤ E ⊢ E ≤ F (R1 0 ) ⊢ n.E ≤ n.F ⊢ E 1 ≤ F ⊢ E 2 ≤ F (R2 0 ) ⊢ E 1 + E 2 ≤ F ⊢ n.E ≤ F 1 (R3 ′ 0 ) ⊢ n.E ≤ F 1 + F 2 ⊢ n.E ≤ F 2 (R3 0 ) ⊢ n.E ≤ F 1 + F 2

Axiomatizations for Regular Weighted Processes

Let N = max{d ❑ (n, m) | n, m ∈ ❑}; by finiteness of ❑, N ∈ |. Let V
be a fixed set of variables, then P R is the set of process expressions generated by the following grammar:

E ::= U | X | n.E | E + E | µX.E n ∈ ❑,X ∈ V
Here we use U to denote the universal process recursively offering any weight in ❑, cf. Proposition 7. Note that we do not incorporate the empty process ❖. Semantically this will mean that all processes in P R are non-terminating, and that the accumulating distance between any pair of processes is finite.

The reason for the exchange of ❖ with U is precisely this last property; specifically, completeness of our axiomatization (Theorem 17) can only be shown if all accumulating distances are finite. The semantics of processes in P R is given as weighted transition systems which are generated by the following standard SOS rules:

n.E n -→ E U n -→ U E 1 n -→ E ′ 1 E 1 + E 2 n -→ E ′ 1 E 2 n -→ E ′ 2 E 1 + E 2 n -→ E ′ 2 E[µX.E/X] n -→ F µX.E n -→ F
As usual we say that a variable X is guarded in an expression E if any occurrence of X in E is within a subexpression n.E ′ . Formally, we define the guarding depth gd(E, X) of variable X in expression E recursively by

gd(U, X) = ∞ gd(X, X) = 0 gd(n.E, X) = 1 + gd(E, X) gd(E 1 + E 2 , X) = min gd(E 1 , X), gd(E 2 , X) gd(µX.E, Y ) = gd(E, X) if X = Y ∞ if X = Y
and we say that X is guarded in E if gd(E, X) ≥ 1. Also as usual, we denote by E[F/X] the expression derived from E by substituting all free occurrences of variable X in E by F , and given tuples

F = (F 1 , . . . , F k ), X = (X 1 , . . . , X k ), we write E[ F / X] = E[F 1 /X 1 , . . . , F k /X k ]
for the simultaneous substitution.

Our inference system for regular processes consists of the set of rules R as shown in Figure 4; whenever [E, F ] ⊲⊳ r may be concluded from R, we write ⊢ R [E, F ] ⊲⊳ r.

Compared to inference system F for finite processes, we note that we have to include the triangle inequality (R4) as an inference rule. Also, the precongruence property of simulation distance is expressed by rules (R1), (R5), and (R6). We will need all those extra rules in the proof of Lemma 14 which again is necessary for showing completeness.

Theorem 13 (Soundness). For closed expressions E, F ∈ P R we have that ⊢ R [E, F ] ⊲⊳ r implies E, F ⊲⊳ r.

Proof. By an easy induction in the proof tree for ⊢ F [E, F ] ⊲⊳ r, using the definition of •, • . In relation to Axiom (A3), we note that

N = max{d ❑ (n, m) | n, m ∈ ❑} implies E, F ≤ ∞ i=0 λ i N = N
1-λ . Our completeness result for regular processes will be based on the following lemmas; here we call an expression E ∈ P R non-recursive if it does not contain any subexpressions µX.E ′ : Lemma 14. For all E ∈ P R and k ∈ ◆ there exist a non-recursive expression F and tuples Ē = (E 1 , . . . , E k ), X = (X 1 , . . . , X k ) for which gd(F, X i ) ≥ k for all i and

⊢ R [E, F [ Ē/ X]] = 0 ⊢ R [F [ Ē/ X], E] = 0
Proof. Repeated use of the unfolding axioms (A6) and (A7), the congruence rules (R1), (R5), and (R6) with r = 0 and of the triangle inequality (R4).

We are now in a position to state our completeness result which enables arbitrary ǫ-close proofs in the sense below. The proof uses unfoldings of recursive expressions as in Lemma 14, and as these unfoldings are finite non-recursive processes, we cannot expect exact completeness.

Theorem 17 (Completeness up to ǫ). Let E and F be closed expressions of P R and ǫ > 0.

Then E, F = r implies ⊢ R [E, F ] ≤ r + ǫ and ⊢ R [E, F ] ≥ r -ǫ.
Proof. Assume E, F = r, and choose k ∈ ◆ such that 2λ k N 1-λ ≤ ǫ. By Lemma 14 we have non-recursive expressions E ′ , F ′ and tuples Ē, F , X, and Ȳ for which gd(E ′ , X i ) ≥ k and gd(F ′ , Y i ) ≥ k for all i, and such that

⊢ R [E, E ′ [ Ē/ X]] = 0 ⊢ R [E ′ [ Ē/ X], E] = 0 ⊢ R [F, F ′ [ F / Ȳ ]] = 0 ⊢ R [F ′ [ F / Ȳ ], F ] = 0 From Lemma 15 it follows that ⊢ R [E ′ [ Ē/ X], E ′ [ Ū/ X]] = 0 ⊢ R [E ′ [ Ū/ X], E ′ [ Ē/ X]] ≤ λ k N 1-λ = ǫ 2 ⊢ R [F ′ [ F / Ȳ ], F ′ [ Ū/ Ȳ ]] = 0 ⊢ R [F ′ [ Ū/ Ȳ ], F ′ [ F / Ȳ ]] ≤ λ k N 1-λ = ǫ 2
Using the triangle inequality and Theorem 13 we now have

E ′ [ Ū/ X], F ′ [ Ū / X] ≤ E ′ [ Ū/ X], E ′ [ Ē/ X] + E ′ [ Ē/ X], E + E, F + F, F ′ [ F / Ȳ ] + F ′ [ F / Ȳ ], F ′ [ Ū/ Ȳ ] ≤ ǫ 2 + 0 + r + 0 + 0 = r + ǫ 2
Only non-recursive expressions are involved here, so that we can invoke Lemma 16 to conclude

⊢ R [E ′ [ Ū/ X], F ′ [ Ū / X]] ≤ r + ǫ 2
Now we can use the triangle inequality axiom (R4) together with the eight equations above to arrive at

⊢ R [E, F ] ≤ r + ǫ Similar arguments show that also ⊢ R [E, F ] ≥ r -ǫ,

Point-wise distance

Again we can easily convert our proof system R into one for point-wise (instead of accumulating) distance. In this case, we obtain R • by replacing inference rule (R1) by (R1 • ) as we did for Proof system F , and (A3) needs to be replaced by

(A3 • ) ⊢ [E, F ] ≤ N
With these replacements we have a sound and ǫ-complete axiomatization of point-wise simulation distance recursive weighted processes.

Theorem 18 (Soundness & Completeness up to ǫ). Let E and F be closed

expressions of P R , then ⊢ R • [E, F ] ⊲⊳ r implies E, F • ⊲⊳ r, and E, F • = r implies ⊢ R • [E, F ] ≤ r + ǫ and ⊢ R • [E, F ] ≥ r -ǫ for any ǫ > 0.
Proof. The proof is similar to that for accumulated distance.

Algorithmic Complexity

In this section we show that, for finite weighted transition systems, computing accumulating distance is polynomial-time equivalent to computing the value of discounted games (DG), hence contained in NP ∩ coNP.We also give a polynomial-time algorithm for computing point-wise distance; hence the conceptually simpler point-wise distance is also computationally easier.

We assume throughout this section that for all weights n, m ∈ ❑, d ❑ (n, m) is polynomial-time computable.

We recall the following definition from [START_REF] Zwick | The complexity of mean payoff games[END_REF]: A two-player game graph G = (S 1 , S 2 , -) over a finite set W ⊆ | of weights is a finite directed bipartite graph with vertices S 1 ∪ S 2 (where states in S i are said to belong to Player i) and edges -∈ S 1 × W × S 2 ∪ S 2 × W × S 1 , in which each vertex has at least one outgoing edge. A memoryless Player-i strategy is a map ξ : S i → W × S 3-i , and it is consistent withif ξ(s) = (a, s ′ ) implies that there exists (s, a, s ′ ) ∈ -(the latter written as s a s ′ ). A pair (ξ, χ) of Player-1 and Player-2 strategies admits a unique sequence (path) of edges t = e 0 e 1 e 2 . . . . We will sometimes write a = w(e) to denote the weight of e = (s, a, s ′ ) ∈ -.

Definition 19. Let G = (S 1 , S 2 , -) be a game graph and 0 ≤ λ < 1 a discounting factor, and s ∈ S 1 ∪ S 2 . The payoff of the discounted game on G from s is defined as v = λ ∞ i=0 λ i w(e i ), where the path t = e 0 e 1 e 2 . . . , starting from s ∈ S 1 , is induced by strategies ξ and χ of Player 1 respectively Player 2 which are such that the objective of Player 1 is to pick ξ maximizing v, whereas Player 2's objective is to pick χ minimizing v.

The decision problem corresponding to discounted games is as follows: Given a game graph G, a starting vertex s 0 , a discount factor 0 ≤ λ < 1, and a threshold v ∈ |, can Player 1 guarantee a payoff of v or more? We use the well-known fact [START_REF] Zwick | The complexity of mean payoff games[END_REF] that this decision problem is in NP ∩ coNP, a fact obtained by reduction to simple stochastic games [START_REF] Condon | The complexity of stochastic games[END_REF].

Theorem 20. Computing accumulating distance is polynomial-time equivalent to computing the payoff for discounted games.

We present two supporting lemmas. Note that our weighted transition systems and discounted games have different weight domains; we use the metric d ❑ to map between them.

Lemma 21. For a given discount factor λ, a WTS (S, T ), and s, t ∈ S, one can construct a game graph G together with a vertex s 0 such that s, t is the payoff of the discounted game on G from s 0 with discount factor √ λ.

Proof. We construct the game G = (S 1 , S 2 , -) with S 1 = S × S, S 2 = S × S × ❑, and edges

(s, t) 0 -(s ′ , t, n) if (s, n, s ′ ) ∈ T (s, t, n) d ❑ (n,m) -----(s, t ′ ) if (t, m, t ′ ) ∈ T
Note that the set of weights of the game is W ⊆ {d ❑ (n, m) | n, m ∈ ❑}, hence finite. Now recall [START_REF] Zwick | The complexity of mean payoff games[END_REF]Thm. 5.1]: For some labeling S 1 ∪ S 2 = {s 1 , . . . , s p }, the payoff x i of G from s i is given as the unique fixed point of the set of equations

x i =        max s i c -s j {c + √ λx j } if s i ∈ V 1 min s i c -s j {c + √ λx j } if s i ∈ V 2
Hence we can let s 0 = (s, t), then s, t = x 0 .

Lemma 22. For any game graph G with start vertex x there exists a WTS with states p x and q x such that the payoff v of the discounted game on G from x is p x , q x .

Proof. Let G = (S 1 , S 2 , -) be the game, with weight set W . Then (S, T ), x models the choice in ξ, ensuring that minimization from q x must match the b-label, acknowledging the move to b ∈ S 2 . Doing so requires taking precisely q x b,0 -→ q ′ y for some y ∈ S 1 . Whenever χ(b) = y ′ minimizes the payoff for Player 2 in vertex b, the corresponding q x b,0 -→ q ′ y ′ ensures 1. the correct cost for the match, i.e. |n -0| = w(xb), and 2. that the maximizing transition from

with (finite) weight set ❑ = (S 1 ∪ S 2 ) × W is defined as follows: S = {p x , q x , q ′ x | x ∈ S 1 } ∪ {p b x | x ∈ S 1 , x -b},
p b x is p b x y,0 --→ p y ′ .
The latter implies that the choice of Player 2 is passed on to Player 1, and that Player 1 must act according to it. This latter property is obtained by the q ′

x S 1 \{y},0 -----→ U transitions, which match any "cheating" maximization and afterwards allow Player 2 to match all possible Player 1 transitions perfectly.

After the Player 1 transition p b x y,0 --→ p y ′ , Player 2 must minimize the game value by choosing a q ′ y ′ y,m --→ q y ′ , thereby adding |0 -m| = w(by ′ ) both to the total value of G and to the total accumulating distance from p x to q x . We have arrived at a new configuration p y ′ , q y ′ which models the vertex y ′ of the game and from which the simulation game can proceed.

Example. We shortly elaborate on the reduction in the preceding proof by considering the game in Figure 5. Assume Player 1 (with diamond-shaped vertices) has an optimal (maximizing) strategy in which A → (n 1 , B) and Player 2 (with square vertices) wants to play B → (n 4 , E) for minimizing the game value. Then the simulation from (p A , q A ) is performed accordingly, Proof of Theorem 20. The result follows directly from Lemmas 21 and 22, together with the additional fact that both reductions are clearly polynomial.

The sought-after property follows:

Corollary 23. The decision problem corresponding to computing accumulating simulation distance of states in a WTS is contained in NP ∩ coNP.

Point-wise distance

To see that point-wise simulation distance is computable in polynomial time, we note that the fixed-point iteration converges in time polynomial in the size of the WTS:

Theorem 24. For a WTS (S, T ) with |S| = p and s 0 , t 0 ∈ S, s 0 , t 0 • may be computed in p 2 steps. The lattice of ≤ k δ relations such defined has at most p 2 elements, hence ≤ δ = ∞ k=0 ≤ k δ can be computed in at most p 2 iterations. To finish the proof, we note that s 0 , t 0 • ≤ δ if and only if s 0 ≤ δ t 0 .

Corollary 25. The decision problem corresponding to computing point-wise simulation distance is contained in P.

[1] Jan A. Bergstra, Alban Ponse, and Scott A. Smolka, editors. Handbook of Process Algebra. North-Holland, 2001.

4 5 15 Figure

 15 Figure 1: Example WTS

3. 4 .

 4 Properties Proposition 6. •, • and •, • • are hemimetrics on S.

(

  A1) follows from ❖,E = 0. (A2) follows from n.E, ❖ = ∞ which is clear by the definition of •, • . (R1) By induction hypothesis, E, F ⊲⊳ r 1 , and as n.E, m.F = d ❑ (n, m)+ λ E, F , it follows that n.E, m.F = d ❑ (n, m) + λr 1 .

Figure 3 :

 3 Figure 3: The F0 proof system

  and for x, y ∈ S 1 and b ∈ S 2 , T is given by: In the last item, U denotes the universal WTS, and the notation q ′ y S 1 \{y},0 -----→ U means that there are transitions q ′ y z,0--→ U for all z ∈ S 1 \ {y}.The metric on ❑ is defined by d ❑ (x, n), (y, m) = |n -m| if x = y and ∞ otherwise. The construction is sketched in Figure5.To see that p x , q x = v, consider a strategy ξ for which ξ(x) = b which maximizes Player 1's payoff from x ∈ S 1 . The transition p x b,n --→ p b

-→ 4 -

 4 i.e. the maximal choice is p A b,n 1 --→ p B A signaling the label b (alternatively c or d could have been chosen). The simulating response is, of course, to match b (since taking e.g. c results in weight d ❑ ((b, n 1 )(c, 0)) = ∞); the transition corresponding to Player 2's minimizing choice B n -E is q A b,0 -→ q ′ E . The next maximizing challenge then has to be p B A e,0 --→ p E , as e.g. p B A f,0 --→ p F would 17

Figure 5 :

 5 Figure 5: DG to WTS Translation (dashed lines represents omitted parts)

Proof.

  For a threshold δ, define iterated δ-simulation relations ≤ n δ , for n ∈ ◆, by ≤ 0 δ = S × S ands ≤ k+1 δ t iff ∀s m -→ s ′ ∃t n -→ t ′ : d ❑ (m, n) ≤ δ and s ′ ≤ k δ t ′

  and as d is a metric, n = m, and we can proceed by induction. The situation if s, t • = 0 is similar.

For the third claim, if (s, t) ∈ R for some unweighted simulation relation R, then any s n -→ s ′ can be matched by a transition t m -→ t ′ with (s ′ , t ′ ) ∈ R, and as d(n, m) < ∞, we can again proceed by induction. Conversely, s, t < ∞ or s, t • < ∞ imply that any transition s n -→ s ′ has a match t m -→ t ′ .

Proof. Repeated use of Axioms (A3) and (A8) together with the congruence rules (R1), (R5), and (R6) with r = 0.

Lemma 16. For closed non-recursive expressions E, F , E, F ⊲⊳ r implies ⊢ R [E, F ] ⊲⊳ r.

Proof. By structural induction similar to the proof of Theorem 10.


Alessio Colombo

Daniele Fontanelli

Dhaval Gandhi

Antonella De Angeli

Luigi Palopoli

Sean Sedwards

Axel Legay

Behavioural Templates Improve Robot Motion Planning with Social Force Model in Human Environments *

An accurate model of human behaviour is crucial when planning robot motion in human environments. The Social Force Model (SFM) is such a model, having parameters that control both deterministic and stochastic elements. We have constructed an efficient motion planning algorithm by embedding the SFM in a control loop that determines higher level objectives and reacts to environmental changes. Low level predictive modelling is provided by the SFM fed by sensors; high level logic is provided by statistical model checking.

To parametrise and improve our motion planning algorithm, we have conducted experiments to consider typical human interactions in crowded environments. We have identified a number of behavioural patterns which may be explicitly incorporated in the SFM to enhance its predictive power. In this paper we describe the results of these experiments and how we parametrise the SFM.

Introduction

Many new applications for robots require them to operate in the proximity of human beings. Examples include robots that execute daily activities in homes or public spaces (e.g., automatic cleaners), robots used to assist older adults or disabled people and a new generation of toys. In contrast to the traditional notion of industrial robots, whose areas of operation must be confined to prevent accidents [START_REF] Asaula | Safety provisions for human/robot interactions using stochastic discrete abstractions[END_REF], these new devices are autonomous and may interact with humans. In this context, the ability to understand human behaviour and social interactions is of primary importance.

Our research is motivated by the activities of the DALi project [START_REF]The DALi Project[END_REF], whose goal is to develop an autonomous device to assist those with reduced ability. In particular, we are developing a 'cognitive engine' (an intelligent motion planner) to help the user reach places of interest in crowded public spaces, such as shopping malls, airports and railway stations. This can be challenging and stressful for those with reduced mental or physical ability.

Our motion planner employs a dynamic model of the local environment to hypothesise trajectories of moving objects. For a given set of alternative actions, these predictions are analysed by higher level logic to decide the user's best course. The predictive model is based on the Social Force Model (SFM) [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF][START_REF] Helbing | Freezing by heating in a driven mesoscopic system[END_REF], using knowledge of prior trajectories and current sensor information. The control loop employs statistical model checking, with the user's objectives and constraints encoded in temporal logic [START_REF] Colombo | Motion Planning in Crowds using Statistical Model Checking to Enhance the Social Force Model[END_REF]. An overview of the motion planner is shown in Fig. 1.

The SFM assumes that human motion respects the physical laws of motion and social conventions that may be modelled as forces. The SFM considers fixed objects (e.g., walls), moving objects (e.g., pedestrians), attractive forces (e.g., between friends or towards points of interest), repulsive forces (e.g., between strangers or away from crowds), frictional forces (when objects collide) and random terms that model uncertainty (unpredictable behaviour). There are alternative models of pedestrian dynamics (e.g., [START_REF] Yu | Centrifugal force model for pedestrian dynamics[END_REF]), however the SFM suits our needs: it is flexible, efficient and is widely adopted in the literature (e.g., [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF][START_REF] Helbing | Freezing by heating in a driven mesoscopic system[END_REF][START_REF] Luber | People tracking with human motion predictions from social forces[END_REF]). Its strengths and weaknesses are well tested and understood [START_REF] Lakoba | Modifications of the helbing-molnar-farkas-vicsek social force model for pedestrian evolution[END_REF][START_REF] Wei-Guo | Evacuation behaviors at exit in CA model with force essentials: A comparison with social force model[END_REF][START_REF] Parisi | A modification of the social force model can reproduce experimental data of pedestrian flows in normal conditions[END_REF]. A mathematical description of the SFM is given in Section 2.

We and our collaborators1 have built prototypes to demonstrate the proof of concept of the scheme shown in Fig. 1. We have verified the efficiency of the algorithm on embedded computing hardware and confirmed that the mobile sensor technology can provide the required information. In this paper we focus on improving the social force model. Our goal is to increase its predictive accu-racy and reduce the need for reactive corrections.

The motion planning challenge

The SFM is a model of continuous interaction (mutual reaction), driven by the desired trajectories of moving agents (denoted v 0 (t), see Section 2). When using the SFM in a generative context (e.g., [START_REF] Helbing | Freezing by heating in a driven mesoscopic system[END_REF]), v 0 (t) may be specified in advance. The desired trajectories of real pedestrians, however, are a function of their objectives and the instantaneous positions and velocities of other pedestrians. In our motion planning application we are able to measure such positions and velocities, but only the objectives of the user are known with any certainty. The intentions of other pedestrians must be inferred from their trajectories. Noting that humans tend to walk in straight lines, we approximate the desired trajectories of other pedestrians piecewise, using short term linear extrapolations of their most recently detected motion. This is sufficient when the motion is smooth and the prediction timescale is short. In reality, the relatively smooth long and medium term trajectories of pedestrians are often interrupted by sudden short term pauses and deviations.

To improve the performance of the SFM, we have therefore investigated in detail the behaviour of human participants in a simulated shopping environment. We have defined an experimental procedure comprising a concurrent verbal protocol (to identify the motivation of the choices made by the participants), video recording and motion tracking using the Kinect platform [START_REF]Kinect for Xbox 360[END_REF]. The reconstructed trajectories of the participants were manually annotated with corresponding motivations, thus identifying a set of behavioural patterns. The data corresponding to each pattern were then used to parametrise the SFM using standard algorithms.

The results of our investigation show that, for some patterns, parametrisation of our existing model is sufficient to produce a good reconstruction of observed behaviours. In other cases, we find that it will be necessary to incorporate the patterns as explicit modifications to the SFM. Since the patterns have recognisable signatures, we conclude that we can construct an improved SFM using behavioural templates.

The social force model

Following [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], our model is constructed in two dimensions, with human agents represented by circular discs and fixed objects represented by lines. In what follows we denote vectors in bold type. Thus, agent i has mass m i centered at position x i ∈ R 2 in the environment, radius r i and velocity v i ∈ R 2 . The SFM is described by a system of linear differential equations

{ ẋi = v i vi = v 0 i -vi τi + fi+ξi mi (1)
v 0 i is the driving (desired) velocity of agent i, represented by a product of speed v 0 i and normalised direction e 0 i . In our algorithm e 0 i is given by the line joining the current position and the next via point. Importantly, since v 0 i is by default set to the user's preferred walking speed, v 0 i is time invariant between via points. τ i is the time taken to react to the difference between desired and actual velocity, while ξ i is a noise term modelling fluctuations not accounted for by the deterministic part of the model. The noise term can also serve to avoid deadlocks and hypothesise alternative trajectories. In our implementation we assume ξ i is normally distributed. In the absence of the exogenous inputs f i and ξ i , the agent's trajectory simply converges to the driving velocity with time constant τ i . f i is the overall force acting on agent i resulting from other objects in the environment and is 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 strangers from getting too close, f att ij is the attractive social force that, e.g., brings friends together, f ph ij is the physical force that exists when two people come into contact. The second term includes the forces acting on agent i as a result of fixed environmental obstacles (e.g., walls): f soc ib is the social force that inhibits agent i from getting too close to the boundaries, f ph ib is the physical force that exists when agent i touches the boundary b. Finally, f att ic is the attractive social force that draws agent i towards fixed objects of incidental interest (shops, cafés, toilets, etc.).

f is principally a function of the distance between an agent and the other objects in the model. d ib is the minimum distance between the circumference of agent i and fixed object b. d ij is the distance between the centres of mass of agents i and j, i.e., the centres of the discs, while r ij = r i + r j is the "touching distance". To aid modelling the different force regimes that exist when agents are not in contact and when they touch (i.e. agents i and j touch if

r ij -d ij ≤ 0) we adopt the function Θ(r ij , d ij ) = max(0, r ij -d ij ).
Using these notions, the various repulsive social and physical forces of (2)) are defined as follows:

f soc ij = {A i exp[(r ij -d ij )/B i ]} n ij Λ(λ i , ϕ ij ) (3)
f ph ij = k 1 Θ(r ij -d ij )n ij + k 2 Θ(r ij -d ij )∆v t ji t ij (4) f soc ib = {A i exp[(r i -d ib )/B i ] + k 1 Θ(r i -d ib )} n ib (5) f ph ib = -k 2 Θ(r i -d ib )(v i • t ib )t ib (6)
n ij (n ib ) is a normalised vector pointing from agent j (fixed object b) to agent i, i.e., the direction of the repulsive force. t ij (t ib ) is a normalised vector tangential to the relative movement of agent i and agent j (fixed obstacle b), i.e., the motion tangential direction. ∆v t ji = (v jv i ) • t ij is the tangential velocity difference. The social forces (3) and ( 5) increase exponentially with reducing distance between objects, with a scale defined by constants A i and B i . In particular, A i is the force acting on agent i at the touching distance; B i is loosely the distance at which the force takes effect.

Λ : R 2 → [0, 1] is a function that gives greater weight to the social force (3) arising from the agents in front of (notionally, seen by) an agent. λ i is a parameter that regulates the effect of Λ on agent i, while ϕ ij is the angle between the directions e 0 i and -n ij , i.e., the field of view of the agent. The physical force (4) between agents comprises a repulsive body compression force (first term) that acts in direction n ij , plus a frictional force (second term) that acts in direction t ij to impede the relative tangential movement of two agents in contact. k 1 and k 2 are constants that define the scale of the physical forces. The physical force ( 6) between an agent and a fixed object is solely described by a frictional term.

Qualitative studies

A qualitative user study was conducted to collect information about how people behave in crowded environments. We considered three main themes: (i) how people interact in confined environments; (ii) how people negotiate shared space and (iii) how people behave with other people in shared space. A simplified shopping context was constructed in a laboratory, as illustrated in Fig. 2. Our experimental procedure involved the use of a concurrent protocol, video-recorded observation and questionnaire administration. These techniques have been found to be particularly effective when conducting experimental investigations [START_REF] Jarke | Scenario management: An interdisciplinary approach[END_REF][START_REF] Rouse | On looking into the black box: Prospects and limits in the search for mental models[END_REF]. The concurrent or talk-aloud protocol -a narration of thought and action during action -was chosen as literature suggests the alternative retrospective protocol (where participants return to view and comment upon their recorded experience) may not accurately reveal participants actual task performance experience. Concurrent protocol participants have been found to go into greater detail and provide more in-depth evaluations [START_REF] Teague | Concurrent vs. post-task usability test ratings[END_REF].

Procedure

The participant sample consisted of 25 University students. Five participants at a time were involved in each experiment, and assigned a different list of four shopping items to collect. They were asked to collect all items in the shortest time possible, verbally explaining all their actions. Participants were specifically arranged around the shopping environment and each of the five shopping lists were unique, to maximise shared space interaction and limit the possibility of processional behaviour. The potential route-behaviour of participants 1 and 4 is illustrated in Fig. 2(a). Participants were asked to collect the items on their list and return to their starting points, while simultaneously verbalising their shopping experience; describing what they saw, where they went, and what they were thinking and doing. The voice-recording functionality of five HUAWEI U8650 Android mobile phones was utilised to capture the speech of the participants. Three Logitech Quickcam Pro 9000 webcams were used to observe and record shared space interaction within the simulated shopping environment. The cameras were coordinated to provide a more encompassing view of the interaction space (Fig. 2(b)). On completion of the experiment, the participants were asked to complete an open questionnaire.

Data analysis

The transcribed concurrent protocols, once synchronized with the video footage and in conjunction with the development of relevant coding schemes, facilitated assessment of a number of variables, including Task Completion Time, Number of Steps taken, and Number of Critical Instances -instances where physical, visual and/or auditory reference was made within the verbal protocols to agent-agent interaction. Detailed video analysis was conducted via the Elan software tool [START_REF] Wittenburg | Elan: a professional framework for multimodality research[END_REF].

Results

Using the Critical Instances markers from the transcribed protocols synchronised to the video footage, we documented and studied the behaviour exhibited by participants when involved in shared space interaction. Critical instances were defined as physical, visual and/or auditory references within the verbal protocols to agent-agent awareness or interaction. Analysis of the identified critical instances revealed a number of common themes and behaviours. The themes that emerged indicated that the behaviour itself was usually employed to either negotiate shared space interaction or avoid collisions within the shared space, and it was possible to categorise the main types of behaviours in two groups: Active and Reactive behaviours.

Active behaviours, in this instance, were considered to be behaviours employed to understand the environment and determine goal strategies toward task completion, environmental awareness and negotiation, mainly focussing around the visual modality: Table 1 provides an example of how the critical instances were interpreted in terms of the taxonomy of behaviours observed during shared space interaction. By examining critical instances within the complete video footage, it was possible to observe the interactional behaviour occurring during such periods according to the

Crit. Inst. Observed Behaviour of Participant (Px)

1 Px aware of other agent in desired space, waits until agent completes task (behaviours (E), (H), (C)) 2

Backward glance locates agent directly behind Px and helps avoid collision (behaviour (B)) 3

Px aware of an agent in front and another approaching agent (behaviours (B), (C)) Px watches the actions of a moving agent and remains still until the agent passes (behaviours (E), (H))

Table 1. Critical instances interpreted in terms of the taxonomy of behaviours. taxonomy, and to determine the frequency with which these behaviours occurred (Fig. 3). As can be observed in Fig. 3, participants utilised the visual modality to engage in either scanning of the overall environment (usually to identify product locations within it), or peripheral vision to monitor the orientational movement of other agents in the environment. As proximity to others in the environment increased, eye-to-eye contact or negotiation between participants was observed to occur. Similarly, under these circumstances, verbal negotiation was also evident in a small number of cases, although this may have been reduced due to the method of recording the concurrent protocols. The most commonly occurring reactive behaviours according to the study were waiting for free space to become available, stationary agents giving priority to moving agents, moving to the left and right, and moving forwards and backwards to create free space, based upon the interactions occurring during the themes defined.

(A) (B) (C) (D) (E) (F) (G) (H) (I) (J) (A) (B) (C) (D) (E) (F) (G) (H) (I) (J)
From the analysis of the videos it was evident that people planned their movement based on their physical distance from others, broadly in line with the assumptions of the SFM. However, we observed modes of behaviour that are not explicitly modelled by the SFM. When somebody came too close to another person, avoidance behaviours were manifested, such as waiting (for the space to be free) or stepping back. These behaviours were clustered in 5 main rules, as reported in Table 2. It is also interesting to note that when two active agents meet, the slower one typically gives way to the faster.

Our study highlights the fundamental role of proxemics [START_REF] Hall | Proxemics -the study of man's spatial relations[END_REF] in human motion. A literature exists concerning psy-chological aspects of spatial behaviour, taking into consideration concepts such as proximity, body orientation, motion in a physical setting, territorial behaviour and privacy. These concepts can provide a theoretical framework to extend the SFM, informing our understanding of people's behaviour in public spaces.

Parametrising the SFM

We conducted a number of motion tracking experiments using the simulated shopping environment described in Section 3. We performed 20 experiments considering two people in the environment and a further 20 experiments using four people. Participants were arranged in specific places around the shopping environment and asked to move according to a shopping list provided in advance. Each shopping list contained a set of places to be visited (via points) in a predefined order, to prompt interactions between participants and thus generate interesting social behaviours. Depending on the number of participants, shopping lists were specifically designed to maximise shared space interaction and limit the possibility of processional behaviour. We recorded video and 3D information of the participants' trajectories using the Kinect sensing apparatus [START_REF]Kinect for Xbox 360[END_REF].

Parameter estimation

Once the experimental data were collected, we decomposed the participants' trajectories according to the SFM, considering situations where there were at least two agents that interacted, having simultaneously crossing trajectories. Since each participant had been asked to travel through a sequence of via points and to pause in each starting and ending position, it was possible to identify the relevant segments of their trajectories in the experimental data. By assuming that the driving velocity of agents is constant along a line joining the initial and desired positions (see Section 2) we were able to infer v 0 .

Due to the constraint imposed by the initial and final configurations, for the purposes of parameter estimation we construct a modified version of the SFM described in Section 2. We first consider the original model (1) and define x 0 i = [x 0 i , y 0 i ] T and v 0 i = [v 0 xi , vy i 0 ] T = [0, 0] T to be the desired final position and velocity, respectively, of agent i. We thus define the Cartesian position error variables xi = x 0 ix i , whose dynamics are given by ẋi = -v i . Furthermore, we define the polar coordinates

ρ i = √ xT i xi , α i = arctan ( ỹi xi ) , (7)
which are respectively the distance and the orientation from the current to the desired position of agent i. We then find that the normalised direction e 0 i turns out to be

e 0 i = [ e 0 xi e 0 yi ] =   xi √ xT i xi ỹi √ xT i xi   = [ c αi s αi ] , (8)
where we adopt the convention c αi = cos(α i ) and s αi = sin(α i ). Hence, the modified SFM is finally given by the dynamical system

   ρi αi vxi vyi    =      -(c αi v xi + s αi v yi ) sα i vx i -cα i vy i ρi v 0 i cα i -vx i τi + fx i +ξx i mi v 0 i sα i -vy i τi + fy i +ξy i mi      . ( 9
)
Using ( 7), ( 8) and ( 9) it is then possible to have an approximate description of the participant accelerations, i.e.,

v i (t k+1 ) -v i (t k ) t k+1 -t k ≈ g i (t k ) + g i (t k+1 ) 2 (10)
where

g i (t k ) = v 0 i e 0 i (t k ) -v i (t k ) τ i + f i (t k ) + ξ i (t k ) m i .
Note that calculating the mean of successive values of g i (•) gives a constant mean acceleration in the sampling period t k+1 -t k .

Having set the desired speed of each agent to v 0 i = max{∥v i (t k )∥}, we estimated values for parameters A i and B i in (3) and ( 5), and for parameter τ i in [START_REF]The DALi Project[END_REF]. We note here here that in our experiments the agents never touch, so the constants k 1 and k 2 are not used and are not estimated. Moreover, the Λ(•) function is for the moment assumed to be equal to 1, without loss of generality, since the trajectories do not take into account interactions from behind. Finally, the mass of the participants m i plays only the role of a factor for the generated forces, hence it is of no relevance to the problem at hand and can be considered as known in advance and removed form the estimation process.

Estimation algorithm

According to the description of the forces and the assumption that the model noise term ξ i is normally distributed, we adopt an iterative Weighted Least Squares (WLS) algorithm to identify the parameters of the dynamical system. Such a choice is justified by the limited number of parameters involved in the estimation and by the relatively small amount of noise in determining the positions of the participants. This method is alternative to [START_REF] Johansson | Specification of the social force pedestrian model by evolutionary adjustment to video tracking data[END_REF], in which a genetic algorithm based solely on video tracking data is used to estimate parameters. In our case, we make use of the high quality output of the Kinect sensor.

Results

Fig. 4 illustrates typical instantaneous output from an experiment: Fig. 4 used in (3), ( 4), ( 5) and ( 6) -and identified by a unique number for tracking purposes.

Fig. 5 illustrates typical global trajectories, reconstructed from the tracked positions of the participants. The global trajectories were divided into local trajectories describing the motion between the pre-defined via points. Even without time information, the locations of the via points are clear in Fig. 5. The local trajectories were then used to estimate parameters of the SFM, to better "predict" the actual trajectories.

In Fig. 6(a) the local trajectories of two participants going in opposite directions are represented by solid lines. The dashed lines represent the trajectories generated by the parametrised SFM. The figure demonstrates the good agreement of the SFM with reality in situations where only relatively small corrections are necessary and there is little conflict.

In Fig. 6(b) the interaction is more conflictual and the negotiation of shared space increases in complexity. In this case our simple linear extrapolation to infer the participants' desired trajectories is not adequate and the SFM's prediction is poor. Interestingly, the point at which the actual and predicted trajectories diverge (identified by a circle in the figure) is the point at which one agent stops to give way to the other. We have found such active behaviour, i.e., moving to facilitate the motion of other agents or waiting for free space, to be most frequent in our experimental context (see Fig. 3).

Conclusions and prospects

Our unmodified social force model is accurate under conditions of continuous flow and when there is little competition for space. Under other circumstances, common in the social environments that are of interest to us, human motion is punctuated by frequent stops and starts, and short term changes of direction. Our experiments have demonstrated that these discrete events are not well captured by a simple social force model. We have nevertheless identified several distinguishable behavioural patterns that may be incorporated as templates into an enhanced SFM. In such a model, these patterns would be recognised and predicted, modifying the desired trajectories of modelled pedestrians. We have shown that this increases the We note that the patterns we have identified might be explained by a more general theory of human interaction. We therefore propose to explore proxemics theory [START_REF] Hall | Proxemics -the study of man's spatial relations[END_REF], which has been developed to explain behaviours of people when they interact in close proximity. We also recognise that our investigations thus far have considered only a limited range of social situations. Key to our ongoing research will be an extensive collection of experimental data that can inform our theory and provide evidence of its applicability and robustness.

Figure 1 .

 1 Figure 1. The DALi motion planner.

  (a) Diagrammatic scheme of participants and shops showing potential trajectories. (b) Video grab of laboratory set-up.

Figure 2 .

 2 Figure 2. Experimental set-up.

  (A) Eye-to-eye negotiation of immediate shared space; (B) Use of peripheral vision in assessment; (C) Visual scanning of environment; (D) Verbal interaction. Reactive behaviours were considered as the reactions of agents in the environment to other agents; the physical reactive movements made to accommodate other agents and successful interaction: (E) Waiting for space/desired location to become clear; (F) Stepping backwards to allow others more room; (G) Moving forwards to allow others more room; (H) Stationary agent yielding to moving agent; (I) Move left; (J) Move right.

  (a) shows a grabbed image of the simulated shopping environment; Fig. 4(b) shows a plan view of the corresponding tracked participants. Each participant is represented by a disk of radius r i -the same value

  (a) Video grab.

Figure 4 .

 4 Figure 4. Instantaneous output of experiment.

Figure 5 .

 5 Figure 5. Trajectories corresponding to experiment shown in Fig. 4.

Distribution of behaviours.

				Eye-to-eye negotiation of immediate space
				Use of peripheral vision to asses environment
				Visual scanning of overall environment
				Verbal negotiation of immediate space
				Waiting for space / desired location to clear
				Step backwards to allow others more space
				Move forwards to allow others more space
				Staionary agent yielding to moving agent
				Move left
				Move right	Active Behaviours
					Reactive Behaviours
	Figure 3. Rule Agent 1 Agent 2	Behaviour
	1	Active	Passive	Agent 2 steps back and allows agent 1
	2	Active	Passive	Agent 1 moves and shares space
	3	Active	Passive	Agent 1 waits for empty space
	4	Active	Active	Agent 1 moves left or agent 2 moves right
	5	Active	Active	Agent 1 waits and gives way to agent 2

Table 2 . Identified behavioural rules.

 2

http://www.ics.forth.gr/cvrl

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


F Anstett-Collin

T Mara

J Goffart
email: jeanne.goffart@univ-savoie.fr

L Denis-Vidal

Nancy -France "Toward circular Economy" Sensitivity analysis for models with dynamic inputs: a case study to control the heat consumption of a real passive house

Keywords: Global Sensitivity Analysis, Karhunen-Loève Expansion, Dynamic Inputs, Building Energy Model

In this communication, we perform the sensitivity analysis of a building energy model. The aim is to assess the impact of the weather data on the performance of a model of a passive house, in order to better control it. The weather data are uncertain dynamic inputs to the model. To evaluate their impact, the problem of generating coherent weather data arises. To solve it, we carry out the Karhunen-Loève decomposition of the uncertain dynamic inputs. We then propose an approach for the sensitivity analysis of this kind of models. The originality for sensitivity analysis purpose is to separate the random variable of the dynamic inputs, propagated to the model response, from the deterministic spatio/temporal function. This analysis highlights the role of the solar gain on a high-insulated passive building, during winter time.

INTRODUCTION

Nowadays, one important feature when designing high-performance buildings is to reduce the energy consumption, for economical but also for environmental purposes. Thus, it appears a need of strategy to better control and optimize the heat consumption of houses. In this framework, softwares for building performance simulation are developed and are constantly updated due to the evolution of materials, techniques, ... This leads to complex, high-dimensional and multi-physics models, presenting uncertain inputs due to measurements, expert judgements. Thus, the problem of their reliability arises. To assess it, it requires some tools to better understand the influence of the parameters responsible for the heat consumption in order to control it. Uncertainty and sensitivity analyses (UASA) can help answering it and evaluating the impact of this lack of knowledge onto the model responses [START_REF] Kurowicka | Uncertainty analysis with high dimensional dependence modeling[END_REF][START_REF] De Rocquigny | Uncertainty in industrial practice, Uncertainty in industrial practice[END_REF][START_REF] Saltelli | Sensitivity Analysis[END_REF]. Uncertainty analysis aims at quantifying the overall uncertainty within a model and sensitivity analysis aims at determining the most influent input parameters onto the model response. Numerous studies have focused on the sensitivity analysis for static non-linear models, i.e. presenting static inputs, for example [START_REF] Saltelli | Sensitivity Analysis[END_REF][START_REF] Blatman | Sparse polynomial chaos expansion and adaptive stochastic finite element using a regression approach[END_REF][START_REF] Borgonovo | A new uncertainty importance measure[END_REF][START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The approaches may be local or global. Local approaches help to determine the impact of a small input variation around a nominal value [START_REF] Turyani | Local methods in sensitivity analysis[END_REF]. Global approaches also allow the determination of the same impact but by varying the input in its entire range of variation. For some applications, as energy building, it can be of great importance to consider the entire uncertainty range of inputs since they can vary within large intervals depending on their meaning. Considering only a single point in that interval as in the case of local study may be not enough informative and robust. Another advantage of global sensitivity analysis is that the sensitivity estimates of individual parameters are evaluated while all the other inputs are varied. In this way, the relative variability of each input is taken into account, thus revealing any existing interactions. Global methods are often based on the analysis of the output variance and are known as ANOVA (ANalysis Of VAriance) techniques [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF][START_REF] Saltelli | A quantitative model independent method for global sensitivity analysis of model output[END_REF]. The present study exclusively focuses on global approaches. Regarding building energy modeling, the literature focuses on sensitivity analysis for building models presenting only static inputs, for instance the thermophysical properties of the materials [START_REF] Tian | A review of sensitivity analysis methods in building energy analysis[END_REF]. However, some phenomena previously neglected for energy-consuming building may become preponderant in the consumption of low energy buildings, as the consideration of weather data, not only the temperature but also the solar radiation, the relative humidity, the wind speed,... But, it is expected that meteorological inputs play a crucial role in designing high-performance buildings. There are very few studies concerning the influence of the weather data which depend on time and are thus seen as dynamic uncertain inputs for the models. This can be explained by the fact that generating coherent weather data randomly is not obvious. One solution is to use experimental measure but it involves a large amount of data, too expensive to manage. Thus, an alternative is to create representative long-term weather data file [START_REF] Crawley | Which weather data should you use for energy simulations of commercial buildings ?[END_REF]. But, while generating static inputs samples is not an issue, it is not straightforward to generate samples that satisfy the desired random fields distribution. The aim of this study is then to assess the impact of the weather data on the heat consumption of a real passive house. This problem is challenging because it involves the one of the sensitivity analysis for models with uncertain dynamic inputs. Such an issue is rarely addressed (except, for instance, in [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF]). The uncertain static inputs can be seen as random variables defined by their marginal distribution and the dynamic ones as random fields, defined by their covariance function. The random variables can then be generated using a random sampling method, for instance the Latin Hypercube random sampling [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. A random field can be represented as a series expansion involving a complete set of deterministic functions with corresponding random coefficients [START_REF] Shinozuka | Simulation of the stochastic process by spectral representation[END_REF][START_REF] Grigoriu | On the spectral representation method in simulation[END_REF][START_REF] Zhang | Orthogonal series expansions of random processes in reliability analysis[END_REF], as the truncated Karhunen-Loève (KL) expansion [START_REF] Huang | Convergence study of the truncated Karhunen-Loeve expansion for simulation of stochastic processes[END_REF][START_REF] Phoon | Simulation of second order processes using Karhunen-Loeve expansion[END_REF]. KL series expansion is based on the eigen-decomposition of the covariance function, involving orthogonal deterministic basis functions and the orthogonal uncorrelated random coefficients. This allows the optimal encapsulation of the information contained in the random process into a set of discrete uncorrelated random variables. The originality of the proposed approach is to separate the random variable of the dynamic inputs, propagated to the model response, from the deterministic spatio/temporal function, using Karhunen-Loève decomposition of the dynamic inputs. Then, for the SA, the influence of the dynamic inputs onto the model response is given by the one of the random coefficients of the Karhunen-Loève decomposition. For the estimation of their sensitivity indices, the Sobol' method is applied [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The proposed approach is then applied to a building energy model to quantify the impact of each weather data input on the performance of a real passive house. The paper is organized as follows. In Section 2, the energy building model which is studied is described. Then, Section 3 is focused on the generation of coherent weather data. In Section 4, sensitivity analysis for static inputs is recalled and then is extended to dynamic inputs. Finally, Section 5 presents the results obtained.

ENERGY BUILDING MODEL

The building energy model studied in this work represents a real building located on an experimental platform named INCAS, introduced by the National Institute of Solar Energy and located in Le Bourget Du Lac, in France. This place presents a temperate continental climate with alpine influence. The building is a single-detached house with low energy consumption (figure 1). The house is double-glazed. The thickness of the wall is 50 cm, made of 15 cm of perpend, 20 cm of insulating material and again 15 cm of perpend. The house contains hundreds of sensors to quantify its thermal behavior. Details about the house can be found in [START_REF] Spitz | Pratical application of uncertainty and sensitivity analysis on an experimental house[END_REF]. The building is divided The output of interest is the heating consumption at the ground and at the first floors. The uncertain dynamic inputs are the outdoor air temperature, the direct and the diffuse solar radiations, the speed and the direction of wind and the relative humidity. A previous study [START_REF] Goffart | Sampling-based method for managing weather inputs in global sensitivity analysis of building energy models, Energy and Buildings[END_REF], not presented here, has focused on the statistical analysis of the weather data. One month of january that is representative of typical winter from 20 years of observations has been used, since the heat consumption is the highest and the most costly, during this period of time, at this place. Moreover, considering only one season allows to respect stationarity of data. The statistical characteristics of a given input are the covariance function and the hourly mean, on one day. The next section describes how to generate coherent weather data randomly.

WEATHER DATA GENERATION

The weather data are considered as the uncertain dynamic inputs of the model. Let denote them as

ω d (θ , x) = {ω d 1 (θ , x), . . . , ω d N d (θ , x)}.
Here N d = 6. The stochastic variable θ is used to indicate the randomness of the input ω d . The variable x represents the time dependence of ω d but for some other applications, it can also represent a space dependence. The random fields ω d i are assumed independent and normally distributed with mean ωi d (x) and a covariance function

C i (x 1 , x 2 ), i = 1, . . . , N d , N d denoting the number of dynamic inputs. The covariance function C i (x 1 , x 2 ) is symmetric, positive definite.
According to [START_REF] Phoon | Simulation of second order processes using Karhunen-Loeve expansion[END_REF], the random fields ω d i , having a mean ωi d (x) and a finite variance σ 2 i , can be approximated using the truncated KL series:

ω d i (x, θ ) ≃ ωi d (x) + M i ∑ ki=1 √ λ ki ξ ki (θ ) f ki (x) (1)
where λ ki and f ki are the deterministic eigenvalues and eigenfunctions of the covariance function C i (x 1 , x 2 ), ξ i is a set of independent standard normal variables and M i is the number of KL-terms. In practice, we retain the first M i eigenmodes that contain the 95% of the variance of the input ω d i . The number of eigenmodes retained depends on the choice of the covariance function and may be very different from one input to another. The key feature for simulating random fields using KL expansion lies on the ability to determine accurately the eigenvalues and eigenfunctions of the covariance function. They are given from the spectral decomposition of the covariance function C i (x 1 , x 2 ), requiring to solve the homogeneous Fredholm integral equation of the second kind given by:

D C i (x 1 , x 2 ) f ki (x 1 )dx 1 = λ ki f ki (x 2 ) (2)
For some covariance functions (first-order Markov or Wiener-Levy processes, for example), the equation ( 2) can be twice differentiable with respect to x 1 . The resulting differential equation can be solved analytically and eigenvalues can be obtained as well. But in most cases, solving equation (2) requires numerical methods, such as the Galerkin one [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF]. The Galerkin approach consists of finding a functional basis for the solution space of the equation, then projecting the solution on the functional basis, and minimizing the residual with respect to the functional basis. Very often, the bases are polynomial or trigonometric. Representation of integral operators using polynomial or trigonometric bases requires approximate integration quadratures. This can be very costly to compute. To avoid tedious quadratures and alleviating computational effort, an alternative is to use a wavelet-Galerkin approach [START_REF] Phoon | Simulation of second order processes using Karhunen-Loeve expansion[END_REF][START_REF] Maleknejada | Using wavelet for numerical solution of Fredholm integral equations[END_REF]. In this case, the representation of integral operators is made in wavelet basis and can be performed without numerical integration, as detailed in the Appendix. The comparison of wavelet-Galerkin method with other available methods in solving the Fredholm integral equation can be found in [START_REF] Huang | Convergence study of the truncated Karhunen-Loeve expansion for simulation of stochastic processes[END_REF].

Once the weather data are generated, the sensitivity analysis can be performed as explained in the next section.

SENSITIVITY ANALYSIS

First, let recall the way of performing sensitivity analysis of a model presenting static inputs. To do so, consider the following model:

y(θ ) = g(ω s (θ )) (3)
where y ∈ R m is the model response of interest, ω s = {ω s 1 , . . . , ω s N s } is a set of N s random variables (static inputs) and g a nonlinear function. In the following, we assume that random variables ω s i are independent and defined by their marginal distribution. They can be generated using a random sampling, Latin Hypercube sampling scheme, for instance, for its ease of implementation, [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. The effect of the uncertain inputs ω s i onto y j (y j ∈ R) can be estimated with sampling-based methods such as Sobol' method [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The first-order sensitivity index of the input ω s i to the model output y j is given by:

S i = V (E(y j |ω s i )) V (y j ) (4)
where V (E(y j |ω s i )) is the variance of the conditional expectation of y j when ω s i is set and V (y j ) is the total variance of the output y j . The first order sensitivity index S i represents the main effect of the input ω s i which corresponds to its contribution alone. The value of S i lies between 0 and 1. The closer to 1 its value is, the more input ω s i contributes to the total variance of the output. Sensitivity indices of higher order can also be computed to assess the influence of the input interactions but it is not presented here, see [START_REF] Saltelli | Sensitivity Analysis[END_REF] for more details. Some models may be complex with a high number of inputs so that analytical calculations of the previous sensitivity index become time consuming or even impossible. It is therefore necessary to estimate them. Some approaches have been proposed in the literature, based on Monte Carlo simulations or spectral expansion of the output [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF][START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF]. The estimator considered in this study is provided by [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF] and the procedure to compute it from few samples to reduce the computational cost is recalled below.

With the method originally proposed in [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], N s + 1 samples are necessary to compute the main effects of the N s inputs. To reduce the simulation cost, a method proposed in [START_REF] Mara | Application of global sensitivity analysis of model output to building thermal simulations[END_REF] has been used, requiring only two samples.

Consider a first sample of inputs with K points, denoted MOSIM14 -November 5-7-2014 -Nancy -France

ω s1 : ω s1 =     ω s1 11 • • • ω s1 1K ω s1 21 • • • ω s1 2K • • • • • • • • • ω s1 N s 1 • • • ω s1 N s K     (5)
The model is evaluated with this first input samples, leading to the output y (1) of dimension m × K, m components of the output, each with K sample points:

y (1) =      y (1) 11 • • • y (1) 1K y (1) 21 • • • y (1) 2K • • • • • • • • • y (1) m1 • • • y (1) mK      (6)
Associate a random permutation function RP i to each input ω s i . A second set, denoted ω s2 , is built from ω s1 by permuting its values.

ω s2 =     RP 1 (ω s1 11 ) • • • RP 1 (ω s1 1K ) RP 2 (ω s1 21 ) • • • RP 2 (ω s1 2K ) • • • • • • • • • RP N s (ω s1 N s 1 ) • • • RP N s (ω s1 N s K )     (7)
Then, the model is evaluated with this second set of input samples, to obtain y (2) :

y (2) =      y (2) 11 • • • y (2) 1K y (2) 21 • • • y (2) 2K • • • • • • • • • y (2) m1 • • • y (2) mK      (8)
where y

(2)

j = RP i (y (1)
j ) and RP i is the random permutation function associated to the input ω s i . To compute the first order sensitivity index (4) corresponding to the effect of ω s i , the values of y (1) are rearranged according to the corresponding permutation RP i . In this way, it is as the output y has been obtained when all the input factors are varying except the one of interest ω s i . Let define y j 2 = y j .y j , where y j is a K-dimensional vector and the symbol (.) denotes the scalar product of two vectors. The following result gives a way for computing the first-order sensitivity index to the output y j .

Ŝi = (y

(1)

j -f (1) u).(RP -1 i (y (2) j ) -f (2) u) y (1) j -f (1) u RP -1 i (y (2) j ) -f (2) u (9)
with

f (1) = 1 K K ∑ k=1 y (1)
jk and f

(2) = 1 K K ∑ k=1 RP -1 i (y (2)
jk ) [START_REF] Goffart | Sampling-based method for managing weather inputs in global sensitivity analysis of building energy models, Energy and Buildings[END_REF] and u the vector of K components equal to 1.

In the case of model presenting dynamic inputs, the problem of determining this sensitivity index is much more complex. In this case, the key point is to generate coherent samples of the dynamic input satisfying the desired random fields distribution. As explained in Section 3, the random fields can be decomposed using the KL expansion [START_REF] Blatman | Sparse polynomial chaos expansion and adaptive stochastic finite element using a regression approach[END_REF]. The influence of the dynamic inputs can then be analyzed through the one of the random coefficients of the KL expansion. The sensitivity analysis is therefore performed as for static inputs, as explained below.

Now consider the following model:

y(x, θ ) = g(ω d (x, θ ), ω s (θ ), x) (11)
where y ∈ R m is the model response of interest, x ∈ D is the spatial/time variable,

ω d = {ω d 1 , . . . , ω d N d
} is a set of N d random fields (dynamic inputs) and ω s = {ω s 1 , . . . , ω s N s } is a set of N s random variables (static inputs). As previously, we assume that random variables ω s i are independent and defined by their marginal distribution. The first step is to generate samples for the dynamic inputs using the KL decomposition of ω d , as detailed in Section 3. Once the eigenmodes are obtained for all the dynamic inputs ω d i (x, θ ) using ( 1), the influence of ω d i (x, θ ) is then given by the one of the random coefficients ξ ki propagated to y(x, θ ). Consider the M i -dimensional random coefficient ξ i , grouping the M i modes of the input ω d i :

ξ i = {ξ 1i , . . . , ξ M i i } (12)
Thus, SA of the model output y(x, θ ) is performed through the random vectors {ξ 1 , . . . , ξ N d , ω s }, with ξ i given by [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. Consequently, the effect of the group of factors ξ i is the one of the dynamic input ω d i and so on. Thus, considering the static and the dynamic inputs of the model, there is

N s + N d ∑ i=1 M i inputs in all to analyze.
The sensitivity indices of the group of factors ξ i = {ξ 1i , . . . , ξ M i i } can be computed as:

S i = V (E(y j |ξ 1i , . . . , ξ M i i )) V (y j ) (13)
where V (E(y j |ξ 1i , . . . , ξ M i i )) is the variance of the conditional expectation of y j when ξ 1i , ξ 2i , . . . and ξ M i i are set. This sensitivity index can be estimated as explained previously.

The sensitivity of ω s can then be computed using ( 4) and the estimator [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF], as presented previously.

The proposed approach is applied to the energy building model presented in Section 2.

RESULTS

Consider the energy building model of Section 2. The model output y = {y 1 , y 2 } is the ideal heating consump-tion at the ground floor (y 1 ) and at the first floor (y 2 ). The daily consumptions are summed over one month, in such a way that the output y does no longer depend on time. However, the inputs ω d are still depending on time. They represent the outdoor air temperature (ω d 1 ), the direct and the diffuse solar radiations (ω d 2 and ω d 3 ), the speed and the direction of wind (ω d 4 and ω d 5 ) and the humidity (ω d 6 ). There is no static inputs for this study. From the hourly mean ωd i and the covariance function C i (x 1 , x 2 ), the sample generation of the dynamic inputs ω d can be carried out, as explained in Section 3, using the Fast Haar wavelet algorithm (see the appendix). To do so, Matlab software has been used. This leads here to 512 modes for each of the 6 inputs, that is M i = 512, ∀i = {1, . . . , 6}. Then, the model is simulated with the generated inputs, using the dedicated software Energy-Plus.

It can be noticed that the heat consumption is higher at the first floor than at the ground floor. The two thermal zones under study have the same volume but there is more glass surface at the ground floor. Thus, it may be assumed that the solar gain is more important at this floor, during winter. The sensitivity indices are computed as explained in Section 4. The 512 ξ i j of a given input are grouped together. Thus, there is 6 sensitivity indices S i to compute, according to [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF]. To do so, the estimator (9) has been used. Moreover, take the advantage of the random permutation trick explained in Section 4, which requires only two samples, here of K = 1000 points, allowing to reduce the computational cost. For this computation, Matlab software has been used.

Figure 2 shows the dispersion of the consumption. For the ground floor, the outdoor temperature and the solar radiation have almost the same sensitivity index (S 1 = 0.46 and S 2 = 0.52). This means that they influence the consumption almost in the same way, but with a slightly prevalence of the solar radiation. Besides, for the first floor, the outdoor temperature has a sensitivity index of S 1 = 0.76, explaining almost 80% of the consumption, compared to low influence of the solar radiation (S 2 = 0.23). On the other hand, the sensitivity indices of the humidity, diffuse radiation, velocity and direction of wind can be neglected showing that they are not influential. Note that small sensitivity indices can be negative because of the numerical approximation of the estimator considered. It is worth noting that these inputs are not influential in the studied model. This conclusion can be different for other models, in the way these variables are taken into account for the calculation of the energy balance. For example, humidity is not relevant here because the ideal heating device considered does not account for the impact of the humidity onto the consumption.

It can be noticed that the sum of the six sensitivity indices is almost one, showing no interaction between the inputs. In conclusion, there is two significant weather inputs, the solar radiation and the outdoor temperature with different influences according to the given floor. These results highlight the capacity of the INCAS house to exploit the solar gains during winter, especially at the ground floor. Indeed, the solar radiation has a greater impact on a high-insulated passive building than on an energy-consuming one. This is relevant to minimize the heat consumption during winter in this low energy house.

CONCLUSION

In this paper, the impact of the weather data on the energy consumption of a passive house has been quantified. This has raised the problem of the generation of coherent weather data. The proposed solution is based on the KL decomposition of the dynamic inputs. An approach for sensitivity analysis of models presenting dynamic inputs has been proposed. The originality is to separate the random coefficients of the input propagated to the output from the spatio/temporal function.

The study has shown that the outdoor temperature and the solar radiation are of prime importance for the INCAS house to reduce the heat consumption. It underlines the capacity of the INCAS house to exploit the solar gains during winter and to confirm its passive strategy to earn energy. This analysis helps to better understand the weather data influence for the heat consumption of a real passive house. It also helps to optimize the capacity of the building to exploit free energy resources. It is worth noting that the results depend on the building of interest (single-detached terraced house, private or industrial building, ...), its orientation, the climate where it is located. SA can be seen here as a methodology to characterize the ability of low efficiency building to exploit free resources of the environment. The next step is to study the influence of both static and dynamic inputs for the heat consumption, the static inputs representing the thermophysical properties of the materials.

Appendix

Using the wavelet-Galerkin approach, the covariance function C i (x 1 , x 2 ) can be expanded as, in the wavelet bases:

C i (x 1 , x 2 ) = N ∑ k=0 N ∑ j=0 A k j ψ k (x 1 )ψ j (x 2 ) = Ψ T (x 1 )AΨ(x 2 ) ( 14
)
where A is an N × N matrix representing the 2D wavelet transform of C i (x 1 , x 2 ). The components A k j of the matrix A are then given by: [START_REF] Kurowicka | Uncertainty analysis with high dimensional dependence modeling[END_REF] where the functions ψ k , k = 1, . . . , N, are Daubechies' wavelets and h k defined by:

A k j = 1 h k h j 1 0 1 0 C i (x 1 , x 2 )ψ k (x 2 )ψ j (x 1 )dx 1 dx 2
h k = 1 0 ψ 2 k (x)dx (16)
We consider the Haar wavelets, the simplest family of Daubechies' wavelets. The Haar mother wavelet is defined by:

ψ(x) =    1 if x ∈ [0; 0.5[ -1 if x ∈ [0.5; 1[ 0 otherwise
A family of orthogonal Haar wavelets over the domain [0; 1] can be generated by shifting and scaling the mother wavelet:

ψ 0 (x) = 1 ψ k (x) = 2 a/2 ψ(2 a x -b), a, b ∈ Z k = 2 a + b, b = 0, 1, . . . , 2 a -1, a = 0, 1, . . . , m -1 (18)
where a and b are the dilation and translational integer constants, m the maximum wavelet level, related to N by N = 2 m . The double integral of (15) corresponds to the 2D wavelet transform of C i (x 1 , x 2 ) and does not require numerical integration. To carry out this 2D wavelet transform, the 1D wavelet transform is applied first on the rows and then on the columns of the matrix containing the values of C i (x 1 , x 2 ) sampled over a N by N grid. The different steps to perform the 1D Haar wavelet transform are recalled below. For further details, see [START_REF] Nievergelt | Wavelets made easy[END_REF][START_REF] Phoon | Simulation of second order processes using Karhunen-Loeve expansion[END_REF].

1. Consider a function F(x) (here the rows of C i (x 1 , x 2 )).

Sample F(x) at N discrete points x i = 2i + 1 2N , i = 0, . . . , N -1. The samples are denoted F i and N = 2 m with m the maximum level of the wavelet, set a priori.

2. Initialize a N × 1 vector denoted a m,k with the samples of F as follows:

a m,k = F k , k = 0, . . . , N -1 (19)
3. The vector is processed using an inverse binary tree, where the topmost layer (mth) contains N nodes with values given by [START_REF] Nievergelt | Wavelets made easy[END_REF]. The nodal values in subsequent layers are computed as:

a j,k = 1 2 (a j+1,2k + a j+1,2k+1 ) (20)
with k = 0, 1, . . . , 2 j -1 and j = m -1, . . . , 1, 0.

4. The wavelet coefficients are evaluated from the nodal values in this binary tree:

c j,k = 1 2 (a j+1,2k -a j+1,2k+1 ) (21)
5. Finally, the wavelet coefficients are given by [a 0,0 , c 1 , . . . , c N-1 ], with:

c i = c j,k , i = 2 j + k (22)
where k = 0, 1, . . . , 2 j -1 and j = 0, 1, . . . , m -1.

6. The function F(x) can then be approximated as:

F(x) = a 0,0 ψ 0 (x) + N-1 ∑ i=1 c i ψ i (x) (23)
with ψ i (x) the Haar wavelets.

To solve (2), the eigenfunctions f ki are approximated by Haar wavelet series:

f ki (x) = N-1 ∑ j=0 d (ki) j ψ j (x) = Ψ T (x)D (ki) (24)
where D (ki) a N × 1 matrix whose components are the wavelet coefficients d (ki) j and Ψ T (x) is the vector of components ψ i (x). Substituting equations ( 14) and (24) into equation [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF] gives the following relation: λ ki Ψ T (x)D (ki) = Ψ T (x)AHD (ki) (25

)
where H is diagonal matrix of constants h k due to the orthogonality condition. This leads to a finite dimensional eigenvalue problem of the form λ ki D (ki) = AHD (ki) , which can be solved easily using an eigenvalue solver.

Figure 1 :

 1 Figure 1: INCAS house

Figure 2 :Figure 3

 23 Figure 2: Dispersion of the heat consumption for the ground floor (1) and for the first floor (2)

Figure 3 :

 3 Figure 3: Sensitivity indices for the 6 inputs, at the ground floor (1) and at the first floor (2)