Alexandre David
email: adavid@cs.aau.dk

Kim G Larsen

Marius Mikučionis
email: marius@cs.aau.dk

Danny Bøgsted Poulsen
email: dannybp@cs.aau.dk

Axel Legay
email: axel.legay@inria.fr

Sean Sedwards
email: sean.sedwards@gmail.com

Dehui Du
email: dehuidu@cs.aau.dk

Statistical Model Checking for Stochastic Hybrid Systems *

This paper presents novel extensions and applications of the UPPAAL-SMC model checker. The extensions allow for statistical model checking of stochastic hybrid systems. We show how our race-based stochastic semantics extends to networks of hybrid systems, and indicate the integration technique applied for implementing this semantics in the UPPAAL-SMC simulation engine. We report on two applications of the resulting tool-set coming from systems biology and energy aware buildings.

Introduction

Statistical Model Checking (SMC) [START_REF] Legay | Statistical Model Checking: An Overview[END_REF][START_REF] Laplante | Probabilistic abstraction for model checking: An approach based on property testing[END_REF][START_REF] Sen | Statistical Model Checking of Black-Box Probabilistic Systems[END_REF][START_REF] Håkan | Probabilistic Verification of Discrete Event Systems Using Acceptance Sampling[END_REF][START_REF] Katoen | The ins and outs of the probabilistic model checker MRMC[END_REF] is an approach that has recently been proposed as new validation technique for large-scale, complex systems. The core idea of SMC is to conduct some simulations of the system, monitor them, and then use statistical methods (including sequential hypothesis testing or Monte Carlo simulation) in order to decide with some degree of confidence whether the system satisfies the property or not. By nature, SMC is a compromise between testing and classical formal method techniques. Simulation-based methods are known to be far less memory and time intensive than exhaustive ones, and are some times the only option. SMC has been implemented in a series of tools that have defeated well-known tools such as PRISM on several case studies. Unlike more "academic" exhaustive (and intractable) techniques, SMC is spreading to various research areas such as systems biology [START_REF] Gong | Computational Modeling and Verification of Signaling Pathways in Cancer[END_REF][START_REF] Kumar Jha | A Bayesian Approach to Model Checking Biological Systems[END_REF] and software engineering [START_REF] Zuliani | Bayesian statistical model checking with application to Simulink/Stateflow verification[END_REF][START_REF] Martins | Statistical Model Checking for Distributed Probabilistic-Control Hybrid Automata with Smart Grid Applications[END_REF], in particular for industrial applications [START_REF] Basu | Statistical Abstraction and Model-Checking of Large Heterogeneous Systems[END_REF][START_REF] Basu | Verification of an AFDX Infrastructure Using Simulations and Probabilities[END_REF][START_REF] Basu | Statistical abstraction and model-checking of large heterogeneous systems[END_REF].

There are several reasons for this success. Firstly, SMC is simple to understand, use and (in principle) to implement. Secondly, no additional modelling or specification effort is needed, provided that the mode-ling formalism used can be given a natural stochastic semantics serving as the basis for interpretation of the specification formalism and as a basis for generating simulation-runs. Thirdly, SMC allows to analyse properties [START_REF] Clarke | On simulation-based probabilistic model checking of mixed-analog circuits[END_REF][START_REF] Basu | Statistical Abstraction and Model-Checking of Large Heterogeneous Systems[END_REF] that cannot be expressed in classical temporal logic, including properties for which classical model checking is undecidable.

In a series of recent works [START_REF] David | Time for statistical model checking of real-time systems[END_REF][START_REF] David | Statistical Model Checking for Networks of Priced Timed Automata[END_REF], we have investigated the problem of Statistical Model Checking for networks of Priced Timed Automata (PTAs). PTAs are timed automata, whose clocks can evolve with different rates, while 1 being used with no restrictions in guards and invariants. In [START_REF] David | Statistical Model Checking for Networks of Priced Timed Automata[END_REF], we have proposed a natural stochastic semantics for such automata, which allows to perform statistical model checking. Our work has later been implemented in UPPAAL-SMC, that is a stochastic and statistical model checking extension of UPPAAL. UPPAAL-SMC relies on a series of extensions of the statistical model checking approach generalised to handle real-time systems and estimate undecidable problems. UPPAAL-SMC comes together with a friendly user interface that allows a user to specify complex problems in an efficient manner as well as to get feedback in the form of probability distributions and compare probabilities to analyse performance aspects of systems. The UPPAAL-SMC model checking has been applied to a wide range of examples from networking and Nash equilibrium [START_REF] Bulychev | Computing Nash Equilibrium in Wireless Ad Hoc Networks: A Simulation-Based Approach[END_REF] through systems biology [START_REF] David | Runtime Verification of Biological Systems[END_REF], real-time scheduling [START_REF] David | Schedulability of Herschel-Planck Revisited Using Statistical Model Checking[END_REF], and energy aware systems [START_REF] David | An Evaluation Framework for Energy Aware Buildings using Statistical Model Checking Submitted[END_REF].

For PTAs clocks evolve with fixed rates depending only on the discrete state of the model, i.e. locations and discrete variables. In this paper, we present and implement an extension of the modelling formalism of UPPAAL-SMC, where clock rates may depend not only on values of discrete variables but also on the value of other clocks, effectively amounting to ordinary differential equations (ODEs). As a first contribution of this paper, we present an extension of our race-based stochastic semantics to networks of stochastic hybrid automata (SHA). One major difficulty is the implementation of this extensions in our engine for generating random runs, in a manner which is correct with respect to the stochastic semantics. In fact UPPAAL-SMC does not solve the equations exactly but currently supports the Euler integration method. A fixed time step (defined by the user) is used by an internal integrator component added to the system. This integrator races with the other processes during its time step and all rates are considered constant and defined by the equations in the model. As a second contribution of the paper we describe this integration in UPPAAL-SMC. It is worth mentioning that while other SMC-based approaches exist for SHAs [START_REF] Zuliani | Bayesian statistical model checking with application to Simulink/Stateflow verification[END_REF][START_REF] Zuliani | Rare-event verification for stochastic hybrid systems[END_REF][START_REF] Ballarini | COSMOS: A Statistical Model Checker for the Hybrid Automata Stochastic Logic[END_REF], none of them do consider ODE-based modelisation for clock rates.

To demonstrate the applicability of this new extension, we apply advanced SMC techniques to two challenging applications coming from systems biology and energy aware buildings. For the case of systems biology, we show how the combination of ODEs and SMC allows us to reason on biological oscillations -a problem that is beyond the scope of most existing formal verification techniques. We model a genetic circadian oscillator, which is used to distil the essence of several real circadian oscillators. For the case of energy aware buildings, we refer to a recently developed framework including components for layout of buildings, availability of heaters, climate and user behaviours allowing to evaluate different strategies for distributing heaters among rooms in terms of the resulting comfort and energy consumption. To indicate central parts of this framework and the clear advantages of modelling the evaluation of room temperatures with ODEs, we illustrate in this paper the framework with a small instance comprising two rooms with a single shared heater.

Structure of the Paper. The remainder of the paper is structured as follows. In the next Section 2 we preview the expressive power of the hybrid extensions of UPPAAL-SMC using an extension of the well-known bouncing ball. Sections 3 and 4 details the semantics of the extended formalism of networks of stochastic hybrid automata. Applications to Energy Aware Buildings and a Biological Oscillator are given in Section 5 and Section 6, respectively. Finally, Section 7 concludes the paper and suggests directions for future research.

Throwing, Bouncing and Hitting Ball

To give an early illustration of the expressive power of the extended modelling formalism of UPPAAL-SMC, we consider a variant of the well-known bouncing ball. In our version 2 the ball is initially thrown against a wall, bounces against it, and then continues its trajectory by falling and bouncing against the floor. In addition, a (inexperienced) player tries to hit the ball randomly according to an exponential distribution. The model is depicted in Fig. 1(a)-(c). The player is modelled as a simple automaton that broadcasts hit! with an exponential distribution of rate 5/2. The x coordinate (Fig. 1a) is initialised to 10 with an uncertain derivative vx uniformly between [-10, -9.5] (the ball is thrown against the wall), after which the ball moves toward the wall (placed at 0). Here, the automaton outputs bounce! on an urgent channel, which forces the transition to take place deterministically at x= 0. After a bounce with a random dampening factor of the velocity vx uniformly between [0.85, 1], the ball continues to move in the opposite direction. The y coordinate (Fig. 1b) is initialised to 10 with an uncertain derivative vy. The model shows the effect of gravitation with vy'= -9.81. The ball bounces with a random dampening factor on the floor (at 0) and when the ball is away from the wall (x≥ 10) then it can be hit by the player provided it is high enough (y≥ 6). Depending on the current direction of the ball, the ball may bounce or it is pushed. One possible trajectory of the ball is shown in Figure 1d. The plot is obtained by checking the query "simulate 1 [x<=40]{y}". The vertical line shows the ball moving to its initial position and should be ignored. The ball bounces as expected against the wall, the floor, and the hitting of the player. UPPAAL-SMC is able to simulate this hybrid system that has a second order ODE, a stochastic controller (the player), and a stochastic environment (random dampening factor).

