Maxime Cordy

Axel Legay
email: axel.legay@inria.fr

Pierre-Yves Schobbens

Louis-Marie Traonouez
email: louis-marie.traonouez@inria.fr

A framework for the rigorous design of highly adaptive timed systems

Keywords: Software Product Lines, Features, Real-time systems, Model-checking, Timed Games

Adaptive systems can be regarded as a set of static programs and transitions between these programs. These transitions allow the system to adapt its behaviour in response to unexpected changes in its environment. Modelling highly dynamic systems is cumbersome, as these may go through a large number of adaptations. Moreover, often they must also satisfy real-time requirements whereas adaptations may not complete instantaneously. In this paper, we propose to model highly adaptive systems as dynamic real-time software product lines, where software products are able to change their features at runtime. Adaptive features allow one to design systems equipped with runtime reconfiguration capabilities and to model changes in their environment, such has failure modes. We define Featured Timed Game Automata, a formalism that combines adaptive features with discrete and real-time behaviour. We also propose a novel logic to express real-time requirements on adaptive systems, as well as algorithms to check a system against them. We implemented our method as part of PyECDAR, a model checker for timed systems.

I. INTRODUCTION

Computers play a central role in modern life and their errors can have dramatic consequences. Proving the correctness of computer systems is therefore an extremely relevant problem for which quality assurance techniques like model checking and testing provide efficient solutions. Testing consists in applying a finite series of test cases to the system. Although it can detect errors, it cannot guarantee their absence. Another of its limitations is that nowadays, systems are embedded and highly configurable, which makes it hard to specify relevant test cases. Model checking [START_REF] Clarke | Design and synthesis of synchronization skeletons using branching-time temporal logic[END_REF] is an automated technique for verifying systems against functional requirements. The approach relies on an exhaustive verification of a bevahioral model of the system against a property expressed in temporal logic. If the system fails to satisfy the property, then the model checking algorithm provides an example of violation. By nature, model checking guarantees the absence of errors. Albeit it suffers from the so-called state-space explosion, it has been widely used and applied on both academic and industry case studies. Model checking was initially intended for closed and static Boolean systems, but has been extended to target increasingly wider classes of systems, including real-time systems.

The recent advances in computer science pose new challenges to model checking. One of the major difficulties is that today's systems often run in open and potentially unsafe environments, which requires them to adapt their behavior in order to accomplish their tasks reliably. In case of highly evolving environment, these adaptations must be performed as quickly as possible, hence the need for self-adaptive systems. These are harder to verify than static, closed systems of which behaviour and environment are known a priori. Applying model checking to such systems requires to represent all its classes of behaviour as well as its capability to transit between them. Moreover, adaptive systems must satisfy multiple goals which may evolve over time and according to changes in the system or its environment [START_REF] Cheng | Software engineering for Self-Adaptive systems: A research roadmap[END_REF].

One way to model an adaptive system is to view it as a set of static programs and transitions between these programs [START_REF] Zhang | Modular verification of dynamically adaptive systems[END_REF].

When the system has to adapt its behaviour, it triggers a transition to one of its other programs. The drawback of this approach is that all these programs must be modelled and verified individually. This approach has huge costs and is intractable. Another difficulty is the need to verify dynamic properties of adaptive systems. Classical logics cannot express them in a proper way. Alternatives to existing model checking techniques of adaptive systems are thus needed.

The static programs composing an adaptive system likely share commonality, as they also have proper parts. An alternative is to organize the variability between these programs into features, a concept borrowed from software product line engineering (SPLE) [START_REF] Cordy | Model checking adaptive software with featured transition systems[END_REF]. In the latter discipline, a feature is an added functionality that meet a requirement of customers. A product of the line is thus obtained by composing desired features together. In the context of adaptive systems, features model differences between the static programs composing the system. Modifications in its behaviour are therefore triggered by changing its features. We name this process reconfiguration. Features constitute an appropriate modelling artifact to reason on runtime variability. Moreover, transposing this concept to adaptive systems permits to benefit from the formal verification techniques currently developed in SPLE.

The behavior of adaptive systems often rely on real-time requirements such as matching deadline or reacting in real-time to fault. For example, a routing protocol must ensure that a data packet must reach the recipient within a certain amount of time (see more in Section II). Unexpected changes in the environment may prevent the satisfaction of these requirements, hence the need for the system to perform adaptations. The reconfiguration process is not always instantaneous, though. The system may require time to change its features, or can have to delay the reconfiguration until it reaches a stable state.

Unfortunately, most of existing model checking techniques for adaptive systems are not capable of handling such constraints.

In this paper, we propose a formal framework to model and verify adaptive systems that must satisfy evolving realtime requirements. We introduce Featured Timed Game Automata (FTGA), a formalism to represent adaptive behaviour, dynamic environment, and real-time. Our model results from the combination of (1) Adaptive Featured Transition Systems [START_REF] Cordy | Model checking adaptive software with featured transition systems[END_REF], a formalism to model dynamic reconfiguration, and evolving environment, and (2) timed automata [START_REF] Alur | A theory of timed automata[END_REF], an established formalism for real-time systems. The semantics of a FTGA is defined as a timed game, where the system plays against the environment. Our formalism differs from existing game-based approaches [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF] in that it concisely models reconfigurations of the systems and evolutions of the environment by exploiting the featured transition approach [START_REF] Classen | Model checking lots of systems: efficient verification of temporal properties in software product lines[END_REF]. This latter provides even more flexibility to our method, which supports not only runtime configuration but also design-time variability. FTGA thus constitute an unified formalism to model the behaviour of real-time adaptive software product lines.

As a second contribution, we propose a new temporal logic to express requirements on FTGA. In [START_REF] Cordy | Model checking adaptive software with featured transition systems[END_REF], we introduced Adaptive Configuration Time Logic (AdaCTL), a variant of the Computational Tree Logic (CTL) to reason on features and reconfigurations. The main differences between AdaCTL and CTL are that (1) the existential and universal quantifiers have a game-based semantics similar to Alternating Tree Logic (ATL) [START_REF] Alur | Alternating-time temporal logic[END_REF], and (2) the satisfaction relation returns a set of configurations rather than a Boolean value. In this paper, we go one step further and introduce Timed Adaptive Configuration Time Logic (T-AdaCTL), a real-time extension of AdaCTL, the semantics of which is inspired from Timed-ATL [START_REF] Henzinger | Timed alternating-time temporal logic[END_REF].

Finally, we design efficient model-checking algorithms to verify an adaptive system modelled as an FTGA against requirements expressed in T-AdaCTL. These algorithms extend efficient timed-game algorithms [START_REF] Cassez | Efficient on-the-fly algorithms for the analysis of timed games[END_REF]. As a proof-of-concept, we implemented our method as part of PyECDAR, a modelchecker for timed systems [START_REF]A tool for timed games and timed specifications[END_REF].

Structure of the paper. In Section II, we introduce our running example and the state of the art is presented in Section III. We define FTGA in Section IV, whereas we introduce T-AdaCTL and our model checking algorithms in Section V. We discuss our implementation in Section VI.

II. INTRODUCTORY EXAMPLE

We present an example inspired by the TCP routing protocol described in [START_REF] Zhang | Modular verification of dynamically adaptive systems[END_REF]. We consider a routing protocol that can work in two different environments: a safe environment, where all the nodes are fully trusted, and an unsafe environment, where some nodes might be corrupted. In an unsafe environment, a message must be encrypted before it is sent. Every operation (routing, sending and encryption) requires time to complete. The behavior of the protocol in the two types of environment are modelled as timed automata in Fig. 1.

The protocol must satisfy safety and liveness properties. When the environment is unsafe, all the messages must be encrypted before they are sent. When the environment is safe, the messages must be sent at most 20 time units after being received. In a changing environment, the protocol must switch between the two configurations in order to adapt itself. This reconfiguration is only possible in state received. We study two different implementations: in the first one, the reconfiguration can occur at most every 25 time units; in the second one, the reconfiguration can always be done but its application requires 5 time units. We want to determine in which implementations the system can satisfy its specifications.

III. STATE OF THE ART

This paper is at the intersection of adaptive systems verification and software product line model checking. We overview relevant related work in both fields.

A. Quality Assurance for Adaptive Systems

In their research roadmap for adaptive systems [START_REF] Cheng | Software engineering for Self-Adaptive systems: A research roadmap[END_REF], Cheng et al. stated that in the context of adaptive systems, the objective of quality assurance is to provide evidence that the system is able to cope with changes in its objectives and its environment. They presented a framework for adaptive systems assurance, in which the system, the goals, and the context are subject to modifications. This results in a succession of models for the system and properties to verify. Based this idea, several verification methods model adaptive systems as a set of programs [START_REF] Allen | Specifying and analyzing dynamic software architectures[END_REF], [START_REF] Kramer | Analysing dynamic change in software architectures: A case study[END_REF], [START_REF] Kulkarni | Correctness of component-based adaptation[END_REF], [START_REF] Zhang | Modular verification of dynamically adaptive systems[END_REF]. To ensure the satisfaction of intended properties in an unstable environment, the system is able to make transitions between those programs, that is, to perform an adaptation to modify its future behaviour. One distinguishes between local properties that specific programs must satisfy, global properties that must be satisfied by any execution of the system, and transitional properties that must hold during an adaptation. To specify the transitional properties, Zhang et al. proposed a new logic called A-LTL [START_REF] Zhang | Using temporal logic to specify adaptive program semantics[END_REF] and algorithms to check a model against properties expressed in this logic [START_REF] Zhang | Modular verification of dynamically adaptive systems[END_REF]. As an alternative, Kulkarni et al. [START_REF] Kulkarni | Correctness of component-based adaptation[END_REF] consider that an adaptive system is a program able to add or remove components during runtime. Instead of traditional model checking, they use proof lattice to verify that all the adaptations satisfy the global properties. Instead of behavioural properties, Filieri et al. [START_REF] Filieri | Run-time efficient probabilistic model checking[END_REF] are interested in verifying non-functional requirements on adaptive systems. They propose novel algorithms to check parametric Markov models. Combining our work with theirs is an interesting perspective, as it could allow us to quantify the impact of adding or removing features at runtime in terms of non-functional properties.

B. Software product line model checking

