domain: info stat math spi

text:
Cyrille Jégourel
email: cyrille.jegourel@inria.fr

Axel Legay
email: axel.legay@inria.fr

Sean Sedwards
email: sean.sedwards@inria.fr

Cross entropy optimisation of importance sampling parameters for statistical model checking

Statistical model checking avoids the exponential growth of states associated with probabilistic model checking by estimating properties from multiple executions of a system and by giving results within confidence bounds. Rare properties are often very important but pose a particular challenge for simulation-based approaches, hence a key objective under these circumstances is to reduce the number and length of simulations necessary to produce a given level of confidence. Importance sampling is a well-established technique that achieves this, however to maintain the advantages of statistical model checking it is necessary to find good importance sampling distributions without considering the entire state space. Motivated by the above, we present a simple algorithm that uses the notion of cross-entropy to find the optimal parameters for an importance sampling distribution. In contrast to previous work, our algorithm uses a low dimensional vector of parameters to define this distribution and thus avoids the often intractable explicit representation of a transition matrix. We show that our parametrisation leads to a unique optimum and can produce many orders of magnitude improvement in simulation efficiency. We demonstrate the efficacy of our methodology by applying it to models from reliability engineering and biochemistry.

Introduction

The need to provide accurate predictions about the behaviour of complex systems is increasingly urgent. With computational power becoming ever-more affordable and compact, computational systems are inevitably becoming increasingly concurrent, distributed and adaptive, creating a correspondingly increased burden to check that they function correctly. At the same time, users expect high performance and reliability, prompting the need for equally high performance analysis tools and techniques.

The most common method to ensure the correctness of a system is by testing it with a number of test cases having predicted outcomes that can highlight specific problems. Testing techniques have been effective discovering bugs in many industrial applications and have been incorporated into sophisticated tools [START_REF] Godefroid | Automated whitebox fuzz testing[END_REF]. Despite this, testing is limited by the need to hypothesise scenarios that may cause failure and the fact that a reasonable set of test cases is unlikely to cover all possible eventualities; errors and modes of failure may remain undetected and quantifying the likelihood of failure using a series of test cases is difficult.

Model checking is a formal technique that verifies whether a system satisfies a property specified in temporal logic under all possible scenarios. In recognition of non-deterministic systems and the fact that a Boolean answer is not always useful, probabilistic model checking quantifies the probability that a system satisfies a property. In particular, 'numerical' (alternatively 'exact') probabilistic model checking offers precise and accurate analysis by exhaustively exploring the state space of non-deterministic systems and has been successfully applied to a wide variety of protocols, algorithms and systems. The result of this technique is the exact (within limits of numerical precision) probability that a system will satisfy a property of interest, however the exponential growth of the state space limits its applicability. The typical 10 8 state limit of exhaustive approaches usually represents an insignificant fraction of the state space of real systems that may have tens of orders of magnitude more states than the number of protons in the universe (∼ 10 80 ).

Under certain circumstances it is possible to guarantee the performance of a system by specifying it in such a way that (particular) faults are impossible. Compositional reasoning and various symmetry reduction techniques can also be used to combat state-space explosion, but in general the size, unpredictability and heterogeneity of real systems [START_REF] Basu | Statistical abstraction and model-checking of large heterogeneous systems[END_REF] make these techniques infeasible. Static analysis has also been highly successful in analysing and debugging software and other systems, although it cannot match the precision of quantitative analysis of dynamic properties achieved using probabilistic and stochastic temporal logic.

While the state space explosion problem is unlikely to ever be adequately solved for all systems, simulation-based approaches are becoming increasingly tractable due to the availability of high performance hardware and algorithms. In particular, statistical model checking (SMC) combines the simplicity of testing with the formality and precision of numerical model checking; the core idea being to create multiple independent execution traces of the system and individually verify whether they satisfy some given property. By modelling the executions as a Bernoulli random variable and using advanced statistical techniques, such as Bayesian inference [START_REF] Sumit | A Bayesian Approach to Model Checking Biological Systems[END_REF] and hypothesis testing [START_REF] Younes | Probabilistic verification of discrete event systems using acceptance sampling[END_REF], the results are combined in an efficient manner to decide whether the system satisfies the property with some level of confidence, or to estimate the probability that it does. Knowing a result with less than 100% confidence is often sufficient in real applications, since the confidence bounds may be made arbitrarily tight. Moreover, statistical model checking may offer the only feasible means of quantifying the performance of many complex systems. Evidence of this is that statistical model checking has been used to find bugs in large, heterogeneous aircraft systems [START_REF] Basu | Statistical abstraction and model-checking of large heterogeneous systems[END_REF]. Notable statistical model checking platforms include APMC [START_REF] Hérault | Approximate probabilistic model checking[END_REF], YMER [START_REF] Younes | YMER: A Statistical Model Checker[END_REF] and VESTA [START_REF] Sen | VESTA: A statistical model-checker and analyzer for probabilistic systems[END_REF]. Well-established numerical model checkers, such as PRISM [START_REF] Kwiatkowska | PRISM: Probabilistic Symbolic Model Checker[END_REF] and UP-PAAL [START_REF] Bengtsson | Uppaal a tool suite for automatic verification of real-time systems[END_REF], are now also including statistical model checking engines to cope with larger models.

A key challenge facing statistical model checking is to reduce the length (steps and cpu time) and number of simulation traces necessary to achieve a result with given confidence. The current proliferation of parallel computer architectures (multiple cpu cores, grids, clusters, clouds and general purpose computing on graphics processors, etc.) favours statistical model checking by making the production of multiple independent simulation runs relatively easy. Despite this, certain models still require a large number of simulation steps to verify a property and it is thus necessary to make simulation as efficient as possible. Rare (unlikely) properties pose a particular problem for simulation-based approaches, since they are not only difficult to observe (by definition) but their probability is difficult to bound [START_REF] Heidelberger | Fast simulation of rare events in queueing and reliability models[END_REF].

The term 'rare event' is ubiquitous in the literature, but here we specifically consider rare properties defined in temporal logic. This distinguishes rare states from rare paths that may or may not contain rare states. The distinction does not significantly alter the mathematical derivation of our algorithm, however it affects the applicability of simple heuristics that are able to find (reasonably) good importance sampling distributions. This is of relevance because our algorithm works by a process of iterative refinement, starting from an initial distribution that must produce at least a few traces that satisfy the property. In what follows we consider discrete space Markov models and present a simple algorithm to find an optimal set of importance sampling parameters, using the concept of minimum cross-entropy [START_REF] Kullback | Information Theory and Statistics[END_REF][START_REF] Shore | Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy[END_REF]. Our parametrisation arises naturally from the syntactic description of the model and thus constitutes a low dimensional vector in comparison to the state space of the model. We show that this parametrisation has a unique optimum and demonstrate its effectiveness on reliability and (bio)chemical models. We describe the advantages and potential pitfalls of our approach and highlight areas for future research.

Importance sampling

Our goal is to estimate the probability of a property by simulation and bound the error of our estimation. When the property is not rare there are standard bounding formulae (e.g., the Chernoff and Hoeffding bounds [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF][START_REF] Hoeffding | Probability Inequalities for Sums of Bounded Random Variables[END_REF]) that relate absolute error, confidence and the required number of simulations to achieve them, independent of the probability of the property. As the property becomes rarer, however, absolute error ceases to be useful and it is necessary to consider relative error, defined as the standard deviation of the estimate divided by its expectation. With Monte Carlo simulation relative error is unbounded with increasing rarity [START_REF]Rare Event Simulation using Monte Carlo Methods[END_REF], but it is possible to bound the error by means of importance sampling [START_REF] Shahabuddin | Importance Sampling for the Simulation of Highly Reliable Markovian Systems[END_REF][START_REF] Heidelberger | Fast simulation of rare events in queueing and reliability models[END_REF].

Importance sampling is a technique that can improve the efficiency of simulating rare events and has been receiving considerable interest of late in the field of statistical model checking (e.g., [START_REF] Clarke | Statistical model checking for cyber-physical systems[END_REF][START_REF] Barbot | Coupling and Importance Sampling for Statistical Model Checking[END_REF]). It works by simulating under an (importance sampling) distribution that makes a property more likely to be seen and then uses the results to calculate the probability under the original distri-bution by compensating for the differences. The concept arose from work on the 'Monte Carlo method' [START_REF] Metropolis | The Monte Carlo Method[END_REF] in the Manhattan project during the 1940s and was originally used to quantify the performance of materials and solve otherwise intractable analytical problems with limited computer power (see, e.g., [START_REF] Kahn | Stochastic (Monte Carlo) Attenuation Analysis[END_REF]). For importance sampling to be effective it is necessary to define a 'good' importance sampling distribution: (i) the property of interest must be seen frequently in simulations and (ii) the distribution of the paths that satisfy the property in the importance sampling distribution must be as close as possible to the distribution of the same paths in the original distribution (up to a normalising factor). The literature in this field sometimes uses the term 'zero variance' to describe an optimal importance sampling distribution, referring to the fact that with an optimum importance sampling distribution all simulated paths satisfy the property and thus there is zero variance in the Bernoulli random variable. It is important to note, however, that a sub-optimal distribution may meet requirement (i) without necessarily meeting requirement (ii). Failure to consider (ii) can result in gross errors and overestimates of confidence (e.g. a distribution that simulates just one path that satisfies the given property). The algorithm we present in Section 3 addresses both (i) and (ii).