In addition, we may perform statistical model-checking in order to estimate the probability that the ball is still bouncing above a height of 4 after 12 time units with the query: Pr[<=20](<> time>=12 and y>=4) which returns the confidence interval [0.44, 0.55] with 95% confidence after having generated 738 runs. We can also test for the hypothesis Pr[<=20](<> time>=12 and y>=4) >= 0.45, which gives a more precise lower bound. The hypothesis holds with a region of indifference ±0.01 and a level of significance of 5% after generating 970 runs.

Networks of Hybrid Automata

UPPAAL-SMC now supports the analysis of stochastic hybrid automata (SHA) that are timed automata whose clock rates can be changed to be constants or expressions depending on other clocks, effectively defining ODEs. This generalizes the model used in our previous work [START_REF] David | Time for statistical model checking of real-time systems[END_REF][START_REF] David | Statistical Model Checking for Networks of Priced Timed Automata[END_REF] where only linear priced automata were handled. Our new release UPPAAL-SMC 4.1.103 supports fully hybrid automata with ODEs and a few built-in complex functions (such as sin, cos, log, exp and sqrt)! Hybrid Automata Intuitively, a hybrid automaton H is a finite-state automaton extended with continuous variables that evolve according to dynamics characterizing each discrete state (called a location). Let X be a finite set of continuous variables. A variable valuation over X is a mapping ν : X → R, where R is the set of reals. We write R X for the set of valuations over X . Valuations over X evolve over time according to delay functions F : R ≥0 × R X → R X , where for a delay d and valuation ν, F(d, ν) provides the new valuation after a delay of d. As is the case for delays in timed automata, delay functions are assumed to be time additive in the sense that

F(d 1 , F(d 2 , ν)) = F(d 1 + d 2 , ν).
To allow for communication between different hybrid automata we assume a set of actions Σ, which is partitioned into disjoint sets of input and output actions, i.e.

Σ = Σ i ⊎ Σ o . Definition 1. A Hybrid Automaton (HA) H is a tuple H = (L, ℓ 0 , X , E, F, I), where: (i) L is a finite set of locations, (ii) ℓ 0 ∈ L is an initial location, (iii) X is a finite set of continuous variables, (iv) Σ = Σ i ⊎ Σ o is

a finite set of actions partitioned into inputs (Σ i ) and outputs (Σ o ), (v) E is a finite set of edges of

the form (ℓ, g, a, ϕ, ℓ ′ ), where ℓ and ℓ ′ are locations, g is a predicate on R X , action label a ∈ Σ and ϕ is a binary relation on R X , (vi) for each location ℓ ∈ L F(ℓ) is a delay function, and (vii) I assigns an invariant predicate I(ℓ) to any location ℓ.

The semantics of a HA H is a timed labeled transition system, whose states are pairs (ℓ, ν) ∈ L × R X with ν |= I(ℓ), and whose transitions are either delay transitions (ℓ, ν)

d -→ (ℓ, ν ′ ) with d ∈ R ≥0 and ν ′ = F(d, ν), or discrete transitions (ℓ, ν) a -→ (ℓ ′ , ν ′ ) if
there is an edge (ℓ, g, a, ϕ, ℓ ′ ) such that ν |= g and ϕ(ν, ν ′ ). We write (ℓ, ν) (ℓ ′ , ν ′ ) if there is a finite sequence of delay and discrete transitions from (ℓ, ν) to (ℓ ′ , ν ′ ).

In the above definition, we have deliberately left open the concrete syntax for the delay function F as well as the invariant I. In UPPAAL-SMC the delay update for a simple clock x -used in (priced) timed automata -is given by an implicit rate x ′ = 1 or an explicit rate x ′ = e appearing in the invariant of ℓ, where e is an expression only depending on the discrete part of the current state. More generally, the effect of the delay function F may be specified by a set of ODEs needing to be solved. It is important to note that in specifying the delay function F and the invariant I, the full syntax of UPPAAL expressionsincluding user-defined functions -is at the disposal. Example 1. Reconsider the extended bouncing ball example from Section 2. Here the automaton for the x coordinate may be initialized to the state (x = 10, vx = -9.8), after which the following transition sequence may occur:

(x = 10, vx = -9.8) 10÷9.8 -→ (x = 0, vx = -9.8) bounce! -→ (x = 0, vx = 9.31)
where in the bounce!-transition the dampening factor has non-deterministically been chosen from the interval [0.85, 1.00] as 0.95. The automaton for the y coordinate may be initialized to the state (y = 10, vy = 5.9) after which the delay function will effectively be given by delay transitions:

(y = 10, vy = 5.9) d -→ (y = -9.81/2d 2 + 5.9d + 10, vy = -9.81d + 5.9)

The ball bounce!'es on the floor, i.e. y = 0 or when (approximately) d = 2.15, and will afterwards -having non-deterministically chosen the dampening factor 0.9 -continue from the state (y = 0, vy = 13.67).

Networks of Hybrid Automata

Following the compositional specification theory for timed systems in [START_REF] David | Timed I/O automata: a complete specification theory for real-time systems[END_REF], we shall assume that NHAs are input-enabled in the sense, that for all states (ℓ, ν) and input actions ι ∈ Σ i , for all HAs j, there is an edge (ℓ j , g, ι, ϕ, ℓ j ′ ) such that ν |= g and ϕ(ν, ν ′ ) for some valuation ν ′ . Also, we shall assume that time always diverges, and that different automata synchronize on matching inputs and outputs as a standard broadcast synchronization [START_REF] Gómez | A Compositional Translation of Timed Automata with Deadlines to UPPAAL Timed Automata[END_REF].

Whenever A j = (L j , X j , Σ j , E j , F j , I j ) ( j = 1 . . . n) are NHA, they are composable into a closed network iff their variable-sets are disjoint (X j ∩ X k = / 0 when j = k), they have the same action set (Σ = Σ j = Σ k for all j, k), and their output action-sets provide a partition of Σ (Σ

j o ∩ Σ k o = / 0 for j = k, and Σ = ∪ j Σ j o ).
For a ∈ Σ we denote by c(a) the unique j with a ∈ Σ j . If ν ∈ R X with X = ∪ j X j , we denote by ν ↓ X j the projection of ν to X j .

Definition 2. Let

A j = (L j , X j , Σ, E j , F j , I j ) (with j = 1 . . . n) be composable NHAs. Their composition (A 1 | . . . | A n ) is the HA A = (L, X , Σ, E, F, I) where (i) L = × j L j , (ii) X = ∪ j X j , (iii) F(ℓ)(d, ν)(x) = F j (ℓ j )(d, ν ↓ X j )(x)
when x ∈ X j , (iv) I(ℓ) = ∩ j I(ℓ j ), and (v) (ℓ, ∩ j g j , a, ∪ j ϕ j , ℓ ′ ) ∈ E whenever (ℓ j , g j , a, ϕ j , ℓ ′ j ) ∈ E j for j = 1 . . . n.

Stochastic Semantics for Networks of Hybrid Automata

Reconsidering again our extended version of the bouncing ball from Section 2, it is clear that there is a constant race between the ball bounce!ing on the floor and the player hit!ing the ball. Whereas the time of bouncing is deterministic -given by the ODE obtained from the (stochastic) effect of the previous bounce! or hit! -the time of hitting is stochastic according to an exponential distribution with rate 5/2. However based on this, a measure on sets of runs of the systems is induced, according to which quantitative properties such as "the probability that the ball with have a height greater than 4 after 12 time-units" become well-defined.

Our early works [START_REF] David | Statistical Model Checking for Networks of Priced Timed Automata[END_REF] -though aimed at stochastic semantics of priced timed automata -is sufficiently general that it also provides the basis for a natural stochastic semantics for networks of HAs, where components associate probability distributions to both the time-delays spent in a given state as well as to the transition between states.

Let A j = (L j , X j , Σ, E j , F j , I j ) ( j = 1 . . . n) be a collection of composable HAs. Under the assumption of input-enabledness, disjointedness of clock sets and output actions, states of the composite NHA A = (A 1 | . . . | A n ) may be seen as tuples s = (s 1 , . . . , s n ) where s j is a state of A j , i.e. of the form (ℓ, ν) where ℓ ∈ L j and ν ∈ R X j . Our probabilistic semantics is based on the principle of independence between components. Repeatedly each component decides on its own -based on a given delay density function and output probability function -how much to delay before outputting and what output to broadcast at that moment. Obviously, in such a race between components the outcome will be determined by the component that has chosen to output after the minimum delay: the output is broadcast and all other components may consequently change state.

Stochastic Semantics of HA Components

The stochastic semantics of HAs refine the non-deterministic choices that may exist with respect to delay, output and next state. We consider the component A j and let St j denote the corresponding set of states. For each state s = (ℓ, ν) of A j , we shall assume that there exist probability distributions for delays, output as well as next-state:

• the delay density function, µ s over delays in R ≥0 , provides stochastic information for when the component will perform an output, thus µ s (t)dt = 1;

• the output probability function γ s assigns probabilities for resolving what output o ∈ Σ j o to generate, i.e. ∑ o γ s (o) = 1;

• the next-state density function η a s provide stochastic information on the next state s ′ = (ℓ ′ , ν ′ ) ∈ R X given an action a, i.e. s ′ η a s (s ′ ) = 1. For outputs happening deterministically at an exact time point d (or deterministic next states s ′ ), µ s (η a s ) becomes a Dirac delta function δ d (δ s ′ ) 4 .