SPL model checking has become a hot topic that received a lot of attention in the last years. Li et al. [START_REF] Li | Interfaces for modular feature verification[END_REF] first proposed to model features as independent automata that are subsequently composed together to form a product. They provide incremental algorithms that can check whether the connection of a new feature preserves the properties of the system. Other work is based on modal transition systems where transition are either mandatory (included in all the configurations) or optional (only part of some products) [START_REF] Fischbein | A foundation for behavioural conformance in software product line architectures[END_REF], [START_REF] Fantechi | Formal modeling for product families engineering[END_REF], [START_REF] Asirelli | A logical framework to deal with variability[END_REF]. Gruler et al. [START_REF] Gruler | Modeling and model checking software product lines[END_REF] extended the CCS process algebra with a variability operator and sketched an algorithm to verify a model against multivalued µ-calculus formulae. In the last years, a family of approaches and tools based on featured transition systems were designed [START_REF] Classen | Model checking lots of systems: efficient verification of temporal properties in software product lines[END_REF], [START_REF] Classen | Symbolic model checking of software product lines[END_REF], [START_REF] Classen | Model checking software product lines with SNIP[END_REF], [START_REF] Cordy | Simulation-based abstractions for software product-line model checking[END_REF], [START_REF] Cordy | Behavioural modelling and verification of real-time software product lines[END_REF], [START_REF] Cordy | Beyond boolean product-line model checking: Dealing with feature attributes and multifeatures[END_REF]. In featured transition systems, transitions are constrained by Boolean formula over the features, such that the system can execute a transition only if its features satisfy the associated formula. This relation does not allow for reconfigurations and is thus a particular case of ours. Cordy et al. [START_REF] Cordy | Model checking adaptive software with featured transition systems[END_REF] extended this formalism with reconfiguration capabilities in order to model adaptive systems. Their work is closely related to this paper but there are also significant differences. First, the authors do not make the distinction between controlled and uncontrolled actions. They assume that non-deterministic transitions are executed by the environment, and are thus always uncontrolled. Second, they do not consider real-time, whereas adaptive systems are often require to adapt their behaviour within a given time frame and these adaptations may not be instantaneous. Third and last, no implementation is available.

IV. FEATURED TIMED GAMES

This section introduces the mathematical model we propose to represent real-time adaptive systems. It includes a representation of an open environment with which the system interacts in real-time. This environment evolves over time, and the system must adapt its behavior to cope with these variations. To make our models concise and facilitate reasoning, we represent both the different functional modes of the system and the state of the environment with adaptive features, i.e., features that can be enabled or disabled at runtime. In standard SPLE, features usually model design-time variability and are thus not meant to be modified at runtime. Our formalism considers these features as a particular case of adaptive features. Therefore, it is flexible enough to support product lines of real-time adaptive systems.

We first introduce the syntax of the model. Then we define its semantics as a timed game. We shall see that timed games are particularly suitable to reason on system's reconfigurations with regard to change in the environment. Beforehand, we recall basic concepts to formally represent runtime variability and real-time.

A. Encoding Variability and Real-Time Constraints

Variability. In SPLE, features usually designate units of difference between software products. We extend this notion to represent the possible adaptations of the system, as well as dynamic characteristics of the environment. Therefore, we distinguish between adaptive and static features, which may or may not change at runtime, respectively. Dependencies between features can be captured in a feature model. In this paper, we define a feature model as tuple d = (F s , F a , F e , d ) where F s is the set of features of the system, F a ⊆ F s contains its adaptive features, and F e denotes the features of the environment. We assume that F s and F e are disjoint and denote their union by F . A configuration of d is any subset of F s ∪ F e . Therefore it denotes a particular variant of the system equipped with specific static and adaptive features, and deployed in a certain type of environment. Finally, d ⊆ P(F ), where P denotes the powerset, is the set of the valid configurations that satisfy the dependencies between the features.

To express that the possible behaviors of the system and the environment may depend on their features, we extend the notion of feature expressions borrowed from featured transition system (FTS) [START_REF] Classen | Model checking lots of systems: efficient verification of temporal properties in software product lines[END_REF]. FTS extends labelled transition systems, such that a transition may only be triggered by a restricted set of configurations. Each transition is labelled with a feature expression, that is a Boolean function exp : P(F ) → {⊤, ⊥} such that exp(p) = ⊤ iff p can execute the transition. We denote by exp ⊆ P(F ) the set of configurations that satisfy exp and by ⊤ the feature expression such that ⊤ = P(F ). Further in this section, we show how we generalize feature expressions to handle reconfiguration and how we combine them with time constraints.

Real-time. Timed Automata are an established formalism to represent real-time behavior. They extend labelled transition systems with real-time clocks of which value evolve as time passes. The clocks evolution and the discrete behaviors of the system are controlled by clock reset added to the transitions, and clock constraints. These constraints are either transitions guards that specifies when the system can execute a transition, or location invariants that defines when the system may remain in a given location. Examples of Timed Automata are shown in Fig. 1 to describe the models of the routing protocol.

Let C be a finite set of clocks.

A clock valuation over C is a function u : C → R ≥0 , that is, u ∈ R C ≥0 .
Given two valuations u and τ , we write u+τ for the valuation defined by (u+τ )(x) = u(x)+τ (x). for all u ∈ Z. We write ϕ to denote the set of valuations that satisfy ϕ.

Then Z ⊆ R C ≥0 is a zone iff Z = ϕ for some ϕ ∈ B(C).
To represent the behavior of system deployed in open environments, a model must distinguish between actions of the system from those of the environment. Timed Game Automata [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF] are Timed Automata where actions are either controllable (actions of the system) or uncontrollable (actions of the environment). In this formalism, the satisfaction of properties is determined by solving a two-player timed game.

B. Featured Timed Game Automata

We are now ready to introduce Featured Timed Game Automata (FTGA) as a formalism to model product lines of real-time adaptive systems. FTGA result from the combination of the encodings presented above. It provides the following modelling facilities:

(1) Open environment. An FTGA distinguishes between controllable and uncontrollable transitions. (2) Real-time. Clock constraints in invariants and transition guards model real-time constraints on the system and its environment.

(3) Variability. Each transition is constrained by a feature expression that defines in which configurations the system or its environment can execute it. It allows one to differentiate between the capabilities of every configuration. (4) Adaptations. The transition relation also encodes which reconfigurations are possible upon the execution of an action by the system or its environment. Formally, FTGA are defined as follows. The adaptation process is encoded as part of function γ. This function is defined such that only adaptive features may only be changed by controllable transitions, and environment features may only be changed by uncontrollable transitions. Formally, let α = (l, a, ϕ, λ, l ′ ) ∈ Trans. For any configurations c, c ′ , e, e ′ , if a ∈ Act c , γ(α)(c∪e, c ′ ∪e ′ ) =⇒ (c\c ′ )∪(c ′ \c) ⊆ F a ∧e ′ = e and if a ∈ Act e , γ(e i )(c ∪ e, c ′ ∪ e ′ ) =⇒ c ′ = c. Moreover, any reconfiguration of the system or the environment must ensure that the new configuration is valid, that is, γ(α)(c ∪ e, c ′ ∪ e ′ ) =⇒ c ′ ∪ e ′ ∈ d . This function provides a flexible encoding to restrict the reconfiguration process. In particular, it is able to specify the minimum and maximum amount of time needed to transit from a given configuration to another one.

To that aim, one may define a self-loop transition constrained by a given clock, and annotated with an action that represents the reconfiguration process.

C. Game semantics

An FTGA specifies the behavior of a set of systems, that is, one per valid configuration. The initial configuration of the system will determine how its behavior may evolve over time. Indeed, static features cannot be changed at runtime and thus fix parts of the system capabilities. Similarly, reconfiguration is not always doable; the initial value of adaptive features may thus impede the system to perform actions early in the execution, which may lead to unavoidable errors.

Accordingly, we define the semantics of an FTGA as a function . : P(F s ) → (Loc × R C ≥0 × P(F s ) × P(F e )) * that associates an initial system configuration with its set of infinite executions. A state of the execution is a tuple s = (l, u, c, e), where l ∈ L is a location, u ∈ R C ≥0 is a clock valuation, c ∈ P(F s ) is a system configuration and e ∈ P(F e ) is an environment configuration such that c ∪ e ∈ d . An initial state is (l 0 , 0, c 0 , e 0 ), where 0 is the valuation that initializes all clocks to zero, and c 0 , e 0 are the initial configuration of the system and the environment, respectively. Whereas the configuration of the system is an input of the semantics function, the initial configuration of the environment is uncontrolled and is thus chosen non-deterministically. Since we consider timed systems, an execution includes two types of transitions: Given that an FTGA considers continuous time, it specifies an infinite number of runs.

Among the transitions executed during a run, some are controlled by the system and others are uncontrolled, i.e. executed by the environment. Also, the system controls how it reconfigure itself, but has no control on the configuration of environment. The achievement of goals can thus be considered as a two-player games where the system plays against the environment. The strategy of one player prescribes a set of moves to perform according to the states previously visited. Each move consists of either delaying or executing an available action. A player can reconfigure itself only after executing an action. Formally, a strategy for the system is a function:

Str C : (Loc×R C ≥0 ×P(F s )×P(F e )) k → (Act c ×P(F s ))∪{τ } with k ≥ 0.
A strategy for the environment is defined symmetrically, except that the environment also selects its initial configuration. A strategy is valid iff (1) it complies with the transition relation and function γ, and (2) it does not lead to time-convergent or zeno runs [START_REF] De Alfaro | The element of surprise in timed games[END_REF]. From now on we consider valid strategies only.

The game proceeds as a concurrent game. In a given state, if one player chooses to delay while the other chooses an action, then this action is performed and the corresponding transition is triggered. If both players select an action then the transition to execute is chosen nondeterministically. Given a system strategy Str C and an environment strategy Str E , the possible outcomes of the game, noted Outcome(Str C , Str E ), are the set of infinite runs