Importance sampling schemes fall into two broad categories: state dependent tilting and state independent tilting [START_REF] De Boer | Adaptive importance sampling simulation of queueing networks[END_REF]. State dependent tilting refers to importance sampling distributions that individually bias ('tilt') every transition probability in the system. State independent tilting refers to importance sampling distributions that change classes of transition probabilities, independent of state. The former offers greatest precision but is infeasible for large models. The latter is more tractable but may not produce good importance sampling distributions. Our approach is a kind of parametrised tilting that potentially affects all transitions differently, but does so according to a set of parameters.

Estimators

Let Ω be a probability space of paths, with f a probability density function over Ω and z(ω) ∈ {0, 1} a function indicating whether a path ω satisfies some property φ. In the present context, z is defined by a formula of an arbitrary temporal logic over execution traces. The expected probability γ that φ occurs in a path is then given by

γ = Ω z(ω)f (ω) dω (1)
and the standard Monte Carlo estimator of γ is given by

γ = 1 N MC NMC i=1 z(ω i )
N MC denotes the number of simulations used by the Monte Carlo estimator and note that z(ω i ) is effectively the realisation of a Bernoulli random variable with parameter γ. Hence Var(γ) = γ(1 -γ) and for γ → 0, Var(γ) ≈ γ. Let f ′ be another probability density function over Ω, absolutely continuous with zf , then Equation ( 1) can be written

γ = Ω z(ω) f (ω) f ′ (ω) f ′ (ω) dω L = f /f ′ is the likelihood ratio function, so γ = Ω L(ω)z(ω)f ′ (ω) dω (2)
We can thus estimate γ by simulating under f ′ and compensating by L:

γ = 1 N IS NIS i=1 L(ω i )z(ω i )
N IS denotes the number of simulations used by the importance sampling estimator. The goal of importance sampling is to reduce the variance of the rare event and so achieve a narrower confidence interval than the Monte Carlo estimator, resulting in N IS ≪ N MC . In general, the importance sampling distribution f ′ is chosen to produce the rare property more frequently, but this is not the only criterion. The optimal importance sampling distribution, denoted f * and defined as f conditioned on the rare event, produces only traces satisfying the rare property:

f * = zf γ (3)
This leads to the term 'zero variance estimator', but only with respect to the Bernoulli variable since, in general, Var(f * ) ≥ 0.

In the context of statistical model checking f usually arises from the specifications of a model described in some relatively high level language. Such models do not, in general, explicitly specify the probabilities of individual transitions, but do so implicitly by parametrised functions over the states. We therefore consider a class of models that can be described by guarded commands [START_REF] Edsger | Guarded commands, nondeterminacy and formal derivation of programs[END_REF] extended with stochastic rates. Our parametrisation is a vector of strictly positive values λ ∈ (R + ) n that multiply the stochastic rates and thus maintain the absolutely continuous property between distributions. Note that this class includes both discrete and continuous time Markov chains and that in the latter case our mathematical treatment works with the embedded discrete time process.

In what follows we are therefore interested in parametrised distributions and write f (•, λ), where λ = {λ 1 , . . . , λ n } is a vector of parameters, and distinguish different density functions by their parameters. In particular, µ is the original vector of the model and f (•, µ) is therefore the original density. We can thus rewrite Equation (2) as

γ = Ω L(ω)z(ω)f (ω, λ) dω
where L(ω) = f (ω, µ)/f (ω, λ). We can also rewrite Equation ( 3)

f * = zf (•, µ) γ
and write for the optimal parametrised density f (•, λ * ). We define the optimum parametrised density function as the density that minimises the cross-entropy [START_REF] Kullback | Information Theory and Statistics[END_REF] between f (•, λ) and f * for a given parametrisation and note that, in general,

f * = f (•, λ * ).

The cross-entropy method

Cross-entropy [START_REF] Kullback | Information Theory and Statistics[END_REF] (alternatively relative entropy or Kullback-Leibler divergence) has been shown to be an effective directed measure of distance between distributions [START_REF] Shore | Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy[END_REF]. With regard to the present context, it has also been shown to be useful in finding optimum distributions for importance sampling [START_REF] Rubinstein | The Cross-Entropy Method for Combinatorial and Continuous Optimization[END_REF][START_REF] De Boer | Adaptive importance sampling simulation of queueing networks[END_REF][START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF]. Given two probability density functions f and f ′ over the same probability space Ω, the cross-entropy from f to f ′ is given by

CE(f, f ′ ) = Ω f (ω) log f (ω) f ′ (ω) dω = Ω f (ω) log f (ω) -f (ω) log f ′ (ω) dω = H(f ) - Ω f (ω) log f ′ (ω) dω (4)
where H(f ) is the entropy of f . To find

λ * we minimise CE(z(ω)f (ω, µ), f (ω, λ)), noting that H(f (ω, µ)) is independent of λ: λ * = arg max λ Ω z(ω)f (ω, µ) log f (ω, λ) dω (5)
Estimating λ * directly using Equation ( 5) is hard, so we re-write it using importance sampling density f (•, λ ′ ) and likelihood ratio function L(ω) = f (ω, µ)/f (ω, λ ′ ):

λ * = arg max λ Ω z(ω)L(ω)f (ω, λ ′ ) log f (ω, λ) dω (6)
Using Equation ( 6) we can construct an unbiased importance sampling estimator of λ * and use it as the basis of an iterative process to obtain successively better estimates:

λ * = λ (j+1) = arg max λ Nj i=1 z(ω (j) i )L (j) (ω (j) i ) log f (ω (j) i , λ) (7)
N j is the number of simulation runs on the j th iteration, λ (j) is the j th set of estimated parameters,

L (j) (ω) = f (ω, µ)/f (ω, λ (j) ) is the j th likelihood ratio function, ω (j) i
is the i th path generated using f (•, λ (j) ) and f (ω

(j)
i , λ) is the probability of path ω (j) i under the distribution f (•, λ (j) ).

A parametrised cross-entropy algorithm

We consider a system of n guarded commands with vector of rate functions K = (K 1 , . . . , K n ) and corresponding vector of parameters λ = (λ 1 , . . . , λ n ). In any given state the probability that command k is chosen is given by

λ k K k K, λ
where the notation •, • denotes a scalar product. For the purposes of simulation we consider a space of finite paths ω ∈ Ω. Let U k (ω) be the number of transitions of type k occurring in ω. We therefore have

f (ω, λ) = n k (λ k ) U k (ω) U k (ω) s=1 K s k K s , λ
Here vector K is indexed by s to emphasise its state-dependence. The likelihood ratios are thus of the form

L (j) (ω) = n k µ k λ (j) k U k (ω) U k (ω) s=1 K s , λ (j) K s , µ
We substitute these expressions in the cross-entropy estimator Equation ( 7) and for compactness substitute

z i = z(ω i ), u i (k) = U k (ω i ) and l i = L (j) (ω i ) to get arg max λ N i=1 l i z i log n k   λ ui(k) k ui(k) s=1 K i,s k K i,s , λ   (8)
= arg max

λ N i=1 n k l i z i u i (k)   log(λ k ) + ui(k) s=1 log(K i,s k ) - ui(k) s=1 log( K i,s , λ )  
We partially differentiate with respect to λ k and get the non-linear system

∂f ∂λ k (λ) = 0 ⇔ N i=1 l i z i   u i (k) λ k - |ωi| s=1 K i,s k K i,s , λ   = 0 (9)
where |ω i | is the length of the path ω i .

Theorem 1. A solution of Equation ( 9) is almost surely a maximum, up to a normalising scalar.

Proof. Using a standard result, it is sufficient to show that the Hessian matrix in λ is negative semi-definite. Consider f i :

f i (λ) = k u i (k)   log(λ k ) + ui(k) s=1 log(K i,s k ) - ui(k) s=1 log( K i,s , λ )  
The Hessian matrix in λ is of the following form with v

(s) k = K (s) k K (s) ,λ and v k = (v (s) k ) 1≤s≤U k (ω) : H i = G -D
where G = (g kk ′ ) 1≤k,k ′ ≤n is the following Gram matrix

g kk ′ = v k , v k ′
and D is a diagonal matrix such that