In UPPAAL-SMC uniform distributions are applied for states where delay is bounded, and exponential distributions (with location-specified rates) are applied for the cases, where a component can remain indefinitely in a location. Also, UPPAAL-SMC provides syntax for assigning discrete probabilities to different outputs as well as specifying stochastic distributions on next-states (using the function random . . a k ) the set of all maximal runs from s with a prefix t 1 a 1 t 2 a 2 . . .t k a k for some t 1 , . . . ,t n ∈ R ≥0 , that is runs where the i'th action a i has been outputted by the component A c(a i ) . Providing the basic elements of a Sigma-algebra, we now inductively define the following measure for such sets of runs:

P A π(s, a 1 . . . a n ) = t≥0 µ s c (t) • ∏ j =c τ>t µ s j (τ)dτ • γ s c t (a 1 ) • s ′ ∏ j η a 1 s t j (s ′ j ) • P A π(s ′ , a 2 . . . a n ) ds ′ dt
This definition requires a few words of explanation: at the outermost level we integrate over all possible initial delays t. For a given delay t, the outputting component c = c(a 1 ) will choose to make the broadcast at time t with the stated density. Independently, the other components will choose to a delay amount, which -in order for c to be the winner -must be larger than t; hence the product of the probabilities that they each make such a choice. Having decided for making the broadcast at time t, the probability of actually outputting a 1 is included. Finally, integrating over all global states s ′ that may result from all components having delayed t time-units and changed state stochastically with respect to the broadcasted action a 1 , the probability of runs according to the remaining actions a 2 . . . a n is taken into account.

Generation of Random Runs

The central component of the SMC engine of UPPAAL-SMC is the efficient generation of random runs according to the stochastic semantics proposed in the preceding Section. UPPAAL-SMC has to integrate the ODEs given in the model while still respecting its stochastic semantics. The tool does not solve the equations exactly and currently supports the Euler integration method. A fixed time step δt (defined by the user) is used by an internal integrator component added to the system. This integrator races with the other processes during its time step and all rates are considered constant and defined by the equations in the model. At the end of every step or if another process wins the race, all the rates are re-evaluated and the resulting values of the clocks are computed according to x new = x old + δt * x ′ old . We plan to implement more robust methods, such as Runge-Kutta's, that give different evaluations of x new . In this case, x ′ old is then derived from the method and not directly from the equation in the model and kept constant during δt. The stochastic semantics and the races between the components are still respected.

The integrator component always picks a delay equal to the given time step δt. The other components pick delays according to their distributions. The winning component delays to the point in time it wanted and tries to take a transition at that point. If it appears that no guard is satisfied, no transition is taken and all components race again. In practice, every component will use the current rates of the clocks to evaluate a lower bound on the delay from which anything can happen to skip delaying for nothing.

Statistical Model Checking

We use SMC [START_REF] Guldstrand | Bisimulation through Probabilistic Testing[END_REF][START_REF] Sen | Statistical Model Checking of Black-Box Probabilistic Systems[END_REF][START_REF] Håkan | Probabilistic Verification of Discrete Event Systems Using Acceptance Sampling[END_REF][START_REF] Legay | Statistical Model Checking: An Overview[END_REF] to estimate and test on the probability that a random run of a network of SHAs will satisfy a given property. Given a model H and a trace property ϕ (e.g. expressed in LTL [START_REF] Pnueli | The Temporal Logic of Programs[END_REF] or MTL [START_REF] Koymans | Specifying Real-Time Properties with Metric Temporal Logic[END_REF]), SMC refers to a series of simulation-based techniques that can be used to answer two questions: (1) Qualitative: is the probability that a random run of H will satisfy ϕ greater or equal to a certain threshold θ (or greater or equal to the probability to satisfy another property ϕ ′ )? and ( 2) Quantitative: what is the probability that a random run of H will satisfy ϕ? In both cases, the answer will be correct up to a user-specified level of confidence, providing a upper bound on the probability that the conclusion made by the algorithm will be wrong. For the quantitative approach, which we will use intensively in this paper, the method computes a confidence interval that is an interval of probabilities that contains the true probability to satisfy the property. Here the confidence level provides the probability that the computed confidence interval indeed contains the unknown probability.

Our UPPAAL-SMC tool-set implements a wide range of SMC algorithms for networks of SHAs not only for reachability and safety properties, but also for general weighted MTL properties [START_REF] Bulychev | Monitor-Based Statistical Model Checking for Weighted Metric Temporal Logic[END_REF][START_REF] Bulychev | Rewrite-Based Statistical Model Checking of WMTL[END_REF]. In addition, the tool offers several features to visualize and reason on the results.

Energy Aware Buildings

UPPAAL-SMC has recently [START_REF] David | An Evaluation Framework for Energy Aware Buildings using Statistical Model Checking Submitted[END_REF] been applied to an evaluation framework for energy aware buildings, where the control of heating is optimized with respect to environmental and user profiles. To indicate central parts of this framework and the clear benefit of modeling with ODEs, we illustrate the framework on a simplified instance and recall results from our case study [START_REF] David | An Evaluation Framework for Energy Aware Buildings using Statistical Model Checking Submitted[END_REF] based on a benchmark for hybrid systems verification [START_REF] Fehnker | Benchmarks for Hybrid Systems Verification[END_REF].

A Simple 2-Room Example

To illustrate the various aspects of the (extended) modeling formalism supported by UPPAAL-SMC, we consider the case of two independent rooms that can be heated by a single heater shared by the two rooms, i.e., at most one room can be heated at a time. Fig. 2(a) shows

on[1]! on[0]! off[1]! off[0]! ON_1 ON_0 3 x<=4 1 OFF 1 x=0 x=0 2 K´'==0 && T[0]´'==K -T[0]/10 +sum(j:int[0,1]) A[0][j]*(T[j]+-T[0])/s on[0]? off[0]? ON T[0]´'==-T[0]/10 +sum(j:int[0,1]) A[0][j]*(T[j]+-T[0])/s Init OFF K=9+random(3) T[0]=T0[0] K´'==0 && T[1]´'==K -T[1]/10 +sum(j:int[0,1]) A[1][j]*(T[j]+-T[1])/s on[1]? off[1]? ON T[1]´'==-T[1]/10 +sum(j:int[0,1]) A[1][j]*(T[j]+-T[1])/s Init OFF K=9+random(3) T[1]=T0[1]
(a) stochastic heater.

(b) room 0. (c) room 1. Figure 2: A simple two room example with an autonomous heater.

the automaton for the heater. It turns itself on with a uniform distribution over time in-between [0, 4] time units. With probability 1/4 room 0 is chosen and with probability 3/4 room 1. The heater stays on for some time given by an exponential distribution (rate 2 for room 0, rate 1 for the room 1). In summary, one may say that the controller is more eager to initiate the heating of room 1 than room 0, as well as less eager to stop heating room 1. The rooms are similar and are modeled by the same template instantiated twice as shown in Fig. 2(b-c). The room is initialized to its initial temperature and then depending on whether the heater is turned on or not, the evolution of the temperature is given by

T ′ i = -T i /10 + ∑ j=0,1 A i, j (T j -T i ) or T ′ i = K -T i /10 + ∑ j=0,1 A i, j (T j -T i )
where i, j = 0, 1 are room identifiers. The sum expression corresponds to an energy flow between rooms and matrix A encodes the energy transfer coefficient between adjacent rooms. Furthermore, when the heater is turned on, its heating is not exact and is picked with a uniform distribution of K ∈ [START_REF] Bulychev | Computing Nash Equilibrium in Wireless Ad Hoc Networks: A Simulation-Based Approach[END_REF][START_REF] David | An Evaluation Framework for Energy Aware Buildings using Statistical Model Checking Submitted[END_REF], realized by the update K=9+random(3).

This example illustrates the support for NSHA in UPPAAL-SMC with extended arithmetics on clocks and generalized clock rates.

Extended Input Language UPPAAL-SMC takes as input NSHA as described above. Additionally, there is support for other features of the UPPAAL model checker's input language such as integer variables, data structures and user-defined functions, which greatly ease modeling. UPPAAL-SMC allows the user to specify arbitrary rates for the clocks, which includes a mix of integer and clock expressions on any location. In addition, the automata support branching edges where weights can be added to give a distribution on discrete transitions. It is important to note that rates and weights may be general expressions that depend on the states and not just simple constants.

Checking Queries