ρ = s 0 τ0 -→s ′ 0 a1 -→s 1 τ1 -→s ′ 1 a2 -→s 2 . . . s n τn -→s ′ n an+1 ---→s n+1 . . . such that • if a i ∈ Act c then Str C (s 0 , . . . , s i ) = (a i , c i+1 ). • if a i ∈ Act e then Str E (s 0 , . . . , s i ) = (a i , e i+1 ). • if τ i ∈ R ≥0 then ∀τ ′ i ∈ [0, τ i [. s i τ -→(l i , u i + τ ′ i , c i , e i ) and Str C (s 0 , . . . , s i , (l i , u i + τ ′ i , c i , e i )) = Str E (s 0 , . . . , s i , s i + τ ′ i ) = {τ }. where s k = (l k , u k , c k , e k ) for any k ∈ N.
Example. Fig. 2 presents an FTGA modelling the routing protocol. The system actions are the plain transitions: route, reconfig, t-reconfig. The environment actions are the dashed transitions: init, receive, encryption, sent. The adaptive feature encrypt determines in which operation modes the system currently is. Two static features p-reconf and t-reconf specifies which of the two configuration methods the system can use (see Section II. Finally, the environment is described with a feature safe that specifies whether the current node in the network can be trusted or not. The function γ is defined in two steps. First, feature expressions are added in the graph to the guard of the transitions, in order to specify which set of features enables the transition. Second, we specify the possible reconfigurations:

• The system may only reconfigure the feature encrypt during the transitions labelled "reconfig". • The environment may only reconfigure the feature safe during the transitions labelled "sent" and "receive".

In consequence, a possible strategy for the environment is to start in a safe configuration, do the init action at y = 25, then the receive action at y = 30, and disable the feature safe during this transition. In reaction, the system strategy can be to start with the system feature p-reconf while the adaptive feature encrypt is disabled, then wait until the environment reaches the location Received. At this point it can do a reconfig action immediately, and enable the feature encrypt during the transition. Finally, at x = 10 it performs the route action to reach the location RoutedUnsafe. The outcome produced by these two strategies is:

" Init, » x = 0 y = 0 - ,  p-reconf safe ff " 25,init ------→ " Ready, » 25 25
- ,  p-reconf safe ff " 5,receive --------→ " Received, » 0 25 - , {p-reconf} " 0,reconf ig ---------→ " Received, » 0 0 - ,  p-reconf encrypt ff " 10,routed --------→ " RoutedUnsafe, » 0

V. TIMED ADACTL MODEL-CHECKING

To express requirements on real-time adaptive systems, we propose T-AdaCTL, a timed extension of the Adaptive Configuration Time Logic (AdaCTL), a logic we recently introduced to reason on reconfigurable systems. We first present its syntax and semantics, and then provide algorithms to check an FTGA against a T-AdaCTL formula.

A. Timed AdaCTL

The formulae of T-AdaCTL are embedded into three levels. The first level is the feature formula, which has the form Ψ ::= [χ]Φ where χ is a feature expression and Φ is a state formula. Intuitively, [χ]Φ defines that if the current configuration of the system and the environment satisfies χ, then the current state must satisfy Φ. Feature formulae can thus define requirements on specific configurations, or even forbid some others. A state formula has the form Φ :

:= ⊤ | a | Ψ 1 ∧ Ψ 2 | ¬Ψ | Aϕ | Eϕ
where a ∈ AP , Ψ, Ψ 1 and Ψ 2 are feature formulae, and ϕ is a path formula. Intuitively, a state satisfies Aϕ (resp. Eϕ) if from this state, the system can come up with a strategy of which the outcome will (resp. may) satisfy ϕ. The path formulae have the form ϕ :

:= Ψ 1 U I Ψ 2 | Ψ 1 W Ψ 2
where Ψ, Ψ 1 and Ψ 2 are feature formulae, I is an interval of R ≥0 with integral bounds, U is called the until operator and W is called the weak until operator. T-AdaCTL extends AdaCTL with a time constraint attached to the until operator, in the same manner as TCTL [START_REF] Alur | A theory of timed automata[END_REF] extends CTL. We omit the next operator of AdaCTL as there is no notion of direct successor in timed systems. Two path operators can be derived from U and W: eventually (♦), such that ♦ I Ψ = ⊤U I Ψ , and forever ( ), such that Ψ = Ψ W⊥. When it comes to state and path formulae, TATL [START_REF] Henzinger | Timed alternating-time temporal logic[END_REF] is a generalisation of T-AdaCTL, as it can express more general time constraints and requirements on the environment too. However, it does not include any notion of features, which makes it inappropriate for expressing properties on our feature-based formalism.

Example. Let us express the properties that the routing protocol must satisfy in T-AdaCTL. The property "If the environment is unsafe, all the messages must be encrypted before they are sent." can be expressed by the formula A ([¬safe]¬RoutedSaf e). This formula specifies that the system can never reach the location RoutedSaf e if the environment is not safe. The property "If the environment is safe, the messages must be sent at most 20 time units after being received." can be expressed by the formula A ([safe]Received ⇒ A♦[0, 20]Ready). It specifies that whenever location Received is reached in a safe environment, the location Ready must be reached within 20 time units.

We provide T-AdaCTL with a formal semantics:

Definition 2 Let G be an FTGA and s = (l, u, c, e) one of its states. Then the satisfiability of a T-AdaCTL feature or state formula by G in state s is determined as follows:

G, s |= [χ]Φ ⇔ c ∪ e ∈ χ ∨ G, s |= Φ G, s |= ⊤ ⇔ ⊤ G, s |= a ⇔ a ∈ L(l) G, s |= Φ 1 ∧ Φ 2 ⇔ G, s |= Φ 1 ∧ G, s |= Φ 2 G, s |= ¬Φ ⇔ G, s |= Φ G, s |= Eϕ ⇔ ∃Str C • ∃Str E • ∃ρ ∈ Outcome(s, Str C , Str E ) • G, ρ |= ϕ G, s |= Aϕ ⇔ ∃Str C • ∀Str E • ∀ρ ∈ Outcome(s, Str C , Str E ) • G, ρ |= ϕ
The semantics of path formulae is similar to that of TCTL path formulae:

G, ρ |= Ψ 1 U I Ψ 2 ⇔ ∃r ∈ I • G, ρ[r] |= Ψ 2 ∧ ∀0 ≤ r ′ < r • G, ρ[r ′ ] |= Ψ 1 G, ρ |= Ψ 1 W Ψ 2 ⇔ (∀r ′ ≥ 0 • G, ρ[r ′ ] |= Ψ 1 ) ∨ (∃r ≥ 0 • G, ρ[r] |= Ψ 2 ∧ ∀0 ≤ r ′ < r • G, ρ[r ′ ] |= Ψ 1 )
where ρ[r] the state reached in ρ at time r.

Note that we assume a continuous-time semantics for timed path operators [START_REF] Bouyer | Model-checking timed temporal logics[END_REF]. We now define the satisfaction of a T-AdaCTL formula by an FTGA. Contrary to classical temporal logics, this relation, noted |= F is not Boolean: it is defined as the set of initial system configurations such that the FTGA satisfies the formula from its initial state.

Definition 3 Let G be an FTGA and Ψ a T-AdaCTL formula.

(G |= F Ψ ) = {c 0 ∈ P(F s ) | ∃e 0 ∈ P(F e ) • c 0 ∪e 0 ∈ d ∧ ∀e 0 ∈ P(F e ) • c 0 ∪ e 0 ∈ d ⇒ G, (l 0 , 0, c 0 , e 0 ) |= Ψ }

B. Model-Checking Algorithms

The semantics of T-AdaCTL is defined over execution paths, of which FTGA contain an infinite number. This means that a model checking procedure for T-AdaCTL must use a symbolic representation to capture this infinite number of runs in a finite data structure. To represent the time domain of symbolic states, we extend the grammar of clock constraints with the negation. Then for a clock constraint, ϕ, ϕ is a federation, i.e. a finite union of zones. In combination with federations, we use feature expressions to encode sets of configurations symbolically, as opposed to representing each configuration individually in separate states. Therefore, our algorithms manipulate symbolic states, i.e. tuples of the form (l, b, ϕ), where l is a location, b is a feature expression, ϕ is a clock constraint. A symbolic state is an abstraction of all the concrete states (l, u, c, e) such that u ∈ ϕ and (c ∪ e) ∈ b .

To model check a T-AdaCTL formula Ψ , we first decompose it into its parse tree, where each node is a subformula. The root is Ψ itself, whereas the leaves are atomic formulae. Then, starting from the leaves, we associate each subformula by the set of symbolic states that satisfy it. This method is similar to the one used to check CTL formulae [START_REF] Pnueli | The temporal logic of programs[END_REF].

We present how to compute the set of symbolic states that satisfy each form of T-AdaCTL formula. For feature and state formulae, the satisfaction rules are the following:

Sat([χ]Φ) = Sat(Φ) ∪ {(l, ¬χ, ⊤) | l ∈ Loc} Sat(⊤) = {(l, ⊤, ⊤) | l ∈ Loc} Sat(a) = {(l, ⊤, ⊤) | a ∈ L(l)} Sat(Ψ 1 ∧ Ψ 2 ) = {(l, b 1 ∧ b 2 , ϕ 1 ∧ ϕ 2 ) | (l, b 1 , ϕ 1 ) ∈ Sat(Ψ 1 ) ∧ (l, b 2 , ϕ 2 ) ∈ Sat(Ψ 2 )} Sat(¬Ψ ) = Sat(Ψ )
where for any S ∈ L × P(F ) × R ≥0 , the complement of S is defined as

S = {(l, b, ϕ)|∄(l, b ′ , ϕ ′ ) ∈ S • b ∩ b ′ = ∅∧ ϕ ∩ ϕ ′ = ∅}.
Computing Sat(Eϕ) and Sat(Aϕ) comes down to solving a two-player game where the system is the verifier and the environment is the spoiler. To that aim, we perform a backward fixed-point computation as it is performed for solving timed games in [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF]. The algorithms are based on discrete predecessors and safe timed predecessors operators. The definition of these operators in FTGA takes into account both variability and real-time, which makes it different from other game-based formalisms. It constitutes the cornerstone and the real novelty of our verification algorithms. Formally, let α = (l, a, ϕ α , λ α , l ′ ) ∈ Trans and (l ′ , b ′ , ϕ ′ ) be a symbolic state. We define the discrete predecessors Pred α (l ′ , b ′ , ϕ ′ ) = (l, b, ϕ) such that:

• b = {c ∪ e | ∃(c ′ ∪ e ′ ) ∈ b ′ • γ(α)(c ∪ e, c ′ ∪ e ′ ) = ⊤} • ϕ = f ree ϕ ′ ∧{x = 0 | x ∈ λ α }, λ α ∧ϕ α ∧Inv(l), where f ree(ϕ, λ) = {u | ∃v ∈ ϕ • ∀x ∈ λ • u(x) = v(x)}.
Observe that the distributivity law applies to this operator:

Pred α i s i = i Pred α (s i )
The discrete predecessors operator can be used to compute the controllable (resp. uncontrollable) moves that allow the system (resp. the environment) to reach (resp. to avoid) a winning state. However, these moves may not be safe as the other player may have concurrent moves. Formally, given a location l and the sets of winning states W in[l ′ ] for each location l ′ , these controllable moves are:

N ext c (l, W in) = α=(l,a∈Actc,ϕ,λ,l ′ ) Pred α (W in[l ′ ])
The uncontrollable moves of the environment are defined symmetrically. The winning moves are obtained through the safe timed predecessors operator. Let s 1 = (l, b 1 , ϕ 1 ) and s 2 = (l, b 2 , ϕ 2 ) be two symbolic states, the safe timed predecessors of s 1 wrt. s 2 are the states that can reach s 1 while avoiding any state from s 2 . They are given by Pred

t (s 1 , s 2 ) = {(l, b 1 ∧ ¬b 2 , Pred t (ϕ 1 , ⊥)), (l, b 1 ∧ b 2 , Pred t (ϕ 1 , ϕ 2 ))} where Pred t (ϕ 1 , ϕ 2 )
is the safe timed predecessors operator for zones as defined in [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF]. This operator has the following property:

Pred t i g i , j b j = i j Pred t (g i , b j ).
In what follows, we denote by Pred(l) the locations from which there is a transition to l.

To enforce that the players' strategies are valid, we compute the deadlock states, which are the states beyond the locations invariant that should not be reached if one player had an urgent action to perform. We denote by DL c (l) (resp. DL e (l)) the deadlock states in location l for which the system (resp. environment) is responsible.

Algorithm 1 computes Sat(AΨ 1 U I Ψ 2 ). The algorithm starts with the winning symbolic states that satisfy Ψ 2 (Lines 3-7), and next performs a backward exploration (Lines 8-19) to discover predecessors of winning states that satisfy Ψ 1 and that the environment cannot impede the system to reach. To check that the time spent to reach the goal in Ψ 2 satisfies the interval constraint I, an additional clock, named clock, is added to the model; this is a standard way to handle timing constraints of logical formula. This extra clock is initialized in I (Line 4) and then decreases during the backward exploration. Sat(AΨ 1 U I Ψ 2 ) is the set of winning states for which the value of the extra clock is zero (Lines 20-25).

Algorithm 2 executes a similar procedure to compute Sat(AΨ 1 W Ψ 2 ). In this case, however, it starts from the states that violate the formula (Line 3) and performs a backward exploration to compute the states from which the system cannot guarantee to avoid losing states (Lines 4-14). Then, the set of winning states for the system is the complement of the set of those states (Line 15).

The algorithms used to compute Sat(EΨ 1 U I Ψ 2 ) and Sat(EΨ 1 W Ψ 2 ) also use similar procedures. The main differences are that the two players now cooperate in order to reach the goal expressed by the path formula. As an example, an algorithm for Sat(EΨ 1 U I Ψ 2 ) can be obtained by adding N ext e (l, W in) to set of states Good in Line 11 of Algorithm 1, and by removing Lined 12-13 (hence setting the set of bad states to empty set in timed predecessors). Deadlock states must be avoided. Indeed, the two players cooperate. In case of a deadlock, they both lose.

VI. IMPLEMENTATION

Our modelling formalism and the associated algorithms have been implemented on top of PyECDAR [START_REF]A tool for timed games and timed specifications[END_REF], a tool for the analysis of timed systems.

In PyECDAR, a model is written in an XML file that follows the format of UPPAAL tool set [START_REF] Behrmann | Developing uppaal over 15 years[END_REF]. This allows us to reuse the intuitive user interface provided by UPPAAL. In our extension, we use two variables for each adaptive features that defines its value before and after the reconfiguration following the transition. This encoding is sufficient to entirely represent the function γ. Additional patterns can be used to facilitate the design of the system. For example, feature expressions can be used in the invariant of a location, which offer another way to specify possible reconfigurations. Similarly, assignments can be used to forbid the reconfiguration of an adaptive feature during a transition.

The original game algorithms of PyECDAR were limited to safety and reachability objectives specific to timed specifications. We have implemented the new game algorithms presented in this paper and the recursive procedure that checks a T-AdaCTL formula. To encode continuous time, we use federations, which are finite unions of DBMs, implemented in UPPAAL DBM Library. We encode feature expressions with BDDs, implemented using PyCUDD library (python bindings for CUDD library [START_REF]Colorado University Decision Diagram Package[END_REF]). Using these two libraries we have implemented a new encoding for symbolic states with both features and time domains, as well as an encoding for finite unions of symbolic states. Therefore, the different operators (union, intersection, negation, discrete predecessors, timed predecessors) are implemented in PyECDAR for unions of symbolic states.

Example. We consider again the routing protocol modelled in Fig. 2. We first use PyECDAR to check it against the T-AdaCTL formula Ψ 1 = A ([¬safe]¬RoutedSaf e). PyEC-DAR computes the satisfaction relation for the formula, which is given by p-reconf ∨ t-reconf ∨ encrypt. It means that the formula is satisfied iff any of the two reconfiguration features are enabled, or feature encrypt is initially enabled. Then, we verify the formula Ψ 2 = A ([safe]Received ⇒ A♦[0, 20]Ready), and we obtain the following result: p-reconf ∨ ¬encrypt. Finally, we consider the formula that combines the previous ones. The satisfaction relation is now restrained to p-reconf, which proves that in order to satisfy both properties at the same time the system requires the p-reconf feature. Note that this result is not the same as the conjunction of the two previous results. Indeed, solving the formula Ψ 1 ∧Ψ 2 comes down to finding configurations in which the system has either strategies to satisfy Ψ 1 or strategies to satisfy Ψ 2 , but there may exist no strategy that satisfies both goals.

VII. CONCLUSION

Figure 1 :

 1 Figure 1: Routing protocol in safe and unsafe environment.

  For λ ∈ P(C), we write u[λ] for a valuation agreeing with u on clocks in C\λ, and setting to 0 the clocks in λ. Let B(C) denote all clock constraints ϕ generated by the grammar ϕ ::= x ≺ k | x-y ≺ k | ϕ∧ϕ, where k ∈ Q, x, y ∈ C and ≺ ∈ {<, ≤, >, ≥}. By U(C) ⊂ B(C), we denote the set of constraints restricted to upper bounds and without clock differences. For ϕ ∈ B(C) and u ∈ R C ≥0 , we write u |= ϕ iff u satisfies ϕ. For Z ⊆ R C ≥0 , we write Z |= ϕ iff u |= ϕ

Definition 1

 1 An FTGA is a tuple G = (Loc, l 0 , C, Act, Inv, Trans, d, γ, AP, L) where Loc is a finite set of locations, l 0 ∈ Loc is the initial location, C is a finite set of clocks, Act = Act c ⊎ Act e , is a finite set of actions partitioned between controllable actions in Act c and uncontrollable actions in Act e , Inv : Loc → U(C) associates an invariant to each location, Trans ⊆ Loc × Act × B(C) × P(C) × Loc is a set of transitions, d = (F s , F a , F e , d ) is a feature model, γ : T rans → (P(F ) × P(F ) → {⊥, ⊤}), specifies for each transition which configurations can execute it, and how the configuration of the system and the environment can evolve, AP is a finite set of atomic propositions, L : Loc → 2 AP labels each location ith atomic propositions it satisfies.

•

  delay transitions: (l, u, c, e) τ -→(l, u + τ, c, e) if τ ∈ R ≥0 and u + τ |= Inv(l). • discrete transitions: (l, u, c, e) a -→(l ′ , u ′ , c ′ , e ′ ) if a ∈ Act and ∃α = (l, a, ϕ, λ, l ′ ) ∈ Trans, such that: u |= ϕ, u ′ = u[λ] and γ(α)(c ∪ e, c ′ ∪ e ′ ) = ⊤.Finally, a run (or execution) in an FTGA is a sequence of states starting from an initial state and alternating delay and discrete transitions: →s n+1 . . .

10 -Figure 2 :

 102 Figure 2: FTGA of the routing protocol.

Algorithm 1 : 1 / 3 ϕ 4 W 5 W 6 end 7 / 8 l 13 N 22 ϕ 27 Algorithm 2 : 1 W ait ← ∅; 2 W 8 N

 113456781322272128 Sat(AΨ 1 U I Ψ 2 )Input: G, Ψ1, Ψ2, I Output: Sat(AΨ1 UI Ψ2) begin * Initialisation * / W ait ← ∅;2 for s = (l, b, ϕ) ∈ Sat(Ψ2) do ← ϕ ∧ clock ∈ I; in[l] ← s; ait ← W ait ∪ Pred(l); * Backward exploration * / while (W ait = ∅) do ← pop(W ait); 9 Good ← W in[l] ∪ (DLe(l) ∩ ¬DLc(l)); 10 Good ← Good ∪ N extc(l, W in); 11 Bad ← Sat(¬Ψ1) ∩ ¬W in[l]; 12 Bad ← Bad ∪ N exte(l, W in); ewW in ← Predt(Good, Bad); 14 if W in[l] N ewW in then 15 W in[l] ← N ewW in; 16 W ait ← W ait ∪ Pred(l); Loc do 20 for s = (l, b, ϕ) ∈ W in[l] do 21 ϕ ← ϕ ∧ clock = 0; ← f ree(ϕ, clock) Sat(AΨ 1 W Ψ 2 ) Input: G, Ψ1, Ψ2 Output: Sat(AΨ1 W Ψ2) begin in ← Sat(¬Ψ1) ∩ Sat(¬Ψ2);3 while (W ait = ∅) do 4 l ← pop(W ait); 5 Bad ← W in[l] ∪ DLc(l); 6 Bad ← Bad ∪ N exte(l, W in); 7 Good ← Ψ2 ∪ N extc(l, W in); ewW in ← Predt(Bad, Good \ Bad); 9 if W in[l] N ewW in then 10 W in[l] ← N ewW in; 11 W ait ← W ait ∪ Pred(l);

Ψ 3 =

 3 A ([¬safe]¬RoutedSaf e) ∧ ([safe]Received ⇒ A♦[0, 20]Ready)

This paper presents a new formal model for highly adaptive real-time systems. Our approach relies on a combination of adaptive feature transition systems with timed automata. The semantics of our model is given as a timed game, which views the system and the environments as concurrent entities. In our setting, requirements are expressed in a new logic called T-AdaCTL, for which we provide a model checking procedure. We have implemented our approach as an extension of the PyECDAR toolset. The new tool has been applied to academic case studies. Our next objective is to evaluate our approach on real-life case studies.


Benoît Boyer

Kevin Corre

Axel Legay

Louis-Marie Traonouez

Statistical Model Checking with Changes and Simulink

Statistical Model Checking (SMC) is a powerful and widely used approach that consists in extracting global information on the system by monitoring some of its executions. In this paper, we add two new stones to the cathedral of results on SMC, that are 1. a new algorithm to detect emergent behaviors at runtime, and 2. an integration of Plasma Lab, a powerful SMC checker, as a library of Simulink. Our results are illustrated on a realistic case study.

Introduction, Motivations

Complex systems such as cyber-physical systems are large-scale distributed systems, often viewed as networked embedded systems, where a large number of computational components are deployed in a physical environment. Each component collects information and offers services to its environment (e.g., environmental monitoring and control, health-care monitoring and traffic control). This information is processed either at the component, in the network or at a remote location (e.g., the base station), or in any combination of these.

Characteristic for nowadays complex systems is that they have to meet a multitude of quantitative constraints, e.g., timing constraints, power consumption, memory usage, communication bandwidth, QoS, and often under uncertainty of the behavior of the environment. There is thus the need for new mathematical foundation and supporting tools allowing to handle the combination of quantitative aspects concerning, for example, time, stochastic behavior, hybrid behavior including energy consumption. The main difficulties being that the state space of nowadays systems is too complex to be analyzed with classical validation technique. Another problem being that the interleaving of information eventually leads to undecidability.

In a series of recent works, the formal methods community has studied Statistical Model Checking techniques (SMC) [START_REF] Hérault | Approximate probabilistic model checking[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF][START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF][START_REF] Legay | Statistical model checking: An overview[END_REF] as a way to reason on quantitative (potentially undecidable) complex problems. SMC can be seen as a trade-off between testing and formal verification. The core idea of the approach is to conduct some simulations of the system and then use results from the statistic area in order to estimate the probability to satisfy a given property. Of course, in contrast with an exhaustive approach, a simulation-based solution does not guarantee a correct result. However, it is possible to bound the probability of making an error. Simulation-based methods are known to be far less memory and time intensive than exhaustive ones, and are sometimes the only option. SMC gets widely accepted in various research areas such as systems biology [START_REF] Clarke | Statistical model checking in biolab: Applications to the automated analysis of t-cell receptor signaling pathway[END_REF][START_REF] Jha | A bayesian approach to mc biological systems[END_REF] or software engineering, in particular for industrial applications. There are several reasons for this success. First, it is very simple to implement, understand and use. Second, it does not require extra modeling or specification effort, but simply an operational model of the system, that can be simulated and checked against state-based properties. Third, it allows to verify properties [START_REF] Basu | Statistical abstraction and model-checking of large heterogeneous systems[END_REF] that cannot be expressed in classical temporal logics. SMC algorithms have been implemented in a series of tools such as Ymer [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF], Prism [START_REF] Parker | The probabilistic model checker prism[END_REF], UPPAAL [START_REF] David | Time for statistical model checking of real-time systems[END_REF], Plasma Lab [START_REF] Boyer | Plasma-lab: A flexible, distributable statistical model checking library[END_REF].

Albeit SMC has already been widely adopted by the research community and by the industry. In this paper we propose two new contributions that have the potential to leverage the applicability of SMC to a broader class of systems.

We first adapt CUSUM [START_REF] Verdier | Adaptive threshold computation for CUSUM-type procedures in change detection and isolation problems[END_REF], an algorithm that can be used to detect changes in signal monitoring. We will show that CUSUM can be used to detect when the probability to satisfy a given property drops below some value. This algorithm offers new possibilities to detect, e.g., emergent behaviors in dynamic systems. Our main contributions will be to extend temporal logic with a change-based operator, discuss the calibration of the algorithm and its integration in a global monitoring procedure, and illustrate the applicability on a concrete case study.

While simulation-based approaches have the potential to be directly integrated in design tools used by industry, most of them still rely on their own input languages. Such an approach is often perceived negatively by the engineers. Indeed, formal verification should behave as a helper and not as yet another language they need to learn. We recently developed Plasma Lab an API-based SMC checker. In this paper, we show that by using the API, Plasma Lab can directly be integrated as a Simulink library, hence offering the first in house tool for the verification of stochastic Simulink models -this tool completes the panoply of validation toolsets already distributed with Simulink (see http://www.mathworks.nl/verification-validation/). Another advantage of our approach is that any advance on SMC that we will implement within Plasma Lab in the future will directly be available to the Simulink users.

Finally, our third contribution is to show the potential of our approach on a realistic case study.

Systems and Problems

Consider a set of states S and a set of states variables SV . Assume that each state variable x ∈ SV is assigned to domain D x , and define the valuation function V , such that V (s, x) ∈ D x is the value of x in state s. Consider also a time domain T ⊆ R. We propose the following definition to capture the behavior of a large class of stochastic systems.

Definition 1 (Stochastic Process).

A stochastic process is a family of random variables X = {X t | t ∈ T }, each random variable X t having range S.

An execution for a stochastic process X is any sequence of observations {x t ∈ S | t ∈ T } of the random variables X t ∈ X . It can be represented as a sequence π = (s 0 , t 0 ), (s 1 , t 1 ), . . . , (s n , t n ), such that s i ∈ S and t i ∈ T , with time stamps monotonically increasing, e.g. t i < t i+1 . Let 0 ≤ i ≤ n, we denote π i = (s i , t i ), . . . , (s n , t n ) the suffix of π starting at position i. Let s ∈ S, we denote P ath(s) to be the set of executions of X that starts in state (s, 0) (also called initial state) and P ath n (s) the set of executions of length n.

In [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF], Youness showed that the executions set of a stochastic process is a measurable space, which defines a probability measure µ over P ath(s). The precise definition of µ depends on the specific probability structure of the stochastic process being studied. We now define the general structure for stochastic discrete event systems.

Definition 2 (Stochastic Discrete Event System (SDES)). A stochastic discrete event system is a tuple Sys = S, I, T, µ, SV, V , where S is a set of states, I ⊆ S is the set of initial states, T is a set of time stamps, µ a probability measure over the executions of Sys, SV is a set of state variables and V is the valuation function.

Properties over traces of Sys are defined via the so-called Bounded Linear Temporal Logic (BLTL) [START_REF] Biere | Linear encodings of bounded LTL model checking[END_REF]. BLTL restricts Linear Temporal Logic by bounding the scope of the temporal operators. Syntactically, we have

ϕ, ϕ ′ := true | false | x ∼ v | ϕ ∧ ϕ ′ | ϕ ∨ ϕ ′ | ϕ ⇒ ϕ ′ | ¬ϕ | (ϕ) | X ≤t ϕ | F ≤t ϕ | G ≤t ϕ | ϕ U ≤t ϕ ′
where ϕ, ϕ ′ are BLTL formulas, x ∈ SV , v ∈ D x and t ∈ Q + and ∼∈ {<, ≤ , =, ≥, >}. As usual, we define F ≤t ϕ ≡ true U ≤t ϕ and G ≤t ϕ ≡ ¬F ≤t ¬ϕ. The semantics of a BLTL is defined with respect to executions π of a SDES using the following rules:

-π |= X ≤t ϕ ≡ ∃i, i = max{j | t 0 ≤ t j ≤ t 0 + t} and π i |= ϕ -π |= ϕ 1 U ≤t ϕ 2 ≡ ∃i, t 0 ≤ t i ≤ t 0 + t and π i |= ϕ 2 and ∀j, 0 ≤ j < i, π j |= ϕ 1 -π |= ϕ 1 ∧ ϕ 2 ≡ π |= ϕ 1 and π |= ϕ 2 -π |= ¬ϕ ≡ π |= ϕ and π |= x ∼ v ≡ V (s 0 , x) ∼ v -π |= true and π |= false
In the rest of the paper, we consider three problems that are 1. the quantitative problem for BLTL, 2. the optimization problem for parametric values and 3. the detection of changes.

Quantitative and Optimization Problems

Given a SDES Sys and a BLTL property ϕ, we use P r[Sys |= ϕ] to denote the probability for an execution of Sys to satisfy ϕ. The quantitative problem consists in computing this probability.

We will also study the optimization problem, that is the one of finding an initial state that maximizes/minimizes the value of a given observation. Consider a set O of observations over Sys. Each observation o ∈ O is a function o : P ath n (s) → D o that associates to each run of length n and starting at s a value in a domain D o . We denote (õ) n the average value of o(π) over all the executions π ∈ P ath n (s). The optimization problem for Sys is to determine an initial state s ∈ I that minimizes or maximizes the value (õ) n , for all o ∈ O.

As an example, an observation can simply be the maximal value of a given parameter along an execution. The average observation then becomes the sum of those observations divided by the number of runs. In this context, the optimization could be to find the initial state that minimizes the value of the parameters.

Change Detection Problem

In this section, we consider the problem of detecting whether the probability to satisfy a given BLTL property ϕ changes at execution time. More precisely, given an execution π = (s 0 , t 0 ), (s 1 , t 1 ), . . . , (s n , t n ) that satisfies ϕ with a probability p 0 <k (resp. p 0 >k) at (s 0 , t 0 ), we are looking for a suffix π i that satisfies ϕ with a probability p≥k (resp. p≤k), with k ∈]0, 1[. As an example, assume that under some traffic conditions, the firemen can extinguish a fire before three hours with a probability smaller than 0.7, it is expected that this probability decreases when the traffic increases. The challenge is to detect the time t when this change happens.

Formally, we consider a sequence of Bernoulli variables X i such that X i = 1 iff π i |= ϕ. An execution π satisfies a change p ≥ k where p = P r[ϕ] at time t, iff P r[X i = 1]<k for t i < t and P r[X i = 1]≥k for t i ≥ t. We note τ ! the fact that an execution is subjected to a change at time t. Using those notations, one can define Boolean propositions over changes and their respective time. One can also combine changes propositions with BLTL formulas, providing that those propositions are not in the scope of temporal operators. We now introduce extended BLTL change-based relations, an extension of BLTL.

:= τ τ ′ + t | τ t | τ ! | ϕ ∈ BLTL | δ ∨ δ ′ | δ ∧ δ ′ | ¬δ | δ ⇒ δ ′ | (δ) with p is a probability identifier, k ∈]0, 1[, t ∈ Q + , ∈ {<, ≤, ≥, >}, ⋆ ∈ {≤, ≥} and ϕ is a BLTL formula.
This extension allows us, e.g., to express conditions such as "if a change occurs at time t, then the system shall reach a state x in less than 10 units of time".

The Algorithms

In this section, we detail our algorithmic solutions to the three problems described in Section 2. Our solutions rely on a simulation-based approach that falls into the family of the so-called Statistical Model Checking (SMC) approach.

Quantitative Verification and Optimization

We first focus on the problem of computing the probability for a SDES Sys to satisfy a given BLTL property. For doing so, we rely on Statistical Model Checking (SMC). SMC estimates the probability that a system satisfies a property using a number of statistically independent simulation traces of an executable model. The idea being to monitor the property on each simulation, and then use an algorithm from the statistic area to compute an estimate of the probability. Those algorithms include Monte Carlo, or importance sampling/splitting [START_REF] Jégourel | Importance splitting for statistical model checking rare properties[END_REF]. Monitoring algorithms for BLTL can be found in [START_REF] Havelund | Synthesizing monitors for safety properties[END_REF]. In our work, the quantitative problem will mainly be solved in the context of calibrating a change algorithm as well as to validate quantitative properties of the model.

Optimization We now show that a simulation approach can also be used to perform an optimization of the model by varying the model parameters and evaluating the observable quantities to optimize. We consider a stochastic state transition system Sys, with a set of initial states I, and a set of observations O.

For each initial state s ∈ I we perform N random simulations of the system Sys(s) and we compute the average value of the observed quantities over the simulations. Therefore, for each observation o ∈ O we compute an estimation

1 N N i=1 o(π i ) of the average value (õ) n after runs of length n.
To solve the optimization problem, we must determine the configurations in I that optimize (minimize or maximize) these quantities. When the problem is defined with several observable quantities, we are faced with a multi-objective problem, and the best configurations are then selected by computing the Pareto frontier of the set of observations.

Observe that the simulation approach allows us to distribute the computation if the architecture permits it.

Change Detection with CUSUM

Assuming a SDES Sys, we consider the change p ≥ k where p = P r[ϕ] with ϕ a BLTL property and k ∈]0, 1[. Let X 1 , . . . X N be a finite set of samples collected during an execution π of Sys. We note p i the probability value at the sample X i . We assume the probability initially satisfies p 0 <k. The problem is stated as:

-H 0 : ∀ 1 ≤ i ≤ N, p i < k i.e. no change occurs -H 1 : ∃ 1 ≤ i ≤ N such that the change occurs at time t: ∀1 ≤ j ≤ N , we have t j < t ⇒ p i < k and t j ≥ t ⇒ p i ≥ k.
In this paper, we will decide between those hypothesis by using the CUSUM algorithm [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]. Like the Sequential Probability Ratio Test (SPRT) [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF][START_REF] Legay | Statistical model checking: An overview[END_REF], the CUSUM comparison is based on a likelihood-ratio test: it consists in computing the cumulative sum of the logarithm of the likelihood-ratio S i over the sample sequence X 1 , . . . X i and detecting the change decision as soon as S i satisfies the stopping rule.

S i = i j=1 s j s j =      ln k pj , if X j = 1 ln 1-k 1-pj , otherwise
The typical behavior of the log-likelihood ratio S i is a global decreasing before the change, and an increasing shape after the change. Then the stopping rule purpose is to detect when the positive drift is sufficiently relevant to detect the change. It consists in saving m i = min 1≤j≤i S j , the minimal value of CUSUM, and compare it with the current value. If the distance is sufficiently great, the stopping decision is taken, i.e., an alarm is raised at time t a = min{t i : S i -m i ≥ λ}, where λ is a sensitivity threshold.

In case there is no detection, we set t a = +∞. Note that we presented CUSUM monitoring for the case p ≥ k, but it could be set up for p ≤ k by defining the stopping rule for the maximum value of CUSUM instead.

CUSUM Calibration: it is important to note that the likelihood-ratio test assumes that the considered samples must be independent. This assumption may be difficult to ensure over a single execution of a system, but several heuristic solutions to guarantee independence exist. One of them consists in finding a state frequently visited during the execution of the system. Collecting exactly one sample each time such a state is visited, ensures independence between samples. Another solution being to introduce delays between the samples.

The CUSUM sensitivity depends on the choice of the threshold λ. A smaller value increases the sensitivity, i.e. the false alarms rate. A false alarm is a change detection at a time when no relevant event actually occurs in the system. Conversely, big values may delay the detection of the changes. The false alarms rate of CUSUM is defined as E[t a ], the expectation of the time before the CUSUM raises an alarm while the system is still running before the change occurs. Ideally, this value must be the biggest as possible E[t a ] → +∞. The detection delay is defined as the expectation time between the actual change of time t and the alarm time t a of the CUSUM: E[t a -t | t < t a ]. Ideally, this value has to be small as possible. In Section 5, we will propose an heuristic that uses the quantitative model checking problem in order to calibrate the algorithm.

The empirical way to choose the stopping rule: Theoretically, the properties of the CUSUM are based on the computation of the Average Run Length function (ARL) [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]. In a very few cases, this function may be computed or approximated using some approximating techniques (Wald or Siegmund) but most of the time, it is too complex to be used and to deduce λ. The alternative [START_REF] Verdier | Adaptive threshold computation for CUSUM-type procedures in change detection and isolation problems[END_REF] consists in simulating Sys 0 , that is a version of the system for which the change does not occur, i.e., the probability will not reach the value k. The idea, which relies on a Monte-Carlo approach, consists in generating the CUSUM values for a huge set of different executions of Sys 0 . The objective being to characterize the duration of false alarms that may be triggered on those situation. In Section 5 we apply such a characterization for change-based BLTL.

Monitoring for Change Relation Satisfiability Let us consider the change relation γ based on τ 1 , . . . , τ n changes. Using the syntax introduced in Section 2.2, it is expressed as let τ 1 and . . . and τ n in γ, where γ contains Boolean operations over changes and BLTL formulas. We use the following monitoring procedure for each atom:

1. For each change τ i , we set a CUSUM monitor that splits the monitoring into one sub-monitor for each random variable. Note that heuristics allows to reuse information between monitoring actions. 2. The proposition τ i ! holds iff t i = +∞. The proposition τ i t holds iff t i t.

Similarly, the proposition τ i τ j + t holds only if t i t j + t but it is undefined if t i = t j = +∞. 3. BLTL formulas can be monitored with classical techniques.

In practice, the tool generates monitors on demand for the given atoms and combine their answers in a Boolean manner.

Simulink and Plasma Lab: an Integration

The results presented in Section 3 have been implemented in the Plasma Lab SMC toolbox available at https://project.inria.fr/plasma-lab/. In this section, we first recap the main features of the tool, and then show how architecture of the implementation can be exploited in order to integrate Plasma Lab within Simulink, hence providing an in shell new verification theory for this widely used language.

On Plasma Lab

Plasma Lab is a compact, efficient and flexible platform for statistical model checking of stochastic models. The tool offers a series of SMC algorithms which includes rare events simulation, distributed SMC, non-determinism, or optimization. The main difference between Plasma Lab and other SMC tools is that Plasma Lab proposes an API abstraction of the concepts of stochastic model simulator, property checker (monitoring) and SMC algorithm. In other words, the tool has been designed to be capable of using external simulators, input languages, or SMC algorithms. This not only reduces the effort of integrating new algorithms, but also allows us to create direct plug-in interfaces with industry used specification tools. The latter being done without using extra compilers.

Figure 1 presents Plasma Lab architecture. More specifically, the relations between model simulators, property checkers, and SMC algorithms components. Simulators features include starting a new trace and simulating a model step by step. Checkers decide a property on a trace by accessing to state values. It also control the simulation, with a state on demand approach. A SMC algorithm component, such as the CUSUM algorithm, is a runnable object. The algorithm process is carried out by the run method. Progress and results are then notified through the Controller API. Samples are obtained by calls to the check function In coordination with this architecture, we use a plugin system to load models and properties components. It is then possible to support new model or property languages. Adding a simulator, checker or algorithm component is pretty straightforward as they share a similar plugin architecture. Thus, it requires only a few classes and methods to get a new component running. Each plugin contains a factory class used by Plasma Lab to instantiate component objects. These components implement the corresponding interface defining their behavior. Some companion objects are also required (results, states, identifiers) to allow communication between components and the Controller API.

On Integrating Plasma Lab within Simulink

We now show how to integrate Plasma Lab within Simulink, hence lifting the power of our simulations approaches directly within the tool. We will focus on those Simulink models with stochastic information. The experienced reader shall observe that our approach is different from the one of Zuliani et al. [START_REF] Zuliani | Bayesian statistical model checking with application to Stateflow/Simulink verification[END_REF] that consists in programming one SMC algorithm within the Matlab toolbox of Simulink. Indeed, the flexibility of our tool will allows us to incrementally add new algorithms to the toolbox without new programming efforts. Moreover, the user will directly use Plasma Lab within the Simulink interface, without tierce party. The approach is also different from the one in [START_REF] David | Smc for stochastic hybrid systems[END_REF] that consists in translating parts of Simulink models into the Uppaal language (which makes it difficult to analyse counter examples). The reader shall observe that Plasma Lab for Simulink offers the first integrated verification tool for Simulink models with stochastic information.

Simulink is a block diagram environment for multidomain simulation and Model-Based Design approach. It supports the design and simulation at the system level, automatic code generation, and the testing and verification of embedded systems continuously . Simulink provides a graphical editor, a customizable set of block libraries and solvers for modeling and simulation of dynamic systems. It is integrated within MATLAB. The Simulink models we considered have special extensions to randomly behave like failures. By default the Simulink library provides some random generators that are not compatible with statistical model checking: they always generate the same random sequence of values at each execution. To overcome this limitation we use some C-function blocks call that generate independent sequences of random draws.

Our objective was to integrate Plasma Lab as a new Simulink library. For doing so, we developed a new simulator plugin whose architecture is showed in Figure 2. One of the key point of our integration has been to exploit MATLAB Control, that is a library 1 allowing to interact with MATLAB from Java. This library uses a proxy object connected to a MATLAB session. MATLAB invokes, e.g. functions eval, feval . . . as well as variables access, that are transmitted and executed on the MATLAB session through the proxy. This allowed us to implement the features of a model component, controlling a Simulink simulation, in MATLAB language. Calls to this implementation are then done in Java from the Plasma Lab plugin.

Regarding the monitoring of properties, we exploit the simulation output of Simulink. More precisely, BLTL properties are checked over the executions of a SDES, i.e. sequences of states and time stamps based on the set of state variables SV . This set must be defined by declaring in Simulink signals as log output. During the simulation these signals are logged in a data structure containing time stamps and are then retrieved as states in Plasma Lab. One important point is that Simulink discretizes the signals trace, its sample frequency being parameterized by each block. In terms of monitoring this means that the sample frequency must be configured to observe any relevant change in the model. In practice, the frequency can be set as a constant value, or, if the model mixes both continuous data flow and state flow, the frequency can be aligned on the transitions, i.e. when a state is newly visited.

An Illustrative Case study

We applied all the three analyses presented on the Simulink model of a temperature controller in a pig shed. This model is inspired by similar studies [START_REF] Jessen | Guided controller synthesis for climate controller using uppaal tiga[END_REF][START_REF] Grabiec | Diagnosis using unfoldings of parametric time petri nets[END_REF][START_REF] David | Optimizing control strategy using smc[END_REF]. The system under control is a pig shed equipped with a fan and a heater to regulate the air temperature. Air temperature in the shed is subjected to random variations due to the variation of external temperature and the variation of the number of pigs that produce heat. The objective of the controller is to counter these variations such that the temperature remains within a given comfort zone. To do so, the controller can activate the heater to increase the temperature, and the fan to bring external air and therefore cool the shed. Then the temperature T of the shed is given by the following differential equation:

T ′ = T ext * Q -T * Q + W heater + W pigs 1 https://code.google.com/p/matlabcontrol/
where T ext is the external temperature, Q = Q min + Q f an is the air flow created by a minimal flow Q min , and an additional flow Q f an when the fan is activated, W heater is the heat produced by the heater, when activated, and W pigs is the heat produced by the pigs. This equation is modeled by the Simulink subsystem of Fig. 3. The number of pigs is computed every 500 time units according to The controller that we study applies a bang-bang strategy that is specified by four temperature thresholds, that is (1) when the temperature goes above TFanOn, the fan is turned on, (2) when the temperature returns below TFanOff, the fan is turned off, (3) when the temperature goes below THeaterOn, the heater is turned on, (4) when the temperature returns above THeaterOff, the heater is turned off. This controller is implemented by Stateflow automata given in Fig. 4. The fan and the heater are subjected to random failures when they are in use. Exponential distributions control the occurrence time of a failure. After a failure a reparation process allows to restart the fan or the heater, but it also takes a random time, exponentially distributed. These failures are modeled by two Stateflow automata, as shown in Fig. 5. In this automaton, rnd is a random number between 0 and 1, and tuse is the duration of use of the fan or heater.

The timings tfail and trepair corresponds respectively to the time of next failure, and the repair time, each chosen according to an exponential distribution with parameter lambdaFail and lambdaRepair, respectively. Additionally, the failure rate increases with usage due to wear and tear. This continue until a replacement is performed, which reset the rate.

Quantitative Verification and Optimization

The controller goal is to maintain the temperature within a comfort zone specified by a minimum and a maximum temperature (resp. T min = 15 • C and T max = 25 • C). We first consider the following values for the controller thresholds:

TFanOn = 22 • C, TFanOff = 20 • C, THeaterOn = 18 • C and THeaterOff = 20 • C.
We apply statistical model-checking to evaluate the efficiency of the controller both in the presence and absence of failures. The first BLTL property that we monitor checks that the system is not in discomfort for an excessive period of time. This is expressed by the following property:

Φ 1 = G ≤t1 F ≤t2 ¬Discomfort
where t 1 is the simulation time, t 2 is the accepted discomfort time, and Discomfort is a predicate that is true when the temperature of the system is outside the comfort zone. A dual of Φ 1 is to check for long periods without discomfort. This is possible with:

Φ 2 = F ≤t1 G ≤t2 Discomfort
Finally, a third BLTL property checks that each period of discomfort is followed by a period without discomfort:

Φ 3 = G ≤t1 G ≤t2 Discomfort ⇒ F ≤t3 (G ≤t4 ¬Discomfort)
Here t 1 and t 2 are as previously, while t 3 ≥ t 2 is the expected time at which the system returns to normal situation, and t 4 is the duration of the period without discomfort.

We use Plasma Lab to estimate the probability to satisfy these properties for different values of the timing constraints, on both models with and without failures. Each property is evaluated over a period of time t 1 = 12000 t.u. with precision ǫ = 0.01 and confidence δ = 0.01. Φ 1 and Φ 2 are evaluated for several values of t 2 . Note that for t 2 = 0, Φ 1 resumes to checking G ≤t1 ¬Discomfort. Φ 3 is evaluated with t 2 = 25 t.u. and several values of t 3 and t 4 .

The results for properties Φ 1 and Φ 2 are presented in Fig. 6. While the probabilities of satisfying Φ 1 show a significant difference between the models with and without failures, the results for Φ 2 are almost identical. This means that discomfort is as frequent in the two models, but it tends to last longer in the presence of failures. The results for Φ 3 are presented in Fig. 7. It shows again that the model without failures recovers quicker from a discomfort period.

Instead of estimating a probability using SMC techniques, we can compute the average value of two quantities in the model, namely the discomfort time, that is the cumulative time when the model is in a discomfort state, and the energy cost, computed with the duration of use of the heater and the fan. We aim at minimizing these two values by choosing adequate values of the model parameters.

Using Plasma Lab we can automatically instantiate the model with a range of values for the four temperature thresholds. We specify the ranges [START_REF] Legay | Statistical model checking: An overview[END_REF][START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF] for THeaterOn and THeaterOff, and [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF]25] for TFanOn and TFanOff, with an increment of 1. We additionally specify the following constraints to select a subset of the possible values of the parameters: TFanOff < TFanOn, THeaterOn < THeaterOff, and THeaterOn < TFanOn.

Using these constraints Plasma Lab generates a set of 225 possible configurations, for each variant of the models, with and without failures. Each configuration is automatically analyzed with 100 simulations. We then plots the average values of the cost and the discomfort in Fig. 8. These graphs helps to select the best values of the parameters by looking at the points that lies on the Pareto frontier of the data.

Change Detection

In our pig shed, the equipment may sometimes fail (heater or fan may break).

In such situation, the shed may be too frequently in the discomfort zone, which may lead to the death of several pigs. As we have seen, the probability of being in the discomfort zone is nominally very low. However, to avoid problems, one should be able to rise a flag as soon as the probability to be in the discomfort zones crosses a given threshold. Our objective is to detect that when such a change happens, there is a maintenance procedure that moves the shed out of the discomfort zone. In our example, this maintenance feature is modeled as a procedure that is regularly applied to the pig shed. Initially, the time between each maintenance is set to a very large value (500000 t.u.). The final objective is to set this time value in order to have an acceptable maintenance delay when the die risk is to heavy for the pigs (emergent behavior). This will be done by detecting changes.

We modeled the property using the change property language we proposed and we used CUSUM algorithm to check it. We first define τ to be the following change: " the probability to be in the discomfort zone more than t 1 = 100 t.u. is greater than 0.35". We are now ready to propose a property that expresses that when the change occurs, then the maintenance must be done in less than t 2 = 1000 t.u. Formally,

φ 4 = let τ = p ≥ 0.35 where p is P r G ≤t1 Discomfort in τ ! ⇒ F ≤τ +t2 Reparation
In order to perform the analysis, the CUSUM algorithm needs a calibration step. We first require an estimate of p 0 , the initial probability of being in the discomfort zone before the change occurs and the stopping sensitivity λ. To estimate p 0 , we disable failures of the temperature regulation system (fans + heaters) in the shed model. The estimate of p 0 is obtained by a Monte Carlo simulation based on the monitoring of G ≤t1 Discomfort. The property was checked every 200 t.u. over execution traces of length 21000 t.u. After 630 sec. of analysis, Plasma Lab returns p 0 ∈ [0, 0.05] with a confidence of 0.9.

Next step is to set the stopping rule, which is again done with a Monte Carlo approach. The objective is to observe several samples (value of the ratio) for several CUSUM. When there is no failure, the curve of samples should decrease. Indeed, it should only increases when failures happen, i.e., when the change happen. In practice, even without failure, the curve may locally increase for a short amount of time, which is due to the uncertainty introduced in the model. The objective is to characterize those local drifts to avoid false alarm.

To do so, we ran 100 executions of the CUSUM and observed 2000 samples per CUSUM (which corresponds to 201000 t.u.). From those experiments, we observed that the mean time (in CUSUM samples) between positive drift is 127.88 t.u. and the mean duration of positive drift is 1.2 samples. The frequency of positive sample is thus 1.2/(127.88 + 1.2), which is in the interval [0, 0.05] as predicted by Monte Carlo algorithm. In order to observe a real alarm one need to push this quotient to 0.35, which is the probability one wants to observe. This amounts to vary the duration of a positive sample, i.e., to replace 1.2 by a higher value in the above quotient. Doing so, we concluded that the probability will become greater than 0.35 when the positive drift is longer that 52 samples. From the definition of Φ 4 , we compute that the drift is ln 0.35 0.05 for each positive sample. We finally set the stopping rule to λ = 52 * ln 0.35 0.05 ≈ 101. We then launched the CUSUM on the model with failures over an execution of 210000 t.u. for the property G ≤t1 Discomfort checked every 200 t.u.. We applied Plasma Lab and observed that the stopping rule was satisfied after the sample 901 that corresponds to the simulation time 103473.34 t.u. We reproduced the experiment several times (20): we observed the change occurred at 102543.23 t.u. in average and in earlier ≈ 101000 t.u. We conclude that to satisfy Property φ 4 the maintenance operation must be scheduled at 100000 t.u.

Conclusion

We have presented extensions of SMC, implemented within the Plasma Lab toolset. We have then shown the flexibility of Plasma by integrating it as a Simulink library. This integration constitutes one of the first proof of concept that SMC can indeed be integrated as feature library in a designed tool used in industry without compiling to an intermediary language. Finally, we have shown the utility of this integration on a large-size case study. Other examples are available at https://project.inria.fr/plasma-lab/. Future work include an integration of Plasma with the FMI standard in order to verify complex heterogeneous systems. Another future work is to extend the power of distributed computing to Plasma-Simulink. The latter is technically challenging as it would require to duplicate compiled code to avoid license duplication and costs.

Definition 3 .

 3 An extended BLTL change relation is defined as: change := p ⋆ k where p = P r[ϕ] prop := let τ = change and τ ′ = change and . . . in δ δ, δ ′

Fig. 1 :Fig. 2 :

 12 Fig. 1: Plasma Lab architecture

Fig. 3 :

 3 Fig. 3: Simulink model of the differential equation controlling the temperature

Fig. 4 :

 4 Fig. 4: Temperature controller Fig. 5: Failure generator

Fig. 6 :Fig. 7 :

 67 Fig. 6: Probability estimation with SMC of satisfying Φ 1 (left) and Φ 2 (right)

Fig. 8 :

 8 Fig. 8: Optimization of the thresholds parameters without failures (left) and with failures (right)


domain: info spi math

text:
Cyrille Jegourel

Axel Legay

Sean Sedwards

A Platform for High Performance Statistical Model Checking -PLASMA

HAL is

Introduction

The need to provide accurate predictions about the behaviour of complex systems is increasingly urgent. With computational power ever-more affordable and compact, man-made systems are inevitably becoming increasingly computerised, distributed and concurrent, creating a correspondingly increased burden to check that they function correctly. At the same time, following the success of the human genome project, there is an increased expectation that computers can provide answers to important questions raised by complex systems in the life sciences.

Complex systems tend to pose two particular challenges to formal verification: the non-determinism caused by concurrency and unpredictable environmental conditions and the size of the state space. Our focus here is model checking, that can verify the most intricate details of a system's dynamical behaviour and where non-determinism may be handled by assigning probabilistic distributions to unknowns and by quantifying results with a probabilityprobabilistic model checking. 'Exact' probabilistic model checking quantifies these probabilities to the limit of numerical precision by an exhaustive exploration of the state space, but is restricted in practise by what can be conveniently stored in memory. Techniques exist to work with a reduced state space (abstraction, lumping, etc.), but the state space of most real natural and man-made systems remain intractable.

Statistical model checking (SMC) avoids an explicit representation of the state space by building a statistical model of the executions of a system and estimating results within confidence bounds. An executable model of the system is run repeatedly and each simulation trace is verified against a property specified ⋆ sean.sedwards@inria.fr in temporal logic. Examples of tools that have successfully applied this approach are [START_REF] Younes | Ymer: A statistical model checker[END_REF][START_REF] Sen | Vesta: A statistical model-checker and analyzer for probabilistic systems[END_REF]. Knowing a result with less than 100% confidence is often sufficient, since the confidence bounds may be made arbitrarily tight, however the key challenges of this approach are to reduce the length (simulation 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.) makes the production of multiple independent simulation runs relatively easy, but it is still necessary to make simulation as efficient as possible. Rare 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].

In what follows we present the prototype of a flexible SMC platform, PLASMA1 , that incorporates an in-built compiler and virtual machine to perform memoryand time-efficient simulations. PLASMA incorporates an efficient discrete event simulation algorithm and uses importance sampling to reduce the necessary number of simulation runs when properties are rare.

2 Software architecture PLASMA adopts a modular architecture to facilitate the extension of its features (Fig. 1). Models can already be specified using the PRISM reactive modules syntax [3] and biochemical syntax of the form A + B → C + D, while the implementation of other modelling formalisms, such as timed automata and procedural programming languages such as C and Java, is in prospect. The input specification is translated into a notional common intermediate language based on elements (referred to as simple commands because they have no explicit synchronisation and no choice of actions) having the structure (guard,rate,actions), where guard, rate and actions are functions over the current state (constants, variables, clocks) of the system. The intermediate language thus expresses the semantics of a system that advances by discrete events: the guard enables the command, the rate resolves non-determinism between enabled commands (and controls the delay in continuous time systems) and the actions update the state of the system. Different input languages may be facilitated by implementing parsers that construct and fill data structures that reflect simple commands. Once the model is represented in the intermediate language it is compiled into an executable form (the model program).

PLASMA uses in-built compilers to create bytecode for execution on its own stack-based virtual machine (VM) that comprises a logic VM (LogicVM in Fig. 1) and a simulation VM (SimVM in Fig. 1). PLASMA's bytecode instructions constitute a domain-specific, low level, platform-independent language designed for efficient statistical model checking. This language contains standard low level instructions, such as push, pop, add, sub, mul, div, etc., as well as non-standard instructions to construct efficient model checking algorithms. The VM is implemented in a high level procedural programming language (currently Java, but the code uses no features that cannot easily be adapted to other languages) and is efficient because it is optimised for its domain of application: high level instructions are efficient sub-parts of model-checking algorithms and all instructions are optimised with respect to the hardware level. The compiler and VM are also sufficiently compact to allow PLASMA to be implemented as a browser application, a distributed component or in an embedded system etc. PLASMA verifies properties specified in bounded temporal logic. Such properties are compiled into bytecode programs (property programs) and then executed on the logic VM. Our current focus is discrete time, however continuous time and other logics may be easily facilitated by implementing additional logic parser-compilers. Overall control of the verification process is maintained by the simulation management kernel (SMK) according to the options specified by the user. In general, the property program executes the model program until it has seen sufficient steps to decide a result and the SMK executes the property program until it has sufficient results to return an answer to the user. In this way, simulation traces contain the minimum number of states necessary to decide the property and the minimum number of simulations are generated. The logic accepts arbitrarily nested path formulae, however formulae that are not nested are particularly memory efficient: by employing a multivalued logic (true, false, undecided) PLASMA need only store the current state of the system. Nested formulae are also handled efficiently. In general, PLASMA stores only a subset of the full trace, having length equal to the maximum sum of the time bounds of any nested formulae.

Stochastic simulation algorithm

PLASMA performs discrete event simulation using the 'method of arbitrary partial propensities' (MAPP [START_REF] Sedwards | A Natural Computation Approach To Biology: Modelling Cellular Processes and Populations of Cells With Stochastic Models of P Systems[END_REF]). The MAPP is based on the Gillespie 'direct method' (DM [START_REF] Daniel | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF]) but performs significantly better in large-scale practical applications than either the DM or the asymptotically better 'next reaction method' (NRM [START_REF] Gibson | Efficient exact stochastic simulation of chemical systems with many species and many channels[END_REF]). In a system of M simple commands, each step of the DM is O(M ) because it iterates through all the commands in the system to find the command to execute and then again to update all the guards rates following the execution of the chosen command. The MAPP divides the M commands into √ M subsets of √ M commands and thus divides choosing a command into two operations of O( √ M ): choosing a subset and choosing a command within the subset. By performing an initial dependency analysis of the system the MAPP avoids updating commands whose guards and rates are not affected during the simulation. Since the average number of dependent commands, D, tends to be smaller than and independent of M , the overall complexity of the MAPP is O( √ M ). The NRM achieves asymptotic complexity of O(log M ) by performing a similar dependency analysis and by arranging the commands in an ordered binary tree whose root always contains the next command to be executed. Choosing a reaction is O(1), but maintaining the invariant property of the tree is proportional to D log 2 M . D is assumed constant, but the fact that it multiplies the complexity of the NRM tends to make the MAPP more efficient in most small to large-scale applications. See Figure 3.

Rare properties and importance sampling

The process of statistical model checking estimates the probability of a property by verifying the execution paths of multiple independent simulation runs. If Ω is a probability space of traces (ω ∈ Ω), f (ω) the probability measure over Ω and z(ω) ∈ {0, 1} is a function indicating whether ω satisfies some property, then the expected probability of the property is given exactly by γ = Ω z(ω)f (ω) dω. This leads directly to the standard Monte Carlo estimator: γ = 1 N N i=1 z(ω i ), where ω i are simulation paths realised under distribution f . As γ → 0, however, it becomes increasingly difficult to bound the relative error in γ and N becomes prohibitively large [START_REF] Heidelberger | Fast simulation of rare events in queueing and reliability models[END_REF]. Importance sampling can be used to reduce N by performing simulations under a 'tilted' (importance sampling) distribution that produces the rare paths more frequently and by compensating for the tilt using the 'likelihood ratio'.

If f ′ is the importance sampling distri- bution then l(ω) = f (ω) f ′ (ω) is the likelihood ratio and γ = Ω z(ω) f (ω) f ′ (ω) f ′ (ω) dω. This leads to the importance sampling estimator γ = 1 N N i=1 l(ω ′ i )z(ω ′ i )
, where ω ′ i are simulation paths realised under the the importance sampling distribution and l(ω ′ i ) is calculated on the fly. By individually optimising all the probabilities in the transition system ('state dependent tilting') it is conceivable to create very good importance sampling distributions, however this is not in general tractable. PLASMA therefore adopts a parameterised (state independent) tilting scheme based on its intermediate language representation of the model. For each simple command in the system, an importance sampling parameter taking strictly positive values is introduced to bias the rates. To test the performance of PLASMA's paramterised importance sampling engine we applied it to repair models from [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] that have previously been considered in the context of state dependent tilting and which may be verified by PRISM's numerical algorithm. We have found that our state independent tilting scheme is nevertheless capable of achieving dramatic increases in performance. For instance, example 1 of [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] considers a property with probability 1.17 × 10 -7 , requiring an expected 10 8 simulation runs to see a few examples. Using just six parameters PLASMA is able to make a 10 6 -fold increase in the frequency of observing the rare event.

Results

Figure 2 illustrates typical performance scaling2 of PLASMA's SMC engine relative to PRISM's numerical algorithm, applying them to increasing instance sizes of a classic probabilistic model checking problem (the randomised dining philosophers protocol of [START_REF] Pnueli | Verication of multiprocess probabilistic protocols[END_REF]). The state space increases exponentially with respect to instance size, hence PRISM's time scaling is exponential and its maximum instance size is here limited by available memory to about 35 philosophers. By estimating a result with (arbitrarily) bounded error, PLASMA can work with much larger models and its performance scales linearly in time proportional to the length of the property formula. Since PLASMA's memory requirement also scales linearly with instance size, its limit is much higher than the maximum shown in Figure 2.

Figure 3 illustrates typical simulation performance scaling of PLASMA's MAPP algorithm in comparison to the direct method of [START_REF] Daniel | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] and PRISM's simulation engine. A stochastic oscillatory model from systems biology [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF] was used as the building block to construct plausible biological models of increasing complexity. Using the DM's O(M ) scaling as reference, the lower order scaling of PLASMA's MAPP algorithm is clear. The performance of PRISM's simulation algorithm is also here limited by memory, but in this case the limit is not related to the state space, which is intractable for even the smallest instance.

Conclusion and future challenges

PLASMA is a compact, efficient and flexible SMC platform that incorporates a novel importance sampling engine. Its broad goal is to take SMC beyond proof of concept and to tackle the analysis of real-world systems. Since such systems are usually not written in abstract modelling languages, we foresee a need to implement other input languages to avoid errors and make the process of verification more convenient. Of particular interest are timed and hybrid systems, since these are commonly used in industrial applications. Importance sampling constitutes a major thread of our research, as it has the potential to dramatically increase the performance of simulation-based techniques. A key challenge is the discovery of good parameterised importance sampling distributions and we are currently developing algorithms to infer these automatically.

Fig. 1 .

 1 Fig. 1. The architecture of PLASMA

Fig. 2 .

 2 Fig. 2.Performance scaling of PLASMA's SMC engine vs. PRISM's numerical algorithm. PLASMA performed 118595 simulations per point.

Fig. 3 .

 3 Fig.3. Simulation algorithm performance scaling using the genetic oscillator of[START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF] as a building block. The DM was implemented in PLASMA.

A demonstration version of PLASMA may be downloaded from https://sites.google.com/site/plasmasmc

All results generated under Windows 7 Enterprise 64-bit and Java 1.6.0 26

32-bit on Intel Core i7 CPU M640 @ 2.80Ghz with 4GB RAM. PRISM 4.0.1 was used.