d kk = u k λ 2 k .
Note that asymptotically

d kk = 1 λ k N s=1 v (s)
k . We write ½ N = (1, . . . , 1) for the vector of N elements 1, hence

d kk = 1 λ k v k , ½ N .
Furthermore, ∀s,

n k=1 λ k v (s) k = 1. So, n k ′ =1 λ k ′ v k ′ = ½ N . Finally, d kk = n k ′ =1 λ k ′ λ k v k , v k ′ . Let x ∈ R n .
To prove the theorem we need to show that -x t Hx ≥ 0.

-x t Hx = x t Dx -x t Gx (10) = k,k ′ λ k ′ λ k v k , v k ′ x 2 k - k,k ′ v k , v k ′ x k x k ′ = k<k ′ λ k ′ λ k x 2 k + λ k λ k ′ x 2 k ′ -2x k x k ′ v k , v k ′ = k<k ′ λ k ′ λ k x k - λ k λ k ′ x k ′ 2 v k , v k ′ ≥ 0
The Hessian matrix H of f is of the general form

H = N i=1 l i z i H i
which is a positively weighted sum of non-positive matrices.

⊓ ⊔

The Hessian is negative semi-definite because if λ is a solution then xλ, x ∈ R + , is also a solution. The fact that there is a unique optimum, however, makes it conceivable to find λ * using standard optimising techniques such as Newton and quasi-Newton methods. To do so would require introducing a suitable normalising constraint in order to force the Hessian to be negative definite. In the case of the cross-entropy algorithm of [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF], this constraint is inherent because it works at the level of individual transition probabilities that sum to 1 in each state.

We note here that in the case that our parameters apply to individual transitions, such that one parameter corresponds to exactly one transition, Equation ( 12) may be transformed to Equation ( 9) of [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] by constraining K, λ = 1. Equation ( 9) of [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] has been shown in [START_REF] Ridder | Asymptotic optimality of the cross-entropy method for markov chain problems[END_REF] to converge to f * , implying that under these circumstances f (•, λ * ) = f * and that it may be possible to improve our parametrised importance sampling distribution by increasing the number of parameters.

The algorithm

Equation ( 9) leads to the following expression for λ k :

λ k = N i=1 l i z i u i (k) N i=1 l i z i |ωi| s=1 K s k K s ,λ (11)
In this form the expression is not useful because the right hand side is dependent on λ k in the scalar product. Hence, in contrast to update formulae based on unbiased estimators, as given by Equation ( 7) and in [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF][START_REF] De Boer | Adaptive importance sampling simulation of queueing networks[END_REF], we construct an iterative process based on a biased estimator but having a fixed point that is the optimum:

λ (j+1) k = Nj i=1 l i z i u i (k) Nj i=1 l i z i |ωi| s=1 K s k K s ,λ (j) (12)
Equation ( 12) can be seen as an implementation of Equation ( 11) that uses the previous estimate of λ in the scalar product, however it works by reducing the distance between successive distributions, rather than by explicitly reducing the distance from the optimum. We offer no proof of convergence here, but assert that if it converges it converges to λ * .

To use the algorithm it is necessary to start with an initial simulation distribution f (•, λ (0) ) and number of simulations N 0 that produce at least a few traces that satisfy the property. The choice of λ (0) and N 0 is highly dependent on the model and the property and must in general be established by heuristics or trial and error. When the number of parameters is small and the property is very rare, an effective strategy is to iterate the algorithm with relatively low N 0 and random parameters until a suitable trace is observed. If the model and property are similar to a previous combination for which parameters were found, those parameters are likely to provide a good initial estimate. Increasing the parameters associated to obviously small rates may help (along the lines of simple failure biasing [START_REF] Shahabuddin | Importance Sampling for the Simulation of Highly Reliable Markovian Systems[END_REF]), however the rareness of a property expressed in temporal logic may not always be related to low probabilities. Finding good initial distributions for arbitrary systems and temporal properties is the subject of ongoing work.

Given a sufficient number of traces from the first iteration, Equation 12should provide a better set of parameters. The expected behaviour is that on successive iterations the number of traces that satisfy the property increases, however it is important to note that the algorithm optimises the quality of the distribution and that the number of traces that satisfy the property is merely emergent of that. As has been noted, in general f (•, λ * ) = f * , hence it is likely that fewer than 100% of traces will satisfy the property when simulating under f (•, λ * ). One consequence of this is that an initial set of parameters may produce more traces that satisfy the property than the final set (see, e.g., Figure 4).

It is conceivable that certain guarded commands play no part in traces that satisfy the property, in which case Equation ( 12) would make the corresponding parameter zero with no adverse effects. It is also conceivable that an important command is not seen on a particular iteration, but making its parameter zero would prevent it being seen on any subsequent iteration. To avoid this it is necessary to adopt a 'smoothing' strategy [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] that reduces the significance of an unseen command without setting it to zero. The strategy adopted for the examples shown below is to divide the parameter of unseen commands by two. The effects of this can be seen in Figure 6. An alternative approach is to add a small fraction of the initial (or previous) parameters to every new parameter estimate. Whatever the strategy, since the parameters are unconstrained it is advisable to normalise them after each iteration (i.e., k λ k = const.) in order to judge progress.

Once the parameters have converged it is then possible to perform a final set of simulations to estimate the probability of the rare property. The usual assumption is that N j ≪ N IS ≪ N MC , however it is often the case that parameters converge fast, so it is expedient to use some of the simulation runs generated during the course of the optimisation as part of the final estimation.

Examples

The following examples are included to illustrate the performance of our algorithm and parametrisation. The first is an example of a chemical system, often used to motivate stochastic simulation, while the second is a standard repair model. All simulations were performed using our statistical model checking platform, PLASMA [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF].

Chemical network

Following the success of the human genome project, with vast repositories of biological pathway data available online, there is an increasing expectation that formal methods can be applied to biological systems. The network of chemical reactions given below is abstract but typical of biochemical systems and demonstrates the potential of statistical model checking to handle the enormous state spaces of biological models. In particular, we demonstrate the efficacy of our algorithm by applying it to quantify two rare dynamical properties of the system.

We consider a well stirred chemically reacting system comprising five reactants (molecules of type A, B, C, D, E), a dimerisation reaction [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF] and two decay reactions [START_REF] Sumit | A Bayesian Approach to Model Checking Biological Systems[END_REF][START_REF] Kahn | Stochastic (Monte Carlo) Attenuation Analysis[END_REF]:

A + B k1 → C (13)
C k2 → D (14)
D k3 → E (15)
Under the assumption that the molecules move randomly and that elastic collisions significantly outnumber unreactive, inelastic collisions, the system may be simulated using mass action kinetics as a continuous time Markov chain [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]. The semantics of Equation ( 13) is that if a molecule of type A encounters a molecule of type B they will combine to form a molecule of type C after a delay drawn from an exponential distribution with mean k 1 . The decay reactions have the semantics that a molecule of type C (D) spontaneously decays to a molecule of type D (E) after a delay drawn from an exponential distribution with mean k 2 (k 3 ). The reactions [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Sumit | A Bayesian Approach to Model Checking Biological Systems[END_REF][START_REF] Kahn | Stochastic (Monte Carlo) Attenuation Analysis[END_REF] are modelled by three guarded commands having importance sampling parameters λ 1 , λ 2 and λ 3 , respectively. A typical simulation run is illustrated in Figure 1, where the x-axis is steps rather than time to aid clarity. A and B combine rapidly to form C which peaks before decaying slowly to D. The production of D also peaks while E rises monotonically.

(i) Pr[♦ C ≥ x] (ii) Pr[♦ D ≥ y]
With an initial vector of molecules (1000, 1000, 0, 0, 0), corresponding to types (A,B,C,D,E), the state space contains ∼ 10 15 states. We know from a static analysis of the reactions that it is possible for the numbers of molecules of C and D to reach the initial number of A and B molecules (i.e., 1000) and that this is unlikely. To find out exactly how unlikely we consider the probabilities of the following rare properties defined in linear temporal logic: (i) ♦ C ≥ x, x ∈ {970, 975, 980, 985, 990, 995} and (ii) ♦ D ≥ y, y ∈ {460, 465, 470, 475, 480, 485}. The results are plotted in Figure 2.

To establish an initial distribution our algorithm (Equation ( 12)) was iterated with random parameters and N 0 = 1000 until one or more traces were observed that satisfied the property in question. No more than 10 such iterations were needed in any case. The algorithm was then iterated 20 times using N j = 1000. Despite the large state space, this value of N j was found to be sufficient to produce reliable results. Starting from randomly chosen initial values, the convergence of parameters can be seen in Figure 3. Figure 4 illustrates that the number of paths satisfying a property can actually decrease as the quality of the distribution improves. Figure 5 illustrates the convergence of the estimate and sample variance using the importance sampling parameters generated during the course of running the algorithm. The initial set of parameters appear to give a very low variance, however this is clearly erroneous with respect to subsequent values. Noting that the variance of standard Monte Carlo simulation of rare events gives a variance approximately equal to the probability and assuming that the sample variance is close to the true variance, Figure 5 suggests that we have made a variance reduction of approximately 10 7 .

Repair model

To a large extent the need to certify system reliability motivates the use of formal methods and thus reliability models are studied extensively in the literature. The following example is taken from [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] and features a moderately large state space of 40,320 states that can be investigated using numerical methods to corroborate our results. The system is modelled as a continuous time Markov chain and comprises six types of subsystems (1, . . . , 6) containing, respectively, (5, 4, 6, 3, 7, 5) components that may fail independently. The system's evolution begins with no failures and with various probabilistic rates the components fail and are repaired. The failure rates are (2.5ǫ, ǫ, 5ǫ, 3ǫ, ǫ, 5ǫ), ǫ = 0.001, and the repair rates are (1.0, 1.5, 1.0, 2.0, 1.0, 1.5), respectively. Each subsystem type is modelled by two guarded commands: one for failure and one for repair. The property under investigation is the probability of a complete failure of a subsystem (i.e., the failure of all components of one type), given an initial condition of no failures. This can be expressed in temporal logic as Pr[X(¬init U failure)].

Figure 6 shows the convergence of parameters (dashed/solid lines) and highlights the effects of the adopted smoothing strategy (circles). Parameters λ 2 and λ 4 (the parameters for repair commands of types 1 and 2, respectively) are attenuated from the outset by the convergence of the other parameters (because of the normalisation). Once their values are small relative to the normalisation constant (12 in this case), their corresponding commands no longer occur and their values experience exponential decay as a result of smoothing (division by two at every subsequent step). Parameters λ 6 and λ 10 (the parameters for repair commands of types 3 and 5, respectively) converge for 12 steps but then also decay. The parameters for the repair commands of types 4 and 6 (solid lines) are the smallest of the parameters that converge. The fact that the repair transitions are made less likely by the algorithm agrees with the intuition that we are interested in direct paths to failure. The fact that they are not necessarily made zero reinforces the point that the algorithm seeks to consider all paths to failure, including those that have intermediate repairs.

Figure 7 plots the number of paths satisfying X(¬init U failure) and suggests that for this model the parametrised distribution is close to the optimum. Figure 8 plots the estimated probability and sample variance during the course of the algorithm and superimposes the true probability calculated by PRISM [26]. The long term average agrees well with the true value (an error of -1.7%, based on an average excluding the first two estimates), justifying our use of the sample variance as an indication of the efficacy of the algorithm: our importance sampling parameters provide a variance reduction of more than 10 5 .

Conclusions and future work

Statistical model checking addresses the state space explosion associated with exact probabilistic model checking by estimating the parameters of an empirical distribution of executions of a system. By constructing an executable model, rather than an explicit representation of the state space, SMC is able to quantify and verify the performance of systems that are intractable to an exhaustive approach. SMC trades certainty for tractability and often offers the only feasible means to certify real-world systems. Rare properties pose a particular problem to Monte Carlo simulation methods because the properties are difficult to observe and the error in their estimated probabilities is difficult to bound. Importance sampling is a well-established means to reduce the variance of rare events but requires the construction of a suitable importance sampling distribution without resorting to the exploration of the entire state space.

We have devised a natural parametrisation for importance sampling and have provided a simple algorithm, based on cross-entropy minimisation, to optimise the parameters for use in statistical model checking. We have shown that our parametrisation leads to a unique optimum and have demonstrated that with very few parameters our algorithm can make significant improvements in the efficiency of statistical model checking. We have shown that our approach is applicable to standard reliability models and to the kind of huge state space models found in systems biology. We therefore anticipate that our methodology has the potential to be applied to many complex natural and man-made systems.

An ongoing challenge is to find ways to accurately bound the error of results obtained by importance sampling. Specifically, the sample variance of the results may be a very poor indicator of the true variance (i.e. with respect to the unknown true probability). Recent work has addressed this problem using Markov chain coupling applied to a restricted class of models and logic [START_REF] Barbot | Coupling and Importance Sampling for Statistical Model Checking[END_REF], but a simple universal solution remains elusive. A related challenge is to find precise means to judge the quality of the importance sampling distributions we create. Our algorithm finds an optimum based on an automatic parametrisation of a model described in terms of guarded commands. This description is usually derived from a higher level syntactic description that is likely optimised for compactness rather than consideration of importance sampling. As such, there may be alternative, equivalent ways of describing the model that produce better importance sampling distributions. Applying existing work on the robustness of estimators (see, e.g., Chapter 4 in [START_REF]Rare Event Simulation using Monte Carlo Methods[END_REF]), we hope to adapt our algorithm to provide hints about improved parametrisation.

Fig. 1 .

 1 Fig. 1. A typical stochastic simulation trace of reactions (13-15).

3 Fig. 3 .Fig. 4 .

 334 Fig. 3. Convergence of parameters for ♦ D ≥ 470 in the chemical model using Nj = 1000.

Fig. 5 . 12 Fig. 6 .

 5126 Fig. 5. Convergence of probability and sample variance for ♦ D ≥ 470 in the chemical model using Nj = 1000.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Convergence of number of paths satisfying X(¬init U failure) in the repair model using Nj = 10000.


S Méance

L Pages

A.-M Gué

P Joseph

PRODUCTION OF 3-10 µm MICROBUBBLES SUITED FOR ULTRASONIC IMAGING BY 2.5D NANOFLUIDIC DEVICES

Keywords: Ultrasound contrast agents, Nanochannels, Nanofluidics

In this paper, we propose a novel route to generate monodisperse microbubbles (MB) in the range size of 3-10 µm. Using deep reactive ion etching, we designed several "terraces" emulsification devices with depths of 200-1000 nm. We performed characterizations and systematical investigations on critical parameters such as gas pressure, liquid flow rate, and chip geometry. We experimentally show here that MB size is mainly dependent on nanochannel height. The robust monodisperse MB production and the biocompatibility of the materials opens an opportunity to integrate such systems in the future diagnostics routines.

INTRODUCTION

For their use in diagnostic echocardiography, ultrasound contrast agents (UCA) have been widely studied for several decades. Injection of UCA MB through the circulatory system imposes a size limitation of about 8 µm, that does not exceed that of pulmonary capillaries. Conversely the ultrasound scattering efficiency is a function of the diameter MB. Hence, commercially produced UCA have a mean size of 2-6 µm. They also have a relatively narrow distribution but they are still polydisperse [START_REF] Lindner | Microbubbles in medical imaging: current applications and future directions[END_REF]. In order to improve the sensitivity of the imaging systems, Hettiarachchi et al. proposed a microfluidic flowfocusing chip for monodisperse UCA generation [START_REF] Hettiarachchi | On-chip generation of microbubbles as a practical technology for manufacturing contrast agents for ultrasonic imaging[END_REF]. The drawback of this approach is that the device depends on a critical dimension: a narrow orifice that can clog with PDMS debris accumulation. Another method based on 2.5D geometries has recently been proposed by Stoffel et al. [START_REF] Stoffel | Bubble production mechanism in a microfluidic foam generator[END_REF]. It is driven by surface tension and relies on the shape instability at a confined-unconfined junction. This mechanism is particularly robust independently on flow rates, and gas pressure. However the MB diameter of 12-27 µm measured in their work was not in the range required for use as ultrasound contrast agents. In this work we experimentally demonstrate the MB generation with size that fits the requirement of contrast echography.

EXPERIMENTAL

The devices were realized by standard microfabrication. Nanochannel (height h between 200 nm and 1 µm, lateral dimension 10 µm x 60 µm) and microchannels (20 µm x 100 µm) were etched in silicon by dry etching. Holes were sand blasted; anodic bonding ensured sealing to a glass substrate. Gas (nitrogen) pressure P g was imposed using a pressure controller. The liquid (deionized water + 2% Tween20) was injected at a flow rate Q w using a syringe pump. Chip design is shown in Figure 1.

RESULTS AND DISCUSSION

In the nanochannel, the curvature introduces a pressure jump, known as the Laplace pressure, between the inside and the outside of a gas-liquid interface. Assuming full wetting, in this case it is ΔP = 2γ/h where γ is the surface tension. When P g > ΔP, the water-gas meniscus advances to the terrace. Once the interface is pinned to the terrace edge, an increase in pressure creates a meniscus unbalanced, causing an overflow of the gas in the microchannel, followed by pinch-off and formation of a new bubble.

Figure 2 shows MB generation for three different gas pressures. When P g is slightly higher than ΔP, only a few MB are formed. Above this pressure, the system is stable and production is monodisperse (Figure 2-b). For P g >5400 mbar, the frequency is too high and the MB are facing an environment increasingly confined. They start to be polydisperse (Figure 2-c).

Figure 3 summarizes all the parameters (gas pressure, flow rate, nanochannel height) involved in this experiment. The first observation relates to the influence of the pressure on the MB size. Indeed when the pressure increases, the diameter increases to its maximum and then decreases. This is consistent with the semi empirical physical model of [START_REF] Stoffel | Bubble production mechanism in a microfluidic foam generator[END_REF] but could also be due to the deformation of the MB in the confined environment. The second point concerns the influence of flowrate Q w . When it increases, the diameter of the MB decreases. In such situation, the geometry can be compared with a T-junction microfluidic device. Then, the higher the flowrate, the greater the shear forces are important to the terrace level causing detachment of the MB. Nevertheless the variations remain weak and the MB size has little dependency on the gas pressure and the liquid flowrate. Finally, the most important parameter is the nanochannel height. As we can observe the size of MB formed by terrace emulsification is roughly proportional to the confining dimension (around 10 times the nanochannel depth).

CONCLUSION

In agreement with existing literature, we show that MB size have little dependency on production frequency. Moreover, we experimentally show that the nanochannel depth is a critical parameter to reach the desired size. Finally, we have implemented a terrace emulsification chip for manufacturing monodisperse MB in the size range desired for ultrasonic imaging.

Figure 1 :

 1 Figure 1: Schematic of the terrace emulsification device.

Figure 2 :

 2 Figure 2: Micrographs of the microbubbles produced at the interface in an h = 400 nm device with Q w =20 µL/min. Experimental conditions: a) P g =4.6 bar; b) P g =5.1 bar; c) P g =5.4 bar.