The fundamental principle in UPPAAL-SMC is to generate runs and evaluate some expression on the states along the obtained run. Runs are always bounded, either by time, by a number of steps, or more generally by cost (when using a clock explicitly). The engine has a built-in heuristic detection of Zeno behaviours to abort the generation of such runs. Examples of the syntax for the different types of bounds are [<=100] for 100 time units since the beginning of the run, [#<=50] for 50 discrete transitions taken from the initial state, and [x<=200] until the clock x reaches 200 5 .

UPPAAL-SMC supports simulations with monitoring custom expressions, probability evaluation, hypothesis testing, and probability comparison. We can simulate and plot the temperatures with the query The query request the checker to provide one simulate run over 600 time units and plot the temperatures of Room(0) and Room(1). The heater in this example is purely stochastic and is not intended to enforce any particular property. Yet, the simulation obtained from this query in Fig. 3 shows that the heater is able to maintain the temperatures within (mostly) distinct intervals.

simulate 1 [<=600]{T[0],T[1]} T[1] T[0]
We can evaluate on a shorter time scale the probability for the temperature of Room(0) to stay below 30 and the temperature of Room(1) to stay above 5 with the queries

Pr[<=100]([] Room(0).Init || T[0] <= 20) Pr[<=100]([] Room(1).Init || T[1] >= 7)
The results are respectively in [0.45, 0.55] and [0.65, 0.75]. The precision and confidence of these socalled confidence intervals are user-defined and influence the number of runs needed to compute the probability. In this example, for having the precision to be ±0.05 with a confidence of 95%, we need 738 runs. In fact if we are only interested in knowing if the second probability is above a threshold it may be more efficient to test the hypothesis Pr[<=100]([] Room [START_REF] Alur | Optimal Paths in Weighted Timed Automata[END_REF].Init || T[1] >= 7) >= 0.69 which is accepted in our case with 902 runs for a level of significance of 95%. To obtain an answer at comparable level of precision with probability evaluation, we would need to use a precision of ±0.005, which would require 73778 runs instead.

The tool can also compare probabilities without needing to compute them individually. We can test the hypothesis that the heater is better at keeping the temperature of Room(1) above 8 than keeping the temperature of Room(0) below 20:

Pr[<=100]([] Room(1).Init || T[1] >= 7) >= Pr[<=100]([] Room(0).Init || T[0] <= 20)
which is accepted in this case with 95% level of significance with just 258 runs.

Results

In [START_REF] David | An Evaluation Framework for Energy Aware Buildings using Statistical Model Checking Submitted[END_REF] we have estimated the comfort time (duration while being in comfortable temperature range) and energy consumption for various weather conditions, user profiles and central controller strategies. Fig. 4 shows six energy consumption estimates in different configurations (a building with 5 rooms and 3 heaters). The energy comparison shows that the dynamic user profile can save more than 33% of energy regardless of the chosen central controller strategy.

Biological Oscillator

One of the key oscillatory behaviours in biology is the circadian rhythm that allows an organism to take advantage of periods of day and night to optimise when to maximise activity and recovery. We show how the genetic circadian oscillator of [START_REF] Barkai | Biological rhythms: Circadian clocks limited by noise[END_REF][START_REF] José | Mechanisms of noiseresistance in genetic oscillators[END_REF] can be modelled and analysed using UPPAAL-SMC. This synthetic model distils the essence of several real circadian oscillators and demonstrates how a reliable system can be constructed in the face of inherent stochasticity. Figure 5a shows a system of differential equations from [START_REF] José | Mechanisms of noiseresistance in genetic oscillators[END_REF]. The equations are typeset in UPPAAL as invariant expression on a location shown in Fig. 5b, where each quantity D A , D R , D ′ A , M A , M R , A, R,C are modelled as continuous clock variables DA, DR, D A, MA etc with rates defined by a corresponding differential equation. The preceding expressions about alphaA, alpha A, alphaR, alpha R, beta A and so on are modelling the constants α A , α ′ A , α R , α ′ R , β A and so on. The assignments on the first transition initialise all the variables with initial conditions. UPPAAL-SMC is then used to simulate the model and provide plots of how the variable values evolve over time, which are displayed in Figure 7a.

dD A /dt = θ A D ′ A -γ A D A A dD R /dt = θ R D ′ R -γ R D R A dD ′ A /dt = γ A D A A -θ A D ′ A dD ′ R /dt = γ R D R A -θ R D ′ R dM A /dt = α ′ A D ′ A + α A D A -δ M A M A dM R /dt = α ′ R D ′ R + α R D R -δ M R M R dA/dt = β A M A + θ A D ′ A + θ R D ′ R -A(γ A D A + γ R D R + γ C R + δ A ) dR/dt = β R M R -γ C AR + δ A C -δ R R dC/dt = γ C AR -δ A C
(a) Ordinary differential equations. The ODE system can be interpreted as a system behavior near the thermodynamic limit (infinite population sizes while maintaining the same concentrations). Alternatively this oscilator can be modeled as a system of stochastic chemical reactions, where each molecule is counted as discrete entity moving according to Brownian motion laws. Using a standard translation between deterministic and stochastic semantics of chemically reacting systems (e.g., Gillespie's algorithm [START_REF] Daniel | Exact Stochastic Simulation of Coupled Chemical Reactions[END_REF]) the coefficients in ODE can be interpreted as reaction rates. The reactions are enumerated in Fig. 6a. Each reaction is then modeled

A + D A γ A -→ D ′ A D ′ A θ A -→ D A + A A + D R α R -→ D ′ R D ′ R θ R -→ D R + A D ′ A α ′ A -→ M A + D ′ A D A α A -→ M A + D A D ′ R α ′ R -→ M R + D ′ R D R α R -→ M R + D R M A β A -→ M A + A M R β R -→ M R + R A + R γ C -→ C C δ A -→ R A δ A -→ / 0 R δ R -→ / 0 M A δ M A -→ / 0 M R δ M R -→ / 0
(a) Reactions. as a separate SHA process shown in Fig. 6b, which can be viewed as an encoding of continuous-time Markov chain process. For example, the first reaction means taking one molecule of each A and D A substance and producing D ′ A with a rate of γ A , which can be interpreted as a transition requiring positive amount of A and DA (modeled as integer variables), consuming one molecule each and producing one D A and the reaction rate is gammaA and proportional to available quantities of A and DA. The resulting trajectories of quantities A, C and R are displayed in Fig. 7b, where the patterns seem to resemble the one from ODE model (Fig. 7a), but the inherent stochasticity result in shaky lines (even saw teeth) and unpredictably fluctuating amplitudes and their periods. Interestingly, our formalism is flexible enough to accommodate a hybrid model combining stochastic aspects as in Fig. 7b and continuous aspects as in Fig. 7a.

The amplitude of each protein quantity can be measured by the queries E[<=75; 2000](max: q), where 75 is the time limit for simulation, 2000 is the number of simulations and q is either A, C, or R. The upper plots of Fig. 8 show the probability density for a range of possible values of amplitude with a vertical line for average value. UPPAAL-SMC can also estimate a distance between peaks by using techniques developed for MITL (Metric Interval Temporal Logic -a more expressive property language than a subset of CTL supported by UPPAAL). The idea of the approach is to translate MITL formula into a monitoring automata which start an auxiliary clock x with a first peak and stop with a second peak [START_REF] Bulychev | Monitor-Based Statistical Model Checking for Weighted Metric Temporal Logic[END_REF]. To detect peaks of A when its amount rises above 1100 and drops below 1000 within 5 time units, we use the formula (in the tool syntax): true U [<=1000] (A>1100 & true U[<=5] A<=1000). Then the distance between peaks can be estimated by measuring maximal value of clock x. The result is shown as logarithm of probability density plots in a second row of Fig. 8. The plots show that in most cases the measured distance between peaks is about 24.2 hours (slightly more than one day-night cycle). Then there are some smaller bumps with several magnitudes lower probability which can be explained by either a) false positive peak as MITL monitor is confused by a sudden stochastic saw tooth in signal A, or b) missing a peak or two, or even three (in C) if the peak is not high enough to be registered, hence the next one is registered instead.

Conclusions

This paper presents an extensions of the UPPAAL-SMC model checker supporting statistical model checking for stochastic hybrid automata, where dynamic of systems can be specified with ODEs. This is a major advance in comparison to existing SMC checkers that can only handle simple derivatives. Currently, the tool applies the Euler integration method. In the future, we plan to implement more robust methods, such as Runge-Kutta's. Another contribution will be to support rare-events as it is the case in [START_REF] Zuliani | Rare-event verification for stochastic hybrid systems[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF].

  Model of y coordinate.

  Trajectory as (x, y) plot.