Figure 3 :

 3 Figure 3: Microbubbles size diameter as a function of P g , Q w , and h.

ACKNOWLEDGEMENTS

We acknowledge CNRS and French National Research Agency for funding (Smart-US program).


Uli Fahrenberg

Kim G Larsen

Axel Legay

Louis-Marie Traonouez

Parametric and Quantitative Extensions of Modal Transition Systems

Modal transition systems provide a behavioral and compositional specification formalism for reactive systems. We survey two extensions of modal transition systems: parametric modal transition systems for specifications with parameters, and weighted modal transition systems for quantitative specifications.

Introduction

Modal transition systems [START_REF] Kim | Modal specifications. In Automatic Verification Methods for Finite State Systems[END_REF][START_REF] Kim | A modal process logic[END_REF] provide a behavioral and compositional specification formalism for reactive systems. They grew out of the notion of relativized bisimulation [START_REF] Kim | A context dependent equivalence between processes[END_REF], which allows for simple specifications of components by allowing the notion of bisimulation to take into account the restricted use that a given component may have in its context.

A modal transition system is essentially a (labeled) transition system, but with two types of transitions: so-called may-transitions which any implementation may (or may not) have, and must-transitions which any implementation is required to have. In fact, ordinary labeled transition systems (or implementations) are modal transition systems where the set of may-and must-transitions coincide. Modal transition systems come equipped with a bisimulation-like notion of (modal) refinement, reflecting that the more must-transitions and the fewer may-transitions a modal specification has the more refined and closer to a final implementation it is.

Example 1. Consider the modal transition system shown in Fig. 1 which models the requirements of a simple email system in which emails are first received and then delivered; must-and may-transitions are represented by solid and dashed arrows, respectively. Before delivering the email, the system may check or process the email, e.g. for encryption or decryption, filtering of spam emails, or generating automatic answers using an auto-reply feature. Any implementation of this email system specification must be able to receive and deliver email, and it may also be able to check arriving email before delivering it. No other behavior is allowed. Such a valid implementation is given in Fig. 2.

The theory of modal transition systems (MTS), or modal specifications as they were called in the paper [START_REF] Kim | Modal specifications. In Automatic Verification Methods for Finite State Systems[END_REF] in the proceedings of the first CAV conference receive deliver check deliver Fig. 1: Modal transition system modeling a simple email system, with an optional behavior: Once an email is received, it may be checked, e.g. be scanned for containing viruses, or automatically decrypted, before it is delivered to the receiver. organized by Joseph Sifakis in Grenoble, 3 was aiming at providing a behavioral compositional specification formalism for reactive systems. At the time of the introduction of MTS, there were two predominant approaches to specifications formalisms and verification methods for reactive and concurrent systems: logical approaches where a specification is a set of properties of implementations (labeled transition systems), and graphical approaches promoted by the various process algebras, where implementations and specifications are systems of the same kind -namely labeled transition systems, and verification amounts to compare such systems with respect to a given behavioral preorder, e.g. bisimilarity.

In search for a complete specification theory, the following properties have been considered desirable (the first three were listed in the early paper [START_REF] Boudol | Graphical versus logical specifications[END_REF]):

expressiveness: the specification formalism should be powerful enough to express all properties of a given implementation. In other words it should be possible to completely specify any labeled transition system, up to bisimulation. modularity: implementations are often made out of several components, and it should be possible to infer satisfaction of an overall specification solely on the basis of sub-specification of the sub-components. refinement: one should have the ability to deal with partial specifications, requiring more and more properties about a system, up to its complete specification.

logical composition: specification should be composable with respect to usual logical operators such as conjunction and (possibly) disjunction. quotienting: given an overall specification S of a composite systems as well as a sub-specification T of a sub-component, the existence of a quotient specification S\T will describe the sufficient and necessary condition of the remaining components in order that S is satisfied by the total systems.

Applying these criteria to the logical and graphical (i.e. bisimulation) framework, as was done in [START_REF] Boudol | Graphical versus logical specifications[END_REF], we see that the logical and graphical frameworks offer complementary advantages: on the graphical side, expressiveness is trivial since a process i a specification of itself. Modularity is usually guaranteed by the fact that bisimulations are compatible with (most) process constructors. On the logical side, expressiveness is achieved if we allow possibly infinite sets of formulae as logical specifications, or admit recursively specified properties. The point of modularity has proved more difficult with early attempts of Sifakis and Graf [START_REF] Graf | A logic for the description of non-deterministic programs and their properties[END_REF] and Holmstrøm [START_REF] Holmström | A refinement calculus for specifications in Hennessy-Milner logic with recursion[END_REF] providing sound and highly usable proof systems for specifications mixing logical and behavioral constructs (as well as fix-point constructs) but lacking accompanying completeness results. Much later the work of Mardare and Policriti [START_REF] Mardare | A complete axiomatic system for a processbased spatial logic[END_REF] provided a first matching completeness result.

In the rest of this paper, we survey two extensions of modal transition systems. The first extension, parametric modal transition systems, is concerned with systems whose behaviors depend on parameters [START_REF] Beneš | Parametric modal transition systems[END_REF]. The second extension, weighted modal transition systems [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF][START_REF] Bauer | General quantitative specification theories with modalities[END_REF] permits to reason on systems whose behaviors depend on quantities. Another paper in this volume [START_REF] Fahrenberg | Specification theories for probabilistic and real-time systems[END_REF] will be concerned with other extensions of modal transition systems which are more closely related to applications.

Acknowledgment. This survey paper presents research which we have conducted with a number of coauthors; in alphabetical order, these are Sebastian S. Bauer, Nikola Beneš, Line Juhl, Jan Křetínský, Mikael H. Møller, Jiří Srba, and Claus Thrane. We acknowledge their cooperation in this work; any errors in this presentation are, however, our own.

Parametric Modal Transition Systems

It is well admitted (see e.g. [START_REF] Raclet | Why are modalities good for interface theories?[END_REF]) that MTS and their extensions like disjunctive MTS (DMTS) [START_REF] Kim | Equation solving using modal transition systems[END_REF], 1-selecting MTS (1MTS) [START_REF] Fecher | Comparing disjunctive modal transition systems with an one-selecting variant[END_REF] and transition systems with obligations (OTS) [START_REF] Beneš | Process algebra for modal transition systemses[END_REF] provide strong support for a specification formalism allowing for step-wise refinement process. Moreover, the MTS formalisms have applications in other contexts, which include verification of product lines [START_REF] Gruler | Modeling and model checking software product lines[END_REF][START_REF] Kim | On modal refinement and consistency[END_REF], interface theories [START_REF] Raclet | Why are modalities good for interface theories?[END_REF][START_REF] Uchitel | Merging partial behavioural models[END_REF] and modal abstractions in program analysis [START_REF] Godefroid | Abstraction-based model checking using modal transition systems[END_REF][START_REF] Huth | Modal transition systems: A foundation for three-valued program analysis[END_REF][START_REF] Nanz | Modal abstractions of concurrent behaviour[END_REF].

Unfortunately, all of these formalisms lack the capability to express some intuitive specification requirements like exclusive, conditional and persistent choices. In [START_REF] Beneš | Parametric modal transition systems[END_REF] the expressive power of MTS and its variants has been extended considerably so it can model model arbitrary Boolean conditions on transitions and also allows to instantiate persistent transitions. The model, called parametric modal transition systems (PMTS), is equipped with a finite set of parameters that are fixed prior to the instantiation of the transitions in the specification. The generalized notion of modal refinement is designed to handle the parametric extension and it specializes to the well-studied modal refinements on all the subclasses of our model like MTS, disjunctive MTS and MTS with obligations.

Motivation

We shall now discuss these limitations on an example as a motivation for the introduction of parametric MTS formalism with general Boolean conditions in specification requirements.

Consider a simple specification of a traffic light controller that can be at any moment in one of the four predefined states: red , green, yellow or yellowRed . The requirements of the specification are: when green is on the traffic light may either change to red or yellow and if it turned yellow it must go to red afterward; when red is on it may either turn to green or yellowRed , and if it turns yellowRed (as it is the case in some countries) it must go to green afterwords.

Fig. 3a shows an obvious MTS specification of the proposed specification. The transitions in the standard MTS formalism are either of type may (optional transitions depicted as dashed lines) or must (required transitions depicted as solid lines). In Fig. 3c, Fig. 3d and Fig. 3e we present three different implementations of the MTS specification where there are no more optional transitions. The implementation I 1 does not implement any may transition as it is a valid possibility to satisfy the specification S 1 . Of course, in our concrete example, this means that the light is constantly green and it is clearly an undesirable behavior that cannot be, however, easily avoided. The second implementation I 2 on the other hand implements all may transitions, again a legal implementation in the MTS methodology but not a desirable implementation of a traffic light as the next action is not always deterministically given. Finally, the implementation I 3 of S 1 illustrates the third problem with the MTS specifications, namely that the choices made in each turn are not persistent and the implementation alternates between entering yellow or not. None of these problems can be avoided when using the MTS formalism.

A more expressive formalism of disjunctive modal transition systems (DMTS) can overcome some of the above mentioned problems. A possible DMTS specification S 2 is depicted in Fig. 3b. Here the ready and stop transitions, as well as ready and go ones, are disjunctive, meaning that it is still optional which one is implemented but at least one of them must be present. Now the system I 1 in Fig. 3c is not a valid implementation of S 2 any more. Nevertheless, the undesirable implementations I 2 and I 3 are still possible and the modeling power of DMTS is insufficient to eliminate them.

Inspired by the recent notion of transition systems with obligations [START_REF] Beneš | Process algebra for modal transition systemses[END_REF], we can model the traffic light using specification as a transition system with (g) PMTS specification S4 Fig. 3: Specifications and implementations of a traffic light controller arbitrary4 obligation formulae. These formulae are Boolean propositions over the outgoing transitions from each state, whose satisfying assignments yield the allowed combinations of outgoing transitions. A possible specification called S 3 is given in Fig. 3f and it uses the operation of exclusive-or. We will follow an agreement that whenever the obligation function for some node is not listed in the system description then it is implicitly understood as requiring all the available outgoing transitions to be present. Due to the use of exclusive-or in the obligation function, the transition systems I 1 and I 2 are not valid implementation any more. Nevertheless, the implementation I 3 in Fig. 3e cannot be avoided in this formalism either.

Finally, the problem with the alternating implementation I 3 is that we cannot enforce in any of the above mentioned formalisms a uniform (persistent) implementation of the same transitions in all its states. In order to overcome this problem, we propose the so-called parametric MTS where we can, moreover, choose persistently whether the transition to yellow is present or not via the use of parameters. The PMTS specification with two parameters reqYfromR and reqYfromG is shown in Fig. 3g. Fixing a priori the (Boolean) values of the parameters makes the choices permanent in the whole implementation, hence we eliminate also the last problematic implementation I 3 .

Definition

We shall now formally capture the intuition behind parametric MTS introduced above. First, we recall the standard propositional logic.

A Boolean formula over a set X of atomic propositions is given by the following abstract syntax

ϕ ::= tt | x | ¬ϕ | ϕ ∧ ψ | ϕ ∨ ψ
where x ranges over X. The set of all Boolean formulae over the set X is denoted by B(X). Let ν ⊆ X be a truth assignment, i.e. a set of variables with value true, then the satisfaction relation ν |= ϕ is given by ν |= tt, ν |= x iff x ∈ ν, and the satisfaction of the remaining Boolean connectives is defined in the standard way. We also use the standard derived operators like exclusive-or ϕ ⊕ ψ = (ϕ ∧ ¬ψ) ∨ (¬ϕ ∧ ψ), implication ϕ ⇒ ψ = ¬ϕ ∨ ψ and equivalence ϕ ⇔ ψ = (¬ϕ ∨ ψ) ∧ (ϕ ∨ ¬ψ).

We can now proceed with the definition of parametric MTS.

Definition 1. A parametric MTS (PMTS) over an action alphabet Σ is a tuple (S, T, P, Φ) where S is a set of states, T ⊆ S × Σ × S is a transition relation, P is a finite set of parameters, and Φ : S → B((Σ × S) ∪ P ) is an obligation function over the atomic propositions containing outgoing transitions and parameters. We implicitly assume that whenever (a, t) ∈ Φ(s) then (s, a, t) ∈ T . By T (s) = {(a, t) | (s, a, t) ∈ T } we denote the set of all outgoing transitions of s.

PMTS has been provided a refinement notion that generalizes the well-studied refinement notions on its subclasses including that of MTS. In the definition, the parameters are fixed first (persistence) followed by all valid choices modulo the fixed parameters that now behave as constants.

First we set the following notation. Let (S, T, P, Φ) be a PMTS and ν ⊆ P be a truth assignment. For s ∈ S, we denote by Tran ν (s) = {E ⊆ T (s) | E ∪ν |= Φ(s)} the set of all admissible sets of transitions from s under the fixed truth values of the parameters.

We can now define the notion of modal refinement between PMTS.

∀M ∈ Tran µ (s) : ∃N ∈ Tran ν (t) : ∀(a, s ′ ) ∈ M : ∃(a, t ′ ) ∈ N : (s ′ , t ′ ) ∈ R ∧ ∀(a, t ′ ) ∈ N : ∃(a, s ′ ) ∈ M : (s ′ , t ′ ) ∈ R .
We say that s modally refines t, denoted by s ≤ m t, if there exists a modal refinement R such that (s, t) ∈ R.

Example 2. Consider the rightmost PMTS in Fig. 4. It has two parameters reqYfromG and reqYfromR whose values can be set independently and it can be refined by the system in the middle of the figure having only one parameter reqY . This single parameter simply binds the two original parameters to the same value.

The PMTS in the middle can be further refined into the implementations where either yellow is always used in both cases, or never at all. Notice that there are in principle infinitely many implementations of the system in the middle, however, they are all bisimilar to either of the two implementations depicted in the left of Fig. 4.

[4] provides an extensive study of the complexity of refinement checking between parametric modal transitions with classification depending on the complexity of obligations as well as the presence or absence of parameters. For each combination the complexity class of the polynomial hierarchy for which modal refinement is complete is provided. In short, the complexities ranges from P-complete to Π p 4 -complete (thus in PSPACE).

Quantitative Modal Transition Systems

Motivated by applications to embedded, real-time and hybrid systems, the modal transition system framework has been extended in order to reason about receive, [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF][START_REF] Bauer | Extending modal transition systems with structured labels[END_REF] deliver, [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF][START_REF] Beneš | Parametric modal transition systems[END_REF] check, [0, 5]

deliver, [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF][START_REF] Bauer | General quantitative specification theories with modalities[END_REF] Fig. 5: Specification of a simple email system, similar to Fig. 1, but extended by integer intervals modeling time units for performing the corresponding actions.

quantitative aspects [START_REF] Bauer | Extending modal transition systems with structured labels[END_REF][START_REF] Juhl | Modal transition systems with weight intervals[END_REF]. With these applications in mind, it is necessary not only to be able to specify quantitative aspects of systems, but also to formalize successive refinement of quantities. To illustrate this extension, consider again the modal transition system of Fig. 1, but this time with quantities, see Fig. 5: Every transition label is extended by integer intervals modeling upper and lower bounds on time required for performing the corresponding actions. For instance, the reception of a new email (action receive) must take between one and three time units, the checking of the email (action check ) is allowed to take up to five time units.

In this quantitative setting, there is a problem with using a Boolean notion of refinement as is done in the preceding section: If one only can decide whether or not an implementation refines a specification, then the quantitative aspects get lost in the refinement process. As an example, consider the email system implementations in Fig. 6. Implementation (a) does not refine the specification, as there is an error in the discrete structure of actions: after receiving an email, the system can check it indefinitely without ever delivering it. Also implementations (b) and (c) do not refine the specification: (b) takes too long to receive email, (c) does not deliver email fast enough after checking it. Implementation (d) on the other hand is a perfect refinement of the specification.