Figure 1 :

 1 Figure 1: Models and a trajectory of a thrown/bouncing ball hit by a player.

  [b] denoting a uniform distribution on [0, b]). Stochastic Semantics of Networks of HA For the stochastic semantics of closed networks of HA consider A = (A 1 | . . . | A n ) with a state space St = St 1 × • • • × St n . For s = (s 1 , . . . , s n ) ∈ St and a 1 a 2 . . . a k ∈ Σ * we denote by π(s, a 1 a 2 .

Figure 3 :

 3 Figure 3: Evolution of the temperatures of the two rooms.

Figure 4 :

 4 Figure 4: Energy consumption estimates for static and dynamic user profiles.

  alphaA'==0 && alpha_A'==0 && alphaR'==0 && alpha_R'==0 && betaA'==0 && betaR'==0 && deltaA'==0 && deltaR'==0 && deltaMA'==0 && deltaMR'==0 && gammaA'==0 && gammaR'==0 && gammaC'==0 && thetaA'==0 && thetaR'==0 && DA'== thetaA*D_A-gammaA*DA*A && DR'== thetaR*D_R-gammaR*DR*A && D_A'== gammaA*DA*A-thetaA*D_A && D_R'== gammaR*DR*A-thetaR*D_R && MA'== alpha_A*D_A+alphaA*DA-deltaMA*MA && MR'== alpha_R*D_R+alphaR*DR-deltaMR*MR && A'== betaA*MA+thetaA*D_A+thetaR*D_R -A*(gammaA*DA+gammaR*DR+gammaC*R+deltaA) && R'== betaR*MR-gammaC*A*R+deltaA*C-deltaR*R && C'== gammaC*A*R-deltaA*C alphaA=50, alpha_A=500, alphaR=0.01, alpha_R=50, betaA=50, betaR=5, deltaMA=10, deltaMR=0.5, deltaA=1, deltaR=0.2, gammaA=1, gammaR=1, gammaC=2, thetaA=50, thetaR=100, DA=1, DR=1, D_A=0, D_R=0, MA=0, MR=0, A=0, R=0, C=0 (b) UPPAAL automaton representation.

Figure 5 :

 5 Figure 5: Dynamics of genetic oscillator.

  UPPAAL model representation.

Figure 6 :

 6 Figure 6: Stochastic model of genetic oscilator.

  Stochastic model simulation plot.

Figure 7 :

 7 Figure 7: UPPAAL-SMC simulations: simulate 1 [<=75] { A, C, R }.

Figure 8 :

 8 Figure 8: Estimated probability density distributions for amplitude and period.

Our version is inspired by the recently announced solution by a 16 year old German school boy (Shouryya Ray), to a

year open problem by Newton concerned with predicting the trajectory of a ball thrown at a wall.

www.uppaal.org.

which should formally be treated as the limit of a sequence of delay density functions with decreasing, non-zero support around d.

It is up to the modeler to ensure that the clock eventually reaches the bound.

* Work partially supported by the VKR Centre of Excellence MT-LAB, the Sino-Danish Basic Research Center IDEA4CPS, and by the CREATIVE project ESTASE.

1 in contrast to the usual restriction of priced timed automata [7,1]


Fabrizio Biondi

Jean Quilbeuf

Axel Legay

Information Leakage by Trace Analysis in QUAIL

Quantitative security techniques have been proven effective to measure the security of systems against various types of attackers. However, such techniques are based on computing exponentially large channel matrices or Markov chains, making them impractical for large programs. We propose a different approach based on abstract trace analysis. By analyzing directly sets of execution traces of the program and computing security measures on the results, we are able to scale down the exponential cost of the problem. Also, we are able to appy statistical simulation techniques, allowing us to obtain significant results even without exploring the full space of traces. We have implemented the resulting algorithms in the QUAIL tool. We compare their effectiveness against the state of the art LeakWatch tool on two case studies: privacy of user consumption in smart grid systems and anonymity of voters in different voting schemes.

Introduction

The protection of privacy and data security is one of the main concerns of computer science. Security often falls down to the impossibility for an attacker to obtain a given secret value. Such an impossibility can be defined by non-interference [START_REF] Goguen | Security policies and security models[END_REF]. However this definition rejects any program which publishes any variable whose value depends on the secret. For instance, publishing the results of an election when each individual vote is secret breaks non-interference. Such a yes/no approach do not consider that an attacker may have a partial information about a secret.

Information-theoretical techniques have the advantage of considering the secret not as an atomic object but as a known number of secret bits, allowing the definition of measures of effectiveness of an attack based on the amount of secret bits that the attack compromises. The amount of secret bits that are compromised by an attack are known as information leakage. Leakage depends on the information about the secret known to the attacker before the attack, known as prior information and usually modeled as a prior probability distribution over the values of the secret. This approach dates back to Denning [START_REF] Denning | A lattice model of secure information flow[END_REF]. Different information leakage measures have been introduced, including Shannon leakage [START_REF] Iii | Toward a mathematical foundation for information flow security[END_REF], min-entropy leakage [START_REF] Smith | On the foundations of quantitative information flow[END_REF] and the g-leakage [START_REF] Alvim | Measuring information leakage using generalized gain functions[END_REF], encoding different security properties of the system. Among the results in the field, Köpf et al. studied leakage of side-channel attacks [START_REF] Backes | Preventing side-channel leaks in web traffic: A formal approach[END_REF][START_REF] Köpf | Automatic quantification of cache side-channels[END_REF], while Boreale has defined leakage for process calculi [START_REF] Boreale | Quantifying information leakage in process calculi[END_REF] and characterized the best attack strategy of an adaptive attacker [START_REF] Boreale | Quantitative information flow under generic leakage functions and adaptive adversaries[END_REF].

In previous work, we have presented a quantitative analysis technique able to precisely quantify information leakage of a system to an arbitrary number of decimal digits after modeling a system-attacker scenario with a Markov chain and analyzing it by applying results in both information theory and graph theory [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF]. The method has been implemented in our QUantitative Analyzer for Information Leakage (QUAIL) tool 1 [START_REF] Biondi | QUAIL: A quantitative security analyzer for imperative code[END_REF].

To our knowledge, QUAIL is the only tool enabling precise and exact information leakage computation, and later tools by multiple authors have used it as comparison [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF][START_REF] Phan | Abstract model counting: a novel approach for quantification of information leaks[END_REF]. Nonetheless, QUAIL needs to produce a full Markov chain model of the systemattacker scenario to produce a meaningful result. Such chain has an exponential number of states in the size of the secrets and observables, making it cumbersome to compute.

The first contribution of this paper is a new algorithm for precise information leakage computation, which we call Exhaustive Trace Exploration (Exp). The algorithm is able to compute information leakage following the same Markovian semantics we introduced previously [4, Section 4] by analyzing the execution traces of the system. This avoids building the full Markov chain model of the system.

The Exp algorithm performs an exhaustive depth-first exploration of all system traces, so that every trace is explored exactly once. The probability and information leakage of each trace is computed, and the information leakage of the system is obtained as the expected value of the information leakage over all traces.

We will show the effectiveness of the algorithm by comparing it with the previous QUAIL implementation and the state of the art LeakWatch tool [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF]. LeakWatch is the most recent of a family of tools for statistical approximation of information leakage developed by Chothia et al. [START_REF] Chothia | A statistical test for information leaks using continuous mutual information[END_REF][START_REF] Chothia | Probabilistic point-to-point information leakage[END_REF]. LeakWatch analyzes Java code, requiring the programmer to annotate the code of the system with secret and observable values, then simulates the system repeatedly using the Java Virtual Machine and estimates the correlation between the secret and observable values. LeakWatch follows a different perspective than QUAIL, since the former simulates the full program a large amount of times to produce an approximated result while the latter performs a single full abstract exploration of the traces of the system to produce the exact value of the leakage.

The second contribution of this paper is the presentation of two scalable case studies for quantitative information leakage analysis and the comparison of LeakWatch and QUAIL on these case studies. The case studies are anonymity of user data in Smart Grids and privacy comparison in voting protocols.

Smart Grids are in the family of interconnected objects and have received a growing interest over the last years. Our case study is based on a real system deployed at fortiss2 labs [START_REF] Koss | Establishing a smart grid node architecture and demonstrator in an office environment using the SOA approach[END_REF]. In our case study, we focus on the negotiation between a set of prosumers and an aggregator. The prosumers (PROducer conSUMERs) consume, store and produce energy. To stabilize the grid, the prosumers negotiate with the aggregator how much energy they will exchange with the grid for the next period of time. This exchange might expose the consumption of one of the prosumers, and, in turn, allow a potential attacker to deduce that a house is empty or that a factory has increased its production. In that example, the difficulty is to decide not only whether the exact information can be deduced or not, but also how well an attack can approximate it. Measuring the leakage with QUAIL indicates how much of the secret is unveiled through the negotiation phase. We show that increasing the number of prosumers also increases security.

In the voting protocols comparison case study, we compare two different voting protocols: the Single Preference where each voter expresses a single vote for his favorite candidate, and the Preference Ranking where each voter ranks all candidates from his most to his least favorite. In both case there are multiple voters and candidates, and the secret is the preference of each voter. Both protocols have a large number of possible secrets and outputs, so they become cumbersome to analyze with both QUAIL and LeakWatch with a small number of voters and candidates. We show that the Preference Ranking protocol is more effective in protecting the privacy of the voters, and that QUAIL can prove it significantly faster and with more precision than LeakWatch. In general, the size of the secret is an important parameter for the comparison: LeakWatch perfoms better on systems with a small secret, while QUAIL scales better with the size of the secret.

The rest of the paper is structured as follows. Section 2 introduces standard concepts and notations. In Section 3 we present the algorithms for leakage computation via trace analysis. Section 4 elaborates on the case studies, and finally Section 5 discusses statistical extensions of the algorithms presented.

Information Leakage with QUAIL

We want to compute the information leakage of a program, i.e. a measure quantifying how much information an attacker infers about the program's secret by observing the program's output. We assume that the attacker has access to the program's source code, unlimited computational power, and some prior information about the secret (e.g. the bit size of the secret). Leakage corresponds to the reduction in the attacker's uncertainty about the secret.

Let h be a random variable with values in a domain D(h) representing the value of the secret and o be a random variable with values in a domain D(o) modeling the value of the output. The information the attacker has on the secret is modeled by a discrete probability distribution, i.e. for a discrete random variable X a function π : D(X) → [0, 1] such that x∈D(X) π(x) = 1. The information that the attacker has on the secret before the attack is modeled by the prior distribution π(h) while the information the attacker has after observing the output is modeled by the posterior distribution π(h|o). We consider the prior distribution as given, since it is part of the model of the attacker. Let U be an uncertainty measure defined on probability distributions, including Shannon entropy, min-entropy, and g-vulnerability. Computing leakage for the measure U reduces to computing the prior and posterior distributions and applying the formula

Leakage U = U (π(h)) -U (π(h|o)) (1)
= U (π(h)) - ō∈D(o) π(o = ō)U (π(h|o = ō)) (2)
We present the syntax of the QUAIL imperative language we use to model programs. We distinguish the variables in public and private variables according to their level of abstraction: public variables have precise values, while private variables have sets of possible values. The observable variable o is public, while the secret variable h is private.

Let v (resp. h) range over names of public (resp. private) variables and x range over reals from [0; 1]. Let L (resp. H) be the set of assignments of values to public variables (resp. sets of values to private variables). Assume that the secret is a private variable h taking values in a known domain D(h) and the observable is a public variable o taking values in a known domain D(o).

Let label denote program points and f (g) pure arithmetic (Boolean) expressions. Assume a standard set of expressions and the following statements:

stmt ::= public intn v | private intn h | v := f(L) | v := rand x | skip | goto label | return | if g(L, H) then goto label_a else goto label_b
The first statement declares a public variable v of size n bits with a given value k, while the second statement similarly declares a private variable h of size n bits with allowed values ranging from 0 to 2 n -1. The third statement assigns to a public variable the value of expression f depending on public variables; assignment to private variables or depending on the value of private variables is not allowed. The fourth statement assigns zero with probability x, and one with probability 1-x, to a public variable. The return statement outputs values of all public variables and terminates. A conditional branch first evaluates an expression g dependent on private and public variables, and it jumps to label label_a if g is true and to label_a otherwise.

Since only a single variable scope exists, loops can be added in a standard way as a syntactic sugar.

Information Leakage by Trace Analysis

We present a method to compute information leakage of a program by analyzing the execution traces of the program. We present in Subsection 3.1 how a single trace is produced according to the markovian semantics of our langage. Subsection 3.2 presents how the leakage is computed from a multiset of final states. Subsections 3.3 and 3.4 present the Exp algorithm and details about its distributed implementation. The Exp algorithm explores all the traces to produce the multiset of final states needed to compute the leakage.

Single Trace Analysis

The Markovity of the semantics allows us to define states containing enough information to determine a probability distribution over all traces originating from any state. Definition 1. A state in a Markovian semantics is a tuple (pc, L, H, p) where pc ∈ N 0 is the program counter, L an assignment function assigning a value to each public variable, H an assignment function assigning a set of values to each private variable, and 0 ≤ p ≤ 1 is the probability of the state.

The initial state of the semantics is (1, ∅, ∅, 1). The successor states of a state (pc, L, H, p) depend on the statement pointed at by the program counter pc, as shown in the semantics in Figure 1. We call a state final if the program counter of the state points to a return statement. The trace analysis terminates when a final state is encountered. This means that the analysis terminates if and only if the program under analysis terminates, so non-terminating programs cannot be analyzed with this technique. We refer to our work in the subject to handle the computation of leakage of non-terminating programs [START_REF] Biondi | Information leakage of non-terminating programs[END_REF].

pc: public intn v := k (pc, L, H, p) → (pc + 1, L ∪ {(L(v) = k, n)}, H, p) pc: private intn h (pc, L, H, p) → (pc + 1, L, H ∪ {H(h) = {0, ..., 2 n -1}, n)}, p) pc: skip (pc, L, H, p) → (pc + 1, L, H, p) pc: v := f(L) (pc, L, H, p) → (pc + 1, L ∪ {(L(v) = f (L), n)}, H, p) pc: v := rand x (pc, L, H, p) → (pc + 1, L[ 0 /v], H, p • x), (pc + 1, L[ 1 /v], H, p • (1 -x))
To analyze a single execution trace, we start from a given Markov state (pc, L, H, p), e.g. the initial state (1, ∅, ∅, 1). The program counter points to the line in the source code to be executed according to the semantics in Figure 1.