Intuitively however, implementations (b) and (c) conform much better to the specification than implementation (a) in Fig. 6: there are no discrepancies in the discrete structure, only the weights are off by 1. Additionally, the quantitative error in implementation (c) occurs later than the one in (b). Hence one may want to say that implementation (d) is in perfect refinement of the specification, (c) is slightly off, (b) is a bit more problematic, whereas implementation (a) is completely unacceptable. A Boolean notion of refinement does not allow to make such distinctions between different negative answers.

To sum up, a Boolean notion of refinement is too fragile for quantitative formalisms. Minor and major modifications in the implementation cannot be distinguished, as both of them may reverse the Boolean answer. As observed e.g. in [START_REF] Luca De Alfaro | Linear and branching system metrics[END_REF], this view is obsolete; engineers need quantitative notions on how modified implementations differ. The introduction of such a quantitative notion of refinement, and its consequences for the specification theory, are the subject of this section, which is based on the papers [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF][START_REF] Bauer | General quantitative specification theories with modalities[END_REF]. Depending on the precise application of our quantitative formalism, there are a few choices which one has to make. One such choice is the precise definition of quantitative refinement, as the way quantitative discrepancies between specifications is measured e.g. depends on whether differences accumulate over time or the interest more lies in the maximal individual differences. Another choice is how to combine quantities during structural composition: when modeling e.g. energy consumption, they should be added; when modeling timing constraints, some form of conjunction should be used.

To facilitate quantitative reasoning on specifications and implementations, we introduce a real-valued distance between specifications such that perfect refinement corresponds to distance 0, small quantitative discrepancies give rise to small distances, and differences in the discrete control structure correspond to distance ∞. For the examples in Figs. 5 and6, we will deduce the following chain of decreasing distances:

∞ = d(I 1 , S) > d(I 2 , S) > d(I 3 , S) > d(I 4 , S) = 0

Weighted modal transition systems

Let Σ be a set of labels with a preorder ⊑ ⊆ Σ ×Σ, and denote by Σ ∞ = Σ * ∪Σ ω the set of finite and infinite traces over Σ. len(σ), for σ ∈ Σ ∞ , denotes the length (finite or infinite) of a trace σ. Let ε ∈ Σ ∞ denote the empty trace, and for a ∈ Σ, σ ∈ Σ ∞ , denote by a.σ their concatenation.

A weighted modal transition system (WMTS) is a tuple S = (S, s 0 , , -→) consisting of a set S of states, an initial state s 0 ∈ S, and must-and maytransitions -→, ⊆ S × Σ × S for which it holds that for all s a -→ s ′ there is

s b s ′ with a ⊑ b.
Intuitively, a may-transition s b t specifies that an implementation I of S is permitted to have a corresponding transition i a -→ j, for any a ⊑ b, whereas a must-transition s b -→ t postulates that I is required to implement at least one corresponding transition i a -→ j for some a ⊑ b. We will make this precise below. An WMTS S = (S, s 0 , , -→) is an implementation if -→ = . Hence in an implementation, all optional behavior has been resolved.

Definition 3. A modal refinement of WMTS S 1 = (S 1 , s 0 1 , 1 , -→ 1 ), S 2 = (S 2 , s 0 2 , 2 , -→ 2 ) is a relation R ⊆ S 1 × S 2 such that for any (s 1 , s 2 ) ∈ R, -whenever s 1 a1 1 t 1 , then also s 2 a2 2 t 2 for some a 1 ⊑ a 2 and (t 1 , t 2 ) ∈ R, -whenever s 2 a2 -→ 2 t 2 , then also s 1 a1 -→ 1 t 1 for some a 1 ⊑ a 2 and (t 1 , t 2 ) ∈ R.
Thus any behavior which is permitted in S 1 is also permitted in S 2 , and any behavior required in S 2 is also required in S 1 . We write

S 1 ≤ m S 2 if there is a modal refinement R ⊆ S 1 × S 2 with (s 0 1 , s 0 2 ) ∈ R.
The implementation semantics of a WMTS S is the set S = {I ≤ m S | I implementation}, and we write

S 1 ≤ t S 1 if S 1 ⊆ S 2 , saying that S 1 thoroughly refines S 2 . It follows by transitivity of ≤ m that S 1 ≤ m S 2 implies S 1 ≤ t S 2
, hence modal refinement is a syntactic over-approximation of thorough refinement.

Distances

Recall that a hemimetric on a set X is a function d : X × X → | ≥0 ∪ {∞} which satisfies d(x, x) = 0 and d(x, y) + d(y, z) ≥ d(x, z) (the triangle inequality) for all x, y, z ∈ X. Note that our hemimetrics are extended in that they can take the value ∞.

We will need to generalize hemimetrics to codomains other than | ≥0 ∪ {∞}. For a partially ordered monoid (▲, ⊑, ⊕, ✵), an ▲-hemimetric on X is a function d : X × X → ▲ which satisfies d(x, x) = ✵ and d(x, y) ⊕ d(y, z) ⊒ d(x, z) for all x, y, z ∈ X. For any set M , let ▲M = (| ≥0 ∪ {∞}) M the set of functions from M to the extended non-negative real line. Then ▲M is a complete lattice with partial order ⊑ ⊆ ▲M × ▲M given by α ⊑ β if and only if α(x) ≤ β(x) for all x ∈ M , and with an addition ⊕ given by (α ⊕ β)(x) = α(x) + β(x). The bottom element of ▲M is also the zero of ⊕ and given by ⊥(x) = 0, and the top element is We have shown in [START_REF] Bauer | General quantitative specification theories with modalities[END_REF][START_REF] Fahrenberg | The quantitative linear-timebranching-time spectrum[END_REF][START_REF] Fahrenberg | Distances for weighted transition systems: Games and properties[END_REF]] that all commonly used trace distances obey a recursive characterization as above. The point-wise distance from [START_REF] Luca De Alfaro | Model checking discounted temporal properties[END_REF], for example,

⊤(x) = ∞.
has ▲ = | ≥0 ∪ {∞}, eval = id and d ▲M m (a.σ, b.τ ) = max(d(a, b), d ▲M m (σ, τ )),
where d : Σ × Σ → | ≥0 ∪ {∞} is a hemimetric on labels. The limit-average distance used in e.g. [START_REF] Černý | Simulation distances[END_REF] has ▲ = (| ≥0 ∪ {∞}) ◆ , the complete lattice of functions

◆ → | ≥0 ∪ {∞}, eval(α) = lim inf j∈◆ α(j) and d ▲M m (a.σ, b.τ )(j) = 1 j+1 d(a, b) + j j+1 d ▲M m (σ, τ ).
For the rest of this section, we fix a recursively specified trace distance. A WMTS (S, s 0 , , -→) is deterministic if it holds for all s ∈ S, s a1 s 1 , s a2 s 2

for which there is a ∈ Σ with td ▲M (a, a 1 ) = ⊤ and td ▲M (a, a 2 ) = ⊤ that a 1 = a 2 and s 1 = s 2 .

Definition 6. The lifted modal refinement distance d ▲M m :

S 1 × S 2 → ▲ between the states of WMTS S 1 = (S 1 , s 0 1 , 1 , -→ 1 ), S 2 = (S 2 , s 0 2 , 2 , -→ 2
) is defined to be the least fixed point to the equations

d ▲M m (s 1 , s 2 ) = max      sup s1 a 1 1 t1 inf s2 a 2 2 t2 F (a 1 , a 2 , d ▲M m (t 1 , t 2 )), sup s2 a 2 -→2t2 inf s1 a 1 -→1t1 F (a 1 , a 2 , d ▲M m (t 1 , t 2 )). We let d ▲M m (S 1 , S 2 ) = d ▲M m (s 0 1 , s 0 2 ). The modal refinement distance is d m = eval • d ▲M m , and we write S 1 ≤ ε m S 2 , for ε ∈ | ≥0 ∪ {∞}, if d ▲M m (S 1 , S 2 ) ≤ ε.
Proposition 1. The modal refinement distance is a well-defined hemimetric, and

S 1 ≤ m S 2 implies S 1 ≤ 0 m S 2 .
The thorough refinement distance between WMTS S 1 , S 2 is

d t (S 1 , S 2 ) = sup I1∈ S1 inf I2∈ S2 d m (I 1 , I 2 ),
and we write