Conditional states and random assignment states have two successors. The successors of a conditional state correspond to the guard being true or false. Since the guard can depend on the secret, both successor states may have positive probability depending on the prior distribution π(h) on the secret, which is available at this time. The successors of a random assignment state correspond to the bit being set to 0 or 1. In both cases the probability of each successor state is computed and one of the successor states with positive probability is chosen to be the next step in the analysis.

Because of the Markovian semantics, each state contains the information to compute the probability distribution over its outgoing transitions. The probability of a trace is computed as the product of the probabilities of the transitions composing the trace. In the successor states of the conditional statement, H|g(L, H) (resp. H|¬g(L, H)) represents the assignment function obtained by removing from the sets of values assigned to the private variables those values that contradict (resp. respect) the guard g(L, H). Similarly, P r(g(L, H)|π(h)) (resp. P r(¬g(L, H)|π(h))) refers to the probability Data: uncertainty measure U , multiset Q of final states Result: posterior uncertainty U (π(h|o)) 1 Initialize π(o) and all π(h, o = ō) to zero;

2 forall the s = (pc, L, H, p) ∈ Q do 3 Let ō = L(o), {k1, ..., kn} = H(h); 4 Set π(o = ō) ← π(o = ō) + p; 5 for i = 1...n do 6 Set π(h = ki, o = ō) ← π(h = ki, o = ō) + p/n; 7 end 8 end 9 For each ō ∈ D(o) let π(h|o = ō) ← π(h, o = ō)/π(o = ō); 10 Return U (π(h|o)) = ō∈D(o) π(o = ō)U (π(h|o = ō))
Algorithm 1: Posterior uncertainty computation that the guard g(L, H) is true (resp. false) considering the prior probability distribution π(h) on the private variables.

When a single trace analysis terminates, a given final state ( pc, L, H, p) is produced, in which pc points to a return statement. The sets of allowed values assigned to the private variables in H have been appropriately reduced to account for the conditional statements visited by the trace.

Posterior Uncertainty Computation Algorithm

We show how to compute the posterior uncertainty U (π(h|o)) of a system with a secret h and an observable o, given the uncertainty measure U and a multiset Q of final states of the system. Q encodes the posterior joint probability of all variables in the system and can be produced by the Exp algorithm presented in Section 3.3.

Let (pc, L, H, p) be a final state in Q, where L represents the assignments of given values to the public variables, H the assignments of sets of values to the private variables, and p the joint probability of such assignments. Since different traces may produce the same final assignments to variables (L, H), the joint probability of these assignments is the sum of the probabilities of all such final states. To apply the formula ( 2

Exhaustive Trace Exploration Algorithm

We present the Exhaustive Trace Exploration Algorithm (Exp) for the efficient computation of information leakage of a program written in an imperative code equipped with Markovian semantics. The algorithm can be used to compute any information-theoretical measure of leakage, including Shannon leakage, min-entropy leakage, guessing entropy and g-leakage. In the algorithm, the prior and posterior distributions of the attacker are computed, and information leakage is computed as the difference between the value of the chosen measure on the prior and posterior distribution, following common practice in the field. The algorithm can be parallelized to take advantage of multicore architectures and is implemented in the current release of the QUAIL tool, available at http://project.inria.fr/quail. Details about the implementation are given in Section 3.4.

The Exp algorithm perform a depth-first exhaustive exploration of the execution traces of the system. Each trace is explored until it gets to a final state, then the final state gets added to the multiset Q of final states. When a state with more than one successor is found during the exploration of a trace, one of the successor state is chosen as the next state to explore and all other successors are put on the stack of the states still to be explored, following a depth-first strategy. We refer back to Section 3.1 for details about the analysis of a single trace. When all traces have been explored, Exp has produced the full multiset of final states Q of the system. The algorithm then computes the information leakage of the system using Algorithm 1. The pseudocode of the algorithm is presented in Algorithm 2. The following theorem proves the correctness of the algorithm: Theorem 2. Algorithm 2 terminates on a terminating program and when it does it outputs the information leakage according to the given uncertainty measure U .

Distributed Implementation

Algorithm Exp is highly parallelizable, and can be implemented to take advantage of multicore and distributed systems. To parallelize it, the stack S and the set Q of final states must be shared between the processes. Different processes can pop states from S, explore the corresponding traces until they find final states, and add those final states to Q. To avoid race conditions, access to both S and Q must be sequential. QUAIL currently implements this strategy on multicore systems.

Algorithm 1 can also be parallelized. Divide Q in subsets Q ō = {s ∈ Q|o = ō in s} for each output ō. It is then possible to utilize a different process for each output ō to analyze each set Q ō and compute the probability and conditional posterior uncertainty of ō independently from each other.

Case Studies

We compare different leakage analysis tools on two case studies. In the first one, we measure information about occupancy of an house based on the information obtained through a Smart Grid. In the second one, we compare the loss of anonymity of the voters due to the publication of the results of an election with different voting protocols.In both case we measure Shannon leakage, setting U (π(h)) = h∈D(h) π(h = h) log 2 π(h = h). The code of all the case studies presented in this paper is available at https://project.inria.fr/quail/tacas15/.

For our both case studies, we compare Exp with LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF]. On the first case study we also compare Exp with the previous algorithm implemented in QUAIL [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF]. LeakWatch evaluates the leakage in Java programs by approximating it from a statistical sample of program executions. In order to run our case studies with LeakWatch, we translate the QUAIL code into Java code. This translation is straightforward except for the declaration of private variables. To handle them, we replace each interval in a declaration by an assignment to a random number within that interval. For instance, the declaration Then we instrument the code to define the secret and observable variables by using the LeakWatch API. In our example, defining the vote as secret and the results as observable is obtained through: ' ,v o t e [ i ] ) ; f o r ( i n t j ; j <C ; j ++) LeakWatchAPI . observe ( r e s u l t [ j ] ) ;

f o r ( i n t i ; i <N ; i ++) LeakWatchAPI . secret ( ' v o t e [ i ]
After compilation of the instrumented code, the analysis is run by the command java -jar leakwatch-0.5.jar SinglePreference.

This executes the SinglePreference class several times and records the secrets and observables to compute the leakage. LeakWatch determines at runtime how much simulations are necessary to obtain a good approximation of the leakage. A Smart Grid is an energy network where every node may produce, store and consume energy. Nodes are called prosumers (PROducer consSUMERS). The Living Lab demonstrator [START_REF] Koss | Establishing a smart grid node architecture and demonstrator in an office environment using the SOA approach[END_REF], is an instance of such a prosumer, whose data can be accessed online 3 . Without coordination between the prosumers, the grid may globally consume or produce too much energy, making it unstable. To stabilize the grid, the prosumers periodically negotiate with an aggregator in charge of balancing the consumption and production among several prosumers. Figure 2 depicts a grid with 3 prosumers. Each prosumer declares its plan, that is, how much it intends to consume or produce during the next period of time. The aggregator sends back a value indicating the excess of energy production or consumption. An excess of 0 indicates that the announced plans are feasible and terminates the negotiation. Otherwise, the prosumers then have to adapt their plan accordingly and send the updated version.

Smart Grid

Switching to a smart grid structure will probably be required in the coming years. Smart grid and smart sensors raise several security and privacy concerns. The platform can ensure the information cannot flow directly between prosumers [START_REF] Bytschkow | Distributed MILS architectural approach for secure smart grids[END_REF]. However, stability requires a feedback from the aggregator that potentially conveys infomation about other prosumer, where only the software can limit information leakage. In general, knowing the consumption of a particular household may reveal some sensitive information about the house (presence of people in the house, type of electrical devices . . . ). Therefore, the consumption of a prosumer should remain secret. The privacy of a prosumer with respect to the aggregator can be ensured in several ways [START_REF] Rottondi | A privacy-friendly framework for vehicle-to-grid interactions[END_REF]. However, each prosumer receives some information about the consumption of other prosumers through the excess value sent back by the aggregator. An attacker might try to use the information obtained through the grid in order to decide whether a given house is occupied or not. In our scenario, we assume different types of houses with different consumptions. Each house is modeled by a private boolean value, which is true if the house is occupied. An occupied house consumes a fixed amount of energy, according to its type. An empty house does not consume anything. Table 1 presents how much a given house consumes, in two different cases that we consider.

For our analyses, we assume that the attacker observes the global consumption, declared by: observable int32 global_consumption;. We consider two variants, depending on whether the attacker focuses on a single house or on all the houses. The secret is either the occupancy of one particular house: secret int1 presence_target:=[0,1]; or the occupancy of all houses: secret array [N] of int1 presence:=[0,1];

Table 2 shows the results of several leakage analyses. The first two columns indicates the case, as presented in Table 1 and the number of houses in the model. For a model with N houses, there is N/3 house of each type. The columns S, M and L indicates the leakage of the variable representing the vacancy of an house of this type. The column "Global leakage" contains the leakage of the whole array of occupancy bits and the column "Global leakage/bit" indicates the average leak for one bit of secret. The average leakage per bit from a global attack is more important that the information obtained by focusing on a single bit. This means than obtaining information about the whole array, for instance the number of occupied houses, is easier than obtaining information about a single bit, i.e. occupancy of a single house. In both cases, the leakage, and thus the anonymity of prosumers, diminishes when the number of houses increases.

In Table 3 we compare the execution time between the previous QUAIL algorithm that builds a global Markov chain from the program, Exp and LeakWatch. The Exp algorithm is much faster as it needs 3s to perform a computation that requires more than one hour (the timeout value) with the previous algorithm. LeakWatch needs less than 1s on the same example. This experiment validates the fact than LeakWatch performs better than Exp when the size of the secret is small (1 bit here).

Voting protocols

In an election, each voter is called to express his preference for the competing candidates. The voting system defines the way the voters express their preference: either on paper in a traditional election, or electronically in e-voting. The voting system also comprehends the additional procedures enforced to guarantee that the voters can vote freely, that they can verify that their vote has been counted and that their vote remains confidential.

After the votes have been cast, the results of the vote are published, usually in an aggregated form to protect the anonymity of the voters. Finally, the winning candidate or candidates is chosen according to a given electoral formula.

In this section we present two different typologies of voting, representing two different ways in which the voters can express their preference: in the Single Preference protocol the voters declare their preference for exactly one of the candidates, while in the Preference Ranking protocol each voter ranks the candidate from his most favorite to his least favorite. Since each protocol we model is concerned only about how the votes are expressed and counted and what results are published, each protocol models a number of electoral formulae. For the same reason, the models are valid both for uninominal and multinominal elections.

Single Preference

The Single Preference protocol typology models all electoral formulae in which each of the N voters expresses one vote for one of the C candidates, including plurality and majority voting systems and single non-transferable vote [START_REF] Norris | Electoral Engineering: Voting Rules and Political Behavior[END_REF]. The votes for each candidate are summed up and only the results are published, thus hiding information about which voter voted for which candidate. The candidate or candidates to be elected are decided according to the electoral formula used.

The secrets and observables are modeled by the following lines of QUAIL code: The secret is an array of integers with a value for each of the N voters. Each value is a number from 0 to C -1, representing a vote for one of the C candidates. The observable is an array of integers with the votes obtained by each of the C candidates. The protocol is simple, and its information leakage can be computed formally, as shown by the following lemma: Lemma 1. The information leakage for the Single Preference protocol with n voters and c candidates corresponds to

- k1+k2+...+kc=n n! c n k 1 !k 2 ! . . . k c ! log 2 n! c n k 1 !k 2 ! . . . k c !
The proof for Lemma 1 is in the Appendix. While the lemma characterizes the solution computed by QUAIL for this case, it is very hard to find such a characterization for any process, so in general QUAIL is the best way to obtain a result. We run QUAIL with the command ./quail single_preference.quail -v 0 -p 5 and obtain 1.8112 showing that the leakage of the Single Preference protocol for 3 voters and 2 candidates is ≈ 1.8112 bits.

Preference Ranking The Preference Ranking protocol typology models all electoral formulae in which each of the n voters expresses an order of preference of the c candidates, including the alternative vote and single transferable vote systems [START_REF] Norris | Electoral Engineering: Voting Rules and Political Behavior[END_REF]. In the Preferential Voting protocol the voter does not express a single vote, but rather a ranking of the candidates; thus if the candidates are A, B, C and D the voter could express the fact that he prefers B, then D, then C and finally A. Then each candidate gets c points for each time he appears as first choice, c -1 points for each time he appears as second choice, and so on. The points of each candidate are summed up and the results are published.

The secrets and observables are modeled by the following lines of QUAIL code: The secret is an array of integers with a value for each of the N voters. Each value is a number from 0 to C! -1, representing one of the possible C! rankings of the C candidates. The observable is an array of integers with the points obtained by each of the C candidates. The full model for this protocol is shown in the Appendix due to space constraints. Table 5 shows the percent error of the leakage value obtained with LeakWatch. Indeed, LeakWatch computes an approximation of the leakage based on simulation, whereas QUAIL computes the exact value. For this reason, the leakage computed by LeakWatch for a given program may change at each invocation of the tool. LeakWatch slightly underestimates the leakage.

We compare the execution time of QUAIL and LeakWatch in Table 6 for Single Preference and in Table 7 for Preference Ranking. These execution times have been obtained on a laptop with a i7 quad-core running at 3.3Ghz and 16G of ram. The results show that QUAIL is significantly faster than LeakWatch on these examples. For Single Proof. It is known that, since the program is deterministic, its Shannon leakage corresponds to the entropy H(π(o)) of the distribution over the output [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF]. We compute

H(π(o)) = - o∈D(o) π(o) log 2 (π(o))
Here o corresponds to a possible output, which is, in our case a vote result. A vote result precises, for each candidate j, the number k j of votes obtained. Furthermore, the total of all votes is the number of voters v.

Fig. 1 :

 1 Fig. 1: Execution rules in Markovian trace semantics.

Theorem 1 . 5 while s is not a final state do 6 Find 8 Assign s := s1; 9 end 10 Add

 1568910 ) U (π(h|o)) = ō∈D(o) π(o = ō)U (π(h|o = ō)), we need to compute the marginal probability distribution π(o) and for each observable output ō ∈ D(o) s.t. π(o = ō) > 0 the corresponding conditional probability distribution on h, i.e. π(h|o = ō). Algorithm 1 computes π(o) and each π(h|o = ō) by analyzing a multiset of final states. For each state (pc, L, H, p) the value of the observable variable o and set of values of the secret variable h are analyzed. The probability of observing the value ō of the observable variable in the state is increased by p, and the probability of observing each of the n values of the secret variable conditioned on ō is increased by p/n. Finally, the probability on each subdistribution π(h, o = ō) is normalized to 1 by dividing it by π(o = ō) to obtain the conditional probability π(h|o = ō). Algorithm 1 terminates, and when it does it outputs the posterior uncertainty U (π(h|o)) of the posterior distribution represented by Q. Data: uncertainty measure U , prior distribution π(h), imperative code of the program Result: leakage of the program according to U 1 Initialize an empty multiset Q; 2 Put the initial state (1, ∅, ∅, 1) on the stack S; 3 while stack S is not empty do 4 Pop a state s = (pc, L, H, p) from the stack S; the successors s1, ..., sn of s that have positive probability; 7 Put s2, ...sn on the stack S; s to Q; 11 end 12 Return U (π(h))-ComputePosteriorUncertainty(U, Q); Algorithm 2: Exp: Exhaustive depth-first trace exploration

  secret array [ N ] of i n t 3 2 v o t e : = [ 0 , C-1]; becomes f o r ( i n t i ; i <N ; i ++) v o t e [ i ] = securerandom . n e x t I n t (C) ;