S 1 ≤ ε t S 2 , for ε ∈ | ≥0 ∪ {∞}, if d t (S 1 , S 2 ) ≤ ε. As for the modal distance, d t is a hemimetric, and S 1 ≤ t S 2 implies S 1 ≤ 0 t S 2 . Theorem 1. For all WMTS S 1 , S 2 , d t (S 1 , S 2 ) ≤ d m (S 1 , S 2 ). If S 2 is determin- istic, then d t (S 1 , S 2 ) = d m (S 1 , S 2 ).

Conjunction

Let : Σ ×Σ ֒→ Σ be a commutative partial label conjunction operator for which it holds, for all b 1 , b 2 ∈ Σ, that there is a ∈ Σ for which both td ▲M (a, b 1 ) = ⊤ and td ▲M (a, b 2 ) = ⊤ iff there exists c ∈ Σ for which both b 1 c and b 2 c are defined. This is to relate determinism (left-hand side of the above) to a similar property for label conjunction which is needed in the proof of Theorem 2. Additionally, we assume that is greatest lower bound on labels, i.e.

- In the definition below, we denote by ρ B (S) the pruning of a WMTS S = (S, s 0 , , -→) with respect to the states in a ("bad") subset B ⊆ S, which is obtained as follows: Define a must-predecessor operator pre : 2 S → 2 S by pre(S ′ ) = {s ∈ S | ∃a ∈ Σ, s ′ ∈ S ′ : s a -→ s ′ } and let pre * be the reflexive, transitive closure of pre. Then ρ B (S) is defined if s 0 / ∈ pre * (B), and in that case,

ρ B (S) = (S ρ , s 0 , ρ , -→ ρ ) with S ρ = S \ pre * (B), ρ = ∩ (S ρ × Σ × S ρ ), and -→ ρ = -→ ∩ (S ρ × Σ × S ρ ). Definition 7. The conjunction of two WMTS S 1 = (S 1 , s 0 1 , 1 , -→ 1 ), S 2 = (S 2 , s 0 2 , 2 , -→ 2 ) is the WMTS S 1 ∧ S 2 = ρ B (S 1 × S 2 , (s 0 1 , s 0 2 ), , -→) given as follows (if it exists): s1 a 1 -→1 t1 s2 a 2 2 t2 a1 a2 defined (s1, s2) a 1 a 2 -→ (t1, t2) s1 a 1 1 t1 s2 a 2 -→2 t2 a1 a2 defined (s1, s2) a 1 a 2 -→ (t1, t2) s1 a 1 1 t1 s2 a 2 2 t2 a1 a2 defined (s1, s2) a 1 a 2 (t1, t2) s1 a 1 -→1 t1 ∀s2 a 2 2 t2 : a1 a2 undef. (s1, s2) ∈ B s2 a 2 -→2 t2 ∀s1 a 1 1 t1 : a1 a2 undef. (s1, s2) ∈ B
Note that conjunction of WMTS may give inconsistent states which need to be pruned away after. As seen in the last two SOS rules above, this is the case when one WMTS specifies a must-transition which the other WMTS cannot synchronize with. Here, the demand on implementations of the conjunction would be that they simultaneously must and cannot have a transition, which of course is unsatisfiable. Theorem 2. Let S 1 , S 2 , S 3 be WMTS.

-If S 1 ∧ S 2 is defined, then S 1 ∧ S 2 ≤ m S 1 and S 1 ∧ S 2 ≤ m S 2 .

-If S 1 ≤ m S 2 , S 1 ≤ m S 3 , and S 2 or S 3 is deterministic, then S 2 ∧ S 3 is defined and S 1 ≤ m S 2 ∧ S 2 .

Structural composition

Let ȅ : Σ × Σ ֒→ Σ be a commutative partial label composition operator which specifies which labels can synchronize. Again we need to relate determinism to an analogous property for label composition, hence we require that it holds, for all b 1 , b 2 ∈ Σ, that there is a ∈ Σ for which both d(a, b 1 ) = ⊤ ▲ and d(a, b 2 ) = ⊤ ▲ iff there exists c ∈ Σ for which both b 1 ȅ c and b 2 ȅ c are defined.

Additionally, we assume that there exists a function P : ▲ × ▲ → ▲ which allows us to infer bounds on distances on synchronized labels. We assume that P is monotone in both coordinates, has P (⊥ ▲ , ⊥ ▲ ) = ⊥ ▲ , P (α, ⊤ ▲ ) = P (⊤ ▲ , α) = ⊤ ▲ for all α ∈ ▲, and that F (a 1 ȅ a 2 , b 1 ȅ b 2 , P (α 1 , α 2 )) ⊑ ▲ P (F (a 1 , b 1 , α 1 ), F (a 2 , b 2 , α 2 )) Intuitively, P gives a uniform bound on label composition: distances between composed labels can be bounded above using P and the individual labels' distances, and [START_REF] Bauer | Quantitative refinement for weighted modal transition systems[END_REF] The next theorem shows that structural composition supports quantitative independent implementability: the distance between structural compositions can bounded above using P and the distances between the individual components.

Theorem 3. For all WMTS S 1 , T 1 , S 2 , T 2 with d m (S 1 S 2 , T 1 T 2 ) = ⊤ ▲ , we have d m (S 1 S 2 , T 1 T 2 ) ⊑ ▲ P (d m (S 1 , T 1 ), d m (S 2 , T 2 )).

Fig. 2 :

 2 Fig.2: An implementation of the simple email system in Fig.1in which we explicitly model two distinct types of email pre-processing.

  Obligation function: Φ(green) = (stop, red ) ⊕ (ready, yellow ) Φ(red ) = (go, green) ⊕ (ready, yellowRed ) Parameters: {reqYfromR, reqYfromG} Obligation function: Φ(green) = ((stop, red ) ⊕ (ready, yellow )) ∧(reqYfromG ⇔ (ready, yellow )) Φ(red ) = ((go, green) ⊕ (ready, yellowRed ))∧(reqYfromR ⇔ (ready, yellowRed ))

Fig. 4 :

 4 Fig. 4: Example of modal refinement

Fig. 6 :

 6 Fig.6: Four implementations of the simple email system in Fig.5.

Definition 4 .

 4 A trace distance is a hemimetric td : Σ ∞ × Σ ∞ → | ≥0 ∪ {∞}for which td(a, b) = 0 for all a, b ∈ Σ with a ⊑ b and td(σ, τ ) = ∞ whenever len(σ) = len(τ ).

Definition 5 .

 5 A recursive specification of a trace distance td consists of a set M with a lattice homomorphism eval : ▲M → | ≥0 ∪ {∞}, -an ▲M-hemimetric td ▲M : Σ ∞ × Σ ∞ → ▲M which satisfies td = eval • td ▲M and td ▲M (a, b) = ⊥ for all a, b ∈ Σ with a ⊑ b, and a function F : Σ × Σ × ▲M → ▲M.

F

  must be monotone in the third coordinate and satisfy, for all a, b ∈ Σ andσ, τ ∈ Σ ∞ , that td ▲M (a.σ, b.τ ) = F (a, b, td ▲M (σ, τ )).Note that the definition implies that for all a, b ∈Σ, td ▲M (a, b) = td ▲M (a.ε, b.ε) = F (a, b, td ▲M (ε, ε)) = F (a, b, ⊥).Hence also F (a, a, ⊥) = td ▲M (a, a) = ⊥ for all a ∈ Σ.

  for all a, b ∈ Σ with a b defined, a b ⊑ a and a b ⊑ b; -for all a, b, c ∈ Σ with a ⊑ b and a ⊑ c, b c is defined and a ⊑ b c.

( 1 )

 1 for all a 1 , b 1 , a 2 , b 2 ∈ Σ and α 1 , α 2 ∈ ▲ for which a 1 ȅ a 2 and b 1 ȅ b 2 are defined. Hence d(a 1 ȅ a 2 , b 1 ȅ b 2 ) ⊑ P (d(a 1 , b 1 ), d(a 2 , b 2 )) for all such a 1 , b 1 , a 2 , b 2 ∈ Σ.

  ensures that this bound holds recursively.Definition 8. The structural composition of two WMTS S 1 = (S 1 ,s 0 1 , 1 ,-→ 1 ), S 2 = (S 2 , s 0 2 , 2 , -→ 2 ) is the WMTS S 1 S 2 = (S 1 × S 2 , (s 1 a2 2 t 2 a 1 ȅ a 2 def. (s 1 , s 2 )

			0 , s 2 0 ),	, -→) with
	transitions defined as follows:
	s 1	a1	1 t 1 s 2

a1ȅa2 (t 1 , t 2 ) s 1 a1 -→ 1 t 1 s 2 a2 -→ 2 t 2 a 1 ȅ a 2 def. (s 1 , s 2 ) a1ȅa2 -→ (t 1 , t 2 )

In fact, the first CAV conference was not called CAV, but had the rather lengthy title "Automatic Verification Methods for Finite State Systems."

In the transition systems with obligations only positive Boolean formulae are allowed.