Fig. 2 :

 2 Fig. 2: Smart grid overview

  secret array [ N ] of i n t 3 2 v o t e : = [ 0 , C-1]; observable array [ C ] of i n t 3 2 r e s u l t ;

secret array [ N ] of i n t 3 2

 2 v o t e : = [ 0 ,C! -1]; observable array [ C ] of i n t 3 2 r e s u l t ;

  Thus the domain of o is D(o) = {k 1 , . . . , k c ∈ N | c j=1 k j = n}. The number of votes with result o = k 1 , . . . , k c correspond to the number of choices of voters for candidate 1 amongst n, i.e. n k1 , times the number of choice of voters for candidate 2 amongst the ones remaining i.e. !...kc! . As the votes are equiprobable, we have p(k 1 , . . . , k c ) = n! c n k1!k2!...kc! , hence the result.

Table 1 :

 1 Consumption of houses wrt size

	Size Case 1 Case 2
	Small	1	1
	Medium 2	3
	Large	3	5

Table 2 :

 2 Leakage of occupancy through the global consumption

						Table 3: Time to compute the leak-
						age for a large house
		Nb of Single house leakage Global Global	Case House Time Time Time
	Case Houses S	M	L leakage leakage/bit		Nb	MC Exp LW
	0	3	0.7500 0.7500 0.7500 2.7500 0.9166	1	3	0.4s 0.1s 0.3s
	0	6	0.0688 0.1466 0.2944 3.4210 0.5701	1	6	1.6s 0.4s 0.3s
	0	9	0.0214 0.0768 0.1771 3.7363 0.4151	1	9	91.4s 1.0s 0.4s
	0	12 0.0135 0.0544 0.1273 3.9479 0.3289	1	12 timeout 2.6s 0.4s
	1	3	1.0000 1.0000 1.0000 3.0000 1.0000	2	3	0.4s 0.1s 0.3s
	1	6	0.1965 0.1965 0.3687 4.0243 0.6707	2	6	1.6s 0.4s 0.5s
	1	9	0.0241 0.0808 0.2062 4.3863 0.4873	2	9	92.1s 1.0s 0.4s
	1	12 0.0074 0.0510 0.1443 4.6064 0.3838	2	12 timeout 2.7s 0.4s

Table 4 :

 4 Voting protocols: percent of secret leaked by Single Preference (on the left) and Preference Ranking (on the right) Experimental Results Table4shows the percentage of the secret leaked by the Single Preference and Preference Ranking protocols for different numbers of voters and candidates. We note that the results for 2 candidates are identical, since in this case in both protocols the voters can vote in only 2 different ways. The results obtained for Single Preference are correct with respect to the formula stated in Lemma 1. The table shows that the Single Preference protocol leaks a greater percentage of its secret than the Preference Ranking protocol.

	SP			Candidates			PR		Candidates
		2	3	4	5	6		2	3	4
	Voters	3 60.4 % 65.7 % 69.0 % 71.3 % 73.0 % 4 50.8 % 56.5 % 60.2 % 62.9 % 64.9 %	Voters	3 60.4 % 61.9 % 62.0 % 4 50.8 % 51.0 % timeout
		5 44.0 % 49.6 % 53.5 % 56.4 % 58.6 %		5 44.0 % 43.4 % timeout
		6 38.9 % 44.4 % 48.3 % 51.2 % 53.5 %		6 38.9 % 37.9 % timeout

Table 5 :

 5 Percent error of the leakage obtained by LeakWatch relatively to the one from QUAIL for Single Preference (on the left) and Preference Ranking (on the right)

	SP			Candidates			PR		Candidates
		2	3	4	5	6		2	3	4
	Voters	3 -3.8 % -3.7 % -3.2 % -2.8 % -2.2 % 4 -5.1 % -3.7 % -2.6 % -2.3 % -2.1 %	Voters	3 -3.8 % -2.6 % timeout 4 -5.1 % -2.6 % timeout
		5 -5.0 % -3.2 % -2.6 % -2.2 % -1.9 %		5 -5.0 % -2.2 % timeout
		6 -5.1 % -3.2 % -2.4 % timeout timeout		6 -5.1 % timeout timeout

http://project.inria.fr/quail

http://fortiss.org

livinglab.fortiss.org

Preference with 6 candidates and 5 voters, QUAIL is 200 times faster than LeakWatch. This shows that QUAIL performs better than LeakWatch with large secrets.

Towards Statistical Estimation: A Perspective

One of the reasons why LeakWatch is faster than the Exp algorithm in some examples is that the former just computes an approximation of the leakage, while the latter computes the precise result. This requires Exp to explore the totality of the execution traces of the system. However, since the space of traces is generally very large, it would be interesting to be able to analyze only a statistically significant sample of the traces and to estimate the leakage of the system based on the sample. Statistical computation of information leakage has been explored by the group of Tom Chothia, leading to the development of the LeakWatch tool we compare against, while Boreale proves that no Monte Carlo estimator exists under very mild conditions for the estimation of channel capacity [START_REF] Boreale | On formally bounding information leakage by statistical estimation[END_REF], a problem closely related to information leakage. Statistical estimation of leakage via trace analysis would follow a different path, as explained in Algorithm 3. We call this algorithm Sim. The Sim algorithm randomly chooses traces to analyze, and terminates when some termination conditions are met. Such conditions can include a certain number or percentage of traces analyzed or a time limit. Contrarily to Exp, Sim may analyze the same trace twice. In that case, the resulting final state should be added only once to the set Q. This can be obtained by enriching the state s with a hash value depending on the trace producing the final state s. After the analysis, leakage is computed from the data available. Compared to Monte Carlo techniques, Sim has the advantage that since we know the probability of each trace explored, the percentage of the probability space actually covered by the simulation is known. In case not all final states are contained in Q, the distribution π(o) computed by Algorithm 1 will be a subdistribution. In this case we increase the posterior uncertainty accordingly, assuming that the unanalyzed part of the posterior distribution behaves like the part that has been analyzed. This assumption is valid only if the final states in Q have been produced by statistically independent analyses, so that they represent a meaningful sample of the full set of final states.

Finally, note that in Algorithm 3 we are choosing the successor state to explore randomly and locally, so the probability of analyzing a trace does not correspond with the probability of the trace. This is useful when some of the traces have a very low probability, as it is often the case with authentication protocols. This kind of rare event simulation is a well established approach in probabilistic model checking. In practice we are simulating for a prior distribution π(h) but choosing the traces to simulate according to a different prior distribution π (h) that generates more interesting traces with a higher probability. This begs for the question of which prior π (s) allows us to explore the trace space more efficiently, thus providing a more meaningful sample in a shorter time. We leave this question open as future work.

A Appendix

Proof (of Theorem 1). Termination of Algorithm 1 is trivial if Q is finite. For the soundness, the algorithm computes posterior uncertainty according to the uncertainty measure U using the formula

The computation reduces to finding the probability distribution π(o) on the observable variable and the conditional probability distribution on the secret π(h|o = ō) for each possible value ō of the observable variable.

Let Q ō be the set {(pc, L, H, p) ∈ Q | L(o) = ō}, and let p(Q ō) be the sum of the probabilities p of the states in

Similarly, let h = {k 1 , . . . , k n } and Q ō, h be the set {(pc, L, H, p) ∈ Q ō | H(h) = h} and let p(Q ō, h) be the sum of the probabilities p of the states in Q ō, h. Then each state in Q ō, h represents a fragment of the the joint distribution on (o, h) that is uniform on h = {k 1 , ..., k n }, and π(o = ō, h ∈ {k 1 , ..., k n }) = p(Q ō, h). The conditional probability of each value of h conditioned by a given ō ∈ D(o) is computed as π(h|o = ō) = π(o=ō,h) π(o=ō) by the normalization in line 9.

Having computed π(o) and all π(h|o = ō), we compute posterior uncertainty with the formula presented above.

Proof (of Theorem 2). We assume that the program under analysis always terminates, i.e. all traces reach a final state in a finite number of steps. The analysis of a single trace follows the Markovian semantics of the language, as explained in Section 3.1.

The algorithm explores the tree of computation traces following a depth-first strategy. Since all traces terminate, eventually all of them are explored and the set Q contains all of the possible final states of the system. Note that the exploration also depends on the prior distribution π(h): values of the secret with a probability zero in the prior distribution are not explored. This behavior is intended, as is avoids unnecessarily exploring traces that have probability zero.

When the full set Q of final states has been produces the algorithm uses Algorithm 1 to compute the posterior uncertainty and applies the formula

so the proof of correctness is a consequence of Theorem 1.