Keywords:

going from Paris to Nantes (150 km/h) and then using this value in the model. In this case, the model would allow us to conclude that the time taken by a train from Paris to Nantes is approximately 2h16. On the other hand, modeling with uncertainties would require to take into account in the model that the speed is not perfectly known. This could be done in several ways (the following description is not exhaustive):

• The first way could be to include in the model a choice between high-speed and low-speed 1. My personal belief -which might be debated, but this is out of the scope of this document -is that probabilities are just a mathematical concept that allows to abstract deterministic, but unknown, phenomena. 4.

ACKNOWLEDGEMENT

First, I would like to thank my reviewers and Jury. I am truly honoured that scientists whom I admire and respect, both for their scientific contributions and for their personality, have accepted to review my work and listen to my presentation. I could not have dreamt of a better jury.

I would also like to thank all my co-authors, without whom this document would be much shorter. In particular, I would like to thank my former Ph.D. students and postdocs: Amine, Eva, Dimitri, Hadrien and Paulin. It has been an honour working with you, and although none of you have (yet) pursued in Academia, it would be my pleasure to work again at your side. Thank you Anicet for the great work we did together. Although you were not officially my Ph.D. student, I would have been proud to consider you as such. There are of course too many names to cite here, but I will insist on a few names: Thank you Etienne for making me part of the ANR PACS crew, and also for the time spent on PIPTAS (what a name...). Thank you Didier for always being available for scientific (or non-scientific) discussions and also for the prank you played at QEST'17 -I will remember. Thank you Kim for taking me as a postdoc and spending a lot of your precious time working with me.

Thank you to all my colleagues at LS2N and beyond (who also might be co-authors, but I hope they can forgive me for thanking them twice). I spend more than half of my days at work and I could not do so without you. In particular, thank you to the AeLoS team for taking me in, thank you to Christian for being the best team leader ever, thank you to Arnaud and Charlotte for being the best officemates ever, thank you to Pascal for so many great discussions (and also for proof-reading this manuscript), and of course thank you to Claude for pushing me to apply in Nantes and supporting me since: best idea ever. Thank you Colin for making me welcome in Nantes (I will remember the first words you ever said to me: "Toi, c'est les automates probabilistes"). Thank you to the ComBi team for giving me so many occasions to drink beers: Audrey, Géraldine, Guillaume (who also proof-read this manuscript), Jérémie, Samuel. A special thank you to Damien E. for so many trips, beers, discussions, coffee-breaks, and often all of them at once. Thank you for your enthusiasm and for convincing me that model-checking and verification can be interesting to others than computer-scientists. I know we don't agree on the nature of probabilities, but it will be my pleasure to try to convince you again over our next beer.

Merci à mes collègues du département informatique, en particulier à ceux qui se sacrifient pour prendre les responsabilités dont personne ne veut : Emmanuel, Eric, Laura, Salima. . . Je tiens aussi à remercier les nombreux et nombreuses collègues qui font en sorte que nous puissions faire notre travail dans les meilleures conditions, qui nous simplifient la vie au maximum et réparent nos bêtises : Anne-Françoise, Annie, Elodie, Caroline, Charlery, Christine, Damien V., Jean-Yves, Karine, Marion, Sophie, Virginie, Vivian et j'en oublie probablement. . . Finalement, un grand merci à mes amis et à ma famille, qui me supportent et me soutiennent en permanence depuis toujours.

TABLE OF CONTENTS

List of Tables and Figures 10

INTRODUCTION Science is a systematic enterprise that builds and organizes knowledge in the form of testable explanations and predictions about the universe. [. . .] Modern science is typically divided into three major branches that consist of the natural sciences (e.g., biology, chemistry, and physics), which study nature in the broadest sense; the social sciences (e.g., economics, psychology, and sociology), which study individuals and societies; and the formal sciences (e.g., logic, mathematics, and theoretical computer science), which study abstract concepts. [. . .]

Science 1 As one can see from the above definition, the aim of Science, in general, is to gather and organize knowledge about given objects of study (which we will call systems in the rest of the document).

As varied as these systems can be, a common concept is used across all modern science branches in order to represent our understanding of these systems: the concept of model. A model is in fact the abstract representation that is carefully built from the gathered knowledge on the system and then used instead of the system itself in order to perform analyses or prediction. When a scientist studies a system through its model, his/her goal is to draw conclusions/predictions about the system itself that are as little biased as possible by the abstraction performed to obtain the model from the "real" system or the perception of it. Although a model is by definition necessarily imperfect, as it is an abstraction that only reflects our understanding (or point of view) of the system, the aim of scientists is to make it as accurate (i.e., faithful to the system) as possible while preserving our ability to use it in our studies. Unfortunately, these two goals are usually antinomic: the more precise the model is, the more difficult it is to use for prediction and analysis for computational reasons.

While natural and social sciences mostly use models as a tool for representing and analyzing given systems, some branches of formal sciences consider the models themselves as their objects of study. For instance, some branches of theoretical computer science are dedicated to developing Introduction languages, formalisms, theories, techniques and tools to build and analyze models, with a priori no restriction on the systems these models represent. This is the broad context in which the contributions presented in this document lie.

Modeling with uncertainties

As explained above, the point of modeling a given system is to represent the knowledge we have of this system and use the model for analysis or prediction. Most of the time, the knowledge we have is incomplete, which is the reason why analysis is required in the first place. There are then two ways to account for this incomplete knowledge: (1) Deterministic modeling and (2) Modeling with uncertainties. Deterministic modeling amounts to only including in the models what is known about the systems, while abstracting away all that is unknown. In this case, the model is an imperfect abstraction of the system and the modeler knows that the results obtained through its analysis might be imperfect themselves, without specific information about these imperfections.

On the other hand, modeling with uncertainties amounts to including the uncertainties we have about the system inside of the model. In this case, the model is an abstraction of the imperfect system perception and its analysis yields information not only about the system itself but also about the approximations that have been performed during abstraction in order to obtain the model.

Example 1. As an illustration, assume that one wants to study the time that a train takes to go from Paris to Nantes. If we deliberately oversimplify the problem, we know that this duration can be computed as the ratio between the distance from Paris to Nantes and the speed of the train. While the distance from Paris to Nantes is perfectly known (342 km), we do not know the exact speed of the train. Indeed, depending on the type of train and the weather conditions, this speed can vary. Deterministic modeling, in this case, would amount to measuring the average speed of trains Introduction trains, thus leading to a non-deterministic model that will yield different outputs depending on the choice of the train.

• Another way of taking this uncertainty into account could be to measure the probability distribution of train speeds and include this probability distribution in the model, leading to a stochastic model. The outputs of such a model would provide detailed information on the resulting probability distribution on possible durations.

• Finally, the speed could be taken into account as a symbolic parameter of the model. The result would be a parametric model, that would then yield a symbolic answer: the duration as a function of the speed of the train.

As one can see in Example 1, modeling with uncertainties can be done in several ways, the most general of which is the use of non-determinism. Indeed, non-determinism allows to represent complete lack of knowledge (but also lack of control) on a given system. Using non-determinism allows to incorporate a notion of choice inside the model, while having no information on how this choice is resolved. On the other hand, probabilities 1 represent some knowledge on the underlying choice: using probabilities allows to specify that there is a choice, and while we do not know what the outcome of this choice is, we have some information on how this outcome is chosen. Finally, the use of parameters allows to perform symbolic analysis, i.e., producing results that are a function of those parameters, despite the uncertainty about their value. In this document, we will propose and study modeling formalisms that account for uncertainties in all these ways. Although modeling with uncertainties yields more informative results than deterministic modeling (as can be seen from the example above), deterministic models are still widespread in most natural and social sciences. The main reasons for this fact are the lack of communication between formal sciences, which have been working on models with uncertainties for a long time, and those who might benefit from such models (see [START_REF] Delahaye | On the power of uncertainties in microbial system modeling: No need to hide them anymore[END_REF] for an argumentation in the context of microbial modeling), but also the lack of efficient methods for analyzing large models with uncertainties. Indeed, since the aim of modeling is to pursue the analysis of the underlying system, some techniques and tools are necessary for performing such analysis.

The most common analysis method used in natural and social sciences is simulation, which consists in running a given model and observing its outputs. Simulation can be used both in the modeling and in the analysis phases: during modeling, simulation is used on intermediate models, and the outputs are compared to existing knowledge on the system (either in the form of Introduction expertise or in the form of experimental data) in order to help the modeler refine its model. In the analysis phase, simulation is used on the final model in order to obtain predictions. In the case of deterministic systems, simulation is an appropriate technique as long as the outputs can be analyzed easily. However, when considering models with uncertainties, simulations need to be performed countless times in order to be able to observe a sufficient portion of the model's potential outputs. In this context, other alternative automated techniques are necessary. A substantial research effort -in particular in theoretical computer science -has been gathered during the last decades in order to develop such techniques. Among them lie automated verification.

Automated verification

Automated verification encompasses a multitude of techniques dedicated to analyzing models. The common ground between all those techniques is the need of having as inputs a formal model, i.e., a model described in an unambiguous language with a well-defined semantics, and a set of formal properties, i.e., queries that are also described in an unambiguous language and that the model/system must (or must not) satisfy. Model verification then consists in an automated (or semi-automated) analysis of the model to determine whether the given property is satisfied. Some verification techniques such as testing [START_REF] Glenford | The art of software testing[END_REF] or statistical model checking [START_REF] Legay | Statistical model checking: An overview[END_REF] are based on extensive simulation; Others, such as model checking [CJGK + 18] instead rely on an exhaustive exploration of the model. Finally, automated theorem-proving [START_REF] Fitting | First-order logic and automated theorem proving[END_REF] uses proof-theory along with the formal semantics of the model to prove (or disprove) that the model satisfies the given property. Whether they rely on simulation, exploration or proof-theory, automated verification techniques always provide more information than simple simulation. Indeed, a simple simulation only represents a single example of what the model can output. In the case of models with uncertainties, the potential outputs are numerous and therefore observing that a single output satisfies (or dissatisfies) the given property is not sufficient to conclude anything about the complete model. Instead, all automated verification techniques offer some guarantees about the obtained results: from percentage of covered test-cases in the context of testing to mathematical proofs for theorem-proving and model checking through error-rate and guaranteed precision for statistical model checking. Moreover, automated verification techniques also provide feedback on the results, mostly in terms of counter-examples when the property is not satisfied.

The unambiguity of the languages used to describe the model and the properties are requirements made necessary by the need for automation. Indeed, in order to develop tools that can

Introduction

automatically analyze models, one needs to assume that the languages used to describe the model and the properties are well-defined, and can be automatically interpreted by a machine.

A plethora of such languages (also called formalisms) have been developed over the years, from programming languages [KR + 88, AGHH00, S + 99, RV98] to more abstract specification formalisms [START_REF] Hopcroft | Introduction to automata theory, languages, and computation[END_REF][START_REF] Milner | Communicating and mobile systems: the pi calculus[END_REF][START_REF] Abrial | The B-Book: Assigning programs to meanings[END_REF][START_REF] Spivey | Understanding Z: a specification language and its formal semantics[END_REF], each of them dedicated to the description of certain types of systems2 . Unfortunately as formalisms become more and more expressive (i.e., capable of representing more and more detailed features of the systems), the accompanying verification techniques also raise more and more complexity, to the point that the automated verification of models expressed in some types of formalisms has been shown to be undecidable [START_REF] Alur | Model-checking for realtime systems[END_REF][START_REF] Esparza | On the decidability of model checking for several µ-calculi and petri nets[END_REF]. More research effort is therefore necessary to propose and study novel modeling formalisms and develop dedicated automated verification techniques. In this document, I will present some such contributions dedicated to the analysis of the behavior of discrete and continuous systems with uncertainty.

Contributions in the past 10 years

Since my Ph.D. thesis in 2010 [START_REF] Delahaye | Modular Specification and Compositional Analysis of Stochastic Systems[END_REF], I have been interested in the complete scope of formal modeling and analysis of systems with uncertainties: from the theoretical study of modeling formalisms and verification techniques to the application of those theories to the analysis of concrete systems. Two main lines of works could in particular be observed in my thesis:

2013-present: Associate professor position

My main interests have not changed since I arrived in Nantes in 2013: my contributions can still be classified in the two lines presented above, i.e., (1) the theoretical study of abstract modeling formalisms and analysis techniques, and (2) the development of practical verification techniques and tools and their application to concrete systems. However, the type of concrete systems I consider has slightly shifted. Indeed, interactions with several colleagues in the context of the ANR project PACS3 have convinced me that parameters play an important role when modeling complex systems as they allow, in particular, symbolic analysis and robust synthesis. My main contributions since then have therefore been aimed at modeling and verifying (timed) probabilistic and parametric systems.

In the theoretical context, I have investigated several formalisms and techniques allowing to describe and analyze parametric systems. In particular, I have participated to the following contributions: Probabilistic Event-B. From 2014 to 2017, I have contributed to the supervision of Mohamed Amine AOUADHI on the topic of Probabilistic Event-B. During this thesis, we have extended the Event-B formalism, which allows for the specification of parametric systems in a refinementbased fashion, to take into account probabilistic information [START_REF] Amine Aouadhi | Moving from Event-B to probabilistic Event-B[END_REF][START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF]. Besides enhancing the Event-B language and semantics to take into account probabilistic choice, we have also investigated the probabilistic refinement of probabilistic Event-B models as well as the automatic proving of some probabilistic properties. Probabilistic time Petri nets. In 2016, we have proposed in [START_REF] Emzivat | Probabilistic time petri nets[END_REF] a new probabilistic syntax and semantics for Petri Nets, where, instead of using the standard stochastic semantics (i.e., where probabilistic choice is resolved using races between concurrent transitions), we offer a discrete probabilistic choice that is independent of the timing information on transitions. Our aim with this formalism was to extend our results to the parametric context, where parameters can range both on the timing information on transitions and on the discrete probabilities. However, meaningful results in the parametric context still elude us.

Introduction

Parametric interval Markov chains. In 2015, I have proposed a first version of Interval Markov

Chains where the interval endpoints can be replaced with parametric functions, along with first attempts to verifying meaningful properties such as consistency [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF]. This line of work has been pursued with many co-authors, first yielding algorithms for parameter synthesis and verification of reachability properties [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]. Between 2016 and 2017, I have then contributed to the supervision of Anicet BART on the adaptation of constraint programming techniques to this context in order to make the verification more efficient [BDL + 17, BDF + 18]. Notably, [BDL + 17] has been awarded the best paper award at the QEST conference in 2017.

In 2016, we have also extended the Parametric Interval Markov Chains formalism to take into account timing information, which has yielded Parametric Interval Probabilistic Timed Automata [START_REF] André | Consistency in parametric interval probabilistic timed automata[END_REF][START_REF] André | Consistency in parametric interval probabilistic timed automata[END_REF]. In this context, we have shown that problems such as consistency and reachability are undecidable. Nevertheless, we have also provided semi-algorithms that allow to solve those problems under certain conditions. In these papers, parameters can only range on the timing information. Again, our intent was to have parameters ranging both on the timing information and probabilities, but the work is still ongoing.

On the other hand, I have have also pursued my work on the verification of concrete systems, by using (parametric) statistical techniques. Parametric statistical model checking. Statistical model checking is a technique that has for a long time been limited to purely probabilistic systems. Indeed, because it relies on statistical sampling, it could not previously be used on systems where the transition probabilities were underspecified. In 2017, during the postdoctoral position of Paulin FOURNIER, we have used an adaptation of a technique called importance sampling in order to extend statistical model checking to parametric systems [START_REF] Delahaye | Statistical model checking for parameterized models[END_REF]. This new technique has been implemented in a prototype tool and applied to a concrete industrial case study targeting the analysis of the flight plan of a UAV 4 [BAD + 19] during the thesis of Ran BAO, which I have contributed to supervise from 2017 to 2020. Statistical model checking for parameter synthesis. Since 2016, I have developed a close collaboration with researchers in Oceanography (Lionel Guidy, LoV, Villefranche -Olivier Aumont, LOCEAN, Paris -Chris Bowler, IBENS, Paris) and Microbiology (Nicholas Bouskill, LBL, Berkeley), where I believe that some of the techniques I develop could be particularly useful. However, the first observation that we made was that those natural sciences mostly use deterministic modeling, despite of the benefits that would be gained by taking uncertainties into account. In 2017, we have therefore proposed an argument in favor of modeling with Introduction uncertainties in the context of Microbiology [START_REF] Delahaye | On the power of uncertainties in microbial system modeling: No need to hide them anymore[END_REF]. We have also adapted a model taken from Oceanography in order to include uncertainties and used statistical model checking for synthesizing the optimal parameter values with respect to experimental data. This result, which has been published in 2020 in Scientific Reports [REG + 20], a reknown generalist journal, has received praises from the Biology community, which is, to my opinion, a major achievement. Automated learning. Finally, I have also considered the problem of automated modeling. Indeed, we would be very interested in some cases in using automated techniques for building models representing the statistical characteristics of some data sets while not necessarily building on the current knowledge of the underlying systems. To this purpose, I have contributed to the supervision of the thesis of Dimitri ANTAKLY from 2017-2020 on the development of automated learning and verification techniques for security assessment. During his thesis, we have developed learning techniques for a recent formalism called Recursive Timescale Graphical Event Models along with verification techniques allowing to select a model that is the most representative of the input data while satisfying some given properties [START_REF] Antakly | Graphical event model learning and verification for security assessment[END_REF]. This new method has for now only been applied to security assessment of software data, but our intent is to port it to the Oceanography context in the near future.

Introduction

-Finally, the last Chapter presents perspectives for my future research topics.

In addition to the chapters presented above, this document contains appendices. In particular, Appendix B gathers the main notations and definitions introduced in each chapter. It is aimed to be detached from the manuscript so the reader can use it while reading the main part of the document.

CHAPTER 1

PARAMETRIC INTERVAL MARKOV CHAINS

The work we present in this chapter was mostly achieved during Anicet Bart's thesis, although it follows a line of work initiated more than 10 years ago on specification theories for probabilistic systems, which has yielded numerous publications on several types of models: Modal This chapter focuses on systems where uncertainty comes from several sources. The first source is the probabilistic nature of the systems themselves (Markov Chains). The second source is the underspecification of those systems: instead of considering probabilistic transitions as usually done in Markov Chains, we consider here intervals of probabilities on the transitions, therefore leading to a specification formalism. In other word, each Interval Markov Chain we analyze is an abstract representation of a familly of systems. Finally, the last source of uncertainty is a second level of underspecification: the endpoints of the probabilistic intervals might not be known precisely, and are therefore replaced with symbolic parameters. While we had already studied parametric Interval Markov Chains in [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF][START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], this chapter presents an alternative and more efficient solution to several problems of interest already solved in [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], using techniques that are adapted from Constraint Programming. For the sake of conciseness, most of the proofs of the results presented in this chapter have been left out of this document. The interested reader can find them in [BDF + 18].

Chapter 1 -Parametric interval Markov chains

Introduction

Discrete time Markov chains (MCs for short) are a standard probabilistic modeling formalism that has been extensively used in the litterature to reason about software [START_REF] James | A Markov chain model for statistical software testing[END_REF] and real-life systems [START_REF] Husmeier | Probabilistic Modeling in Bioinformatics and Medical Informatics[END_REF]. However, when modeling real-life systems, the exact value of transition probabilities may not be known precisely. Several formalisms abstracting MCs have therefore been developed. Parametric Markov chains [START_REF] Alur | Parametric real-time reasoning[END_REF] (pMCs for short) extend MCs by allowing parameters to appear in transition probabilities. In this formalism, parameters are variables and transition probabilities may be expressed as polynomials or rational functions over these variables. A given pMC represents a potentially infinite set of MCs, obtained by replacing each parameter by a given value. pMCs are particularly useful to represent systems where dependencies between transition probabilities are required. Indeed, a given parameter may appear in several distinct transition probabilities, which requires that the same value is given to all its occurences. Interval Markov chains [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF] (IMCs for short) extend MCs by allowing precise transition probabilities to be replaced by intervals, but cannot represent dependencies between distinct transitions. IMCs have mainly been studied with three distinct semantics interpretations. Under the once-and-for-all semantics, a given IMC represents a potentially infinite number of MCs where transition probabilities are chosen inside the specified intervals while keeping the same underlying graph structure. The interval-Markov-decision-process semantics (IMDP for short), such as presented in [START_REF] Chatterjee | Model-checking omega-regular properties of interval Markov chains[END_REF][START_REF] Sen | Model-checking Markov chains in the presence of uncertainties[END_REF], does not require MCs to preserve the underlying graph structure of the original IMC but instead allows a finite "unfolding" of the original graph structure: new probability values inside the intervals can be chosen each time a state is visited. Finally, the at-every-step semantics, which was the original semantics given to IMCs in [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF], does not preserve the underlying graph structure too while allowing to "aggregate" and "split" states of the original IMC in the manner of probabilistic simulation.

Model checking algorithms and tools have been developed in the context of pMCs [DJJ + 15, HHWZ10, KNP11] and IMCs with the once-and-for-all and the IMDP semantics [START_REF] Chakraborty | Model checking of open interval Markov chains[END_REF][START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF]. State of the art tools [DJJ + 15] for pMC verification compute a rational function on the parameters that characterizes the probability of satisfying a given property, and then use external tools such as SMT solving [DJJ + 15] for computing the satisfying parameter values. For these methods to be viable in practice, the allowed number of parameters is quite limited. On the other hand, the model checking procedure for IMCs presented in [START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF] is adapted from machine learning and builds successive refinements of the original IMCs that optimize the probability of satisfying the given property. This algorithm converges, but not necessarily to a global optimum. It is worth noticing that existing model checking procedures for pMCs and IMCs strongly rely on their underlying graph structure (i.e., respect the once-and-for-all semantics). However, in [START_REF] Chatterjee | Model-checking omega-regular properties of interval Markov chains[END_REF] the authors perform model checking of ω-PCTL formulas on IMCs w.r.t. the IMDP semantics and they show that model checking of LTL formulas can be solved for the IMDP semantics by reduction to the model checking problem of ω-PCTL on IMCs with the IMDP semantics. For all that, to the best of our knowledge, no solutions for model checking IMCs with the at-every-step semantics have been proposed yet.

In this chapter, we focus on Parametric interval Markov chains [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] (pIMCs for short), that generalize both IMCs and pMCs by allowing parameters to appear in the endpoints of the intervals specifying transition probabilities, and summarize the contributions obtained during Anicet Bart's thesis. Details and proofs of these contributions can be found in [BDF + 18, BDL + 17]. The main results are the following.

First, we formally compare abstraction formalisms for MCs in terms of succinctness: we show in particular in Proposition 2 that pIMCs are (strictly) more succinct than both pMCs and IMCs when equipped with the right semantics. In other words, everything that can be expressed using pMCs or IMCs can also be expressed using pIMCs while the reverse does not always hold.

Second, we prove in Theorem 1 that the once-and-for-all, the IMDP, and the at-every-step semantics are equivalent w.r.t. reachability properties, both in the IMC and in the pIMC settings. Notably, this result gives theoretical backing to the generalization of existing works on the verification of IMCs to the at-every-step semantics.

Third, we study the parametric verification of fundamental properties at the pIMC level: consistency, qualitative reachability, and quantitative reachability. Given the expressivity of the pIMC formalism, the risk of producing a pIMC specification that is incoherent and therefore does not model any concrete MC is high. Consistency is therefore of paramount importance. In order to provide solutions to consistency and reachability, we model these problems using constraint encodings and use state of the art constraint solvers for solving them. Constraints are first order logic predicates used for modeling and solving combinatorial problems [START_REF] Rossi | Handbook of Constraint Programming[END_REF]. A problem is described with a list of variables, each in a given domain of possible values, together with a list of constraints over these variables. Such problems are then sent to solvers which decide whether the problem is satisfiable, i.e., if there exists a valuation of the variables satisfying all constraints, and in this case computes a solution. We therefore propose, in Section 1.4.1 and Proposition 3, constraint encodings for deciding whether a given pIMC is consistent and, if so, synthesizing parameter values ensuring consistency. We then extend, in Section 1.4.2 and Proposition 4, these encodings to qualitative reachability, i.e., ensuring that given state labels are reachable in all (resp. none) of the MCs modeled by a given pIMC. Finally, in Section 1.5.2 and Theorem 2, we focus on the quantitative reachability problem, i.e., synthesizing parameter values such that the probability of reaching given state labels satisfies fixed bounds in at least one (resp. all) MCs modeled by a given pIMC. While consistency and qualitative reachability for pIMCs have already been studied in [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], the constraint encodings we propose are significantly smaller (linear instead of exponential). To the best of our knowledge, our results have provided the first solution to the quantitative reachability problem for pIMCs.

The last contribution, presented in Section 1.6, is the implementation of all our verification algorithms in a prototype tool that generates the required constraint encodings and can be plugged to any SMT solver for their resolution.

Background

We start by introducing notions and notations that will be used throughout this chapter. Given a finite set of variables X = { x 1 , . . . , x k } we call domain a set of values associated to a variable in X. We write D x for the domain associated to the variable x ∈ X and D X for the set of domains associated to the variables in

X. A valuation v over X is a set v = { (x, d) | x ∈ X, d ∈ D x } of elementary valuations (x, d)
where for each x ∈ X there exists a unique pair of the form (x, d) in v. When clear from the context, we write v(x) = d for the value given to variable x according to valuation v. A rational function f over X is a division of two (multivariate) polynomials g 1 and g 2 over X with rational coefficients, i.e., f = g 1 /g 2 . We write Q for the set of rational numbers and Q X for the set of rational functions over X. The evaluation v(g) of a polynomial g under the valuation v replaces each variable x ∈ X by its value v(x).

An atomic constraint over X is a Boolean expression of the form f (X) g(X), with ∈ { ≤, ≥, <, >, = } and f and g two functions over variables in X. An atomic constraint is linear if the functions f and g are linear. A constraint over X is a Boolean combination of atomic constraints over X.

Given a finite set of states S, we write Dist(S) for the set of probability distributions over S, i.e., the set of functions µ : S → [0, 1] such that ∑ s∈S µ(s) = 1. We write I for the set containing all interval subsets of [0, 1]. In the following, we consider a universal set of symbols A that we use for labelling the states of our structures. We call these symbols atomic propositions. We will use Latin alphabet in state context and Greek alphabet in atomic proposition context.

Constraints. Constraints are first order logic predicates over a given set of variables. Informally,

1.2. Background
a Constraint Satisfaction Problem (CSP) consists of a set of variables associated to given domains and subject to a set of constraints. A CSP is satisfiable if there exists a valuation of the variables that satisfies all constraints. Checking satisfiability of constraint problems is difficult in general, as the space of all possible valuations has a size exponential in the number of variables.

Definition 1 (Constraint Satisfaction Problem). A Constraint Satisfaction Problem (CSP) is a tuple Ω = (X, D,C)
where X is a finite set of variables, D = D X is the set of all domains associated to the variables from X, and C is a set of constraints over X.

We say that a valuation over X satisfies Ω if and only if it satisfies all constraints in C. We write v(C) for the satisfaction result of the valuation of the constraints C according to v (i.e., true or false). In the following we call CSP encoding a scheme for formulating a given problem into a CSP. The size of a CSP corresponds to the number of variables and atomic constraints appearing in the problem. Note that, in constraint programming, having less variables or less constraints during the encoding does not necessarily imply faster solving time of the problems.

Discrete Time Markov Chains.

Definition 2 (Discrete Time Markov Chain). A Discrete Time Markov Chain (MC for short) is a tuple M = (S, s 0 , p, L), where S is a finite set of states containing the initial state s 0 , L : S → 2 A is a labelling function, and p : S × S → [0, 1] is a probabilistic transition function such that for all s ∈ S, the function p(s, .) is a distribution (i.e., s → p(s, s ) ∈ Dist(S)).

The labelling function L of a MC M = (S, s 0 , p, L) associates to each state s ∈ S a set of atomic propositions. Given a state s ∈ S, L(s) is called the label of s. We write MC for the set containing all discrete time Markov chains.

A Markov Chain can be represented as a directed graph where the nodes correspond to the states of the MC and the edges are labelled with the probabilities given by the transition function of the MC. In this representation, a missing transition between two states represents a transition probability of zero. As usual, given a MC M , we call a run of M a sequence of states obtained from executing M , i.e., a sequence w = s 1 , s 2 , . . . such that the probability of taking the transition from s i to s i+1 is strictly positive: p(s i , s i+1 ) > 0, for all i. A run w is finite if and only if it belongs to S * , i.e., it represents a finite sequence of transitions from M .

Example 2. Figure 1.1 illustrates the MC M 1 = (S, s 0 , p, L) ∈ MC where the set of states S is given by { s 0 , s 1 , s 2 , s 3 , s 4 }, the atomic propositions are restricted to { α, β }, the initial state is s 0 , and the labelling function L is as follows: The sequences of states w = s 0 , s 1 , s 2 , w = s 0 , s 2 , and w = s 0 , s 2 , s 2 , s 2 are three (finite) runs from the initial state s 0 to the state s 2 .

{(s 0 , ∅), (s 1 , { α }), (s 2 , { β }), (s 3 , { α, β }), (s 4 , { α })}.
In order to prove some of our results, we need a notion of bisimulation between Markov chains. We therefore use the classical notion of probabilistic bisimulation.

Definition 3 (Probabilistic bisimulation for Markov chains [BK08]). Let M = (S, s 0 , p, L) be a Markov chain. A probabilistic bisimulation on M is an equivalence relation R on S such that for all pair of states (s 1 , s 2 ) ∈ R , we have -L(s 1 ) = L(s 2 ), and p(s 1 , T ) = p(s 2 , T ) for each equivalence class T ∈ S/R .

Given two MCs M 1 and M 2 , we say that M 1 and M 2 are bisimilar if and only if there exists a probabilistic bisimulation R on M 1 ∪ M 2 with (s 1 0 , s 2 0 ) ∈ R , where M 1 ∪ M 2 is a MC consisting of all states, transitions and labels of M 1 and M 2 1 .

Reachability.

A Markov chain M defines a unique probability measure µ M over the set of infinite sequences s 0 , s 1 , . . . ∈ S ω (see [START_REF] Baier | Principles of model checking[END_REF] for details). According to this measure, the probability of a finite sequence of states w = s 0 , s 1 , . . . , s n , written P M (w) is the product of the probabilities of the transitions transitions involved in this sequence, i.e., P M (w

) = p(s 0 , s 1 )• p(s 1 , s 2 )•. . .• p(s n-1 , s n ).
Naturally, if w is not a run of M , we have P M (w) = 0. Given a MC M , the overall probability of reaching a given state s from the initial state s 0 is called the reachability probability and written P s 0 M (♦s) or P M (♦s) when clear from the context. This probability is computed as the sum of the probabilities of all finite runs starting in the initial state and reaching the state s for the first time. Formally, let reach s 0 (s) =

Markov chain abstractions

{ w ∈ S * | w = s 0 , . . . s n with s n = s and s i = s ∀0 ≤ i < n } be the set of such runs. We then define P M (♦s) = ∑ w∈reach s 0 (s) P M (w) if s = s 0 and 1 otherwise. This notation naturally extends to the reachability probability of a state s from a state t that is not s 0 , written P t M (♦s) and to the probability of reaching a state label Γ ⊆ A written P s 0 M (♦Γ). One could also extend it to the reachability probability of an atomic proposition α by considering all state labels containing α.

In the following, we say that a state s (resp. a label Γ ⊆ A) is reachable in M if and only if the reachability probability of this state (resp. label) from the initial state is strictly positive.

Example 3 (Example 2 continued). In Figure 1.1 the probability of the run w = s 0 , s 2 , s 1 , s 1 , s 3 is 0.3 • 0.5 • 0.5 • 0.5 = 0.0375 and the probability of reaching the state s 1 from s 0 is

P s 0 M 1 (♦s 1 ) = p(s 0 , s 1 )+Σ +∞ i=0 p(s 0 , s 2 )•p(s 2 , s 2 ) i •p(s 2 , s 1 ) = p(s 0 , s 1 )+ p(s 0 , s 2 )•p(s 2 , s 1 )• (1/(1- p(s 2 , s 2 ))) = 1.
Furthermore, the probability of reaching { β } corresponds to the probability of reaching the state s 2 , which is 0.3 here.

Markov chain abstractions

Modeling an application as a Markov Chain requires knowing the exact probability for each possible transition of the system. However, this can be difficult to compute or to measure in the case of a real-life application (e.g., because of precision errors or limited knowledge). In this section, we start with a generic definition of Markov chain abstraction models. Then we recall three abstraction models from the literature, respectively pMC, IMC, and pIMC, and finally we present a comparison of these existing models in terms of succinctness. A Markov chain abstraction model is a specification theory for MCs. It consists of a set of abstract models, called specifications, each of which representing a (potentially infinite) set of MCs -implementations -together with a satisfaction relation defining the link between implementations and specifications. As an example, consider the powerset of MC (i.e., the set containing all possible sets of Markov chains). Clearly, (2 MC , ∈) is a Markov chain abstraction model, which we call the canonical abstraction model. This abstraction model has the advantage of representing all possible sets of Markov chains but it also has the disadvantage that some Chapter 1 -Parametric interval Markov chains

s 0 s 1 α s 2 β s 3 α, β s 4 α 0.7 0.3 1 -p p p 1 -p 1 1 Figure 1.2 -pMC I s 0 s 1 α s 2 β s 3 α, β s 4 α [0, 1] [0, 1] [0.5, 1] [0.3, 0.5] [0, 0.6] [0.2, 0.6] [0, 0.5] 1 [0, 0.5] [0.5, 0.6] Figure 1.3 -IMC I
Markov chain abstractions are only representable by an infinite extension representation. Indeed, recall that there exist subsets of [0, 1] ⊆ R which cannot be represented in a finite space (e.g., the Cantor set [START_REF] Cantor | Überunendliche, lineare punktmannig faltigkeiten v [on infinite, linear point-manifolds (sets)[END_REF]). We now present existing MC abstraction models from the literature.

Existing MC abstraction models

Parametric Markov Chain is a MC abstraction model from [START_REF] Alur | Parametric real-time reasoning[END_REF] where a transition can be annotated by a rational function over parameters. We write pMC for the set containing all parametric Markov chains.

Definition 5 (Parametric Markov Chain). A Parametric Markov Chain (pMC for short) is a tuple I = (S, s 0 , P, L, X) where S, s 0 , and L are defined as for MCs, X is a set of variables (parameters) ranging over [0, 1], and P : S × S → Q X associates with each potential transition a parameterized probability.

Let M = (S, s 0 , p, L) be a MC and I = (S , s 0 , P, L , X) be a pMC. The satisfaction relation |= p between MC and pMC is defined by M |= p I if and only if S = S , s 0 = s 0 , L = L , and there exists a valuation v of X such that p(s, s ) equals v(P(s, s )) for all s, s in S.

Example 4. Figure 1.2 shows a pMC I = (S, s 0 , P, L, X) where S, s 0 , and L are identical to those of the MC M from Figure 1.1, the set X contains only one variable p, and the parametric transitions in P are given by the edge labelling (e.g., P(s 0 , s 1 ) = 0.7, P(s 1 , s 3 ) = p, and P(s 2 , s 2 ) = 1p). Note that the pMC I is a specification containing the MC M 1 from Figure 1.1 (using the valuation v such that v(p) = 0.5).

Interval Markov Chains extend MCs by allowing to label transitions with intervals of possible probabilities instead of precise probabilities. We write IMC for the set containing all interval Markov chains.

Markov chain abstractions

Definition 6 (Interval Markov Chain [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF]). An Interval Markov Chain (IMC for short) is a tuple I = (S, s 0 , P, L), where S, s 0 , and L are defined as for MCs, and P : S × S → I associates with each potential transition an interval of probabilities.

Example 5. Figure 1.3 illustrates an IMC I = (S, s 0 , P, L) where S, s 0 , and L are similar to the MC given in Figure 1.1. By observing the edge labelling we see that P(s 0 , s 1 ) = [0, 1], P(s 1 , s 1 ) = [0.5, 1], and P(s 3 , s 3 ) = [1, 1]. On the other hand, the intervals of probability for missing transitions are reduced to [0, 0], e.g., P(s 0 , s

0 ) = [0, 0], P(s 0 , s 3 ) = [0, 0], P(s 1 , s 4 ) = [0, 0].
In the literature, IMCs have been mainly used with three distinct semantics: at-every-step, interval-Markov-decision-process, and once-and-for-all. All these semantics are associated with distinct satisfaction relations which we now introduce.

The once-and-for-all IMC semantics [DJJ + 15, WTO + 11, PLSVS13] is alike to the semantics for pMC, as introduced above. The associated satisfaction relation |= o I is defined as follows.

Definition 7 (Once-and-for-all satisfaction). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the once-and-for-all semantics, written M |= o I I , if and only if (T,t 0 , L M ) = (S, s 0 , L I ) and for all reachable states s in M and all state s ∈ S, p(s, s ) ∈ P(s, s ).

In this sense, we say that MC implementations using the once-and-for-all semantics need to have the same structure as the IMC specification. Remark that this definition only targets reachable states in the MC. Indeed, some states of the IMC could have inconsistent outgoing transitions. While these states cannot be "satisfied" in MC implementations, satisfaction should still be possible if they are unreachable.

Next, the interval-Markov-decision-process IMC semantics (IMDP for short) [START_REF] Chatterjee | Model-checking omega-regular properties of interval Markov chains[END_REF][START_REF] Sen | Model-checking Markov chains in the presence of uncertainties[END_REF] operates as an "unfolding" of the original IMC by picking each time a state is visited a possibly new probability distribution which respects the intervals of probabilities. Thus, this semantics allows one to produce MCs satisfying IMCs with a different structure. Formally, the associated satisfaction relation |= d I is defined as follows.

Definition 8 (Interval-Markov-decision-process satisfaction). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the IMDP semantics, written for the identity mapping π). Note that in [START_REF] Chatterjee | Model-checking omega-regular properties of interval Markov chains[END_REF][START_REF] Sen | Model-checking Markov chains in the presence of uncertainties[END_REF] the authors allow the Markov chains satisfying the IMCs w.r.t. |= d I to have an infinite state space. In this work we consider Markov chains with a finite state space only, therefore MC implementations with respect to the IMDP semantics consist in a finite unfolding of the specification IMC, which corresponds to a transient IMDP semantics followed by a once-and-for-all semantics. Finally, the at-every-step IMC semantics, first introduced in [JL91], operates as a simulation relation based on the transition probabilities and state labels, and therefore allows MC implementations to have a different structure than the IMC specification. Compared to the previous semantics, in addition to the unfoldings this one allows to "aggregate" and "split" states from the original IMC. Formally, the associated satisfaction relation |= a I is defined as follows.

Definition 9 (At-every-step satisfaction). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the at-every-step semantics, written M |= a I I if and only if there exists a relation R ⊆ T × S such that (t 0 , s 0 ) ∈ R , and whenever (t, s) ∈ R , we have 1. the labels of s and t correspond:

L M (t) = L I (s), 2. there exists a correspondence function δ (t,s) : T → (S → [0, 1]) such that (a) ∀t ∈ T if p(t,t ) > 0 then δ (t,s) (t ) is a distribution on S (b) For all s ∈ S, (Σ t ∈T p(t,t ) • δ (t,s) (t )(s )) ∈ P(s, s ), and (c) For all (t , s ) ∈ T × S, if δ (t,s) (t )(s ) > 0, then (t , s ) ∈ R .
Remark that reachability of the states from M is not required for the at-every-step satisfaction relation. Indeed, in this case, states that are not reachable in M can be absent from the satisfaction relation, and therefore do not impact at-every-step satisfaction. Example 6 illustrates the three IMC semantics and Proposition 1 compares them. We say that an Regarding M 2 , note that two probability distributions have been chosen for the state s 1 from I . This produces two states t 1 and t 1 in M 2 and changes the structure of the graph. Thus, M 2 |= o I I . I with the mapping π(t i ) = s i for all t i and π(t 1 ) = s 1 . Remark that there is no state t 4 here, but having an unreachable state t 4 and defining π(t 4 ) = s 4 would not prevent satisfaction. Trivially, we also have M 2 |= a I I with the relation

R 2 = { (t, s) | π(t) = s }.
Finally, in the MC M 3 with state space T the state s 3 from I has been "split" into two states t 3 and t 3 and the state t 14 "aggregates" the states s 1 and s 4 from I . The relation R 3 ⊆ T × S containing the pairs (t 0 , s 0 ), (t 14 , s 1 ), (t 14 , s 4 ), (t 2 , s 2 ), (t 3 , s 3 ), and (t 3 , s 3 ) is a satisfaction relation between M 3 and I such as defined by |= a I . For instance, consider the pair (t 2 , s 2 ) ∈ R . For this pair, define the correspondence function δ (t 2 ,s 2 ) such that δ (t 2 ,s 2 ) (t 2 )(s 2 ) = 1, δ (t 2 ,s 2 ) (t 14 )(s 1 ) = 0.5, δ (t 2 ,s 2 ) (t 14 )(s 4 ) = 0.5, and δ (t 2 ,s 2 ) (t )(s ) = 0 otherwise. We can verify the following:

1. For all t such that p 2 (t 2 ,t ) > 0, δ (t 2 ,s 2 ) (t ) is a distribution on S. This is the case for t 2 and t 14 .

2. For all s ∈ S,

(Σ t ∈T 2 p 2 (t 2 ,t ) • δ (t 2 ,s 2 ) (t )(s )) ∈ P(s 2 , s ): [s 1 ] : (Σ t ∈T 2 p 2 (t 2 ,t ) • δ (t 2 ,s 2 ) (t )(s 1 )) = 0.8 * 0.5 = 0.4 ∈ [0, 0.6] = P(s 2 , s 1 ) [s 2 ] : (Σ t ∈T 2 p 2 (t 2 ,t ) • δ (t 2 ,s 2 ) (t )(s 2 )) = 0.2 * 1 = 0.2 ∈ [0.2, 0.6] = P(s 2 , s 2 ) [s 4 ] : (Σ t ∈T 2 p 2 (t 2 ,t ) • δ (t 2 ,s 2 ) (t )(s 4 )) = 0.8 * 0.5 = 0.4 ∈ [0, 0.5] = P(s 2 , s 4 ) 3. For all (t , s ) ∈ T 2 × S, if δ (t 2 ,s 2 ) (t )(s ) > 0, then (t , s ) ∈ R 2 .
s 0 ∅ s 1 α s 2 α [0, 0.4] [0.6, 1] [0, 1] [0, 1] [0, 1] [0, 1] IMC I t 0 ∅ t 1 α t 2 α 0.4 0.6 0.5 0.5 0.5 0.5 MC M 1 t 0 ∅ t 1 α t 2 α t 2 α 0.4 0.6 1 1 1 MC M 2 t 0 ∅ t 1 α t 2 α 0.5 0.5 0.5 0.5 0.5 0.5 MC M 3 Figure 1.6 -IMC I , MCs M 1 , M 2 , and M 3 such that M 1 |= a I I , M 1 |= d I I and M 1 |= o I I ; M 2 |= d I I and M 2 |= o I I ; M 3 |= a I I and M 3 |= d I I .
[0, 0.6] = P(s 2 , s 1 )).

Proposition 1. The at-every-step satisfaction relation is (strictly) more general than the interval-Markov-decision-process satisfaction relation which is (strictly) more general than the once-andfor-all satisfaction relation.

Examples illustrating this result are given in Figure 1.6 and a complete proof is provided in [BDF + 18].

Parametric Interval Markov Chains, as introduced in [DLP16], abstract IMCs by allowing (combinations of) parameters to be used as interval endpoints in IMCs. Under a given parameter valuation the pIMC yields an IMC as introduced above. pIMCs therefore allow the representation, in a compact way and with a finite structure, of a potentially infinite number of IMCs. Note that one parameter can appear in several transitions at once, requiring the associated transition probabilities to depend on one another. Let X be a finite set of parameters and v be a valuation over X. By combining notations used for IMCs and pMCs the set I(Q X ) contains all parametrized intervals over [0, 1], and for all

J = [ f 1 , f 2 ] ∈ I(Q X ), v(J) denotes the interval [v( f 1 ), v( f 2 )] if 1.3. Markov chain abstractions s 0 s 1 α s 2 β s 3 α, β s 4 α [0, 1] [0, 1] [q, 1] [0.3, q] [0, p] [0.2, p] [0, 0.5] 1 [0, 0.5] [0.5, p] Figure 1.7 -pIMC P 0 ≤ v( f 1 ) ≤ v( f 2
) ≤ 1 and the empty set otherwise2 . We write pIMC for the set containing all parametric interval Markov chains.

Definition 10 (Parametric Interval Markov Chain [DLP16]). A Parametric Interval Markov

Chain (pIMC for short) is a tuple P = (S, s 0 , P, L, X), where S, s 0 , L, and X are defined as for pMCs, and P : S × S → I(Q X ) associates with each potential transition a (parametric) interval.

In [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], we introduced pIMCs where parametric interval endpoints are limited to linear combination of parameters. In this chapter, we extend the pIMC model by allowing rational functions over parameters as endpoints of parametric intervals. Given a pIMC P = (S, s 0 , P, L, X)

and a valuation v, we write v(P ) for the IMC (S, s 0 , P v , L) obtained by replacing the transition function P from P with the function P v : S × S → I defined by P v (s, s ) = v(P(s, s )) for all s, s ∈ S. The IMC v(P ) is called an instance of pIMC P . Finally, depending on the semantics chosen for IMCs, three satisfaction relations can be defined between MCs and pIMCs. Example 7. Consider the pIMC P = (S, s 0 , P, L, X) given in Figure 1.7. The set of states S and the labelling function are the same as in the MC and the IMC presented in Figures 1.1 and 1.3 respectively. The set of parameters X has two elements p and q. Finally, the parametric intervals from the transition function P are given by the edge labelling (e.g., P(s 1 , s 3 ) = [0.3, q], P(s 2 , s 4 ) = [0, 0.5], and P(s 3 , s 3 ) = [1, 1]). Note that the IMC I from Figure 1.3 is an instance of P (by assigning the value 0.6 to the parameter p and 0.5 to q). Furthermore, as said in Example 6, the In the following, we consider that the size of a pMC, IMC, or pIMC L, written |L|, corresponds to its number of states plus its number of transitions not reduced to 0, [0, 0] or ∅. We will also often need to consider the predecessors (Pred), and the successors (Succ) of some given states. Given a pIMC with a set of states S, a state s in S, and a subset S of S, we write: ). On the other hand, IMCs with the IMDP and at-every-step semantics allow unbounded unfolding of the state space with different probabilistic choices, which implies that pMCs would need infinitely many states to represent the same set of MC implementations. Finally, we show that the set of MC implementations of an IMC with the once-and-for-all semantics can also be represented with a pMC of the same size as the original IMC. As a consequence, pMCs are strictly more succinct than IMCs with the once-and-for-all semantics. This is formalized in the following lemma. This result is illustrated in Figures 1.8 and 1.9. For the sake of conciseness, we only give some a proof-sketch. A complete proof is provided in [BDF + 18]. Let (L 1 , |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models. Recall that according to the succinctness definition (cf.

-Pred(s) = { s ∈ S | P(s , s) / ∈ { ∅, [0, 0] } } -Succ(s) = { s ∈ S | P(s, s ) / ∈ { ∅, [0, 0] } } -Pred(S ) = s ∈S Pred(s ) -Succ(S ) = s ∈S Succ(s ) 1.3.2 Abstraction model comparisons

Markov chain abstractions s

1 α s 0 β [0, 1] [0, 1] 1 IMC I s 1 α s 0 β p 1-p 1 pMC P 1 s 1 α s 2 γ s 0 β p 1 q 1 1-p 1 -q 1 1-q 1 q 1 1 pMC P 2 s 1 α s 2 α . . . . . . s n α s 0 β p 1 q 1 1-p 1 -q 1 p 2 q 2 1-p 2 -q 2 p n 1-p n 1 pMC P n
Definition 11) L 1 ≤ L 2 if there exists L 2 ∈ L 2 such that either L 1 ≡ L 2 for all L 1 ∈ L 1 , or the only way to have L 1 ≡ L 2 is with L 1 exponentially larger than L 2 .

(1-2) It is easy to see that all of the MCs given in Figure 1.9 satisfy I , given in Figure 1.8 with

|= a I and |= d I .
Because State s 1 can be "unfolded" as many times as we want, with different output probabilities, we would need pMCs such as P n from Figure 1.8 with infinitely many states to represent the same set of MC implementations.

(3-4) pMC P 2 from Figure 1.8 enforces the transitions from s 1 to s 0 and the transition from s 2 to s 0 to have the same probability. Since IMCs do not allow such dependencies between transitions, the set of MCs satisfying P 2 cannot be matched by any IMC semantics.

Chapter 1 -Parametric interval Markov chains

(5) For the same argument as above, we have (IMC, |= o I ) ≤ (pMC, |= p ). On the other hand, any transition equipped with an interval of the form [l, u] (IMC), with 0 ≤ l ≤ u ≤ 1 can be matched with a transition with probability l + (ul) • p (pMC), where p is a fresh parameter. Since this transformation preserves the size of the original IMC, we have (pMC, |= p ) < IMC, |= o I .

We now compare pMCs and IMCs to pIMCs. We show that (pIMC, |= o pI ) is equivalent to (pMC, |= p ) and pIMC is strictly more succinct than IMC for the three semantics.

Our comparison results are presented in the following proposition.

Proposition 2. The Markov chain abstraction models can be ordered as follows w.r. The proofs of items (1-5) are straightforward. Item 6 is more involved but can be obtained using a construction that derives, from any given pIMC P , a pMC I that represents the same set of MC implementations up to bisimilarity. The intuition of this construction is to use fresh parameters and several copies of the transitions in P to encode the constraints on the interval endpoints.

Details of this construction are provided in [BDF + 18].

Qualitative properties

As seen above, pIMCs are a succinct abstraction formalism for MCs. The aim of this section is to investigate qualitative properties for pIMCs, i.e., properties that can be evaluated at the specification (pIMC) level, but that entail properties on its MC implementations. pIMC specifications are very expressive as they allow the abstraction of transition probabilities using both intervals and parameters. Unfortunately, as it is the case for IMCs, this allows the expression of incorrect specifications. In the IMC setting, this is the case either when some intervals are ill-formed or when there is no probability distribution matching the interval constraints of the outgoing transitions of some reachable state. In this case, no MC implementation exists that satisfies the IMC specification. Deciding whether an implementation that satisfies a given specification exists is called the consistency problem. In the pIMC setting, the consistency problem is made more complex because of the parameters which can also induce inconsistencies in some cases. One could also be interested in verifying whether there exists an implementation that reaches some target states/labels, and if so, propose a parameter valuation ensuring this property. This problem is called the consistent reachability problem. Let P = (S,s 0 ,P,L,X) be a pIMC, we write C ∃c (P ) for the CSP produced by C ∃c according to P . Any solution of C ∃c (P ) will correspond to a MC satisfying P . In C ∃c (P ) = (X c , D c ,C c ), the sets of variables X c and domains D c are as follows: we use one variable π p with domain [0, 1] per parameter p in X; one variable θ s s with domain [0, 1] per transition (s, s ) in { { s } × Succ(s) | s ∈ S }; and one Boolean variable ρ s per state s in S. These Boolean variables will indicate for each state whether it appears in the MC solution of the CSP (i.e., in the MC satisfying the pIMC P ). Finally, we use the following constraints, for each state s ∈ S:

(1) ρ s 0 , (only for s 0 )

(2) ¬ρ s ⇔ Σ s ∈Pred(s)\{ s } θ s s = 0, (only for s = s 0 )

(3) ¬ρ s ⇔ Σ s ∈Succ(s) θ s s = 0 (4) ρ s ⇔ Σ s ∈Succ(s) θ s s = 1 (5) ρ s ⇒ θ s s ∈ P(s, s ), for all s ∈ Succ(s)
Recall that given a pIMC P the objective of the CSP C ∃c (P ) is to construct a MC M satisfying P . Constraint (1) states that the initial state s 0 appears in M . Constraint (2) ensures that for each non-initial state s, variable ρ s is set to false if and only if s is not reachable from its predecessors. Constraint (4) ensures that if a state s appears in M , then its outgoing transitions form a probability distribution. On the contrary, constraint (3) propagates non-appearing states (i.e., if a state s does not appear in M then all its outgoing transitions are set to zero). Finally, constraint (5) states that, for all appearing states, the outgoing transition probabilities must be selected inside the specified intervals. Example 8. Consider the pIMC P given in Figure 1.7. Figure 1.10 describes the variables in C ∃c (P ): one variable per transition (e.g., θ 1 0 , θ 2 0 , θ 1 1 ,. . . ), one Boolean variable per state (e.g., ρ 0 , ρ 1 ,. . . ), and one variable per parameter (π p and π q ). The following constraints correspond to the constraints (2), (3), (4), and (5) generated by our encoding C ∃c for the state 2 of P :

¬ρ 2 ⇔ θ 2 0 = 0 ¬ρ 2 ⇔ θ 1 2 + θ 2 2 + θ 4 2 = 0 ρ 2 ⇔ θ 1 2 + θ 2 2 + θ 4 2 = 1 ρ 2 ⇒ 0 ≤ θ 1 2 ≤ π p ρ 2 ⇒ 0.2 ≤ θ 2 2 ≤ π p ρ 2 ⇒ 0 ≤ θ 4 2 ≤ 0.5
Finally, Figure 1.11 describes a solution for the CSP C ∃c (P ). Note that given a solution of a pIMC encoded by C ∃c , one can construct a MC satisfying the given pIMC w.r.t. |= o I by keeping all states s such that ρ s is equal to true and considering the transition function given by the probabilities in the θ s s variables. The following proposition shows that our encoding works as expected. The proof is provided in [BDF + 18]. Proposition 3. A pIMC P is existentially consistent if and only if C ∃c (P ) is satisfiable.

We note that this existential consistency encoding is linear in the size of the pIMC, which is better than the exponential encoding we proposed in [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], which enumerates the powerset of the states in the pIMC resulting in deep nesting of conjunctions and disjunctions.

Qualitative reachability

Let P = (S, s 0 , P, L, X) be a pIMC and Γ ⊆ A be a state label. We say that Γ is existentially reachable in P if and only if there exists an implementation M of P where Γ is reachable (i.e., P M (♦Γ) > 0). In a dual way, we say that Γ is universally reachable in P if and only if Γ is reachable in any implementation M of P . As for existential consistency, we use a result from [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF] (6) ω s 0 = 1, (only for s 0 )

(P ) = (X c ∪ X , D c ∪ D ,C c ∪C ) is such that (X c , D c ,C c ) = C ∃c (P ), X
(7) ω s = 1, (only for s = s 0 )

(8) ρ s ⇔ (ω s = 0) (9) ω s > 1 ⇒ s ∈Pred(s)\{ s } (ω s = ω s + 1) ∧ (θ s s > 0), (only for s = s 0 ) (10) ω s = 0 ⇔ s ∈Pred(s)\{ s } (ω s = 0) ∨ (θ s s = 0)

Quantitative properties

We now move to the verification of quantitative reachability properties in pIMCs. Quantitative reachability has already been investigated in the context of pMCs and IMCs with the onceand-for-all semantics. In this section, we propose a major theoretical contribution: a theorem showing that the three IMC semantics are equivalent with respect to quantitative reachability, which allows the extension of all results from [WTO + 11, BLW13] to the at-every-step semantics.

Based on this result, we also extend the CSP encodings introduced in Section 1.4 in order to solve quantitative reachability properties on pIMCs regardless of their semantics.

Equivalence of all semantics w.r.t quantitative reachability

Given an IMC I = (S, s 0 , P, L) and a state label Γ ⊆ A, a quantitative reachability property on I is a property of the type P I (♦Γ) ∼ p, where 0 ≤ p ≤ 1 and ∼ ∈ { ≤, <, >, ≥ }. Such a property is verified if and only if there exists a MC M satisfying I (with the chosen semantics) such that

P M (♦Γ) ∼ p.
While the techniques we propose here work for all values of p, the techniques for qualitative reachability properties are usually more efficient when p = 0 or 1.

As explained above, existing techniques and tools for verifying quantitative reachability properties on IMCs only focus on the once-and-for-all and the IMDP semantics. However, to the 1.5. Quantitative properties

s 0 s 1 Γ s 2 Λ [0, 1] [0, 1] 1 1 IMC I s 0 s 1 Γ s 2 Λ 1 0 1 1 MC M 1 s 0 s 1 Γ s 2 Λ s 0 s 1 Γ s 2 Λ 1 0 1 1 0 1 1 1 MC M 2 Figure 1.12 -IMC I , MC M 1 ,
P M ≤ (♦Γ) ≤ P M (♦Γ) ≤ P M ≥ (♦Γ).
Informally, MC implementations for the at-every-step semantics may contain several copies of the original IMC states, each with a different probability of eventually reaching Γ. We therefore build M ≤ (resp. M ≥ ) by choosing as sole representative of a given IMC state the one from M that has the lowest (resp. highest) probability of eventually reaching Γ. This way we ensure that M ≤ (resp. M ≥ ) satisfies the original IMC with the once-and-for-all semantics and has a lower (resp. higher) probability of eventually reaching Γ than M . This is formalised in the following theorem. The proof, which is particularly intricate, can be found in [BDF + 18].

Theorem 1. Let I = (S, s 0 , P, L) be an IMC, Γ ⊆ A be a state label, ∼ ∈ { ≤, <, >, ≥ }, and 0 < p < 1. I satisfies P I (♦Γ) ∼ p with the once-and-for-all semantics if and only if I satisfies P I (♦Γ) ∼ p with the IMDP semantics if and only if I satisfies P I (♦Γ) ∼ p with the at-every-step semantics.

Thus, the three semantics |= o I , |= d I , and |= a I are equivalent with respect to the minimal and the maximal quantitative reachability properties. However, note that the equivalence for reachability properties does not trivially extend to more general properties. For instance, consider the formula P =1 ( X (Γ ∧ P =1 (♦Λ))), where Γ and Λ are two state labels. This property states that the probability of visiting a state labeled with Γ in the next step and almost certainly reaching a state labeled with Λ afterwards is 1. Figure 1.12 contains an IMC I such that this property does semantics, enforcing that the probability of visiting a state labeled with Γ as the next state is 1 imposes that the branch leading to Λ is given a probability 0, hence preventing from almost certainly visiting a state labeled with Λ afterwards (cf. the MC M 1 in Figure 1.12). On the contrary, the |= a I semantics allows one to unfold the structure, therefore to first assign a probability 1 to the transition leading to Γ and then give a probability 1 to the transition leading to Λ (cf. the MC M 2 in Figure 1.12). Note that the same would be possible under the |= d I semantics.

Constraint encodings

Note that the result from Theorem 1 naturally extends to pIMCs. We therefore exploit this result to construct a CSP encoding for verifying quantitative reachability properties in pIMCs. As in Section 1.4, we extend the CSP C ∃c , that produces a correct MC implementation for the given pIMC, by imposing that this MC implementation satisfies the given quantitative reachability property. In order to compute the probability of reaching state label Γ at the MC level, we use standard techniques from [START_REF] Baier | Principles of model checking[END_REF] that require the partitioning of the state space into three sets S =1 , S =0 , and S ? that correspond to a subset of states reaching Γ with probability 1, a subset of states from which Γ cannot be reached, and the remaining states, respectively. Once this partition is chosen, the reachability probabilities of all states in S ? are computed as the unique solution of a linear equation system (see [START_REF] Baier | Principles of model checking[END_REF], Theorem 10.19, p.766). We now explain how we identify states from S =0 , S =1 and S ? , and how we encode the linear equation system, which leads to the resolution of quantitative reachability.

Let P = (S, s 0 , P, L, X) be a pIMC and Γ ⊆ A be a state label. We start by setting

S =1 = { s | L(s) = Γ }.
We then extend C ∃r (P ) in order to identify the set S =0 . Let C ∃r (P , Γ) = (X r ∪ X , D r ∪ D ,C r ∪ C ) be such that (X r , D r ,C r ) = C ∃r (P ), X contains one Boolean variable λ s and one integer variable α s with domain [0, |S|] per state s in S, D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

(11) α s = 1, if Γ = L(s) (12) α s = 1, if Γ = L(s) (13) λ s ⇔ (ρ s ∧ (α s = 0)) (14) α s > 1 ⇒ s ∈Succ(s)\{ s } (α s = α s + 1) ∧ (θ s s > 0), if Γ = L(s) (15) α s = 0 ⇔ s ∈Succ(s)\{ s } (α s = 0) ∨ (θ s s = 0), if Γ = L(s)
Note that variables α s play a symmetric role to variables ω s from C ∃r : instead of indicating the existence of a run from s 0 to s, they characterize the existence of a run from s to a state labeled with Γ. In addition, due to constraint (13), variables λ s are set to true if and only if there exists a run with non zero probability from the initial state s 0 to a state labeled with Γ visiting s. Thus, Γ cannot be reached from states such that λ s = false. Therefore, S =0 = { s | λ s = false }, which is formalised in Proposition 5.

Proposition 5. Let P = (S, s 0 , P, L, X) be a pIMC and Γ ⊆ A be a state label.

(16) ¬λ s ⇒ γ s = 0 (17) λ s ⇒ γ s = 1, if Γ = L(s) (18) λ s ⇒ γ s = Σ s ∈Succ(s) γ s θ s s , if Γ = L(s)
As a consequence, variables γ s encode the probability with which state s eventually reaches Γ when s is reachable from the initial state and, 0 otherwise. Proposition 6. Let P = (S, s 0 , P, L, X) be a pIMC and Γ ⊆ A be a state label. Formally, let ∼ ∈ { ≤, <, ≥, > } be a comparison operator, we write ∼ for its negation (e.g., ≤ is >). This leads to the following theorem.

Theorem 2. Let P = (S, s 0 , P, L, X) be a pIMC, Γ ⊆ A be a label, p ∈ [0, 1], ∼ ∈ { ≤, <, ≥, > } be a comparison operator, and (X r , D r ,C r ) be C ∃r (P , Γ): Example 9. Due to the number of constraints in C ∃r , we illustrate our approach on a very small pIMC P given in the top left of Figure 1.13. In this example, we are interested in two properties:

-CSP (X r , D r ,C r ∪ (γ s 0 ∼ p))
"Does there exist an implementation such that P(♦{ α }) ≥ 0.5?" and "Are all implementations such that P(♦{ α }) ≤ 0.5?". These properties are not trivial because the parameter q appears both on the right side of the interval leading to 1 and on the left side of the interval leading to 2.

In the bottom of Figure 1.13, we report the constraints obtained in C ∃r (P , { α }). Some of the most trivial constraints have been removed for readability.

We first consider the CSP obtained by adding the constraint ensuring the first quantitative reachability property: (QR 1 ) : γ 0 ≥ 0.5. Solving this CSP yields a solution that allows to build a MC M , given in the top right of Figure 1.13, satisfying the quantitative reachability property.

The valuations obtained by solving the given CSP are such that π p = 0, π q = 0.5, θ 1 0 = θ 2 0 = 0.5, and, more importantly, γ 0 = 0.5.

We now consider the CSP obtained by adding the constraint checking the second quantitative reachability property: (QR 2 ) : γ 0 > 0.5. In this case, the CSP is unsatisfiable, which proves that all implementations are such that P(♦{ α }) ≤ 0.5. This is not really surprising as the probability of going to 2 is necessarily greater than the probability of going to 1.

Prototype implementation and experiments

Our results have been implemented in a prototype tool3 which generates the above CSP encodings, and CSP encodings from [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] as well. In this section, we first present our benchmark, then we evaluate our tool for the qualitative properties, and we conclude with the quantitative properties.

Benchmark

MCs have been used for many decades to model real-life applications. PRISM [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF] 4 is a reference for the verification of probabilistic systems. In particular, it is able to verify properties for MCs. As said in Section 1.2, pIMCs correspond to abstractions of MCs. PRISM references several benchmarks based on MCs 5 . Note first that we only consider in this section pIMCs with linear parametric expressions. In this context all CSPs encodings for verifying the 1.6. Prototype implementation and experiments For the first benchmark used for verifying qualitative properties, we constructed the pIMCs from existing MCs by randomly replacing some exact probabilities on transitions by (parametric) intervals of probabilities. Our pIMC generator takes 4 arguments: the MC transition function; the number of parameters for the generated pIMC; the ratio of the number of intervals over the number of transitions in the generated pIMC; the ratio of the number of parameters over the number of interval endpoints for the generated pIMC. The benchmarks used are presented in Table 1.14. We selected 5 applications from PRISM [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF]: HERMAN -the self-stabilisation protocol of Herman from [START_REF] Kwiatkowska | Probabilistic verification of herman's self-stabilisation algorithm[END_REF]; EGL -the contract signing protocol of Even, Goldreich & Lempel from [START_REF] Norman | Analysis of probabilistic contract signing[END_REF]; BRP -the bounded retransmission protocol from [START_REF] D'argenio | Reachability analysis of probabilistic systems by successive refinements[END_REF]; CROWDSthe crowds protocol from [START_REF] Shmatikov | Probabilistic model checking of an anonymity system[END_REF]; and NAND -the nand multiplexing from [START_REF] Norman | Evaluating the reliability of NAND multiplexing with PRISM[END_REF]. Each one is instantiated by setting global constants (e.g., N for the application HERMAN, L and N for the application EGL) leading to more or less complex MCs. We used our pIMC generator to generate a heterogeneous set of benchmarks: 459 pIMCs with 8 to 15 102 states, and 28 to 21 567 transitions not reduced to [0, 0]. The pIMCs contain from 2 to 250 parameters over 0 to For the second benchmark used for verifying quantitative properties we extended the NAND model from [START_REF] Norman | Evaluating the reliability of NAND multiplexing with PRISM[END_REF]. The original MC NAND model has already been extended as a pMC in [DJJ + 15], where the authors consider a single parameter p for the probability that each of the N nand gates fails during the multiplexing. We extend this model to pIMC by considering one parameter for the probability that the initial inputs are stimulated, and we have one parameter per nand gate to represent the probability that it fails. We consider 4 pIMCs with 104 to 7 392 states, and 147 to 11 207 transitions not reduced to [0, 0]. The pIMCs contain from 4 to 12 parameters appearing over 82 to 5 698 transitions.

0 ∅ 1 α 2 β [0, 1] [p, q] [q, 1] [1, 1] [1, 1] pIMC P 0 ∅ 1 α 2 β π p = 0 π q = 0.5 0 0.5 0.5 1 1 MC M satisfying P(♦{ α }) ≥ 0.5 (1) ρ 0 (2) (¬ρ 1 ⇔ (θ 1 0 = 0)) ∧ (¬ρ 2 ⇔ (θ 2 0 = 0)) (3) (¬ρ 1 ⇔ (θ 1 1 = 0)) ∧ (¬ρ 2 ⇔ (θ 2 2 = 0)) (4) (ρ 0 ⇔ (θ 0 0 + θ 1 0 + θ 2 0 = 1)) ∧(ρ 1 ⇔ (θ 1 1 = 1)) ∧ (ρ 2 ⇔ (θ 2 2 = 1)) (5) (ρ 0 ⇒ (0 ≤ θ 0 0 ≤ 1)) ∧ (ρ 0 ⇒ (π p ≤ θ 1 0 ≤ π q )) ∧ (ρ 0 ⇒ (π q ≤ θ 2 0 ≤ 1)) (6) ω 0 = 1 (7) (ω 1 = 1) ∧ (ω 2 = 1) (8) (ρ 1 ⇔ ω 1 = 0) ∧ (ρ 2 ⇔ ω 2 = 0) (9) ((ω 1 > 1) ⇒ (ω 1 = ω 0 + 1) ∧ θ 1 0 > 0) ∧ ((ω 2 > 1) ⇒ (ω 2 = ω 0 + 1) ∧ θ 2 0 > 0) (10) ((ω 1 = 0) ⇔ (ω 0 = 0) ∨ (θ 1 0 = 0)) ∧ ((ω 2 = 0) ⇔ (ω 0 = 0) ∨ (θ 2 0 = 0)) (11) α 1 = 1 (12) (α 0 = 1) ∧ (α 2 = 1) (13) (λ 0 ⇔ (ρ 0 ∧ (α 0 = 0))) ∧ (λ 1 ⇔ (ρ 1 ∧ (α 1 = 0))) ∧ (λ 2 ⇔ (ρ 2 ∧ (α 2 = 0))) (14) ((α 0 > 1) ⇒ [(α 0 = α 1 + 1) ∧ (θ 1 0 > 0)] ∨ [(α 0 = α 2 + 1) ∧ (θ 2 0 > 0)]) ∧((α 2 > 1) ⇒ [(α 2 = α 2 + 1) ∧ (θ 2 2 > 0)]) (15) (α 0 = 0) ⇔ [(α 1 = 0) ∨ (θ 1 0 = 0)] ∧ [(α 2 = 0) ∨ (θ 2 0 = 0)] (16) (¬λ 0 ⇒ (γ 0 = 0)) ∧ (¬λ 1 ⇒ (γ 1 = 0)) ∧ (¬λ 2 ⇒ (γ 2 = 0)) (17) λ 1 ⇒ (γ 1 = 1) (18) λ 0 ⇒ (γ 0 = γ 1 • θ 1 0 + γ 2 • θ 2 0 + γ 0 • θ 0 0 ) (QR 1 ) γ 0 ≥ 0.5 (QR 2 ) γ 0 > 0.

Constraint modeling

Given a pIMC in a text format our tool produces the desired CSP according to the selected encoding (i.e., one from [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], C ∃c , C ∃r , or C ∃r ). Recall that our benchmark only consider linear parametric expressions on transitions. The choice of the contraint programming language for implementing a CSP encoding depends on its nature (e.g., the type of the variables: integer vs. continuous, the kind of the contraints: linear vs. non-linear). Table 1.16 summarizes the nature of the encodings where SotA stands for the encoding from [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] answering the existential consistency problem. Thus, SotA, C ∃c , and C ∃r can be implemented as Mixed Integer Linear Programs (MILP) [START_REF] Pablo | Mixed integer linear programming formulation techniques[END_REF] and as Satisfiability Modulo Theory (SMT) programs [START_REF] Barrett | Satisfiability modulo theories[END_REF] with QF_LRA logic (Quantifier Free Linear Real-number Arithmetic). This logic deals with Boolean combinations of inequations between linear polynomials over real variables. Note that, QF_LRA does not deal with integer variables. Indeed logics mixing integers and reals are harder than those over reals only. However, all integer variables in our encodings can be replaced by real-number variables 6 . Each integer variable x can be declared as a real variable whose real domain bounds are its original integer domain bounds; we also add the constraint x < 1 ⇒ x = 0. Since we only perform incrementation of x this preserves the same set of solutions (i.e., ensures integer integrity constraints). Finally, due to the non-linear constraints in C ∃r , these encodings are implemented as SMT programs [START_REF] Barrett | Satisfiability modulo theories[END_REF] with the QF_NRA logic (Quantifier Free Non linear Real-number Arithmetic). We use the same technique as for C ∃c and C ∃r for replacing integer variables by real-number variables. We chose the programming language Python for implementing our CSP modeler. We do not evaluate any arithmetic expression while generating CSPs, and numbers in the interval endpoints of the pIMCs are read as strings and no trivial simplification is performed while modeling. We do so to avoid any rounding of the interval endpoints when using floating point numbers.

Experiments have been realized on a 2.4 GHz Intel Core i5 processor. Time out has been set to 10 minutes. Memory out has been set to 2Gb. Table 1.17 presents the average size (i.e., the number of variables and the number of constraints) of the instances considered for each set of benchmarks introduced in Table 1.14, as well as the average encoding and solving time for the existential consistency problem using (1) SMT SotA encoding, (2) SMT C ∃c encoding, and (3) MILP C ∃c encoding. First, note that all pIMCs are successfully compiled when using our C ∃c encoding while the SotA encoding produces out of memory errors for 4 sets of benchmarks: more than 20% of the instances (see OM cells in Table 1.17). We recall that the SotA encoding is defined inductively, and that it iterates over the power set of the states. In practice, this implies deep recursions joined with enumeration over the power set of the states. The exponential gain exposed in Section 1.4 is visible in terms of number of variables and constraints in Table 1.17, and in terms of encoding time in Figure 1.18. Each dot in Figure 1.18 corresponds to one instance of our benchmark. While the encoding time ranges between 0 and 1s when using the C ∃c encoding, it varies between 0 and 500s when using the SotA encoding (if it does not run out of memory).

MILP

Solving

We chose Z3 [START_REF] De | Z3: An efficient smt solver[END_REF] in its last version (v. 4.4.2) as SMT solver. We chose CPLEX [cpl10] in its last version (v. 12.6.3.0) as MILP solver. Both solvers have not been tuned and we use their default strategies. Experiments have been realized on a 2.4 GHz Intel Core i5 processor. Time out has been set to 10 minutes. Memory out has been set to 2Gb.

Table 1.17 presents the resolution time for the existential consistency problem on our first benchmark using (1) SMT SotA encoding, (2) SMT C ∃c encoding, and (3) MILP C ∃c encoding. While the SotA CSPs are larger than the C ∃c CSPs, the solving time for the SotA CSPs appears to be competitive compared to the solving time for the C ∃c CSPs. The scatter plot in Figure 1.19 (logarithmic scale) compares solving times for the SMT SotA encoding and SMT C ∃c encoding. However when considering the resolution time of the problem (i.e., the encoding time plus the solving time) the C ∃c encoding clearly answers faster than the SotA encoding. Finally, the comparison between the solving time using SMT C ∃c encoding and MILP C ∃c encoding is illustrated in Figure 1.20. It shows that Both the loss of safety (passing from real numbers with Z3 SMT resolution to floating point numbers with CPLEX MILP resolution) and the fact that CPLEX is highly optimized for MILP problems whereas Z3 is a "generic" solver lead to a non negligible gain in terms of resolution time (near to an exponential gain in our benchmark). Indeed the SMT C ∃c encoding requires 50 seconds to complete the solving process while the MILP C ∃c encoding needs less than 5 seconds for the same instances.

Table 1.21 summarizes the results w.r.t. our second benchmark: the pIMC sizes (in terms of states, transitions, and parameters), the CSP sizes (in terms of number of variables and constraints), and the resolution time using the Z3 solver. Note first that we perform pre-processing when verifying reachability properties: using a simple graph analysis, we eliminate some states that cannot, trivially, reach the goal states. This explains why C ∃r has less variables and constraints than C ∃c . Finally, note the order of magnitude between the resolution time required for solving the qualitative properties vs. the quantitative properties w.r.t. our encodings. Indeed, we did not succeed in solving pIMCs with more than 300 states and 400 transitions for quantitative properties while we verified pIMCs with more than 10 000 states and 20 000 transitions in the qualitative context.

Conclusion

In this chapter, we have compared several Markov Chain abstractions in terms of succinctness and we have shown that Parametric Interval Markov Chain is a strictly more succinct abstraction formalism than other existing formalisms such as Parametric Markov Chains and Interval Markov Chains. In addition, we have proposed constraint encodings for checking several properties over pIMC. In the context of qualitative properties such as existential consistency or consistent reachability, the size of our encodings is significantly smaller than other existing solutions. In the quantitative setting, we have compared the three semantics for IMCs and pIMCs and showed that the semantics are equivalent with respect to quantitative reachability properties. As a side effect, this result ensures that all existing tools and algorithms solving reachability problems in IMCs under the once-and-for-all semantics can safely be extended to the IMDP and at-every-step semantics with no changes. Based on this result, we have then proposed CSP encodings addressing quantitative reachability in the context of pIMCs regardless of the chosen semantics. Finally, we have developed a prototype tool that automatically generates our CSP encodings and that can be plugged to any constraint solver accepting the SMT-LIB format as input.

Perspectives

The results presented in this chapter could be extended in several manners. First, our tool for pIMC verification could be extended in order to manage other, more complex, properties (e.g., supporting the LTL language in the spirit of what Tulip [WTO + 11] does). However, as shown in Section 1.5.1, the techniques we use, based on the equivalence of the three IMC semantics for reachability properties, cannot be easily extended to more general properties. Some work is therefore needed to find alternative techniques when considering other properties than reachability.

In this line of work, we have investigated the reward-bounded reachability problem in the context of Parametric Markov Reward Models (i.e., pIMC where each state is associated to a reward). We have shown that in the context of model checking rPCTL, an extension of PCTL where each path formula is equipped with a specification of a bound on the accumulated reward, the three semantics (once-and-for-all, IMDP and at-every-step) are not equivalent. However, we have also shown that the at-every-step and IMDP semantics are equivalent when restricted to reward-bounded reachability properties. Although this does not provide the same simplification 1.8. Perspectives as if they were equivalent to the once-and-for-all semantics (allowing to preserve the structure), we have still proposed an algorithm for the verification of reward-bounded reachability properties using a reduction to the model checking problem for parametric Markov Chains. These results have been presented in a paper that is under evaluation at the moment. The consistency and reachability problems could also be extended to other type of interval-based specifications. In particular, we have also studied these problems in the context of parametric Interval Probabilistic Timed Automata. In this context, parameters can range both on timing constraints (as is the case in parametric Timed Automata) or on interval endpoints (as is the case here). In [AD16, ADF20], we have restricted ourselves to parameters only appearing in timing constraints, leaving the intervals free of parameters. With this restriction, we used the results presented in this chapter to propose an algorithm solving the problem in the non-parametric case. We have shown that the problem is undecidable in the parametric setting, but we have proposed a syntactic condition on the use of parameters that ensures decidability along with a construction that solves the problem in this case. Considering the full class where parameters can range both on timing parameters and interval endpoints remains to be done. Another line of extension could be to consider a parametric extension of a more expressive formalism than Interval Markov Chains. For instance, it would be natural to consider a parametric extension of Constraint Markov Chains [CDL + 11]. Our intuition is that the results presented here would naturally extend to this setting. The only difficulty would be to make sure that the equivalence result on the three semantics is still valid in this context. Finally, we have only focused until now on existential problems, i.e., determining whether there exists parameter values ensuring certain properties. We have also proposed algorithms and techniques for synthesizing solutions to these problems. However, a technique for representing and analyzing the set of all solutions to these problems still eludes us, and would be of particular importance in the context of understanding or optimizing the systems under study. CHAPTER 2

A PROBABILISTIC EXTENSION FOR EVENT-B

The work we present in this chapter was achieved during Mohamed Amine Aouadhi's thesis. This work was initiated from the observation that many software systems are developed using specification and refinement techniques using proof-based modeling formalisms such as B [START_REF] Abrial | The B-Book: Assigning programs to meanings[END_REF], Event-B [Abr10], Z [START_REF] Spivey | Understanding Z: a specification language and its formal semantics[END_REF] and their associated toolsets. Although uncertainty can be taken into account in such models, the only way of expressing it is through non-determinism. As explained in the previous chapter, uncertainty is a very important notion when modeling complex systems and needs to be taken into account inside the models in order to make verification and analysis of those systems more accurate. As also explained before, uncertainty can be expressed in other, more informative ways, than non-determinism. In particular, the use of probabilities inside the modeling formalism is a way of being more precise when modeling certain kinds of uncertainties than limiting ourselves to non-determinism.

To the best of our knowledge, few attempts had been made at extending the Event-B formalism with probabilities when we started this line of work. In particular, the existing attempts had been limited in the manner in which probabilistic information could be included in Event-B models. In this work, we instead propose a probabilistic extension of Event-B where probabilities can appear in any place where choices are available, and can be introduced in the models in a progressive manner through probabilistic refinement. This chapter is mostly based on [START_REF] Amine Aouadhi | Moving from Event-B to probabilistic Event-B[END_REF][START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF]. For the sake of conciseness, some of the results obtained in [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF] as well as detailed proofs have been left out of this document. However, the interested reader can refer to [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF] for details.

Introduction

As systems become more and more complex, with randomised algorithms [START_REF] Motwani | Randomized algorithms[END_REF], probabilistic protocols [START_REF] Abrial | A mechanically proved and incremental development of ieee 1394 tree identify protocol[END_REF] or failing components, it is necessary to add new modeling features in order to take into account complex system properties such as reliability [START_REF] Villemeur | Reliability, Availability, Maintainability and Safety Assessment[END_REF], responsiveness [START_REF] Wesley | Estimating task response time with contentions for real-time distributed systems[END_REF][START_REF] Kishor S Trivedi | Recent advances in modeling response-time distributions in real-time systems[END_REF], continuous evolution, energy consumption etc. One of these features is probabilistic reasoning to introduce uncertainty in a model or to mimic randomised behaviour. Probabilistic modeling formalisms have therefore been developed in the past, mainly extending automata-based formalisms [START_REF] Stoelinga | An introduction to probabilistic automata[END_REF][START_REF] Martin L Puterman | Markov decision processes: discrete stochastic dynamic programming[END_REF]. One such probabilistic modeling formalism has been presented in Chapter 1. Abstraction [START_REF] Katoen | Abstraction of probabilistic systems[END_REF][START_REF] Dehnert | On abstraction of probabilistic systems[END_REF], refinement [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF] and model checking algorithms [START_REF] Bianco | Model checking of probabilistic and nondeterministic systems[END_REF][START_REF] Baier | Principles of model checking[END_REF] have been successfully studied in this context. However, the introduction of probabilistic reasoning in proof-based modeling formalisms has been, to the best of our knowledge, quite limited [HA12, ST96, Hoa05, Gol98, BFG + 14, HMM05, APM09, Hur03]. Translations from proof-based models are always possible. However, the use of automata-based verification in this context is inconvenient due to the state-space explosion in the translation.

Unfortunately, other sources of non-determinism than assignments have been left untouched, although the authors argue that probabilistic choice between events or parameter values can be achieved by transformations of the models that embed these choices inside probabilistic assignments. While this is unarguably true, such transformations are not trivial and greatly impede the understanding of Event-B models. Moreover, these transformations would need to be included in the refinement chain when designers need it, which would certainly be counterintuitive to engineers. In this chapter, we extend these works by proposing a probabilistic extension of Event-B and presenting some ways of introducing probabilistic reasoning within Event-B. As the design process within Event-B is based on refinement, we propose to provide a standard description of a system by a set of models related by refinement. According to this modeling process, probabilities can be introduced in several manners that are illustrated in Figure 2.1. -Models on the left side of the picture are standard Event-B models that only contain non-deterministic choices but no probabilistic information;

Non-determinism

-Models on the right side of the picture are fully probabilistic Event-B models where all choices are probabilistic;

-Models in the centre of the picture are mixed Event-B models, where both non-deterministic and probabilistic choices are present.

Depending on the system that is being developed, the development process could always stay on the left side (when the system does not have any probabilistic aspect), and therefore end in the bottom left of the picture; or the development process could move to the right side and either end in the bottom right when the system is fully probabilistic or in the middle when the system contains both probabilistic and non-deterministic aspects.

In any case, there are many ways both to add standard (non-deterministic) details in the model and to add probabilistic information. Fig. 2.1 provides three generic development processes (green, red and blue), which we detail below.

Assuming the model under development is fully probabilistic, one could consider starting with an abstract non-deterministic version of the model, then progressively refining it in a standard way until a satisfying level of details is achieved. Once enough details have been introduced, all non-deterministic choices can be refined into probabilistic choices in one shot (this last step is called probabilisation). This whole process is depicted in blue in Fig. 2.1. Obviously, one could also consider starting with the probabilisation step, therefore obtaining an abstract fully probabilistic model. From this model, details can then be introduced through probabilistic refinement. This whole process is depicted in red in Fig. 2.1. While the fully probabilistic counterpart of the standard Event-B refinement still eludes us at this point, we nevertheless propose some restricted refinement steps for fully probabilistic systems through context refinement and the introduction of new probabilistic events. Finally, the designer could decide to interleave the introduction of new details in the model with the introduction of probabilistic information. In this context, intermediate models are mixed models and one has to consider the standard refinement of mixed models, the partial probabilisation operation (that only turns some of the non-deterministic choices into probabilistic choices), the introduction of probabilistic events in a standard (non-deterministic) model and the introduction of standard (non-deterministic) events in a probabilistic or mixed model. That last development process is depicted in green in Fig. 2.1. This chapter summarizes the results obtained during Mohamed Amine Aouadhi's thesis and is based on [START_REF] Amine Aouadhi | Moving from Event-B to probabilistic Event-B[END_REF][START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF]. In these papers, we have provided the scientific foundations in order to allow all the design possibilities presented above in the Event-B framework. In particular, we have proposed some new syntactic elements for writing probabilistic and mixed Event-B models in the Event-B framework. The consistency of such models is then expressed, as in standard Event-B, in terms of proof obligations. In order to prove the correctness of our approach, we have shown that the operational semantics of such models can be expressed in terms of (potentially infinite-state) Markov Chains -for fully probabilistic models -and (potentially infinite-state) Markov Decision Processes -for mixed Event-B models -therefore resembling the LTS operational semantics of standard Event-B models.

In this chapter, we propose several operations for introducing details and probabilistic aspects in probabilistic Event-B models. 1 In particular, we focus on the introduction of new probabilistic events in a given model. In the standard Event-B setting, convergence is a required property for proving a refinement steps as soon as new events are introduced in the model. The counterpart property in the probabilistic setting is almost-certain convergence, which has already been studied in [START_REF] Mciver | Abstraction, refinement and proof for probabilistic systems[END_REF][START_REF] Thai | The development of a probabilistic B-method and a supporting toolkit[END_REF] in the context of probabilistic programs and the standard B method, and in [START_REF] Thai | Reasoning about almost-certain convergence properties using Event-B[END_REF][START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] in the context of non-deterministic Event-B models with only probabilistic assignments. While the authors of [START_REF] Mciver | Abstraction, refinement and proof for probabilistic systems[END_REF][START_REF] Thai | The development of a probabilistic B-method and a supporting toolkit[END_REF] propose hypotheses under which probabilistic while loops almost-certainly converge, these hypotheses cannot be directly applied to our setting as they would require a translation from the probabilistic Event-B setting to the standard probabilistic B setting which is not trivial. In addition, some new conditions would need to be exhibited in the probabilistic Event-B setting that ensure that the hypotheses on the standard probabilistic B setting are met. Instead, have chosen to exhibit conditions at the probabilistic Event-B level and we show that these conditions ensure almost-certain convergence of the operational semantics of the model. On the other hand, [Hoa14, HH07] focus on almost-certain convergence at the probabilistic Event-B level for probabilistic Event-B models where probabilities only appear inside probabilistic assignments, but cannot appear in the choice between enabled events or in the choice of parameter values. However, we show that the proof obligations developed in this context are not sufficient for our models. We therefore propose new sufficient conditions, expressed in terms of proof obligations, for the almost-certain convergence of a set of fully probabilistic events. While the conditions we exhibit are more constrained than those from [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] concerning events and parameters, they are also less restrictive concerning probabilistic assignments.

Finally, some of the results presented in this chapter have been implemented in a prototype plugin for Rodin, which we briefly present at the end of the chapter.

Background

We start by introducing notations for (probabilistic) transition systems and basic elements of the Event-B method that will be used throughout the chapter.

Transition systems

In the following, recall that Dist(S) denotes the set of distributions over a given set S, i.e., the set of functions δ : S → [0, 1] such that ∑ s∈S δ(s) = 1. Remark that the definitions we give here are slightly different from the definitions given in Chapter 1. Indeed, in the context of Event-B, the set of atomic propositions that label states in our probabilistic models may depend on the corresponding Event-B model and therefore have to be included in the following definitions. Moreover, probabilistic transition systems in this context (such as Markov Chains), need to take into account action names that label transitions.

Labelled Transition System [BK08]

A labelled transition system (LTS for short) is a tuple

M =(S,

Event-B

Event-B [Abr10] is a formal method used for the development of complex discrete systems. Systems are described in Event-B by means of models. For the sake of simplicity, we assume in the rest of the chapter that an Event-B model is expressed by a tuple M=( v,I( v ), V( v),Evts, Init ) where v= {v 1 . . . -Predicate (non-deterministic) assignment:

v n } is a set of variables, I ( v) is an invariant, V( v)
x :| Q( t, v,x,x ') means that the variable x is assigned a new value x' such that the predicate Q( t, v,x,x ') is satisfied.
-Enumerated (non-deterministic) assignment: Before-after predicate and semantics The formal semantics of an assignment is described by means of a before-after predicate (BA) Q( t, v,x,x ') . This BA describes the relationship between the values of the variable before (x) and after (x') the execution of an assignment. Before-after predicates are as follows:

x :∈ {E 1 ( t, v) . . . E n ( t , v)} means
-the BA of a deterministic assignment is

x'= E( t , v),
-the BA of a predicate assignment is Q( t, v,x,x ') , and

-the BA of an enumerated assignment is x'∈{E 1 ( t, v) . . . E n ( t , v)}.
Recall that the action S j ( t , v) of a given event e j may contain several assignments that are executed in parallel. Assume that v 1 . . . v i are the variables assigned in S j ( t , v) -variables

v i+1 . . . v n are thus not modified -and let Q( t, v,v 1 ,v' 1 ) . . . Q( t, v,v i ,v' i )
be their corresponding BA. Then the BA S j ( t , v, v ') of the event action S j ( t , v) is:

S j ( t , v, v ') = Q( t, v,v 1 ,v' 1 ) ∧ . . . ∧ Q( t, v,v i ,v' i ) ∧ (v' i+1 =v i+1 ) ∧ . . . (v' n =v n )
Proof obligations The consistency of a standard Event-B model is characterized by means of proof obligations (POs) which must be discharged. Discharging all the necessary POs allows to prove that the model is sound with respect to some behavioural semantics. Formal definitions of standard Event-B POs are given in [START_REF] Abrial | Modeling in Event-B: system and software engineering[END_REF]. In the following, we only recall the most important of them: (event/INV) for invariant preservation, which states that the invariant still holds after the execution of each event in the Event-B model M. Given an event e i with guard G i ( t , v) and action S i ( t , v), this PO is expressed as follows:

I ( v) ∧ G i ( t , v) ∧ S i ( t , v, v ') I ( v ') (event/INV)
Operational semantics As established in [START_REF] Bert | Construction of finite labelled transistion systems from B abstract systems[END_REF], the semantics of an Event-B model

M=( v,I( v ), V( v),Evts, Init )
is expressed in terms of a labelled transition system (LTS) M =(S, s 0 , Acts, T , AP, L) where S is a set of states, each state in S being uniquely identified by its label; s 0 ∈ S is the initial state obtained by executing the Init event; Acts is the set of actions (event names); AP is the set of atomic propositions: a set of predicates that correspond to the valuations of v and satisfy the invariant I ( v); L: S → AP is a labelling function that provides the valuations of the variables v in a given state; and T ⊆ S × Acts × S is the transition relation corresponding to the actions of the events of M.

Refinement in Event-B

In Event-B, the process of modeling systems is based on the theory of refinement. Many research work have focused in the development of the theory of refinement [BvW89, Bac89, MM06]: it appears from the literature that refinement is used in two related concepts in computer science.

The first one considers refinement as a top-down program development methodology when the system is firstly described by an abstract specification and progressively refined by other ones. Each abstract specification will be more detailed and new details can be introduced during refinement. The second concept concerns the preservation of correctness between abstract and refined specifications.

Event-B refinement supports both concepts: it is a mechanism for introducing details about the static and dynamic properties of a model while preserving correctness. For the static part, refinement in Event-B allows a detailed description of the state space by the introduction of new variables (i.e., data/context refinement [START_REF] Jifeng He | Data refinement refined[END_REF][START_REF] De Roever | Data refinement: model-oriented proof methods and their comparison[END_REF]). Concerning the dynamic aspects, refinement in Event-B also allows a more detailed description of the execution of the system by adding new events processing the new introduced variables or by giving more details on the

Background

events of an abstract model. For more details on the theory of refinement in Event-B, we refer the reader to the Action Systems formalism [START_REF] Back | Refinement calculus, part II: Parallel and reactive programs[END_REF] which has inspired the development of Event-B. In Event-B, under a number of conditions expressed as proof obligations, a concrete model N=( w, J( v, w),Evts c , Init c ) may refine an abstract model M=( v,I( v ), V( v),Evts a , Init ). In this case, w is a set of variables containing some (preserved) variables of the abstract model and some new variables introduced by refinement, J( v, w) is an invariant that must provide a relation between the (removed) abstract variables ( v) and the new concrete variables ( w), Evts c is the set of concrete events that contains both the refined events (Evts c ∩ Evts a ) and the new events introduced by refinement (Evts c \ Evts a ), and Init c is the concrete initialisation event.

Events Each abstract event from the set Evts a can be refined by one or more concrete events.

Moreover, several events from the abstract set Evts a can be refined a single one. The first case is called splitting while the second one is called merging. In this chapter, we do not address the cases of event splitting and merging, we only consider the refinement of an abstract event e a by only one concrete event e c as follows. We remark that the event e c may contain one more component Wi c ( t , ū, v, v ', w, w') which denotes a witness. A witness links the abstract parameters t and the abstract variables v' to concrete parameters ū and concrete variables w'.

As the Event-B refinement process allows a more detailed description of the execution of the system, it is necessary to introduce new events (Evts c \ Evts a ) which characterize the evolution of the new added variables during refinement.

Proof obligations In Event-B refinement, the behaviour of the concrete model must be compatible with the behaviour of the abstract one. This constraint is verified and maintained by some proofs obligations dedicated to refinement. All the refinement POs are presented in [START_REF] Abrial | Modeling in Event-B: system and software engineering[END_REF]. In the following, we recall only the most important ones.

For the refinement of an abstract event e a by a concrete event e c , the two following POs must be satisfied:

(1) Guard strengthening. The guard of e c is as least as strong as the guard of e a :

I ( v) ∧ J( v, w) ∧ G c ( ū, w) ∧ Wi c ( t, ū, v, v ', w, w') G a ( t, v) (grd/STRENGTH)
(2) Simulation. The action of e c simulates the action of e a :

I ( v) ∧ J( v, w) ∧ G c ( ū, w) ∧ S c ( ū, w, w') S a ( t , v, v ') (act/SIM)
The last ensures that when a concrete event is executed, what it does is not contradictory with what the abstract event does.

When introducing new events during refinement, it is then necessary to show that their introduction cannot prevent the system from behaving as specified in the abstract model. In particular, it is necessary to show that such new events are convergent, in the sense that they cannot keep control indefinitely: at some point the system has to stop executing new events in order to follow the behaviour specified in its abstract model.

In order to prove that a set of events is convergent in Event-B, we have to introduce a natural number expression V( v), called a variant, and show that all convergent events strictly decrease the value of this variant. As a consequence, when the variant hits zero, it is guaranteed that no convergent event can be performed. That leads to two POs to be discharged:

(1) Numeric variant. Under the guard G i ( t , v) of each convergent event e i , the variant V( v) is greater or equal to 0.

I ( v) ∧ G i ( t , v) V( v)∈NAT (event/var/NAT)
(2) Convergence. The action S i ( t , v) of each convergent event e i must always decrease the variant V( v).

I ( v) ∧ G i ( t , v) ∀ v '. S i ( t , v, v ') ⇒V( v')<V( v) (event/VAR)

Running example: Simple peer-to-peer protocol

We now introduce a running example, based on a simplified scenario of a peer-to-peer protocol inspired from Bittorent [bit]. The description of the complete case study can be found in [pri]. The model considers a set of N clients trying to download a file that has been partitioned into K blocks. Initially, no block has been downloaded by any client. The protocol ends when all the clients have successfully downloaded all the blocks. At this level of abstraction, we only consider one event (AllDL) describing in one statement the whole execution of the protocol. magicDB is a parameter chosen in such a way that for all client c ∈ 1..N and block b ∈ 1..K, we have magicDB(c →b)=finished. The substitution DB := magicDB corresponds to the (somehow magical) download of all the blocks by all the clients in one shot. Notice that in reality, the download of all the blocks of the file by all the clients is not done in one shot. It is instead made gradually by successive attempts. Introducing these attempts is the purpose of the first refinement, which we present hereafter.

Initial protocol model

Step-by-step download We now present a first refinement of the protocol. For this purpose, we enlarge the set of variables and events. The resulting model is presented in Fig 2 .3. We introduce a new variable DBin that contains a matrix which represents the state of download of each block at each iteration of the model. For each client and each block, the corresponding state could be finished -indicating that the client has successfully downloaded the block; incomingindicating that the client is currently trying to download the block; or empty -indicating that the download of the block has not yet started. Initially, as in the abstract model, no block has been downloaded by any client, therefore DBin:= ( 1..N × 1..K) × {empty}. Furthermore, we impose that each client c is not currently trying to download more than one block c, as indicated in the invariant in Fig. 2.3.

To model the download process in a step by step manner, we introduce two new events: Start1DL and Finish1DL. In the event Start1DL, a client c and a block b are chosen in a non-deterministic manner in such a way that the download of the block b by client c has not yet started; furthermore, the considered client c is not currently trying to download another block k; then, the client starts to download the block (DBin(c →b):=incoming). The event Finish1DL models that a client c terminates the downloading of a block b in a similar manner. The events Start1DL and Finish1DL As we introduce two new events, we have to show their convergence by introducing the variant given in Fig. 2.3. Each time Start1DL and Finish1DL are activated, their actions increase the numbers of incoming or finished in DBin, therefore the variant clearly decreases.

MODEL P2P 2 REFINES P2P 1 VARIABLES DB DBin INVARIANTS DBin ∈ 1..N × 1..K → {empty,incoming,finished} ∧ ∀ c . (c ∈ 1..N ⇒ card({b | b ∈ 1..K ∧ DBin(c →b)=incoming}) ≤ 1) VARIANT 2 × N * K -2 × card({c →b |n∈1..N ∧ b ∈ 1..K ∧ DBin(c →b) = finished}) -card({c →b |n∈1..N ∧ b ∈ 1..K ∧ DBin(c →b) = incoming}) EVENTS Init = begin DB:=(1..N × 1..K) × {empty} DBin:=(1..N × 1..K) × {empty} end DLFinished = refines AllDL when DB=(1..N × 1..K) × {empty} ∧ DBin=(1..N × 1..K) × {finished}
Figure 2.4 gives an extract of the operational semantics of this first refinement in terms of transition system, with N=2 and K=2. In this figure, a small dot in the matrix indicates that the block has not been downloaded yet, while an empty (resp. filled) bullet represents that the block is incoming (resp. successfully downloaded). For readability purposes, transitions in this TS have been annotated with the corresponding parameter values. The indicated v corresponds to the value of the variant in the corresponding state. As expected, the value of v decreases each time a new event is performed.

Remark that, at this point, any download attempt is ultimately successful. The purpose of the Introducing failures In this refinement P2P 3 , we want to take into consideration some possible failures during the download process. More precisely, two possible failures can occur when a download is incoming: a failure can be critical (in this case, the download must be aborted) or not (in this case, the download continues). We therefore simply add to the previous Event-B model a new event FailureDL, given in Fig 2 .5. Its action non-deterministically chooses a value from {empty, incoming}: empty is chosen when the failure aborts the download and incoming is chosen otherwise. This new event has the same guard as the event FinishDB. As a consequence, both events are enabled at the same time and the choice between them is non-deterministic,

DB = . . . . , DBin = . • . . DB = . . . . , DBin = . • • . DB = . . . . , DBin = . • . • Finish1DL(1,2) Start1DL(2,1) Start1DL(2,2) DB = . . . . , DBin = . • • . DB = . . . . , DBin = . • • . Finish1DL(1,2) Finish1DL(2,1) DB = . . . . , DBin = • • • . DB = . . . . , DBin = . • • . Start1DL(1,1) Finish1DL(2,1) DB = . . . . , DBin = • • • . DB = . . . . , DBin = • • • . Finish1DL(2,1) Finish1DL(1,1) DB = . . . . , DBin = • • • . DB = . . . . , DBin = • • • • Finish1DL(1,1) Start1DL(2,2) DB = . . . . , DBin = • • • • Start1DL(2,2) DB = . . . . , DBin = • • • • Finish1DL(2,2) DB = • • • • , DBin = • • • • DLFinished v = 8 v = 7 v = 7 v = 7 v = 7 v = 6 v = 6 v = 6 v = 5 v = 5 v = 4 v = 4 v = 3 v = 3 v = 2 v = 2 v = 1 v = 0
DB = . . . . , DBin = • . . . DB = . . . . , DBin = . • . . DB = . . . . , DBin = . . • . DB = . . . . , DBin = . . . • Start1DL(1,1) Start1DL(1,2) Start1DL(2,1) Start1DL(2,2) DB = . . . . , DBin = . • . . DB = . . . . , DBin = . • • . DB = . . . . , DBin = . • . • FailureDL(1,2) FailureDL(1,2) Finish1DL(1,2) Start1DL(2,1) Start1DL(2,2) DB = . . . . , DBin = . • • . DB = . . . . , DBin = . • • . Finish1DL(1,2) Finish1DL(2,1) FailureDL(2,1) FailureDL(1,2) FailureDL(1,2)+(2,1) v = 8 v = 7 v = 7 v = 7 v = 7 v = 6 v = 6 v = 6 v = 5 v = 5
Figure 2.6 -Extract of the transition system of the second refinement of the simple P2P protocol, with N=2 and K=2

As we have introduced a new event, we need to prove its convergence. Unfortunately, due to the non-deterministic nature of the assignment in the event FailureDL, it is impossible to provide a variant that decreases regardless of its outcome. The refinement between this new model and the previous one is therefore not correct.

We will address the same problem in the probabilistic setting and prove that, in this case, the probabilistic version of FailureDL almost-certainly converges. As a consequence, unlike here, the refinement is correct in the probabilistic setting.

Fully probabilistic Event-B

Fully probabilistic Event-B

We recall that our goal is to introduce probabilistic reasoning within Event-B. In this section, we present the basic elements of syntax and semantics for fully probabilistic Event-B models. We begin by presenting the sources of non-determinism in Event-B and explaining how they can be replaced by probabilistic choices in the context of fully probabilistic models. We then present the syntax of such fully probabilistic Event-B models. In order to ensure the consistency of these models, we present the set of new POs specific to the introduced elements and how standard POs can be adapted in this context. Finally, in order to prove the correctness of our approach, we propose operational semantics of such models in terms of Markov chains.

Introducing probabilistic choices

In Event-B, non-determinism can appear in three places: the choice of the enabled event to be executed, the choice of the parameter value to be taken and the choice of the value to be assigned to a given variable in a non-deterministic assignment. To obtain a fully probabilistic Event-B model, we propose to replace all these non-deterministic choices with probabilistic ones. In the following, we go through these three sources of non-determinism and explain how to turn them into probabilistic choices.

Choice of the enabled event In standard Event-B, when several events are enabled in a given configuration, the event to be executed is chosen non-deterministically. In order to resolve this non-deterministic choice, we propose to equip each probabilistic event with a weight. In configurations where several probabilistic events are enabled, the probability of choosing one of them will therefore be computed as the ratio of its weight against the total value of the weights of all enabled events in this state. Using weights instead of actual probability values is convenient as the set of enabled events evolves with the configuration of the system. Using probability values instead would require to normalise them in all configurations. Moreover, for the sake of expressivity, we propose to express the weight W i ( v) of a probabilistic event e i as an expression over the variables v of the fully probabilistic Event-B model. The probability of executing a given event can therefore evolve as the system progresses. A probabilistic event is therefore allowed to be executed only if i) its guards is fulfilled and ii) its weight is strictly positive.

Choice of the parameter values

In standard Event-B, events can be equipped with parameters.

In each configuration where this is possible, a valuation of the parameters is chosen such that the guard G i ( t , v) of the event is satisfied. When there are several such parameter valuations, one of them is selected non-deterministically. We therefore propose to replace this non-deterministic choice by a uniform choice over all parameter valuations ensuring that the guard of the event is satisfied. The uniform distribution is a default choice but our results can be extended to any other discrete distribution.

Non-deterministic assignments Recall that non-deterministic assignments in Event-B are expressed in two forms: predicate non-deterministic assignments and enumerated non-deterministic assignments.

We propose to replace predicate non-deterministic assignments by predicate probabilistic assignments written

x:⊕Q x ( t, v,x')
Instead of choosing non-deterministically among the values of x' such that the predicate

Q x ( t , v,x '
) is true as in standard predicate non-deterministic assignments, we propose to choose this new value using an uniform distribution. For simplicity reasons, we enforce that this uniform distribution must be discrete, and therefore that the set of values x' such that Q x ( t , v,x ') is true must always be finite. As above, the uniform distribution we propose by default could be replaced by any other discrete distribution.

We propose to replace enumerated non-deterministic assignments by enumerated probabilistic assignments written

x:= E 1 ( t, v) @ p 1 ⊕ . . . ⊕ E m ( t, v) @ p m
In this structure, the variable x is assigned the expression E i with probability p i2 . In order to define a correct probability distribution, each p i must be strictly positive and smaller or equal to 1, and they must sum up to 1. Although rational numbers are not natively handled in Event-B, we assume that the context of the model allows the use of rational numbers. In practice, that can be done by defining a "Rational" theory in Rodin using the theory plug-in providing capabilities to define and use mathematical extensions to the Event-B language and the proving infrastructure [START_REF] Butler | Practical theory extension in event-b[END_REF].

Remark that standard deterministic assignments are retained, but can also be considered as enumerated probabilistic assignments where m=1.

e i = weight W i ( v) any t where G i (t , v) then SP i (t , v)
end Refining all non-deterministic choices into probabilistic choices has side effects on the syntax of events and models. In probabilistic Event-B, we therefore propose to use the syntax on the side for a probabilistic event e i where W i ( v) is the weight of the event, G i ( t , v) is the guard of the event and SP i ( t , v) is a probabilistic action, i.e., an action consisting only of deterministic and probabilistic assignments which are executed in parallel. Remark that the Before-after predicate SP i ( t , v, v ') of such a probabilistic event will be identical to the BA of its standard (non-deterministic) counterpart.

For simplicity reasons we impose, as in standard Event-B, that the initialisation event must be deterministic. The results we present in the rest of the chapter can nevertheless easily be extended to probabilistic initialisation events.

Definition 12 (Fully Probabilistic Event-B Model). A Fully probabilistic Event-B model is a tuple M=( v,I( v ), PEvts, Init ) where v= {v 1 . . . v n } is a set of variables, I ( v) is the invariant, PEvts is a set of probabilistic events and Init is the initialisation event.

Running Example A probabilistic version of the P2P model from Section 2.3 is given in Fig. 2.7. This model has the same variables, the same invariants and the same events as the Event-B model from Figs. 2.3 and 2.5.

The events of the model P2P P are annotated with specific weights. The risk of download failures decreases with the number of successful downloads: each time a block is successfully downloaded, the weight of Finish1DL increases whereas the weight of FailureDL decreases. The weight of Start1DL models that the probability of starting a new download decreases with the number of blocks being currently downloaded.

In case of failure, we fix the probability of aborting the download to 40%. This probability is introduced in the event FailureDL by using an enumerated probabilistic assignment instead of a non-deterministic one: the variable DBin(c →b) is assigned the value empty with a probability 4 10 (the download aborts) and the value incoming with a probability 6 10 (the download continues).

Consistency

As in standard Event-B, the consistency of a fully probabilistic Event-B model is defined by means of proof obligations (POs). In this section, we therefore introduce new POs specific to fully probabilistic Event-B and explain how we adapt standard Event-B POs in order to prove the consistency of fully probabilistic Event-B models. Numeric weight For simplicity reasons, we impose that the expression W i ( v) representing the weight of a given probabilistic event must evaluate to natural numbers.

MODEL P2P P VARIABLES DB DBin INVARIANTS DB∈ 1..N × 1..K → {empty,finished} ∧ DBin ∈ 1..N × 1..K → {empty,incoming,finished} ∧ ∀ c . (c ∈ 1..N ⇒ card({b | b ∈ 1..K ∧ DBin(c →b)=incoming}) ≤ 1) VARIANT 2 × N * K -2 × card({c →b |n∈1..N ∧ b ∈ 1..K ∧ DBin(c →b) = finished}) -card({c →b |n∈1..N ∧ b ∈ 1..K ∧ DBin(c →b) = incoming}) EVENTS Init = begin DB:=(1..N × 1..K) × {empty} DBin:=(1..N × 1..K) × {empty} end DLFinished = weight N * K when DB=(1..N × 1..K) × {empty} ∧ DBin=(1..N × 1..K) × {finished}
I ( v) ∧ G i ( t , v) W i ( v) ∈ NAT (event/WGHT/NAT)
Parameter values finiteness In order to be able to use a discrete uniform distribution over the set of parameter valuations ensuring that the guard of a probabilistic event is satisfied, we impose that this set must be finite.

I ( v) finite ({ t | G i ( t , v)})
Enumerated probabilistic assignments well-definedness and feasibility In all enumerated probabilistic assignments, it is necessary to ensure that the discrete probability values p 1 . . . p n define a correct probability distribution. Formally, this leads to two POs:

(1) Probability values p i in enumerated probabilistic assignments are strictly positive and smaller or equal to 1.

0 < p i ≤ 1 (event/assign/pWD1)
(2) The sum of the probability values p 1 . . . p n in enumerated probabilistic assignments must be equal to 1.

p 1 + . . .+ p n = 1 (event/assign/pWD2)
Feasibility of enumerated probabilistic assignments is trivial: as soon as at least one expression

E i ( t , v
) is present and well-defined, it always returns a value.

Predicate probabilistic assignment well-definedness and feasibility In order to define a discrete uniform distribution over the set of values of a variable x making the predicate Q x ( t , v,x ') of the corresponding assignment satisfied, we impose that this set must be finite.

I ( v) ∧ G i ( t , v) ∧ W i ( v)>0 finite ({x' | Q x ( t , v,x ') }) (event/assign/pWD3)
Feasibility of predicate probabilistic assignments is ensured by the standard feasibility PO [START_REF] Abrial | Modeling in Event-B: system and software engineering[END_REF] inherited from Event-B. It ensures that the set {x' | Q x ( t , v,x ') } is not empty.

Modifications to standard POs

In standard Event-B, if we want to prove that an event is enabled, we need to prove that its guard is satisfied. However, in fully probabilistic Event-B, we additionally need to prove that its weight is strictly positive. We therefore modify standard and optional Event-B POs as follows.

Invariant preservation The invariant must be preserved by all enabled probabilistic events.

I ( v) ∧ G i ( t , v) ∧ W i ( v) > 0 ∧ SP i ( t, v, v ') I ( v ') (event/pINV)
Deadlock freedom In all acceptable configurations, there must exist at least one enabled probabilistic event.

I ( v) (G 1 ( t, v) ∧ W 1 ( v) > 0) ∨ . . . ∨ (G n ( t, v) ∧ W n ( v) > 0 ) (model/pDLF)
For the sake of understanding, we hereby insist on the separation between the guard of an event, which reflects the classical notion of enabledness, and the fact that its weight must be strictly positive. Obviously, one could also automatically re-write the guard of all probabilistic events in order to include the condition on its weight. This solution would allow conserving most of the standard Event-B consistency POs without modifications in the probabilistic setting.

Running example

Consider the fully probabilistic Event-B model P2P P given in Fig. 2.7: all the weight expressions are natural numbers (event/WGHT/NAT) and, for each event, the number of acceptable parameter valuations is finite (event/param/pWD). For the probabilistic enumerated assignment on the event FailureDL, each given probability is a rational p such that 0 < p ≤ 1 (event/assign/pWD1) and their sum is clearly equal to 1 (event/assign/pWD2). The invariant is always preserved by each probabilistic version of the events (event/pINV). The model P2P P is therefore consistent.

Semantics

As explained in Section 2.2.2, the operational semantics of standard Event-B models is expressed in terms of Labelled Transition Systems. In the following, we extend this work by presenting the operational semantics of fully probabilistic Event-B models in terms of Discrete Time Markov Chains (MC).

Remark that our goal, unlike in [START_REF] Tarasyuk | Reliability assessment in Event-B development[END_REF][START_REF] Tarasyuk | Integrating stochastic reasoning into Event-B development[END_REF] is not to translate our models into MCs and use standard model checking techniques to verify them. Instead, we aim at reasoning directly on fully probabilistic Event-B models and benefiting from the symbolic proof mechanism that is the signature of the Event-B approach. The following MC semantics are nevertheless introduced as a demonstration of the correctness of our approach and results.

Notations

Let M=( v, I( v ), PEvts, Init ) be a fully probabilistic Event-B model and σ be a valuation of its variables. Given a variable x ∈ v, we write [σ]x for the value of x in σ. Given an expression E( v) over variables in v, we write [σ]E( v) (or [σ]E when clear from the context) for the evaluation of E( v) in the context of σ. Given an expression E( t, v) over variables and parameters, we write [σ, θ]E( t, v) for the evaluation of E( t, v) under parameter valuation θ and variable valuation σ.

Given a probabilistic event e i with a set of parameters t and a valuation σ of the variables, we write T e i σ for the set of parameter valuations θ such that the guard of e i evaluates to true in the context of σ and θ. Formally,

T e i σ = {θ | [σ, θ]G i ( t, v) = true}.
Recall that parameter valuations are chosen uniformly on this set. We write P T e i σ for the uniform distribution on the set T e i σ .

Fully probabilistic Event-B

Given a valuation σ of the variables and a probabilistic event e i , we say that e i is enabled in σ iff (a) the weight of e i evaluates to a strictly positive value in σ and (b) either e i has no parameter and its guard evaluates to true in σ or there exists at least one parameter valuation θ such that the guard of e i evaluates to true in the context of σ and θ, i.e., T e i σ = ∅. Let e i be a probabilistic event in PEvts and let x be a variable modified by e i . Recall that x can be modified only by one assignment within the action of e i . If x is modified by an enumerated probabilistic assignment

(x := E 1 ( t, v)@ p 1 ⊕ . . . ⊕ E m ( t, v)@ p m (m ≥ 1)), then we write E e i (x)
for the set of all expressions that can be assigned to the variable x by this assignment.

E e i (x) = {E 1 ( t, v), ..., E m ( t, v)}
The probability of choosing an expression E i among all others expressions is written P e i

x (E i ) = p i . Given a probabilistic event e i , we write Var(e i ) for the set of variables in v that are modified by the action of e i , i.e., the variables that appear on the left side of an assignment in SP i ( t, v). Recall that a variable x ∈ Var(e i ) must be on the left side of either a predicate probabilistic assignment or an enumerated probabilistic assignment. Let e i ∈ PEvts be a probabilistic event, x ∈ Var(e i ) be a variable, σ, σ two valuations of the variables v and θ a valuation of the parameter values associated to the event e i such that e i is enabled in the context of σ and θ, and leads the system to σ . If x is modified by an enumerated probabilistic assignment of e i , then we write E e i (x)| σ σ,θ for the set of expressions in E e i (x) such that their evaluation in the context of σ and θ returns the value of x in the valuation σ .

Formally,

E e i (x)| σ σ,θ = {E ∈ E e i (x) | [σ, θ](E( t, v)) = [σ ]
x} If e i is not equipped with parameters, then this subset is written

E e i (x)| σ σ .
If x is modified by a predicate probabilistic assignment (x : ⊕Q x ( t, v, x )), then we write V e i θ,σ (x) for the set of values x that make the predicate Q x ( t, v, x ) true when evaluated in σ and θ.

V e i θ,σ (x) = {x | [σ, θ]Q x ( t, v, x ) = true}
If e i is not equipped with parameters, then this subset is written V e i σ (x). Let e i be a probabilistic event and let x be a variable in Var(e i ), given an original valuation σ of the variables, a valuation θ of the parameters of e i and a target valuation σ of the variables, we write P e i σ,θ (x, σ ) for the probability that x is assigned the new value [σ ]x when executing e i from the valuation σ and with parameter valuation θ. If e i is not equipped with parameters, this is written P e i σ (x, σ ). In the following, we always use the more general notation and assume that it is replaced with the specific one when there are no parameters. Formally, this probability is given by:

(1) if x is modified by an enumerated probabilistic assignment, then:

P e i σ,θ (x, σ ) = ∑ E∈E e i (x)| σ σ,θ P e i x (E)
(2) if x is modified by a predicate probabilistic assignment, then:

P e i σ,θ (x, σ ) = 1 card(V e i θ,σ (x)) if [σ ]x ∈ V e i
θ,σ (x) and 0 otherwise.

MC operational semantics

Informally, the operational semantics of a fully probabilistic Event-B model

M=( v, I( v ), PEvts, Init ) is a Probabilistic LTS [[M]] = (S,
Acts, P, s 0 , AP, L) where the states, labels, actions, atomic propositions and initial state are similarly obtained as for the standard LTS semantics of Event-B. The only difference with the standard LTS semantics is that the transitions are equipped with probabilities, which we explain below. In the following, we identify the states with the valuations of the variables defined in their labels.

Intuitively, the transition probabilities are obtained as follows: Let e i ∈ PEvts be a probabilistic event, x ∈ v be a variable and s, s be two states of [[M]] such that (s, e i , s ) is a transition in the standard LTS semantics, i.e., where e i is enabled in s and there exists a parameter valuation θ ∈ T e i s , if any, such that the action of e i may take the system from s to s . The probability assigned to transition (s, e i , s ) is then equal to the product of (1) the probability that the event e i is chosen from the set of enabled events in state s, (2) the probability of choosing each parameter valuation θ, and (3) the overall probability that each modified variable is assigned the value given in s under parameter valuation θ.

Definition 13 (Fully Probabilistic Event-B operational Semantics). The operational semantics of a fully probabilistic Event-B model M=( v, I( v ), PEvts, Init ) is a PLTS [[M]] = (S, Acts, P, s 0 , AP, L) where S, Acts, s 0 , AP, and L are defined as in the standard LTS semantics of Event-B models (see Section 2.2.2), and P : S × Acts × S → [0, 1] is the transition probability function such that for a given state s, for all e i , s ∈ Acts × S, we have P(s, e i , s ) = 0 if e i / ∈ Acts(s) or ∃x ∈ X\{Var(e i )} st [s]x = [s ]x and otherwise 2.4. Fully probabilistic Event-B P(s, e i , s ) =

[s]W i ( v)

∑ e j ∈Acts(s) [s]W j ( v) (1)
× ∑ θ∈T e i s P T e i s (θ)

(2) × ∏

x∈Var(e i ) P e i s,θ (x, s )

As expected, the following proposition shows that the semantics of a fully probabilistic Event-B model as defined above is indeed a MC. The proof is available in [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF].

Proposition 7. The operational semantics of a fully probabilistic Event-B model M satisfying the POs given in Section 2.4.2 is a MC. We only present this detailed construction to illustrate the operational semantics of our model as defined above. This MC will not be used within the design process in probabilistic Event-B.

P2PP DB = . . . . , DBin = . . . . 1|Start1DL 1 4 |(1,1) 1 4 |(1,2) 1 4 |(2,1) 1 4 |(2,2) DB = . . . . , DBin = • . . . DB = . . . . , DBin = . • . . DB = . . . . , DBin = . . • . DB = . . . . , DBin = . . . • 1 1 1 1 1 2 |FailureDL 1 8 |Finish1DL 3 8 |Start1DL 1|(1,2) 1|(1,2) 1 2 |(2,1) 1 2 |(2,2) DB = . . . . , DBin = . • . . DB = . . . . , DBin = . • • . DB = . . . . , DBin = . • . • 1 1 1 6 10 4 10 1 5 |Finish1DL 4 5 |FailureDL 1 2 |(1,2) 1 2 |(2,1) 1 2 |(1,2) 1 2 |(2,1)
We now explain how some of the probability values in the MC from Fig. 2.8 are computed. We only focus on interesting (and complex) examples and leave out the rest of the computation to the reader. From the state referenced (s 1 ) on Fig. 2.8, three events are enabled:

-FailureDL with a weight value of 4. The probability of choosing this event is therefore 4 4+3+1 = 1 2 ; -Start1DL with a weight value of 3. The probability of choosing this event is therefore 3 4+3+1 = 3 8 ; -Finish1DL with a weight value of 1. The probability of choosing this event is therefore 1 4+3+1 = 1 8 . When choosing the event Finish1DL, only one valuation for the parameters is possible with a probability of 1. The corresponding action is deterministic, and therefore executed with probability 1. The global probability of leaving the state (s 1 ) using the event Finish1DL is therefore 1 8 × 1 × 1 = 1 8 . When choosing the event Start1DL, two possible valuations for the parameters are possible, with a probability 1 2 for each of them. The action corresponding to event Start1DL is deterministic, and the parameter valuation (2, 1) allows to reach state (s 2 ). As a consequence, the global probability of reaching (s 2 ) from (s 1 ) using the event Start1DL is 3 8 × 1 2 × 1 = 3 16 . When choosing the event FailureDL, only one valuation for the parameters is possible. The corresponding action is probabilistic, leading to two different states (with probabilities 6 10 of going back to (s 1 ) and 6 10 of going to another state). As a consequence, the global probability of looping on (s 1 ) using the event FailureDL is 1 2 × 1 × 6 10 = 3 10 . From the state referenced as (s 2 ) on Fig. 2.8, we only focus on one interesting transition. Among the two events that can be enabled, we consider the event FailureDL: the probability of choosing this event is 4 5 . Two possible valuations for the parameters are then possible, with a probability of 1 2 for each of them. Then, for each parameter valuation, the corresponding action is probabilistic, leading to different states. What makes this transition interesting is that for different parameter valuations, some actions lead to the same state: 2.5 Introducing probabilities in event-B models through refinement

Our main goal is to enable modeling probabilistic behaviours within Event-B. As explained earlier (and illustrated in Fig. 2.1), this can be done in several ways while preserving the refinement-based approach inherent to Event-B. The first way is to transform a standard (group of) event(s) into a probabilistic (group of) event(s) in a refinement step. This transformation therefore allows to turn a standard or mixed Event-B model into a mixed or fully probabilistic Event-B models. This transformation is straightforward but subject to certain conditions that we call probabilisation feasibility POs. These conditions ensures that the resulting model will be consistent, e.g., that the set of parameter values such that the guard of an event subject to probabilisation is finite. The study of this transformation has been left out of this document for the sake of conciseness, but can be found in [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF].

A second way of introducing probabilistic information in the design process is by introducing new probabilistic events through (probabilistic) refinement. One principal aspect of refinement in Event-B is the addition within a refinement step of new variables and new events acting on those variables. In this section, we explain how to introduce new probabilistic events in a given (abstract) model. Regardless of the type of model (non-deterministic, mixed or fully probabilistic), it is necessary to show that the introduction of these new events cannot prevent the system from behaving as specified in the abstract model. Recall that this is usually done by proving that the set of new events introduced in the refinement step converges, i.e., that events from this set cannot keep control indefinitely. As a consequence, at some point, the system must stop the execution of new events in order to execute the behaviour proposed in the abstract model.

We therefore propose a solution in order to prove that a given set of probabilistic events almostcertainly converges. In the following, we only focus on fully probabilistic Event-B models, but all the presented results can be generalised to standard or mixed Event-B models (see [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF] for details).

In standard Event-B refinement, it is required to show that a given set of events always converges.

On the contrary, in probabilistic Event-B, it is only required to prove that a given set of probabilistic events almost-certainly converges. In other words, we are interested in showing that, in all states of the system where convergent probabilistic events can be executed, the probability of eventually taking a non-convergent event or reaching a deadlock is 1 (i.e., the probability of infinitely executing convergent events is 0). This property has already been investigated in [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] and [START_REF] Thai | Reasoning about almost-certain convergence properties using Event-B[END_REF], in the context of events having probabilistic actions but where non-determinism is still present between events. In this context, Hallerstede and Hoang propose in [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] sufficient conditions for a set of events to almost-certainly converge. These conditions can be summarised as follows: As in standard Event-B, one needs to exhibit a natural number expression V( v) called a variant. Unlike in the standard setting, only one resulting valuation of the execution of each convergent event needs to decrease this variant. Indeed, in this case, the probability of decreasing the variant is strictly positive. Unfortunately, using such a permissive condition is not sufficient in our context: there might also be a strictly positive probability of increasing the variant. Therefore, Hallerstede and Hoang require the introduction of another natural number expression U( v) which must maximise the variant V( v) and never increase. The proposition from [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] is refined in [START_REF] Thai | Reasoning about almost-certain convergence properties using Event-B[END_REF], where

Introducing probabilities in event-B models through refinement

Hoang requires in addition that the probabilities considered in probabilistic assignments are bounded away from 0. This is ensured by requiring that the set of values that can be returned by a probabilistic assignment is finite.

Adaptation to probabilistic events

We now show how to adapt the results proposed in [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF] and [START_REF] Thai | Reasoning about almost-certain convergence properties using Event-B[END_REF] to new probabilistic events introduced in a fully probabilistic Event-B model. Since there are no non-deterministic choices between enabled events, it is not anymore necessary to require that all enabled events in a given configuration may decrease the variant. We therefore start by relaxing the condition proposed in [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF]: we only require that, in all configurations where a convergent event is enabled, there is at least one convergent event for which at least one resulting valuation decreases the variant.

(1) Almost-certain convergence. In all configurations where at least one convergent event is enabled, there must exist at least one valuation v' obtained after the execution of one of these enabled events which decreases the variant.

I ( v) ∧ ((G i ( t , v) ∧ W i ( v) > 0) ∨ . . .∨(G n ( t, v) ∧ W n ( v) > 0)) (model/pVar) (∃ v '. G i ( t , v) ∧ W i ( v)>0∧ SP i ( t, v, v') ∧ V( v')<V( v)) ∨ . . . ∨ (∃ v '. G n ( t, v) ∧ W n ( v)>0∧ SP n ( t, v, v') ∧ V( v')<V( v))
As in [START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF], we also require that convergent events can only be enabled when the variant is positive and that the variant is bounded above. In order to simplify the reasoning, we propose to use a constant bound U, as in [START_REF] Thai | Reasoning about almost-certain convergence properties using Event-B[END_REF].

(2) Numeric variant. Convergent events can only be enabled when the variant is greater or equal to 0.

I ( v) ∧ G i ( t , v) ∧ W i ( v)>0 V( v)∈NAT
(event/var/pNAT)

(3) Bounded variant. Convergent events can only be enabled when the variant is less or equal to U.

I ( v) ∧ G i ( t , v) ∧ W i ( v)>0 V( v)≤ U (event/pBOUND)
Finally, the finiteness of the set of values that can be returned by a probabilistic assignment is already ensured by the syntax for enumerated probabilistic assignments and by PO (event/as-Inadequacy of adapted POs Unfortunately, as we deal with potentially infinite-state systems, POs 1-3 presented above are not anymore sufficient for proving that the probability of eventually executing a non-convergent event or reaching a deadlock is 1. Indeed, although the probability of decreasing the variant is always strictly positive because of PO (model/pVar) and although the number of values that can be returned by a given probabilistic assignment is always finite, the combination of event weights and parameter choice can make this value infinitely small in some cases. In this case, it is well known that almost-certain reachability/convergence is not ensured. This problem is a direct consequence of the unboundedness of the weights of convergent events as well as of the number of acceptable parameter values, which, by getting arbitrarily big, cause the probability of decreasing the variant to get arbitrarily small. Two examples illustrating this fact are given below.

Example 10 (Necessity of bounding event weights). In this example, we show by means of an example of a probabilistic Event-B model the necessity of bounding the weights of new probabilistic events in order to ensure almost-certain convergence.

Consider the probabilistic Event-B model M1 and the corresponding MC semantics given in Fig. 2.10. This model has two variables: x and y and three events evt1, evt2 and evt3, two of which (evt1 and evt2) are convergent. The variant of this model is x and the bound on the variant is clearly U = 2.

In states where x = 1, only convergent events evt1 and evt2 are enabled and the local probability of choosing evt1 is 1 y while the local probability of choosing evt2 is y-1 y . In states where x = 2, only evt1 can be chosen with probability 1. In states where x = 0, the only enabled event is the (non-convergent) event evt3.

Clearly, the model M1 satisfies proof obligations (model/pVar), (event/var/pNAT) and (event/pBOUND). However, as we show below, the probability of eventually taking a non-convergent event is strictly smaller than 1 from all states where x > 0 because the probability of decreasing the variant, although strictly positive in all states, gets infinitely small from states where x = 1 as y increases.

Without loss of generality, we compute the probability of eventually taking evt3 from the initial state where x = 1 and y = 2. The reasoning starting from other states is similar. This probability is equal to the sum of (1) the probability of directly taking evt1 from (1, 2),

(2) the probability of reaching (1, 4) and taking evt1 from (1, 4),

(3) the probability of reaching (1, 8) and taking evt1 from (1, 8), . . . probability of reaching state (1, 2 i ) with i > 2 and taking evt1 from this state is strictly smaller than 1 2 i+1 . As a consequence, the probability of eventually taking evt3 from the initial state is strictly smaller than

Clearly, (1) is equal to 1 2 , (2) is equal to 1 2 • 1 4 = 1 8 , (3) is equal to 1 2 • 3 4 • 1 8 < 1
1 2 + ∞ ∑ i=2 1 2 i+1 = 3 4
Therefore, M1 does not almost-certainly converge.

The behaviour we expose here is a direct consequence of the unboundedness of the weights of convergent events, which, by getting arbitrarily big, cause the probability of decreasing the variant to get arbitrarily small.

Example 11 (Necessity of bounding event parameter values). We now use a similar example to show the necessity of bounding the number of admissible parameter values in new probabilistic events in order to prove their almost-certain convergence. The probabilistic Event-B model M2 and its corresponding MC semantics, given in Fig. 2.11 are similar to the ones presented in Fig. 2.10.

In this case also, we observe that the probability of eventually executing non-convergent event Additional Proof Obligations. We therefore adapt classical results from infinite-state MC to our setting and propose sufficient conditions in terms of proof obligations to prove the almostcertain convergence of the set of new introduced events. Informally, the following POs ensure that the probability of decreasing the variant cannot get infinitely small by requiring that both the weights of convergent events and the number of potential values given to parameters in convergent events are bounded.

(4) Bounded weight. The weight of all convergent events must be bounded above by a constant upper bound BW.

I ( v) ∧ G i ( t , v) W i ( v) ≤ BW (event/wght/BOUND)
(5) Bounded parameters. The number of potential values for parameters in convergent events must be bounded above by a constant upper bound BP.

I ( v) card({ t | G i ( t , v)}) ≤ BP (event/param/BOUND)
The following theorem shows that the conditions presented above are sufficient for guaranteeing the almost-certain convergence of a given set of events in a probabilistic Event-B model. The full proof of this theorem is particularly intricate and therefore left out of this document, but a detailed proof can be found in [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF]. In the following, we provide the intuition of the proof.

Proof Sketch

Almost certain convergence of the Event-B model is proven using its operational semantics, i.e., showing the following convergence property: the probability measure of the set of runs that eventually reach a deadlock or use a non-convergent event is 1. In order to take into account the difference between convergent and non-convergent events, we start by defining a slightly extended version of the MC semantics of M. In this version, all the states are replicated in order to "remember" the last event executed. We observe that the MC semantics of M has a potentially infinite set of states. Therefore, showing the convergence property is not trivial. In order to prove it, we use the global coarseness property introduced in [MHA07], which is a sufficient condition for the "decisiveness" of infinite-state Markov Chains. Formally, given a Markov Chain M = (S , s 0 , P ) and a target set of states F ⊆ S, we say that M is globally coarse w.r.t. F iff there exists some minimal bound α > 0 such that for all state s ∈ S, the probability of eventually reaching F from s is either 0 or greater or equal to α. It is then shown in [START_REF] Mayr | Decisive markov chains[END_REF] that whenever a Markov Chain M is globally coarse w.r.t. the set F , the probability of eventually reaching either F or a set of states F from which F cannot be reached is 1 from any state of M .

We therefore use this result on the semantics of M using a "clever" partition of the state space.

Remark. The additionnal POs imposing boundedness of weights and event parameters are therefore sufficient (in addition to the adaptation of standard convergence POs presented earlier) for proving the convergence of a given set of events. While these POs are certainly restrictive, they are easily provable and seem consistent with the requirement on the boundedness of the variant. Identifying less restrictive conditions in general is still an open question.

Recall that, as explained at the end of Sec. 2.3, we cannot prove that the peer-to-peer protocol (with failures) always terminates, because we cannot prove the convergence of the nondeterministic Event-B model P2P 3 . Indeed, the event FailureDL fails to decrease the corresponding variant, therefore preventing convergence. Now consider the fully probabilistic Event-B model P2P P given in Fig. 2.7 from Sec. 2.4. To show that the protocol always terminates we have to show that the set of events {Start1DL, Finish1DL, FailureDL} almost-certainly converges (i.e., converges with probability 1).

-The variant is numeric (event/var/pNAT) and we can take the expression 2 × N × K as a possible upper bound for the variant (event/pBOUND);

-The weights of the convergent events are bounded by N × K + 1 (event/wght/BOUND);

-The possible parameter values are bounded by N × K (event/param/BOUND); -Finally, in each possible configuration, the event Finish1DL decreases the variant with a positive probability (model/pVar).

Since the model P2P P satisfies all the required POs, Theorem 3 ensures that the set {Start1DL,

|FailureDL v = 8 v = 7 v = 7 v = 7 v = 7 v = 6 v = 6 v = 6 v = 5 v = 5
Figure 2.12 -Extract of the MC of the simple P2P protocol, with N=2 and K=2, illustrating the decrease of the variant

Conclusion

In order to illustrate that the variant indeed decreases with a positive probability from all states using probabilistic event Finish1DL, we provide in Fig. 2.12 an extract of the MC semantics of P2P P (for N=2 and K=2), where all states are labelled with the value of the variant.

Conclusion

In this chapter, we have presented an extension to the Event-B formalism that allows describing models with probabilistic aspects. We have focused on fully probabilistic models but the complete theory including also mixed Event-B models has been developped and can be found in [START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF].

In this context, we have provided proof obligations for the consistency of fully probabilistic models and expressed their operational semantics in terms of probabilistic transition systems. Moreover, we have also explained how the addition of probabilistic information can be done either as a standalone artefact (probabilisation of a standard model) or as a part of the design process that can be interleaved with standard refinement steps. In particular, we have focused on the addition of new probabilistic events in fully probabilistic models and developed sufficient conditions in terms of proof obligations in order to show that a given set of probabilistic events is almost-certainly convergent, which is a required property in this context in standard Event-B.

Most of our results have been implemented in a probabilistic plugin for the Rodin Platform. The Rodin tool [ABH + 10] is an Eclipse-based IDE for designing models in Event-B. It allows the creation of Event-B models, the automatic generation of POs and it incorporates some provers for discharging the necessary POs. Rodin is based on a set of plugins, that facilitate its extension to support new functionalities.

Our plugin is still under development, but it already supports the specification of fully probabilistic Event-B models and the generation of some dedicated POs presented in this chapter. The plugin also allows the probabilisation of a non-deterministic Event-B model: it automatically generates the corresponding probabilistic Event-B model. Once this latter is generated, the developer must complete the weight of each probabilistic event and probability values for each quantitative probabilistic assignment.

Perspectives

standstill and is still missing some essential features. We plan on taking over this implementation in the future.

On the theoretical side, although we have considered the addition of new probabilistic events in a probabilistic model, a complete counterpart to standard refinement for the probabilistic setting still eludes us. The problem mainly lies in two operations that are allowed in standard Event-B refinement: the split and merge operations.

-The split operation allows, in one refinement step, to transform one abstract event into multiple concrete events, while allowing the guards of these concrete events to be more restrictive than the guard of the abstract event. In the probabilistic setting, the problem mainly lies in the repartition of the weights of concrete events w.r.t the weight of the abstract event depending on which of the guards are satisfied. We have not found yet a satisfying solution that does not restrict in a too strict way the original split operation. Indeed, a simple but restrictive solution is to impose that the guards of the concrete events must be identical to the guard of the original abstract event. In this case, we only have to impose that the sum of the weights of the concrete events is equal to the weight of the abstract event.

-The merge operation, on the contrary, allows for a single concrete event to refine several abstract events at once. The same problem regarding the repartition of the weights of the original events arises in this setting, which has prevented us from finding a satisfying solution yet.

The verification of probabilistic properties has also been left aside for the moment, as their expression and meaning in the Event-B formalism is not trivial. The only way of proving probabilistic properties on a probabilistic Event-B model at the moment is to use its MC/MDP operational semantics, which is far from ideal. In the future, we plan on studying those two aspects (expression and verification of probabilistic properties) while limiting ourselves to the Event-B syntax and proof mechanism.

Finally, we have not yet taken advantage of the parametric information present in (probabilistic)

Event-B models. The main reason is that parameters in Event-B are used as state and/or transition variables. As a consequence, the expertise we have in the probabilistic parametric domain (as presentend in Chapters 1 and 3 for example) cannot be transfered in a trivial manner to the Event-B setting. However, two lines of future work seem promising.

-Although this would slightly restrict the expressive power of probabilistic Event-B, using distinct sets of variables for parameters used in event weights and parameters only used 2.7. Perspectives in state variable assignments would allow expressing the operational semantics of probabilistic Event-B models as parametric Markov chains, therefore enabling their verification using standard techniques in this context. The next step will be to adapt those techniques to the Event-B syntax and proof mechanism.

-Parametric statistical model checking, as will be presented in Chapter 3, is a verification technique that can be applied to any parametric system that can be simulated or sampled.

In particular, this could be used directly on systems expressed in the Event-B formalism, or on the translation of such systems to a programming language. Although Automated techniques to translate an Event-B specification to a C program already exist, they have not yet been extended to the probabilistic setting.

CHAPTER 3

STATISTICAL MODEL CHECKING FOR

PARAMETRIC SYSTEMS

The work we present in this chapter was achieved during Paulin Fournier's postdoc and Ran Bao's thesis. The contributions we present here follow up on a long series of works on statistical model checking [LDB10, BBB + 10a, BBB + 10b, BBD + 12, BBB + 12, NBB + 15] initiated more than 10 years ago. This line of work was initiated with the observation that, in many cases, standard probabilistic model checking techniques fail to scale up to the size of industrial systems and are inadequate when automata-based models of the systems are not easily available. As a consequence, more efficient and flexible alternative techniques have been developed. Statistical Model Checking (SMC) in particular is a simulation-based technique that offers many advantages in this context. Indeed, while its aim is not to compute the exact probabilities with which a system satisfies a given property, SMC allows to estimate this probability (or compare it with a fixed threshold) with a formally guaranteed precision and error-rate. Since it is simulation-based, SMC is not limited to automata-based models: it can be used on any executable model regardless of the formalism used to describe it. Moreover, because of its simulation-based nature, SMC is more dependant on the complexity of the property it verifies than on the complexity of the system itself. As a consequence, SMC can easily scale up to industrial-systems' size as long as they can be simulated efficiently.

However, SMC also suffers from limitations: in particular, it is mostly restricted to purely probabilistic systems, which prevents its use in the context of non-deterministic and/or parametric systems as considered in Chapters 1 and 2. In this chapter, we try to solve this problem by developing a version of SMC for parametric systems. In this context, instead of yielding a numeric approximation of the probability that a purely probabilistic system satisfies a property, pSMC yields a polynomial function, whose variables are the parameters of the system, that adequately approximates the probability that a parametric system satisfies the given property for all values of the parameters. Along with this probability, we provide parametric confidence intervals that formally guarantee a parametric precision. The technique is implemented in a prototype tool with two versions: one that allows the parametric verification of automata-based systems described in a standard language, and a second version that allows to verify parametric Python models. While the first version is mainly used for comparison to state-of-the-art parametric probabilistic model checkers such as PRISM [START_REF] Kwiatkowska | PRISM: Probabilistic Model Checking for Performance and Reliability Analysis[END_REF], PROPhESY [DJJ + 15] and Storm [START_REF] Dehnert | A storm is coming: A modern probabilistic model checker[END_REF], the second one has been used on a concrete case-study where we analyze the security of a UAV flight plan while taking into account the parametric precision of its sensors and parametric wind force. For the sake of conciseness some details of the case-study and some proofs of the theoretical results have been left out of the document, but the interested reader can find them in [DFL19, BAD + 19].

Introduction

Nowadays, modeling and abstracting are widely accepted as crucial steps in the understanding and study of real-life systems. In many cases, it is necessary to incorporate probabilities in the models to cope with uncertainty, to abstract complex behaviour, or to introduce randomness. Markov chains and Markov decision processes, in particular, have been widely studied as shown in Chapters 1 and 2.

Statistical model checking. Though exact verification methods are known for such models they usually require solving huge equation systems [START_REF] Baier | Principles of model checking[END_REF], and therefore have scalability issues with the biggest models. A way to avoid this complexity is to consider approximation techniques through simulation. In particular, Monte-Carlo simulation techniques used in the context of statistical model cheking [START_REF] Legay | Statistical model checking: An overview[END_REF] allow to infer the real behaviour of the system via independent simulations up to a computable precision.

Parametric Markov Chains. The values given to probabilistic transitions can have a huge impact on the behaviour of the system. In the early stages of development, it may therefore be useful to have an insight on how the values of transition probabilities affect the system in order to be able to set the best value in terms of convergence speed for example. To this purpose, parametric Markov chains have been introduced in [START_REF] Alur | Parametric real-time reasoning[END_REF]. They allow to replace the probability values given to transitions by parameter variables, and therefore to be able to give guarantees on the system for all possible values of the parameters.

Results. The aim of this chapter is extend the application of Monte-Carlo simulation to parametric Markov chains in order to approximate the probability of the considered property as a polynomial function of the parameters. In addition, we also derive a confidence interval on the obtained probabilities as a polynomial function of the parameters. Aside from using parameterized models, which comes in particular with better flexibility in the modeling, robustness of the results, and usability at the earliest stages of conception, the expected benefits of our new approach are largely those of such simulation techniques for non-parameterized Markov chains: better scalability through a reduced memory footprint, a complexity that is largely independent of the model complexity (be it in terms of size of the state-space, or of used features as long as they are executable). More specific to our approach, since we derive polynomial function approximations, where exact methods lead to rational functions, these results should be easier to post-process. Finally the complexity of our approach is largely independent of the number of parameters.

In order to experimentally confirm the interest of our approach we have implemented it in a (fairly crude) prototype in Python, and we report very encouraging results on case-studies from the literature as well as a complete analysis of a case-study provided by our industrial partner PIXIEL GROUP1 .

Related work. Model checking and parameter synthesis for parametric Markov chains have been widely studied in the last decade [DJJ + 15, HHWZ10, BDL + 17, DLP16], as previously mentioned in Chapter 1. To the best of our knowledge, all existing works on this topic focus on exact techniques which either produce constraints on the parameter values [BDL + 17, DLP16, BDF + 18] (as presented in Chapter 1) or compute the probability of satisfying a given property as a rational function of the parameters [DJJ + 15, HHWZ10]. While these techniques have the advantage of precision, they only scale to models having few parameters. We conjecture that the approximation technique we propose in this chapter will be advantageous in the context on models with a larger number of parameters because it allows to produce polynomial instead of rational functions.

On the other hand, statistical model checking [START_REF] Legay | Statistical model checking: An overview[END_REF] has, to the best of our knowledge, never been applied to parametric models as such. The closer existing techniques are reinforcement learning algorithms combined with statistical model checking, which have been applied in the context of non-deterministic and probabilistic models such as Markov decision processes [HMZ + 12] (where non-determinism could be replaced with parameters). While these techniques allow to compute the best (or worst) probability of satisfying a given property, they do not provide error precision or confidence intervals. Moreover, since they only compute an approximation 3.2. Background measure is written P M in this chapter. This gives us that for any finite run ω = s 0 , s 1 , . . . s n , P M (ω) = Π n i=1 P(s i-1 , s i ). In the rest of the chapter, we only consider finite runs.

Example 12 (Examples of properties). In this chapter, we consider properties on bounded runs and we aim at computing approximations for the following values:

Reachability P M (♦ ≤l s). Recall from Chapter 1 that a run ω = s 0 , s 1 , . . . is said to reach a state s in less than l steps, written ω |= ♦ ≤l s, if there exists i ≤ l such that s i = s.

Safety P M ( =l E). A run ω = s 0 , s 1 , . . . is said to be safe for a set of states E ⊆ S during l steps,

written ω |= =l E, if for all i ≤ l, s i ∈ E.
Expected reward E l M (r). Given a reward function r : Γ(l) → R we write E l M (r) = ∑ ω∈Γ(l) P M (ω)r(ω) for the expected value of r on the runs of length l.

Notice that for any property ϕ ⊆ Γ(l), P M (ϕ) = E l M (1 ϕ ) where 1 ϕ is the reward function such that 1 ϕ (ω) = 1 if ω ∈ ϕ and 0 otherwise. In the following of the chapter we will thus focus on properties of the form E l M (r).

Given two Markov chains M 1 = (S 1 , s 1 0 , P 1 ) and M 2 = (S 2 , s 2 0 , P 2 ) we say that M 1 and M 2 have the same structure if (S 1 , s 1 0 ) = (S 2 , s 2 0 ) and for all state s, s ∈ S 1 , P 1 (s, s ) > 0 if and only if P 2 (s, s ) > 0.

Parametric Markov chains

We now recall the context of Parametric Markov Chains, in a slightly modified version of the one given in Chapter 1.

Given a finite set of parameters X we write Poly(X) for the set of all real (multivariate) polynomials on X. Given a parameter valuation v ∈ R X and a polynomial f ∈ Poly(X), the evaluation of f under valuation v is written v( f ).

As for Markov Chains, we omit in this chapter the labeling function of Parametric Markov Chains.

In the following, a Parametric Markov chain is a tuple M = (S, s 0 , P, X) such that S is a finite set of states, s 0 ∈ S is the initial state, X is a finite set of parameters, and P : S × S → Poly(X) is a parametric transition probability function.

Remark (Rational function). Notice that in the definition of pMC given above, in contrast to the one given in Chapter 1, we restrict ourselves to real (multivariate) polynomials on X. However, all our results could naturally be extended to rational functions instead, i.e., functions of the form d/q where d, q ∈ Poly(X).

Let M be a pMC and v ∈ R X be a valuation of the parameters of M . Let P v be the transition probability function obtained under valuation v, i.e., P v (s, s ) = v(P(s, s )) for all s, s ∈ S. We say that v is a valid parameter valuation with respect to M if the tuple (S, s 0 , P v ) is a Markov chain, i.e., v defines valid probability distributions for the transitions of M . If v is a valid parameter valuation w.r.t M , the resulting Markov chain is written M v .

Remark (Consistency). Notice that one can obtain the set of all valid parameter valuations as the result of a set of constraints stating that each transition has a probability between 0 and 1 and that the sum of outgoing transition probabilities is 1 for all states. The problem asking whether a pMC admits valid parameter valuations, and computing them is called the consistency problem (see Chapter 1). In this chapter, we do not address the consistency problem further.

Given a pMC M , a run ω of M is a sequence of states s 0 , s 1 , . . . such that for all i ≥ 0, P(s i , s i+1 ) = 0 (i.e., the probability is either a strictly positive real constant or a function of the parameters). As for Markov chains we write Γ M (l) for the set of all finite runs of length l and Γ M for the set of all finite runs.

Observe that for any valid parameter valuation v, Γ M v (l) ⊆ Γ M (l) since v may assign 0 to some transition probabilities.

Example 13. A pMC M 1 = (S, s 0 , P, X) is given as an example in Figure 3.1. In this figure, we have S = { 0 . . . 4 }, s 0 = 0, and X = { p, q, r }. As depicted, some of the transitions are parametric. For instance, P(1, 0) = p. Let v be the parameter valuation such that v(p) = v(q) = 0.5 and v(r) = 0. According to the definition above, v is a valid parameter valuation for M 1 . Indeed, under this parameter valuation, all the transitions have values between 0 and 1 and the probabilities of the outgoing transitions of all states sum up to 1. The resulting Markov chain M v 1 is given in Figure 3.2.

An example of run in M 1 is ω = 0, 1, 2, 2. The length of ω is |ω| = 3. As explained above, remark that ω is not a run of M v 1 because the probability of the transition going from 1 to 2 has been set to 0 due to the parameter valuation v. On the other hand, all runs of M v 1 are also runs of M .

Approximation in parametric Markov chains

We now move to the main result of this chapter: a simulation-based method for approximate verification of parametric Markov chains based on Monte-Carlo. We start by recalling the central limit theorem, which is at the heart of our approach. 1 for parameter valuation v such that v(p) = v(q) = 0.5 and v(r) = 0 Theorem 4 (The central limit theorem (see e.g. [START_REF] Ross | A First Course In Probability[END_REF])). Let X 1 , X 2 , . . . be a sequence of independent and identically distributed random variables, each having mean γ and variance σ 2 . Then the distribution of

X 1 + • • • + X n -nγ σ √ n tends to the standard normal distribution as n → ∞. That is, for -∞ < a < ∞, lim n→∞ P X 1 + • • • + X n -nγ σ √ n ≤ a = 1 √ 2π a -∞
e -x 2 /2 dx.

We now recall standard Monte-Carlo before presenting our contribution.

Standard Monte-Carlo analysis

The aim of this chapter is to propose a statistical verification method, based on Monte-Carlo, for approximating the expected value of a given reward function r on the runs Γ(l) of a given pMC M . In order to provide some intuition, we briefly recall how standard Monte-Carlo analysis works in the context of statistical model checking of Markov chains. In this context, a set of n samples of the runs of the MC is produced. Each of these samples is evaluated, yielding a value 1 if it satisfies the desired property and 0 otherwise. According to the central limit theorem, the mean value of the samples provides a good estimator of the probability that a random run satisfies the desired property. Moreover, the central limit theorem provides a confidence interval that only depends on the number of samples (provided this number is large enough).

Unfortunately, as the transition probabilities are not known a priori in the context of pMCs, this technique cannot be applied directly (since we cannot produce samples according to the parametric transition probabilities). The method we propose in the following is in line with a technique called importance sampling (see [START_REF] Reuven | Simulation and the Monte Carlo method[END_REF] for a description). The purpose of this technique is to sample a stochastic system using a chosen probability distribution (which is not the original distribution present in this system) and "compensate" the results using a likelihood ratio in order to estimate a measure according to the original distribution. In the context of SMC, importance sampling has mainly been used in order to estimate the probability of rare events [START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF] and/or to reduce the number of required samples in order to obtain a given level of guarantee [START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF]. It has also been used in the context of parametric continuoustime Markov chains in order to estimate the value of a given objective function on the whole parameter space while using a reduced number of samples [START_REF] Bortolussi | Smoothed model checking for uncertain continuous-time markov chains[END_REF]. However, to the best of our knowledge, importance sampling has never been used in order to produce symbolic functions of the parameters as we do here. The intuition of this method is to fix the transition probabilities to an arbitrary function f , which we call normalization function, and to use these transition probabilities in order to produce samples of the pMC M . However, instead of evaluating the obtained runs by directly using the desired reward function r, we define a new (parametric) reward function r that takes into account the parametric transition probabilities. We show that, under any parameter valuation v, the evaluation of the mean value of r on the set of samples is a good estimator for the expected value of the reward r on M v . The central limit theorem also allows to produce parametric confidence intervals. Contrary to standard Monte-Carlo applied to Markov chains, the precision of the obtained approximation not only depends on the size of the sample, but also on the choice of parameter valuation as well as on the chosen normalization function. These issues are further discussed in Section 3.6.

Parametric estimation

In the following, we introduce some notations, then move to the random variable corresponding to our new parametric reward function, and finally show the correctness of our approach. Along the rest of the section, we consider a pMC M = (S, s 0 , P, X) and a reward function r :

Γ M → R. Given a function f : S×S → [0, 1] we say that f is valid w.r.t M if for all states s, ∑ s ∈S f (s, s ) = 1.
Given a valid function w.r.t M , let M f be the MC obtained from M by replacing P by f .

Examples of such valid functions include the evaluation of the parametric transition probability by a valid valuation, or the function u M such that, for each state, the distribution on its successors is uniform for all successors allowed in M (it respects the structure of M ). In the following, f is called a normalization function and, in particular, u M is called the uniform normalization function. The choice of a good normalization function is discussed in Section 3.6. Let Pa : Γ M → Poly(X) be a parametric reward function defined inductively as follows: Pa(s 0 ) = 1 and for all run ω • s • s ∈ Γ M , Pa(ω • s • s ) = Pa(ω • s)P(s, s ). Note that Pa can be seen as the counterpart of P for parametric Markov chains. Indeed for any valid valuation v and any run ω ∈ Γ M v we have

P M v (ω) = v(Pa(ω)). (3.1)
We now define the parametric reward function r that will allow us to estimate the expectation of r. Given any valid normalization function f and any run ω ∈ Γ M , let r be such that r (ω) = Pa(ω)

P M f (ω) r(ω).
We now prove our main result. Let ω ∈ Γ M f (l) be a random sample of M f and let Y be the random variable defined as follows Y = r (ω) = Pa(ω)

P M f (ω) r(ω).
The following computation shows that, under any valid parameter valuation v such that M f and M v have the same structure, we have

v(E(Y )) = E l M v (r). v(E(Y )) = v   ∑ ω∈Γ M f (l) P M f (ω)Y   (3.2) = v   ∑ ω∈Γ M f (l) P M f (ω)
Pa(ω)

P M f (ω) r(ω)   (3.3) = ∑ ω∈Γ M f (l) v(Pa(ω))r(ω) (3.4) = ∑ ω∈Γ M f (l) P M v (ω)r(ω) (3.5) = ∑ ω∈Γ M v (l) P M v (ω)r(ω) (3.6) = E l M v (r) (3.7) (3.
P M f (ω) = 0; (3.5
) is a consequence of (3.1); and finally, since Γ M v (l) = Γ M f (l) because M v has the same structure as M f we obtain (3.6).

Our adaptation of the Monte-Carlo technique for pMC is thus to estimate the expected value of Y in order to obtain a good estimator for the expectation of r.

Let ω 1 , . . . , ω n be a set of n runs of length l of M f . Let Y i be the random variable with values in Poly(X) such that Y i = r (ω i ). Notice that the Y i are independent copies of the random variable Y . The random variables Y i are therefore independent and identically distributed. Let γ be the parametric function giving their mean value and σ 2 be the parametric function giving their variance.

By the results above, for all valid parameter valuation v such that M v and M f have the same

structure, E l M v (r) = v(E(Y )) = v(E(∑ n i=1 Y i /n)) = v(γ).
Our parametric approximation of the expected value is therefore

γ = n ∑ i=1 Y i /n.

Parametric confidence intervals

We now use the central limit theorem in order to compute the confidence intervals associated to this estimation. As for the estimation γ itself, the obtained confidence intervals will be given by parametric functions.

Since the random variables (Y i ) are independent and identically distributed, each having mean γ and variance σ 2 , the expected value and variance of their sum Y 1 + • • • + Y n is as follows:

E(Y 1 + • • • + Y n ) = nγ and Var(Y 1 + • • • + Y n ) = nσ 2 .
Recall that both γ and σ are parametric functions. Let Z be the (parametric) random variable such that

Z = (Y 1 + • • • +Y n -nγ) √ nσ 2 .
By the central limit theorem, v(Z) tends toward the standard normal distribution N (0, 1), for all valid parameter valuation v such that M v and M f have the same structure, as n → ∞. Therefore, for all valid parameter valuation v such that M v and M f have the same structure, assuming that n is large enough, we obtain that

P(-z ≤ v(Z) ≤ z) ≈ ϕ ϕ ϕ(z) -ϕ ϕ ϕ(-z) = 2ϕ ϕ ϕ(z) -1, where ϕ ϕ ϕ(z) = z -∞ exp(-x 2 /2)dx/( √ 2π
) is the cumulative distribution function of the standard normal distribution N (0, 1). As a consequence, by definition of Z, we obtain that for all valid parameter valuation v such that M v and M f have the same structure, and for n large enough,

3.4. Implementation 2ϕ ϕ ϕ(z) -1 ≈ P v(γ) -z v(σ) √ n ≤ v( γ) ≤ v(γ) + z v(σ) √ n .
The (parametric) random interval

I = ( γ -zσ/ √ n, γ + zσ/ √ n) is called a confidence interval
for γ with level 2ϕ ϕ ϕ(z) -1. Typically we use z = 1.96 since 2ϕ ϕ ϕ(1.96) -1 = 0.95 (or z = 2.56 for which the level is 0.99). According to our hypothesis, this (parametric) confidence interval is only valid for valid parameter valuations v such that M v and M f have the same structure.

At this point, it is important to recall that the size of the confidence interval for γ is also parametric: it is equal to 2 • z σ √ n . Indeed, the value of σ depends both on the valuation of the parameters and on the choice of the normalization function. We now explain how to compute a parametric estimation of σ. From classical probability theory, we know that an unbiased estimator for the variance σ 2 is:

σ 2 = 1 (n -1) n ∑ i=1 Y 2 i - n (n -1) γ 2 .
For n large enough, the (parametric) confidence interval associated to our estimation γ of γ can therefore be estimated by the (parametric

) interval I = ( γ -z σ/ √ n, γ + z σ/ √ n). The (parametric)
size of this interval is therefore given by 2z σ/ √ n. Thus, for an estimation of γ with a confidence interval of level 0.95 (for a large enough value of n), the (parametric) size of the confidence interval for the estimation of the variance is 3.92 σ/ √ n. Recall that σ is a parametric function that also depends on the choice of the normalization function.

Implementation

We implemented our technique in Python to validate the approach. All the following experiments have been realized on a 2,5 GHz Intel Core i7 processor. The prototype is still in the early development stages, thus we only experimented with the uniform normalization function as well as with valid valuation normalization functions. Moreover, the size of the samples is set manually. This implementation aims at validating our approach and not at competing with tools such as PARAM [START_REF] Moritz Hahn | PARAM: A model checker for parametric Markov models[END_REF], PROPhESY [DJJ + 15] or Storm [START_REF] Dehnert | A storm is coming: A modern probabilistic model checker[END_REF].

No optimizations have been implemented and we believe that this prototype could obtain much better results with simple optimizations. The code (still in development) is available at https://github.com/paulinfournier/MCpMC.

Toy example

We first tested our program on the pMC M 1 given in example 13. The considered property is the probability of reaching state 4. The results of our program on this example are presented in Figure 3.3a. The number of simulations was set to 10 000 and the length of the run bounded to 100. The normalization function used is the uniform normalization. We also give in Figure 3.3b the graphical representation of the theoretical probability of reaching state 4. The parsing of the model took around 3ms and the simulations took around 17 seconds. The memory consumption of the program is depicted in Figure 3.4. Note that even if the program uses around 80MB of memory, most of it is due to the loading of python libraries. A careful analysis of the memory consumption shows that the actual model only uses 0.7MB, the simulations use 6.5MB and the figure uses around 4MB. The memory consumption could thus be easily reduced by using a lighter programming language with a better handling of memory such as C.

To show that our method is transparent with regard to the number of parameters we extended this model with 100 parameters {p 0 , ..., p 99 }. Since q remains the same and r is always 1 -(p i + q), this new model is equipped with a total of 101 parameters. We therefore consider an unfolding of the pMC given in Figure 3.1 of depth 101. Each time we enter a copy of state 1 or 3, new parameters are used. For example, the nth copy of states 1 and 3 use parameter q n mod 100 instead of q. The result of the experiment is shown in Figure 3.5. Since there are too many dimensions to plot in 2D, we only plot the result with respect to parameter q for random valuations of the parameters in (p i ). The number of simulations and simulation length are set as before. Note that the size and shade of the confidence intervals are not what one could expect but this is due to the fact that each point is evaluated for a different random valuation of the parameters in (p i ).

The important point of this experiment is that the results were obtained in 21 seconds, which is approximately the same as with only 2 parameters. Note that this increase in time is due to a slight increase in the complexity of the model (the states are not really duplicated but we have to keep track of the current depth in order to consider the right parameters in each simulation).

Figure 3.5 -Extension of M 1 with 101 parameters

Zeroconf

The Zeroconf model, taken from the PARAM website2 , models the management of a network.

When a new host joins the network, it randomly selects an address among K potential candidates.

If there are already m hosts in the network the probability of collision is q = m/K, which is a parameter of the model. In order to detect collisions, the new host asks the others whether the address is free. If he receives a positive answer, then the new hosts considers that his address is valid. The second parameter of the model is the probability p that the new hosts do not receive an answer. In this case, the new host retries at most n times (here n = 140). In case the new host does not receive an answer after n attempts, it decides that its address is valid. We consider the expected number of attempts until the address is considered valid (either because of positive answer or because of the absence of answers after n attempts). In Figure 3.6 we present the results of our approach for 10 000 simulations, the uniform normalization function and for a simulation length of 500. Using our prototype, the experiment took around 60s. In Figure 3.6, we present as well the result obtained with PARAM (taken from their website). Remark that the shape of the distribution we obtain is similar to the one obtained with PARAM. The size of the confidence intervals are given as the color of the points.

The crowds protocol

The Crowds protocol [START_REF] Michael | Crowds: Anonymity for web transactions[END_REF] aims at preserving anonymity of Internet users. To do so, each message is sent via random routes, with the assumption that a corrupted router can only see the local sender of a message. This protocol guarantees that the probability of a corrupted router observing the real sender (and not just routing another user's message) is small. A model of this protocol as been proposed in Prism [START_REF] Shmatikov | Probabilistic model checking of an anonymity system[END_REF] and later extended with two parameters in [START_REF] Moritz Hahn | PARAM: A model checker for parametric Markov models[END_REF]. One parameter PF represents the probability of forwarding a message to a random selected member (and therefore the probability of delivering the message directly to the receiver is 1 -PF). The other parameter badC represents the probability of a router being corrupted. The reward function we are interested in corresponds to the probability of dishonest router observing the real sender more than other crowd members. Figure 3.7 represents the output of our approach on this model for 10 000 simulations, using the uniform normalization function. This experiment was done with the same values as used on the PARAM website, that is 5 honest crowd members and 7 different path reformulates3 . The simulation time took around 6 minute and used around 80MB (again most of it due to loading of python libraries). Again, one can observe that the shape of the obtained distribution is similar to the one obtained with PARAM. The only zones where the two distributions are slightly different are those where the confidence intervals we produced are the largest.

Application to UAV flight plan analysis

One of the advantages of SMC is that the model under study does not need to be represented in a fully detailed manner as a MC or pMC: any symbolic representation that allows runs to be sampled is sufficient. As a consequence, the model does not have to be written in the prism language for instance, which has limitations on the use of real-valued variables and constants.

We have therefore developed a version of the prototype tool presented in the previous section that can take as an input any program written in Python that follows conditions ensuring that the underlying model is a pMC. In particular, this allows taking advantage of all the algorithmic power of python libraries in the description of a model, which is not easily done in standard model checking tools such as PRISM [START_REF] Kwiatkowska | PRISM: Probabilistic Model Checking for Performance and Reliability Analysis[END_REF], PARAM [START_REF] Moritz Hahn | PARAM: A model checker for parametric Markov models[END_REF], Prophesy [DJJ + 15] or Storm [START_REF] Dehnert | A storm is coming: A modern probabilistic model checker[END_REF].

In this section, we describe the application of the parametric statistical model checking technique presented in the previous sections to a complete case-study targeting the analysis of the trajectory of a UAV while taking into account parameters representing the precision of some of its sensors as well as some environment conditions.

In the following, we start with a presentation of the context of the case-study. We then move to a description of the process we followed in order to build the model describing the UAV in this context. Finally, we present experimental results obtained using the Python version of the tool presented in the previous section.

Context of the case-study

Unmanned Aerial Vehicles (UAV) are more and more present in our lives through entertainment or industrial activities. They can be dangerous for their environment, for instance in case of a failure when an UAV (aka a drone) is flying above a crowd. Unfortunately until today, there does not exist any kind of UAV regulation around the world. Only some recommendations are used; for instance in order to avoid accidents in case of malfunctioning, a drone should never fly above a crowd. Guaranteeing that a drone never flies above a crowd requires to pay close attention to the drone trajectory computation as well as to the accuracy of the measurements concerning its immediate position in space and its movements. However, a rigorous study is necessary to ensure reliability of the drone control system, for instance by decreasing the risks of failure using the appropriate tuning of the drone flying parameters which impact the computation of its trajectory. Accordingly, the questions are how to prove that the UAV failure probability is low and which parameters have to be taken into account to ensure human safety during performances including UAVs.

High-quality aircrafts such as Hexarotors can easily avoid the majority of minor failures related to hardware because they can fly with only five motors and the probability of concurrent failure of more than two motors is in general insignificant. In the same way, in case of battery failure, the UAV is able to land down on a specified area without any safety issue for the environment as long as it is situated in a safe zone where humans are not endangered. However, software failure may be a lot more problematic and complex to study. In this case, the UAV behavior might become unpredictable. One critical issue in this context is the potential inaccuracy of position estimation in drone systems, either as a result of inaccurate sensor measurements or of misinterpretation of data coming from those sensors. Moreover, besides aircraft system failure consideration, there is also a far more critical aspect to take into account: the weather environment. Therefore, a general approach to improve UAV safety is to study the impact of inaccuracy in position measurements on the resulting flight path compared to a given, fixed, flight plan while taking into account weather conditions. The purpose of this case-study is therefore to provide means to study the reliability of UAVs in the context of a given flight plan. In order to do that, we have to build a formal model which (2) to tune those parameters in order to reduce the system failure probability. To this intent, we thoroughly study the UAV system, formalize it and analyze it using parametric statistical model checking. Among the components of a drone system, we particularly focus on the Flight Controller (FC), which is responsible for computing estimations of the UAV position during its flight in order to adapt its trajectory to a given predefined flight plan. We therefore build a formal model of the flight controller in terms of parametric probabilistic models that takes into account the potential inaccuracy of the position estimation. Since UAVs are particularly sensitive to the weather environment (and in particular to wind conditions), we also enhance our model in order to take into account potential wind perturbations. Since wind force can drastically vary from one point of a given flight plan to another, we also use parameters to encode the wind force and allow our model to adapt to particular weather conditions.

Building a formal model of a UAV

We now present our method to build the UAV model. Recall that we are interested in studying human safety w.r.t. UAV flights. In order to do so, we use inputs such as the intended flight plan of the UAV in order to partition the spatial environment into zones that can either be safe or unsafe for the public. The aim of our model is then to evaluate the probability for a UAV to enter a "forbidden" (i.e., unsafe) zone.

We start by explaining how the zones are computed with respect to the given flight plan. We then show how the UAV software can be decomposed into components and focus on the most important ones. Finally, we detail how the formal models for the important components are built and present the resulting global model. The size of each safety zone is not definitely fixed; it can be defined for a specific requirement or for a given application. In practice the safety zones are specifically defined for a flight environment and for a given flight plan. The main principle is that no human should be present in Zones 1 to 3, while a few people can be present in Zone 4 and most people can be present in Zone 5. As a consequence, the probability that the UAV endangers humans is directly proportional to the probability that it enters Zones 4 or 5. In the following of this section, our target will therefore be to estimate this probability.

Safety zones

Drone components

We now move to the decomposition of the UAV hardware and UAV software into components and introduce the most important component in the UAV system: the flight controller (FC). The FC is responsible for collecting data from various sensors, using this data to compute the precise position and attitude of the drone and adjust the attitude in order to follow the given flight plan to the best of its ability.

Notice the difference between position and attitude: while the position of the UAV is defined by 3-dimensional coordinates x, y and z, its attitude is the collection of yaw, pitch and roll measurements for the UAV compared to the vertical (see Figure 3.9). The attitude allows to control the movement of the UAV: by controlling the speed of each motor, one can control which motor will be the highest, and hence control the direction the UAV will fly to.

Flight controller. As explained above, the FC is the central component in any UAV as it is responsible for collecting data from sensors and translating them to the UAV attitude. An overview of the FC of an UAV is given in Figure 3.10. Remark that the FC can be linked to 3.5. Application to UAV flight plan analysis components responsible for communicating with a remote control. While these components are necessary in order to allow a pilot to take over when the automatic flight mode of the UAV fails, we will consider in the following that this is not the case and that the UAV we study are always in automatic flight mode (i.e., following the pre-defined flight plan).

As one can see from Figure 3.10, the intuitive behavior of the FC is as follows. The filter uses sensors measurements in order to compute the current drone position and attitude. Since the data can be noisy and inaccurate, the filter uses complex algorithms in order to clean the noises in the measurements and compute a realistic position and attitude. Remark that in some cases, the filter can itself introduce inaccuracy in the computed position and attitude, which can be problematic. Once the estimated current position and attitude are computed, the Proportional Integral Derivative (PID) uses this information to compute the local trajectory that the drone has to follow in order to be as close as possible to its intended flight plan. This local trajectory is then transformed into a new value for the attitude of the drone. Finally, Modulation transforms this attitude into signal to the Electronic Speed Controller (ESC) which is responsible for controlling each motor's speed. Recall that we are interested in computing the probability that a UAV enters a forbidden zone while following its flight plan. By construction, as long as the position and attitude measurements are perfect, there is no reason why the UAV should deviate from its intended trajectory, and therefore the probability that it enters a forbidden zone is null. However, as explained above, the data gathered from sensors can be noisy and inaccuracy can sometimes be introduced through filtering. In this case, the estimated position and attitude of the UAV can be faulty, resulting in a deviation from the intended flight plan and potentially leading to a forbidden zone. It is therefore of paramount importance to study how the filters work and to take into account in our formal model the potential inaccuracy of position and attitude measurement. This is illustrated in Figure 3.11, where we assume that the effective position of the drone is in A, but the estimated position that takes into account sensor and filter noise is in A . In this case, the role of the FC is to ensure that the drone goes to B despite the current deviation from its flight plan. As a consequence, the orders sent to the ESC aim at following the trajectory from A to B. However, since the effective position of the drone is in A, this will lead the drone to follow the trajectory from A to B instead, therefore leading to an effective deviation from the flight plan.

For the sake of conciseness, we skip the details of the trajectory computation (they can be found in The resulting probabilistic choice depicted in the box labelled Filter Computation therefore dictates the distance between the exact and estimated position of the UAV. This choice is followed by a computation in the box labelled Safety Zone Computation that computes the exact coordinates of the next position of the drone and allows to decide the safety zone to which this position belongs. When the wind is not taken into account, the result of this computation is enough to decide whether the model should pursue its execution. When the wind is taken into account, another step follows, depicted in the box labelled Wind Computation, where other probabilistic parameters are used in order to decide the wind strength (we assume that the direction is constant) and a new position taking into account these perturbations is computed. Finally, the zone to which this last position belongs is computed and, depending on whether this zone is safe, the model goes on to another position estimation. Remark that the position estimation frequency and the position and distance of checkpoints in the flight plan are given as inputs to the model. The position of checkpoints in the flight plan allows to compute the required UAV speed, while the frequency of the position estimation allows to fix the number of position estimations that will happen in a given flight plan (i.e. the number of loops the model goes through, at most).

Implementation, experimentations and results

The model presented above has been implemented in a Python program, and takes into account deviations on axes x and y, potential change of intermediate checkpoints in case of delay along a complex (more than 2 points) trajectory as well as wind perturbations. For the sake of simplicity, we assumed here that the wind direction was constant with respect to the trajectory (i.e., wind angle is constant in the local Cartesian system associated to the drone), but the generalisation to a varying wind direction is straightforward.

Remark. Preliminary versions of our model (that do not take all details into account) have been implemented in the Prism input format and run on Prism and PARAM. In all cases but the most simplistic versions (discretized position with only 5 possible values), this has resulted in Out of Memory or Timeout errors.

As in the previous section, running our prototype tool on this model returns polynomials representing (1) the estimations of the probability that the drone enters an unsafe region during its flight and (2) the associated confidence interval. Instead of reporting the resulting polynomials, which consist in thousands of terms, we will only present the evaluation of these polynomials using realistic values for the parameters. We defined two scenarios (Scenario 1, Scenario 2) with one set of values of parameters for each scenario. For these two scenarios, ProbaF0 (resp F1, F2, F3, F4) models the probability that the estimated position is from 0 to 2m (resp. 2 -4m, 4 -6m, 6 -8m, 8 -10m) from the real position. In the first (resp. second) scenario, we have set these values to 0.15/0.3/0.4/0.1/0.05 (resp. 0.1/0.25/0.35/0.2/0.1). According to field experiments, the first scenario is more realistic than the second one. Similarly, the wind parameters correspond to the probability of having a wind force of 0 -20km/h,20 -30km/h, 30 -50km/h and 50 -70km/h respectively and have been set to 0.55/0.43/0.01/0.01 (which corresponds to typical weather conditions in Nantes, France) for the numerical evaluation. In both scenarios, Zone 4 (resp. 5) is situated 8m (resp. 50m) from the flight plan.

Remark. Although we only provide the evaluation of the polynomials for two scenarios, the advantage of our technique is that we could evaluate the polynomials for any given parameter valuation. Indeed, once the polynomials are obtained, their evaluation with a given parameter value is very efficient: no further simulations of the model are given. On the contrary, using standard SMC techniques would require performing again the whole set of simulations every time a new parameter valuation is chosen.

In Table 3.13, we gather the results for running the simulation for the two considered scenarios; the simulation with pSMC is performed with 10k, 20k and 50k samples. Each time, the two polynomials are computed and then evaluated using the parameter values given above. In order to illustrate the stability of our results despite their statistical nature, each complete scenario was performed two times (labelled V1 and V2 in the table). The value reported in the table represents 3.6. Potential improvements the probability of the UAV eventually reaching Zones 4 or 5 during its flight. Experiments were performed using three different versions of the model. M 1 is a simple version that does not take wind perturbations into account. M 2 is an intermediate version where wind force is probabilistic but not parametric (the values are set to the scenario given above). Finally, M 3 is the complete version where wind force is parametric. For each model, the number of parameters involved is reported in the column #p. The flight plan consists of only 3 points, resembling the one shown in Figure 3.8, with a total flight duration of 5s and a position estimation frequency of 1Hz. Remark that the probabilities of entering the forbidden zones are quite high. This is not surprising as Zone 4 is situated 8m from the intended trajectory and the precision of position estimation can be up to 10m. These values have been made deliberately high for the purpose of this study but can be chosen more realistically when verifying the real model.

Although the analysis of the obtained polynomial is not an easy task, the user can use those polynomials to get an estimation of the probability that the drone enters the forbidden zones for any chosen parameter values without having to run the model again. More than that, this case-study has confirmed the potential of parametric statistical model checking and shown that our prototype tool, even in its current crude form, can scale up to systems well out of reach of state-of-the-art (parametric) model checkers.

However, there are still many open questions on how to analyze and take advantage of the results offered above (besides the evaluation of the polynomials). Moreover, many choices in our technique have been made by default, but could have been made otherwise. In the next section, we discuss some promising lines of work to enhance the parametric statistical model checking technique itself.

Potential improvements

In this section we explore some variants of our technique in order to tackle some of its inherent problems. As a reminder of the notations and definitions given in Section 3.3, the random variable representing the binary output of each of simulations is written Y and the normalization function f is the one that assigns values to transition probabilities in order to sample the considered pMC M .

Choice of normalization function

The choice of the normalization function f has a huge impact on the convergence speed of our method. Figure 3.15 illustrates this on a simple example where the property tested is the reachability of state win in the three state pMC given Figure 3.14. In Figure 3.15, where the number of simulations is set to 500, we can see that with a normalization function assigning 0.99 to p the confidence intervals are much larger (in particular for small values of p) than with a normalization function assigning 0.01 to p.

This can be explained by the following computation that bounds the variance of Y .

v(Var(Y )) = v E l M f (Y 2 ) -E l M f (Y ) 2 (3.8) = v   ∑ ω∈Γ M f (l) P M f (ω)y(ω) 2   -E l M v (r) 2 (3.9) = ∑ ω∈Γ M f (l) (v(Pa(ω))) 2 r(ω) 2 /P M f (ω) -E l M v (r) 2 (3.10) = ∑ ω∈Γ M f (l) P M v (ω)r(ω) 2 P M v (ω)/P M f (ω) -E l M v (r) 2 (3.11) ≤ m ∑ ω∈Γ M f (l) P M v (ω)r(ω) 2 -E l M v (r) 2 (3.12) = mE l M v (r 2 ) -E l M v (r) 2 (3.13) = (m -1)E l M v (r 2 ) +Var l M v (r) (3.14)
Where m = max ω∈Γ M f (l),r(ω) =0 P M v (ω)/P M f (ω).

In the example, with the normalization function assigning 0.99 to p we have m = (1p)/0.01 = 100 -100p, thus for small values of p the obtained variance is great. On the contrary, considering the normalization function assigning 0.01 to p, we have m = (1p)/0.99, which is always quite small.

Following this remark, a good direction to improve the implemented prototype would be to automatically choose the normalization function minimizing m. Note that such a normalization function can be hard to compute, but one could first use the uniform normalization function then find which valid valuation maximizes the size of the confidence interval and restart an evaluation for this valid valuation. This technique may produce good results since the size of the confidence intervals decreases only with the square root of the number of runs. Thus decreasing m may be faster than increasing the number of simulations.

Modification of the structure

Notice that our technique is valid only for valuations such that M f and M v have the same structure. Indeed, to obtain the approximation of the parametric expected value we divide by the number of runs. However, when evaluating this parametric expected value on valuations modifying the structure of the MC M f , it may be the case that some of the runs are impossible for this valuation. For example, any run that takes a transition with a parametric probability p is impossible to obtain for a valuation v such that v(p) = 0, although it may appear in the

simulations of M f if f (p) = 0.
To address this problem, instead of dividing by the total number of runs obtained when simulating under the normalization function f , we can divide by the number of runs that are possible for the considered valuation v. Formally, instead of the estimator (∑ n i=1 Y i )/n, we can use (∑ n i=1 Y i )/N where N is the function of v defined by:

N(v) = |{i|P M v (ω i ) > 0}|
This can be computed on-the-fly. Defining the estimated variance and the random confidence interval with the same technique would give us a parametric approximation technique which is valid for any valid valuation (regardless of structure preservation). Note however that this is valid only if the structure of M f allows more runs than the structure of

M v i.e., if Γ M v (l) ⊆ Γ M f (l).

Otherwise, this technique only gives an approximation of E

l M v (r(ω)|ω ∈ Γ M f (l)).
Notice also that this is useful only if N(v) is large enough. Indeed if N(v) << n it may be more relevant to estimate Y for a normalization function giving the same structure as M v .

Complement of the property

In classic Monte-Carlo, there is no need to consider the complement of the property. Indeed, when one approaches both the probability of ϕ and that of ¬ϕ by γ and γ ¬ respectively, it is always the case that γ = 1γ ¬ . In our approach however, since the probabilities of the runs are normalized, this does not always hold. Consider for example the pMC given figure 3.14 with a normalization function such that p = 1ε with ε close to 0. With high probability, our approach would give that the probability of reaching win is γ = 0 since it would miss the run reaching win which has a really low probability for this normalization function (note that this problem only appears if the number of simulations is not big enough). However, if we estimate the negation as well, we have γ ¬ = (n * p/(1ε))/n = p/(1ε), which is much closer to the truth for any valuation with p ≈ 1.

It is thus relevant to consider both the property and its negation when considering a parametric approximation of a probability. Note that this holds only for the approximation of a probability and not for the general expected value of a reward. Note also that this is a good example of the importance of a good normalization function.

To implement this in practice, one could approach both the probabilities of ϕ and ¬ϕ and whenever v( γ ¬ ) ≈ 1v(γ) one could either increase the number of simulations or restart the approximation for the normalization function v.

Conclusion

In this chapter, we have presented a new technique for statistical model checking of parametric Markov chains. This technique is based on a parametric adaptation of the standard Monte-Carlo analysis. We have shown that our technique allows to estimate the expected value associated to any reward function by a polynomial function of the parameters, and proposed as well a parametric confidence interval for the estimation. Compared to exact model checking techniques, the technique we have proposed here offers the same benefits as standard simulation techniques for non-parametric models: better scalability and a complexity which is largely independent of the model complexity (be it in the size of the state space, in the type of used features, or in the number of parameters). Contrary to existing statistical model checking techniques for non-deterministic systems (such as the one presented in [HMZ + 12]), the one we have proposed in this chapter allows to compute parametric confidence intervals. Finally, our technique has been implemented in a prototype tool using the Python language. This prototype tool accepts inputs both in the PRISM language and in the form of Python programs. It has been tested on a set of benchmarks and on an industrial case-study with encouraging results.

Perspectives

As we have shown in this chapter, parametric statistical model checking is a technique that offers a lot of promises, but is, in its current state, far from being used at its full potential. Besides the potential improvements presented in Section 3.6, a lot of work remains in order to compare this technique (and our prototype tool) to state-of-the-art tools such as PARAM [START_REF] Moritz Hahn | PARAM: A model checker for parametric Markov models[END_REF], Prophesy [DJJ + 15] or Storm [START_REF] Dehnert | A storm is coming: A modern probabilistic model checker[END_REF]. In particular, as soon as our tool is in a more optimized state, we plan on running a set of benchmark to compare its scalability in terms of model size with the one of Storm [START_REF] Dehnert | A storm is coming: A modern probabilistic model checker[END_REF]. We conjecture that such an optimized tool as Storm will be orders of magnitude more efficient than our tool on small models. However, we expect that the resources needed by Storm will grow linearly with the model size, while those needed by our tool should remain fairly constant. On the other hand, the resources needed by our tool should grow linearly with the size of the property that we verify, whereas those needed by Storm should not. Still, when considering the results presented in Section 3.5, we are confident that regardless of the property, our tool should be able to handle models that Storm cannot.

The parametric statistical model checking technique we have presented can also be adapted to verify models with non-determinism, such as Markov Decision Processes. Indeed, the nondeterministic choices in this context could be replaced by parametric choices and our technique could then be used in order to study optimal schedulers. Moreover, we could also envision using the version of our tool that takes as input a Python program to study the optimization of controllers: by considering a non-deterministic/probabilistic/parametric program, one could build a parametric controller that could then be analyzed and optimized in an automated manner by using our technique (i.e., executing the controller in parallel with the system in order to be able to analyze its effect on the system). This could be done for instance on programs generated using the probabilistic Event-B formalism presented in Chapter 2.

Finally, the most exciting perspective, to our point of view, is the ability of this technique and the associated tools to analyze models and systems that are described in any programming language. In the past years, we have indeed been in contact with researchers from other domains, such as Oceanography and we have learnt from those interactions that, although the verification techniques and the modeling formalisms we develop in our domain are very powerful and could be very useful to them, the effort that needs to be done in order to master both the modeling formalisms and the techniques is too important for potential users out of our domain. However, because of its simplicity and because it can be used on existing models written in any language without much effort, we have been able to convince some of them that statistical model checking can be integrated in the parameterization process of some of their models. Notably, we have developed a technique, called the Statistical Model Checking Engine (SMCE), that can be used almost automatically on any probabilistic and parametric model in order to tune the model to make it fit with experimental data. Compared to other existing parameterization techniques, the SMCE not only yields the set of optimal parameter values, but also provides the user with a detailed correlation analysis of the whole set of parameters -which is outside the scope of any other parameterization technique used in this context, to the best of our knowledge. This work has been published in a famous broad-audience scientific journal [REG + 20], which is, to our opinion, a major achievement. Our aim now is therefore to pursue the development of both sides of our tool, in order to propose state-of-the-art automated parametric analyses that could be both used by researchers from our domain but also by any modeler able to write an executable program representing their system. Finally, we also aim in the long term at linking the analysis we do here with automated abstraction techniques based on (deep-)learning in order to reduce the cost of sampling very complex models. Succeeding in this challenge would then allow us to use our techniques to verify large-scale models such as those developed for the analysis of climate change.

Introduction

In most of nowadays real life jobs or applications, employees, workers, consumers and users are using connected devices. Hence, they are prone to not only external exploits but also internal misuse that can become dangerous. Consequently, it is important to monitor the behavior of users and workers in their environment to ensure a "secure" work flow. For instance, monitoring the behavior of truck drivers can allow the supervisor to verify if they are getting enough rest time while on route. Computer scientists have been putting a lot of effort in researches in order to build a safe perimeter in which connected devices can be used [WHA + 16, A + 15]. Contrarily to previous chapters, we consider here that a precise model of the systems under study is not available. More than that, we also consider that the systems can only be analyzed through a given set of observations that cannot be extended. The uncertainty here therefore comes from the lack of information we have on the behavior of the system/user.

In this chapter, we use behavior analytics in order to classify behaviors, from "standard" to "dangerous", by rating their level of dangerousness. The key to establishing an adequate behavior analytic is a set of faithful data recorded from the concerned system. In 2020 we generate in ten minutes more data than the entire data that is recorded throughout history until 2003. Thus, clearly we have no real problems with generating data nor storing it, hence the problem is our capability of extracting useful information from this data. This is where data mining techniques [START_REF] Van Der Aalst | Process Mining: Discovery, Conformance and Enhancement of Business Processes[END_REF] come into play.

In order to build a secure access to data in a real world system and to ensure its safeness from any upcoming potential threat, one should learn the dependencies and behaviors of the different components of the system. One must then identify malicious behaviors and act at the right moment to intercept them. This is where model-based Machine Learning, that allows to build a model of the concerned system/user based on observed data, comes into play. Many types of modeling formalisms exist in the literature, each developed for its own purpose. The ones that were considered in the previous chapters, such as Event-B or (parametric) (Interval) Markov Chains, are better tailored to the verification of systems whose models are available. Other families of modeling formalisms, such as Hidden Markov models [RVVDA08, KA98, VM98] or (Dynamic) Bayesian Networks [START_REF] Murphy | Dynamic bayesian networks: representation, inference and learning[END_REF], are better tailored to the learning of uncertain discrete probabilistic models. In this work, we insist on the importance of timing information in the dependencies between variables of the system, and therefore focus on continuous-time graphical probabilistic models. Several such famillies of models exist, e.g., continuous time Bayesian networks [START_REF] Nodelman | Continuous time Bayesian networks[END_REF], Markov jump processes [START_REF] Rao | Fast MCMC sampling for Markov jump processes and extensions[END_REF], Poisson networks [START_REF] Rajaram | Poisson-networks: A model for structured point processes[END_REF] and Graphical Event Models (GEMs) [START_REF] Gunawardana | A model for temporal dependencies in event streams[END_REF]. In this work, we chose to use Recursive Timescale Graphical Models (RTGEMs) [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] a sub-family of GEMs, that present advantages compared to the other formalisms. In particular, they are designed to universally approximate any smooth, nonexplosive, stationary, multivariate temporal marked point process [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF].

In this chapter, we therefore propose a strategy for analyzing the level of dangerousness of a given system only observed through pre-generated data by using RTGEM models. The strategy unfolds as follows:

Background

(1) We start by learning an RTGEM model that is representative of the given data;

(2) We analyze this model and decide whether it satisfies given security properties;

(3) When it does not, we explore the neighborhood of this model to find an RTGEM that does satisfy the property and, if one is found, measure the distance between this safe model and the original one. This distance gives us the level of dangerousness of the original data.

The chapter is dedicated to developping and experimenting the necessary techniques to achieve the above strategy. In Section 4.2, we start with properly introducing the formalism and methods used later on. In particular, we recall statistical model checking (already introduced in Chapter 3), which will be used in the verification process of RTGEM, and provide a progressive definition of the RTGEM formalism. Section 4.3 introduces the main techniques required for our strategy, such as the learning and sampling algorithms for RTGEM, and defines the notion of distance used for determining the level of dangerousness of a given RTGEM. Finally, Section 4.4 details the experiments performed to evaluate our strategy and Sections 4.5 and 4.6 conclude the chapter.

Background

We start with some background definitions that will then be used in the rest of the chapter. As seen in the previous chapters, model checking suffers from the state explosion problem when the models are large. In the context of this Chapter, our aim is to build models that represent complex statistical dependencies between several variables. As a consequence, the models we will build might be large. Moreover, standard models used for representing statistical dependencies (e.g., Bayesian models [START_REF] Maxwell | A bayesian approach to learning bayesian networks with local structure[END_REF][START_REF] Murphy | Dynamic bayesian networks: representation, inference and learning[END_REF] or probabilistic graphical models [TEF + 05, GM16]) are not adapted to standard model checking techniques. As a consequence, we have chosen to build our strategy on statistical model checking [START_REF] Legay | Statistical model checking: An overview[END_REF], which we now recall.

Statistical model checking

As intoduced in Chapter 3, SMC is a simulation-based technique that allows estimating the probability that a model satisfies a given linear property while giving formal guarantees on the precision and the error-rate of this estimation. Moreover, one of the advantages of SMC is that it can be applied to any type of stochastic model that can be executed. SMC comes in two flavors: (1) quantitative SMC that allows to verify quantitative properties (i.e., compute the probability with which the model satisfies a property); and (2) qualitative SMC that allows to verify qualitative property (i.e., check whether the probability with which a model satisfies a property is greater than a given threshold). Quantitative SMC, and in particular Monte-Carlo estimation, has been presented in details in Section 3.3.1. We now briefly present the qualitative version of SMC.

In the following, we consider a stochastic system S and a linear property ϕ. As previously mentioned, SMC is based on simulations, and each simulation is represented by a trace that has a binary outcome for satisfying the property. Therefore, we consider a discrete random variable B i following a Bernoulli distribution of parameter p, that is associated to each trace and can take two values: 0 (if the property is not satisfied) or 1 (if the property is satisfied). Assume that Pr[B i = 1] = p is the probability of satisfying the property and Pr[B i = 0] = 1p is the probability of rejecting the property. Qualitative SMC allows to decide if p is greater (or smaller) than a given threshold θ.

Qualitative approach. The original qualitative approaches are proposed in [START_REF] Hakan | Verification and planning for stochastic processes with asynchronous events[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF],

and are based on hypothesis testing. The idea is to test two hypothesis, H 0 : p ≤ θ against H 1 : p > θ, while bounding the probability of making an error. A Type-I error is when we accept H 0 while H 1 holds and a Type-II error is when we accept H 1 while H 0 holds. Therefore, we define two parameters α and β, with (α, β) the "strength" of the test, being the pair of bounding errors. An ideal performance of the test is when the Type-I error is equal to α and the Type-II error is equal to β, and both α and β are small. However, in practice it is impossible to ensure a low α and β simultaneously. In order to avoid this problem, an indifference region [p 0 , p 1 ] is defined.

Let p = P S (ϕ) the probability of ϕ being satisfied on the system S. In order to determine whether

p > θ (qualitative property, θ being in [p 0 , p 1 ]), one can test H 0 : p ≤ p 0 (instead of H 0 : p ≤ θ)
against H 1 : p > p 1 (instead of H 1 : p > θ. However, the tuning of the indifference region in practice is constrained by the loss of precision (when the size of the region increases) and by the fact that we cannot conclude about the result if p is inside [p 0 , p 1 ]. We present two hypothesis testing algorithms: Single Sampling Plan (SSP) and Sequential Probability Ratio Test (SPRT).

The Single Sampling Plan consists in specifying a constant c and a number of simulations n to test whether ∑ n i=1 b i > c, where b i is the outcome of the Bernoulli random variable, in order to see which hypothesis is accepted. Thus, the hardest part in this algorithm is to compute values for the pair (n, c) with respect to the pair (α, β) and the indifference region. The number n increases when we minimize the size of the indifference region and the parameters α and β. An optimization algorithm was proposed in [START_REF] Hakan | Verification and planning for stochastic processes with asynchronous events[END_REF] in order to determine a pair (n, c) where n is minimal.

Background

The Sequential Probability Ratio Test is an approach proposed by Wald [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF], consisting in choosing two values A and B (A > B) ensuring the strength of the test and that are computed using α, β and the indifference region. Unlike SSP, this approach does not fix a number of simulations a priori but instead dynamically checks whether the number of simulations done so far allows to conclude. As soon as a conclusion can be obtained, the algorithm stops. If m is the number of already performed observations, the test is based on the following metric :

p 1m p 0m = m ∏ i=1 P(B i = b i | p = p 1 ) P(B i = b i | p = p 0 ) = p d m 1 (1 -p 1 ) m-d m p d m 0 (1 -p 0 ) m-d m
,

where d m = ∑ m i=1 b i . It is shown in [You05]
that we can accept H 1 after m samples if p 1m p 0m ≥ A, whereas we can accept H 0 if p 1m p 0m ≤ B. Using this approach we can reduce the number of simulations in certain scenarios compared to the Single Sampling Plan approach (where the number of samples is fixed to n), and avoid doing nm useless samples. However, it is also demonstrated that convergence can be very slow when p is too close to θ. In this case, SSP might be more advantageous.

We now progressively introduce the types of models that will be used in the rest of this chapter: Recursive Timescale Graphical Event Models (RTGEM).

Graphical event models

We start with a formal definition of the structure we use to represent our input data.

Marked point processes

Let S be a totally ordered set. A point process is a set of points that are positioned (following the order) in S. A marked point process is a point process that contains additional features at each point. In other words, a marked point process (m.p.p.) is composed of a point process and marks associated witch each point. If we consider a set of marks M, an m.p.p. can be expressed as the pair:

{(s i , m i ) : i = 1, ..., n} where s i ∈ S, m i ∈ M and s i ≤ s i+1 for all 1 ≤ i ≤ n.

Henceforth, we assume that the space S that is used is the one-dimensional space of time, and the marks are labels describing each point. Marked point processes are used for expressing event streams in this framework, i.e. random labeled (marked) events (points) that arrive at different times {t i } (unidimensional space). The word point is used for describing events as being instantaneous in the time dimension and for the sake of simplicity, the notation l for labels is going to be used instead of m for marks.

In practice, we tend to use a timeline to represent data that arrive at irregular intervals. A timeline can be defined as a sequence of pairs {time, event} capturing the relative frequency and ordering of events [START_REF] Jeremy | Forest-based point process for event prediction from electronic health records[END_REF]. In figure 4.1, an example is shown of how a timeline can be seen as a multivariate marked point process. In the following, we define event streams as well as stochastic models that are capable of representing such an evolutionary process.

If we take a step back, a point process can be described by unrolling a stochastic model that defines inter arrival times between different events, i.e. defining the time of the next event based on all the times of previous events. In the following, we will explore a family of stochastic models based on conditional intensity functions, called Conditional Intensity Models (CIM). We start by defining the data type that is generally used for learning these models.

Event streams. An event stream consists in a timed sequence of events with strictly increasing timestamps. An event stream can be written as: (t 1 , l 1 ), . . . , (t n , l n ), with 0 < t i < t i+1 < t * (t * = t n+1 ) for all 1 ≤ i ≤ n -1 and where l i are labels chosen from a finite set of labels L. Hence, a single event is defined by the pair (t i , l i ) but for the sake of simplicity we sometimes refer to an event only by its label l when the context is clear enough. We note that in the following, t 0 = 0 and t * = t n+1 are used as conventions, as well as the fact that two events cannot occur at the 126 4.2. Background exact same time. An event stream can be considered as an m.p.p. and the data that can be derived from it can be written as x t * . We write |x t * | for the size of our data x t * (the number of events in the sequence, here it is n). The history at time t is the set of all the events that occurred before t: h i denotes the ith history h i = (t 1 , l 1 ), ..., (t i-1 , l i-1 ). With this in mind, we define a conditional intensity function. In the following, and for the sake of simplicity, data will sometimes be written as x, history as h and time as t, without indexing them explicitly at every instance, but providing enough context to understand their use.

Definition 14. A conditional intensity function, defines the risk of observing a given event at a certain time, depending on the observed history in some data x. We recall that E[X] is the expected value of a random variable X. Mathematically a conditional intensity function λ is written as:

λ(t) = lim ∆ t →0 E[N(t,t + ∆ t )|h t ] ∆ t
where N(T ) denotes the number of points (events) occurring in a time interval T in x.

In other words, the conditional intensity function allows to specify the mean inter arrival times between two events in a process. Furthermore, we recall an important assumption for the following about this approach: there cannot be more than one event in an infinitesimal interval of time.

Before moving forward to the definition of a CIM, we start with an introductory example about a common probability distribution based on an intensity function. Consider for instance the Poisson distribution [START_REF] Frank | Handbook of the Poisson Distribution[END_REF], that is a probability distribution representing a number of events occurring in a fixed interval of time under the assumption that any event is independent with regards to the duration from the previous events. A Poisson distribution is used for modeling processes such as the number of cars arriving to a garage between two given times. The probability function of a random variable X that follows a Poisson distribution of parameter λ (where also the expectation and the variance are equal to λ) is given by: P(X = k) = λ k e -λ k! . We note that in order to sample from such distribution, the inter arrival times between every two events should be computed. The inter arrival time follows the exponential distribution in the case of a Poisson distribution, with a mean of 1/λ, i.e. the mean time between every two events in a Poisson distribution is 1/λ.

A CIM that models a marked point process must treat each event individually, thus each label l ∈ L must be associated with its own conditional intensity function λ l (t | h). The conditional intensity function, describes the risk of having the event l at time t given a certain history h of events (the totality of the events or a part of them depending on the CIM). For instance, in a Markovian CIM where an event can be particularly dependent from a set of events (called parents), the conditional intensity functions satisfy the following property:

λ l (t | h) = λ l (t | [h] Pa(l) )
where Pa(l) is the set of parents of l and [h] Pa(l) is the history of the parents only.

Definition 15. [DVJ07] Formally, a CIM θ is a set of indexed conditional intensity functions {λ l (t | x)} l∈L . The data likelihood function is the joint density function of all the points in the m.p.p. that can be factorized into a product of all conditional intensity functions and can be written as:

p(x | θ) = ∏ l∈L n ∏ i=1 λ l (t i | h i ; θ) 1 l (l i ) e -Λ l (t i |h i ;θ)
where Λ l (t | h; θ) = t -∞ λ l (τ | x; θ)dτ for the data x and the indicator function 1 l (l ) is one if l = l and zero otherwise.

We now introduce Piecewise-Constant Conditional Intensity Models (PCIMs).

Piecewise-constant conditional intensity models

Piecewise-Constant Conditional Intensity Models (PCIMs) are based on the assumption that intensity functions are constant in given time intervals. In other words, there is a mapping associated with these models that assigns a parameter λ for every label in L according to an active state. The active state is determined in function of the time t and the data x.

These models are important because they form the basis of graphical event models and they ensure fast and efficient learning and inference. These models are related to Poisson networks [TEF + 05], since Poisson networks also contain piecewise constant parameters in their mapping (see more detail in [START_REF] Gunawardana | A model for temporal dependencies in event streams[END_REF]). However, it was experimentally shown in [START_REF] Gunawardana | A model for temporal dependencies in event streams[END_REF] that they are, by orders of magnitude, faster to train than Poisson networks.

Let L be a set of labels. Each label l ∈ L is associated to a set of discrete states Σ l . For each label l and each state s ∈ Σ l we have a parameter λ ls . An active state s for each label (correspondingly 4.2. Background the active parameter λ ls ) has to be determined by a mapping σ l : T × X -→ Σ l (with T the set of all possible times and X the set of all possible data). Definition 16. A PCIM is defined by local structures S l = (Σ l , σ l (t, x)) and local piecewise constant parameters λ ls . Let S = {S l } l∈L be the set of all discrete states and θ = {λ ls } l∈L,s∈Σ l be the set of associated parameters. The PCIM data likelihood knowing the structure and the parameters in this case can be written as:

p(x | S, θ) = ∏ l∈L ∏ s∈Σ l λ M ls (x) ls e -λ ls T ls (x)
where M ls (x) is the number of occurrence of events of type l while s is active in the event sequence x, and T ls (x) is the total duration while s was active for event type l.

In practice, the learning of a PCIM is divided into two processes: the local structure learning process that could be done using a decision tree to alternate between active states and define adequate mapping from the data to the states set similarly to [START_REF] Maxwell | A bayesian approach to learning bayesian networks with local structure[END_REF]; and the parameters learning process that uses statistical estimates depending on the type of the model and knowing the structure. Commonly, a greedy search is used where for a fixed structure (for instance a given decision tree and a given mapping for states) the parameters are calculated, and the structure is continuously modified until we cannot find a better model based on a selection criterion or a certain "gain".

The theory of CIMs is a root to many continuous time formalisms. It is an expressive tool that can help build graphical models such as Poisson networks and GEMs, that will be detailed later on. Moreover, PCIMs favor the discretization of states instead of the discretization of time, creating a set of constant parameters associated with each "scenario" in time and history of past events in which an event can occur. In the following, we define a family of graphical models that is directly used in the rest of the chapter.

Graphical event models

Definition 17. A Graphical Event Model (GEM) is a tupple G = ((L, E), θ) where (L, E) is a directed graph and θ a set of parameters. A GEM can represent event streams of the type x as defined previously, as well as the dependencies between the different labels (or events) in time.

In this case, the likelihood of the data knowing the graph and its parameters is written as:

p(x t * | t * ) = |x t * | ∏ i=1 λ l i (t i | h i ) |x t * |+1 ∏ i=1 e -∑ l∈L t i t i -1 λ l (τ|h i )dτ , (4.1)
with h |x t * |+1 the entire history of the event stream, including x t * .

In figure 4.2, an example of a GEM is shown with four labels (nodes) A, B, C and D, each with its corresponding parameter (conditional intensity function). The structure of a GEM (directed graph) has no constraints over cycles or loops, giving it a superior degree of expressive power compared to other graphical models based on DAGs, like Dynamic Bayesian Networks for instance, that cannot directly allow cycles or loops in a single time-slice. In this formalism, one should note that the conditional intensity functions λ l (t | h) are not piecewise-constant, which means that they do not take constant values for a certain period of time so the models are a family of CIMs and not PCIMs.

A λ A (t | h) B λ B (t | h) C λ C (t | h) D λ D (t | h)
In this work we are only interested in Markov m.p.p.s with respect to a certain GEM. The conditional intensity functions λ l (t | h) in this case satisfy the following property for all t and h:

λ l (t | h) = λ l (t | [h] Pa(l) )
where Pa(l) are the parents of l in G. This means that the conditional intensity of a certain label l at time t only depends on the history of the parents of l and not the entire history of the process.

The dependencies between the events can be easily depicted on the graph of the GEM. In A fully connected GEM can represent any m.p.p. with labels in L [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF].

We now introduce Timescale GEMs, a subfamily of GEMs that is more granular and can be more easily used in practice.

Timescale graphical event models

A Timescale GEM (TGEM) is a GEM where each dependency between two events (edge in the graph) is defined for a given finite timescale which specifies the temporal horizon and the granularity of the dependency represented by that edge. Formally, a timescale is a set T of half-open intervals (a, b] (with a ≥ 0 and b > a) that form a partition of some interval (c,t h ] (with c ≥ 0), where t h is the highest value of T and is called the horizon of the corresponding edge.

Definition 18. A TGEM M = (G, T ) consists of a GEM G = ((L, E), θ) and a set of timescales T = T e (e∈E) corresponding to the edges E of the graph of G.

We write c l (h,t) as the parent count vector of bounded counts (of occurrences) over the intervals in the timescales of the parents of l (in order to have bounded counts, a maximum threshold should be fixed). We provide an example to explicitly show parent count vectors and parameters (see Example 14).

The conditional intensity of a node (or a variable) labeled l in the GEM only depends on the history of the number of occurrences of its parents within the corresponding timescale. As a consequence, a given parent count vector yields a unique conditional intensity function for the corresponding parameter. The global conditional intensity functions are therefore piecewiseconstant:

λ l (t | h) = {λ l,c l (h,t) } c l (h,t)∈C l .
We use C l to denote the set of all possible parent count vectors of label l. For the following example (Figure 4.3), we consider that all TGEMs are bounded by 1 (making the parent count vectors binary), thus only the fact that a parent has occurred (or not) within the corresponding timescale is important and not the number of times a parent occurs. . The same applies to the rest of the variables, except the variable D that is independent of all other variables, so its conditional intensity function is written as λ D . All conditional intensity functions are equal to constants making them piecewise-constant depending on the corresponding combination of parents.

In the case of TGEMs, the likelihood of the data in equation 4.1 is simplified and written as follows:

p(x t * | t * ) = ∏ l∈L ∏ j∈C l λ n t * ,l, j (x t * ) l, j
e -λ l, j d t * l, j (x t * ) , (4.2)

where the sufficient statistics n t * ,l, j (x t * ) and d t * l, j (x t * ) are the count of l-events and the durations, respectively, when the parent count vector was equal to j (a certain combination of parents).

In [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF], the sufficient statistics are formally written as:

n t * ,l, j (x t * ) = |x t * | ∑ i=1 1(l i = l)1(c l (h i ,t i ) = j) d t * ,l, j (x t * ) = |x t * |+1 ∑ i=1 t i t i-1 1(c l (h i , τ) = j)dτ
where t |x t * |+1 represents the duration d between the last event occurence and the final time t * in 4.2. Background the data.

Recursive timescale graphical event models

In practice, the learning process of GEMs is not an easy task, especially because they are not piecewise-constant intensity models, thus the learning of the conditional intensity functions can become hard (since they are not constants). The difficulty of learning GEMs, along with the necessity of have a "universal" Graphical Event Model that is easy to learn and to handle, led to the creation of a subclass of TGEMs called Recursive Timescale Graphical Event Models (RTGEMs).

Definition 19. The family of RTGEMs is defined recursively to be the finite closure of the empty model M 0 = ((L, {}), {}) under a set of allowed operators O F = {add, split, extend} (discussed in the following).

RTGEMs as described in [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] can universally approximate any smooth multivariate temporal point process. RTGEMs are learned recursively using the set of operators introduced above. More details on the learning procedure are given in the following.

Furthermore, it is proven in the works of [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] that using a greedy search algorithm to learn this class of models always shows structural and parametric consistency. In other words, the learned model converges in probability to the optimal model when the learning data size is increased, unlike the more general cases of TGEMs and GEMs for which, to the best of our knowledge, consistency has not been proven.

A finite consistent RTGEM learning procedure [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] is to do a forward greedy search (for model construction) followed by a backward greedy search (for model refinement), both based on model selection, and using data that is faithful. Bayesian Information Criterion (BIC) [S + 78] is a very general model criterion that has been adapted to select the "best" (or a "better") RTGEM when doing a greedy search, and is written as follows for a model M:

S t * (M) = log(p(x t * | t * ; M, λ t * ,l, j (x t * ))) -∑ l∈L | C l | •log(t * ),
where the left term of the difference is the likelihood of the data knowing the history, the model and the calculated maximum likelihood estimates of the λ functions (written λ). The right term represents the complexity of the model, with Dim(M) = ∑ l∈L | C l | as the dimension of an RTGEM, multiplied by the size of the sample used for learning to define the overall complexity.

Given a model M, the maximum likelihood estimate λ of the λ parameters for M (knowing its structure) is as follows [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] :

λ t * ,l, j (x t * ) = n t * ,l, j (x t * ) d t * l, j (x t * )
Now that RTGEMs are introduced, we explain how such models can be learned, sampled and compared.

Learning, sampling and comparing RTGEMs

The purpose of GEMs is to represent statistical dependencies between variables that one can observe on a data stream. However, the process of building a GEM that faithfully represents the statistical dependencies present in a given data stream can be arduous. RTGEMs have been introduced in [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] in order to solve this problem. By restricting the forms that the dependencies can take to timescales and the set of operations that can be done to build these timescales, the authors have shown that the learning time is drastically reduced. In the following, we present the set of operators that can be used in order to build an RTGEM and briefly explain how the learning process is performed. Sampling the model is also very important in our context in order to be able to perform SMC. We therefore explain next how RTGEMs can be sampled.

Finally, we introduce one of our contributions: the definition of a distance measure that allows to compare two given RTGEMs.

Learning RTGEM

As explained in Sextion 4.2.2, the learning procedure defined in [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] consists of a forward greedy search for model construction followed by a backward greedy search for model refinement. Both the forward and the backward greedy search rely on a model selection criterion. As explained in Section 4.2.2, we usually use the Bayesian Information Criterion adapted to RTGEMs for this purpose.

The inital model consists of a graph with no edges whose vertices are the variables present in the given data stream. During the forward greedy search, the edges of the model (representing the dependencies between the variables) are progressively built by using a limited number of forward elementary operators. The set of forward elementary operators is the following:

O F = {add, split, extend}. The "add" operator adds a non-existing edge to a model and its corresponding timescale T = (0, c], with c a constant (also called horizon). The "split" operator splits one interval (a, b] in the timescale of a chosen edge into two intervals (a, a+b 2 ], ( a+b 2 , b]. The extend operator extends the horizon of a chosen edge by adding the interval (t h , 2t h ], with t h being the previous horizon.

After a number of forward operations (either a fixed number or when the improvement of the model selection criterion is sufficiently low), the forward search ends and a backward search starts. In this backward search, symmetric backward operators are used. For the sake of convenience, and since there are no further information in the literature about these operators, one can write O -1 F = {reverse_add, reverse_split, reverse_extend} for these symmetric backward operators.

The "reverse_add" operator removes a chosen edge with only one interval in its timescale (to make it the exact inverse of the "add" operator). The "reverse_split" operator is used for merging two consecutive intervals in a timescale on a chosen edge that have been initially split. The "reverse_extend" operator removes the highest (the last) interval in a timescale on a chosen edge only if the upper bound of this interval was initially created by an extend. Again, once a sufficient number of operators have been applied, or when the improvement of the model selection criterion is sufficiently low, the backward search ends and a definitive model is produced.

Sampling RTGEM

Once the model is learned, our aim is to verify if it satisfies given security properties. Since we use SMC for this verification, we need to produce samples of the learned RTGEMs. The sampling of such models is similar to the sampling of Poisson Networks [START_REF] Rajaram | Poisson-networks: A model for structured point processes[END_REF]. The needed "elementary" sampling techniques are demonstrated using the example in Figure 4.4 where different simple RTGEMs are shown. In Figure 4.4a the sampling is straightforward since there is only a single constant rate λ A . In Figure 4.4b, the sampling of A is straightforward but the sampling of B will depend on the sampling of A. As a consequence, A must be sampled before B, then the active λ B should be chosen respectively in regards to the occurrences of A along the timeline, and finally B can be sampled using the corresponding λ B value in each section of the timeline. Finally, in Figure 4.4c, things become more complicated: none of the variables can be sampled straightforwardly because of the cycle between A and B. Therefore, some sort of "competition" is

A λ A (a) A λ A B λ B,0 λ B,1 (0, 5] (b) A λ A,0 λ A,1 B λ B,0 λ B,1 (0, 5] (0, 4] (c)
Figure 4.4 -Illustration of elementary sampling techniques for TGEMs created between the variables and two time values τ A and τ B are sampled using the corresponding rates λ A,c l and λ B,c l (c l is the corresponding parent configuration, c l = 0 initially). The highest value (between τ A and τ B ) is rejected because the one that comes before is supposed to change the corresponding conditional intensity function of the other (making the latter one wrongly sampled). Note that the corresponding rates are the conditional intensity functions that are active (with regards to the parents configuration) at the moment of the sampling. Afterwards, the lowest value (between τ A and τ B ) is accepted if it is within the firing period. The firing period is the lowest interval on the timescale of the edge entering the corresponding node (for A it is between 0 and 5), because usually after that interval there is a switch of parents configuration and subsequently a switch in the conditional intensity function.

Using these elementary techniques, sampling can be generalized to a larger number of nodes, but two operations (illustrated in Figure 4.5) need to be performed beforehand. there exists a path between every single pair of nodes. The first operation is illustrated in Figure 4.5b where strongly connected components (nodes of the same color) are brought together forming a new "node" locally and the directed graph is transformed to a DAG at the end of the operation. The second operation (illustrated in Figure 4.5d) is topological ordering [START_REF] Thomas | Introduction to algorithms second edition[END_REF], and it defines an order that can be used for sampling the new nodes of the DAG using the three elementary operations. Note that the topological order is not unique, for instance 1 and 1 could be sampled in random order. A detailed description of the sampling procedure introduced above is given in [START_REF] Antakly | [END_REF].

We now introduce one of our theoretical contributions: a notion of distance that allows to compare two given RTGEMs.

Distance between RTGEMs

To the best of our knowledge, there is no existing metric between RTGEMs. The most intuitive distance measure that comes to mind is the minimal number of operations needed to move from one RTGEM to another since they are built using a recursive procedure. However, such a distance is not accurate in our context because not all operators add (or remove) the same information every time.

In the literature, the popular Hamming distance [START_REF] Hamming | Error detecting and error correcting codes[END_REF] has been adapted for some probabilistic graphical models such as Bayesian networks [START_REF] Tsamardinos | The max-min hill-climbing bayesian network structure learning algorithm[END_REF]. In the following, we propose an extension of the Structural Hamming Distance (SHD), adapted to RTGEMs, where we evaluate the amount of differing information on two different edges.

A timescale in an RTGEM can be represented by a vector v = [0, a, b, c, . . .] containing the successive endpoints values. Given two RTGEMs with the same set of labels L, G 1 = ((L, E 1 ), θ 1 ) and )∈V e id \(0,0)

G 2 = ((L, E 2 ), θ
|v e 1 i -v e 2 j | min(v e 1 i , v e 2 j )   + v e nid v e nid + v e id .
Remark that for each pair of matched endpoints, the relative difference (scaled by its minimum) is taken into account in order to penalize the fact that they may not be perfectly identical.

Obviously, this factor is smaller than 1 (otherwise either v e 1 i or v e2 j would be closer to 0 than to its counterpart, which would imply that (v e 1 i , v e 2 j ) / ∈ V e id ). As a consequence, pairs of imperfectly matching endpoints are less penalized than endpoints that have no match, but still more penalized than pairs of identical endpoints.

Finally, if there are no detected imperfect matches (i.e., V e id = {(v e 1,i , v e 2, j ) : v e 1,i = v e 2, j }), then we have

d * (v e 1 , v e 2 ) = d(v e 1 , v e 2 ).
The proximal distance measure, that will be used in the rest of the chapter, takes into account this similarity between endpoints and is defined similarly to the SHD but using the elementary distance d * defined above instead of d:

SHD * (G 1 , G 2 ) = |(E 1 \ E 2 ) ∪ (E 2 \ E 1 )| + ∑ e∈E 1 ∩E 2 d * (v e 1 , v e 2 ). (4.5)
Although this is not a problem in our context, we remark that the proximal distance defined above is not, strictly speaking, a distance metric but only a semimetric (it does not satisfy the triangular inequality). Counter examples are given in [START_REF] Antakly | [END_REF]. Because it suits our purpose, we will nevertheless call it a distance in the following, but this is an abuse of language.

Experimentations

As explained before, the aim of this work to measure how secure a system, which we can only observe through pre-generated outputs, is w.r.t. a given qualitative security property of the form

P(ϕ | M) > c
where ϕ is a linear property and c is a constant 2 To this aim, we have developped a strategy that can be summed up by Algorithm 1.

The only inputs we have in this context are the query ϕ and some data D. The algorithm is in three Otherwise, the third step (lines 5-6) consists in searching the neighbourhood of M o for a model M * that satisfies the security property and, if one is found, compute the distance between M o and M * . If none is found then the empty model ⊥ is returned and the distance is considered infinite.

Remark that the algorithm presented above is generic, i.e., it does not depend on the choice of the modeling formalism, verification technique or distance notion. In the following, we present experiments performed for each step of this algorithm using the techniques and notions presented in the previous sections, but other techniques and notions could be used instead.

To the best of our knowledge, there is no existing tool in the literature implementing the learning, sampling or verification of RTGEM. We have therefore developed a branch of the PILGRIM library3 dedicated to the manipulation of RTGEM, which we use in the rest of the chapter. All the experiments presented hereafter have been performed on a Windows 10 Enterprise i5 laptop, with a 2.3 GHz CPU and 8GB of RAM.

RTGEM learning performance

We start with an experiment designed to evaluate the performance of the learning algorithm for RTGEM presented in Section 4.3.1. To this purpose, we have considered a benchmark consisting in several random RTGEM (which we call references) generated with a prefixed complexity. From these RTGEM, we generate data samples of different length and use the learning technique of Section 4.3.1 on those samples to learn new RTGEMs. We then measure the proximal distance between the learned RTGEMs and their corresponding reference RTGEM to evaluate the accuracy of the learning algorithm w.r.t. the size of the input data, and compute the mean learning times to evaluate the performance of the learning algorithm.

For this experiment, we used RTGEMs of different complexity (number of vertices and timescales), written |C l |. For each RTGEM complexity, we generated 50 random models4 from which we sampled data of varying length, written t * . We then use each sampled data to perform the RTGEM learning algorithm.

The mean values of learning times are presented in Table 4.6, where we can see the variations in learning duration depending on the size of the sample and the complexity of the RTGEM. Note that we fix the same number of nodes (labels) for every batch of graphs for the sake of convenience. In addition, the parameters λ are drawn randomly from a set of parameters, whose values are fixed beforehand in a way not to bias the tests (not very small nor very large values). From these results we can conclude that the bigger the complexity and the sample size, the more time consuming the learning is, which is not surprising. Still, we can see that learning of RTGEM can be performed on a standard computer with acceptable performance even with large data representing large RTGEM.

H H H H H H t * Dim(M) |C l | = 8 |C l | = 12 |C l | = 32 |C l | = 44 |C l | =
In Figure 4.7, we show plots representing the variation of the mean of the normalized distances between the learned and original RTGEM in function of the size of the sample for different complexities. The standard deviation is indicated by the shadowed curve around the plot for each complexity. The normalized distance is the SHD * (re f erence, learned) divided by the maximum possible distance between two graphs of the given complexity, i.e., the square of the total number The main observations that can be made based on this plot is the common trend between all graphs, from every complexity, to get closer to the reference model as the sample size increases. Moreover, we can see that the more complex the reference model is the harder it is to learn a model that is close to it (more data is needed). However, we can observe (mainly between t * = 10000 and t * = 25000) some fluctuations that are due to the learning of inaccurate models.

We conjecture that this is happening when there is enough data to select an RTGEM that models inaccurate behaviors, but still not enough data to learn the accurate dependencies. Some explanations and more detailed experiments are given in [START_REF] Antakly | [END_REF].

Statistical model checking RTGEM

We now move to an evaluation of the second step of Algorithm 1: verifying whether the fittest model satisfies the given security property. Note that sampling performance is an important part of SMC. In [Ant20], we discuss and evaluate an optimization of the sampling strategy that is

Experimentations

A λ A = 0.1 B λ B,0 = 1 λ B,1 = 0.7 C λ C,0 = 0.5 λ C,1 = 0.3 D λ D,0 = 0.2 λ D,1 = 0.1 (0, 2] (0, 5] (0, 2]
Figure 4.8 -The RTGEM M 1 on which we apply statistical model checking to verify the query ϕ 1 used in the following experiments, which we call early stopping. We start with the verification of a simple property in Example 15 by using a quantitative approach (with Monte-Carlo simulations) and a qualitative approach using the Sequential Probability Ratio Test (SPRT).

Example 15. Consider the graph of the RTGEM M 1 in Figure 4.8 and the query ϕ 1 = 100 (B[0, 5]), meaning that event B must occur at least every 5 time units during a 100 time units period. We want the model to satisfy the security property P(ϕ 1 | M 1 ) > 0.80. For the sake of coherence we always use the name security query to refer to the linear property ϕ, and the name security property for the qualitative inequality P(ϕ | M) > c with M the model and c a threshold (c ∈ [0, 1]).

The first experiments targets the quantitative approach based on Monte-Carlo simulations first intoduced in Section 3.3.1. A precision (δ) and an error rate (ε) must be fixed beforehand. The number of simulations (N) is also computed beforehand via the formula:

N = ln(2) -ln(ε) 2δ 2
This technique allows to compute an estimate p(ϕ 1 |M 1 ) to P(ϕ 1 | M 1 ) with guaranteed precision. Table 4.9 shows the results for this quantitative test. Unfortunately, to the best of our knowledge, there are no exact techniques that are adapted to this kind of formalism in order to compute the real probability P(ϕ 1 | M 1 ) and compare it with the obtained estimations.

Remark that although the property we are verifying is qualitative, it is also possible to apply a quantitative approach, have a numerical estimate (with a certain precision) and afterwards compare the estimate with the desired threshold. The results for this approach confirm that, as expected, for a high precision (δ) and a low error rate (ε) the number of simulations becomes high.

The second experiment targets the qualitative SPRT approach. In this context, we use a centered indifference region of size δ, centered in 0.8 (the target value) and strength parameters (α, β) (see Section 4.2.1). Contrary to the quantitative approach, the required number of simulations is not fixed beforehand, but rather determined dynamically, depending on the outcome of each simulation. This technique computes a decision about the security property P(ϕ 1 | M 1 ) > 0.8. Remark that the comparison of the performance of the quantitative and qualitative approach is situational 5 . Their relative performance can vary depending on the model and the property, although it is common knowledge that the qualitative approach is usually orders of magnitude faster if the estimated probability is not too close to the fixed threshold. Since we do not know beforehand whether this condition is satisfied, we will use both techniques in parallel in the following experiments, allowing the quantitative technique to conclude when the qualitative technique does not converge fast enough.

5. Further experiments can be found in [START_REF] Antakly | [END_REF] Chapter 4 -Learning and Verification of Graphical Event Models 4.13. For the sake of simplicity we do not represent the parameters of the label "Logout" for all of the remaining figures, since it is not a target node with regards to ϕ. does not satisfy the security property. Hence, the space exploration technique is repeated twenty separate times in this test, in order to check whether every time we find a model that satisfies the property. Similarly, we check whether it is the same proximal secure model that is found every time or if there are different possible ones.

The space exploration is executed on the learned model M o KO 1 and only one proximal secure model (M * KO 1 ) is found in its neighborhood (see Figure 4.16). As previously mentioned, we use the Random Walk technique to explore the neighborhood of the learned model. The maximum number of allowed modifications (chosen operators) in this test is set to 25. Remark that the proximal secure model does not make sense w.r.t. the real-life application. Indeed, our aim is only to find a mathematical object, a combination of parameters, that are close to the fittest model M o KO 1 and provide samples of behaviors which are "secure" with regards to the property.

In Table 4.17, we show some results of the algorithm for this test. We write o for the average number of applied modifications before finding a proximal secure model, τ for the average duration of the space exploration process (including the formal verification at each step), and n for the number of distinct proximal secure models found. The success rate is the number of times we found a proximal secure model over the total number of executed Random Walks. We use the same notations for the remaining test. In Table 4.22, we show some results of the algorithm for this test. We use the same notations as in the previous case. The last column shows the distances for the two proximal models M * KO 1 and M * KO 2 in this order. Notice that the success rate is not perfect in this case: three times out of twenty we did not find a proximal secure model (failure rate of 3/20). The proximal secure model

Discussion

These experiments have confirmed that our strategy for evaluating the security of a dataset w.r.t. a given security property is usable in practice for models of a medium complexity. However, we also remark from the performance evaluation on the KO 2 case that this strategy might not scale for models of a higher complexity.

Nevertheless, this experiment has allowed us to confirm that the main bottleneck of our strategy is the verification of the intermediate models, in particular due to the complexity of sampling data from RTGEMs. Indeed, we can see in Table 4.23 that the mean duration for verifying a secure model representative of the given data is around 10 minutes. Since an average of 88 operators were necessary to obtain the secure model, we can conclude that a consequent part of the total time (1h 15min) has been dedicated to running SMC on intermediate models.

Moreover, the experiment also confirms our intuition that finding a secure proximal model and computing the distance to the fittest model is more efficient when the input data corresponds to a model that is closest to satisfying the security property. Indeed, we can see that the whole procedure was performed in an average of 14 min in the KO 1 case while it took an average of 1h15min for the KO 2 case. In order to verify that KO 1 and KO 2 correspond to our intuition (i.e., KO 1 closest to verifying P(ϕ | M o KO 1 ) > 0.8, we have applied quantitative statistical model checking to the models used for generating both datasets. We obtained the following estimations: p(ϕ | M o KO 1 ) = 0.728 and p(ϕ | M o KO 2 ) = 0.585. Of course, the above observations have to be taken lightly as they only reflect a few experiments, and would need a lot more for confirmation.

Conclusion

In this chapter, we have presented a strategy for measuring the security of a model observable through a given set of data w.r.t. a security property. To do that, we have used the RTGEM formalism along with a learning algorithm based on a forward and backward greedy search. Our models have been verified against the property using SMC techniques, and we have used a random walk technique to explore the neighbourhood of the learned model to find a secure one when this was needed.

Besides the immediate benefits of being able to rate the dangerousness of a given set of data, this work is, to the best of our knowledge, one of the first attempts at putting together learning and formal verification techniques without using automata-based models as an in-between (as is 4.6. Perspectives done in [START_REF] Leucker | Learning meets verification[END_REF][START_REF] Neider | Applications of automata learning in verification and synthesis[END_REF]). This confirms that SMC is an appropriate verification technique even in the context of data-driven systems, where automata-based modeling formalisms are often inconvenient. This confirmation opens the way to many interesting perspectives, which we now present.

Perspectives

The technique presented in this chapter is a first attempt at verifying systems based on a given set of observations. It mainly serves as a proof of concept showing that SMC verification and RTGEM-based learning can be combined toward this end. However, there are many ways in which this technique can be improved, and many other domains to which it can be applied.

Potential improvements. We have observed that the bottleneck of our strategy is related to RTGEM sampling during the verification phase. In particular, we have observed that an improvement to sampling performance can be obtained by limiting the sampling to nodes related to the property under investigation (see [START_REF] Antakly | [END_REF] for details). We conjecture that other techniques could be used to drastically enhance the sampling/verification performance such as taking advantage of partial samples drawn from unmodified parts of a previous RTGEM while exploring its neighborhood.

Similarly, we have chosen a random walk technique for exploring the neighborhood of the fittest model. While this technique has been slightly optimized by starting the random search with a reduced number of operators, we conjecture that one could take a better advantage of the given property in order to produce heuristics that would considerably reduce the required number of intermediate models that are verified before finding a secure proximal model.

Another potential improvement concerns the notion of distance between RTGEMS. Indeed, the notion of distance we have defined, which allows us to quantify the security level of a model, is not directly related to security. Although it provides a measure of how far a given set of data is from a secure model, this measure highly depends on the data itself and the complexity of the learned model. As a consequence, the absolute value obtained as an output for one set of data is difficult to compare to the one obtained for a different set. In the future, we plan to investigate other notions of distance that will allow us to compare the security levels obtained for different sets of data.

Other application domains. The technique developped in this chapter could be applied to other goals than security. In particular, RTGEM are a formalism that could be very useful in other scientific domains, such as oceanography. In the context of the TARA Ocean community, in particular, a profusion of heterogenous timed data is collected and studied through statistical analysis. Using graphical models such as RTGEMs to represent this data would be particularly interesting as it would allow analyzing timing dependencies between the observed variables (which is mostly not done today). Moreover, the performances observed during the experiments performed in this chapter makes us confident that RTGEMs would be able to handle sufficiently large datasets to be usable in this context. While we observed that sampling of RTGEMs is hard, it is nothing compared to the sampling of complex natural system models (e.g., ocean biogeochemical models [AÉT + 15]). While some scientists have proposed using neural networks to simulate the outputs of climate models with high speed and precision [START_REF] Hutson | Ai shortcuts speed up simulations by billions of times[END_REF], the lack of understanding that still lies behind neural networks is a problem. An intermediate solution would be to use graphical models such as RTGEMs, which are well understood and give precise insight on the phenomena they represent, to enhance the speed of complex models sampling while preserving understandability of the model itself.

PERSPECTIVES

As can be seen from this document, my contributions in the past 10 years have followed two (somewhat separate) main lines of work. On the theoretical side, I have developed and studied modeling formalisms and verification techniques dedicated to the analysis of discrete and continuous parametric probabilistic systems. On the practical side, I have developed statistical methods and tools, and applied them to the learning and verification of complex systems. Unfortunately, in most cases, the modeling formalisms we develop on the theoretical side are not the ones we use on the practical side. Three main reasons stand out for this gap between theory and practice:

• The verification methods developed in the context of abstract modeling formalisms such as the ones presented in Chapters 1 and 2 are computationally expensive and often fail to scale up to models of industrial-size.

• The act of modeling a complex system requires deep expertise both on the system itself and on the chosen modeling formalism. Unfortunately, those who develop modeling formalisms and analysis techniques (who we will call model developers) rarely possess the domain expertise required to use these formalisms on concrete complex systems. On the other hand, those who use models in practice (who we will call modeling practitioners) rarely possess the theoretical expertise required to use abstract modeling formalisms.

• The existing barrier between natural/social sciences and formal sciences is such, that modeling formalisms are often not perfectly tailored to the use of modeling practitioners, and modeling practitioners are even rarely aware of the existence of such modeling formalisms.

While theoretical advances on modeling formalisms and verification techniques are without doubt important, I believe that the long-term aim of such research should be to provide better tools for the analysis of concrete systems. In order to do that, researchers in formal sciences need to make verification techniques more efficient and lessen the effort required from concrete modeling practitioners for using our modeling formalisms. During the past years, numerous interactions with Biologists and Oceanographers, in particular in the context of the Tara Ocean GO-SEE research federation 6 , have allowed me to refine the following research axes, which I believe could contribute to work toward this aim.

Enhancing parametric modeling and analysis for natural sciences

From the practical point of view, (parametric) statistical model checking stands out as a very promising technique for analyzing complex natural models. One of the main advantages of SMC in this context is the limited number of requirements on the chosen modeling formalism. Indeed, as we have shown in [REG + 20], SMC can be applied with very light adaptations to any stochastic system modeled as a computer program, which is common in natural sciences. More requirements are needed in the context of parametric SMC, but the prototype tool we presented in [BAD + 19] shows that Python programs where parameters are expressed in a certain form (i.e., inheriting from predefined Python classes) can also be analyzed. Again, this can be generalized to other programming languages.

However, the types of parameters that can be taken into account in these models is not completely satisfactory. Indeed, the parametric SMC theory we developed is based on parametric Markov chains, where parameters can only range over transition probabilities. On the other hand, models of natural systems are often based on sets of parameterized differential equations, where the nature of the parameters is substantially different. Even when the differential equations are discretized, translating the equation parameters into transition probabilities is not an easy task. This is the main reason why [REG + 20] is based on standard and not parametric SMC. Nevertheless, the only requirements about the applicability of (parametric) SMC are that (1) the model is executable, and (2) once the parameter values are fixed, the resulting model is purely stochastic. These requirements are met by models consisting of sets of parameterized differential equations, therefore there is no theoretical reason why (parametric) SMC methods could not be adapted to this context. In order to do so, we need to study: (i) new modeling formalisms where differential equation parameters can be taken into account as such, (ii) the form of the probability distribution on their traces, and (iii) how the results presented in Chapter 3 can be adapted to this new context.

Developing new efficient abstraction techniques

Another limitation of (parametric) SMC is the number of simulations that are required in order to produce accurate results. Indeed, in the context of complex natural systems, producing simulations of a complete model can be computationally expensive. While a sufficient reduction of the number of required simulations is out of the picture 7 , we believe that the solution lies in the development of abstraction methods that will reduce the computational complexity of performing simulations. Several such abstraction techniques can be envisioned. In particular, we plan to pursue the work initiated during my thesis on stochastic abstraction [BBB + 12], a technique that allows to drastically reduce the size of the model in order to perform fast simulations. While the technique shows promise, the difficulty in this context is to be able to link the formal precision guaranteed by SMC algorithms with the loss of details resulting of the abstraction. Another lead, promoted in recent works [START_REF] Hutson | Ai shortcuts speed up simulations by billions of times[END_REF] is to use AI techniques based on neural networks in order to build a new model of the system whose outputs are statistically equivalent to those of the original one but that is orders of magnitude faster to simulate. However, the lack of understanding that still lies behind neural networks is a problem. Instead, we plan on developing such abstraction techniques based on other graphical probabilistic models such as the ones presented in Chapter 4, which yield more understandable models. Beyond that, we will use our expertise on compositional interface theories [Del10, CDL + 10, CDL + 11, DFH + 12, DFLL13b, BDF + 13, DKL + 13, FLDL18] to study compositional versions of this technique, where only parts of the model are abstracted, in the spirit of what we did with stochastic abstraction [BBB + 10a, BBB + 12].

Building a verification platform for practitioners

The echos we have received from the Biological community for the results presented in [REG + 20] have convinced us of the interest that a tool for automatic analysis and parameterization of natural systems expressed in standard programming languages would raise. Our intent in the near future is therefore to automatize the processes we have followed in [REG + 20] in order to adapt the original model for SMC analysis and to build a platform that could be used by any modeling practitioner interested in such analysis, without further requirements. Our intent, as we did for the prototype tool presented in [BAD + 19] is to propose two front-ends for this platform: one dedicated to the formal analysis of parametric "computer-science" models (i.e., for the use of formal science researchers) where state-of-the-art verification techniques can be implemented 7. There exist statistical techniques for reducing this number of simulations, but they are limited and their effect is not sufficient compared to the computational complexity of producing a meaningful number of simulations. and tested; and the second one dedicated to the concrete analysis/verification of models of natural systems (i.e., for the use of modeling practitioners). Besides the development of new verification techniques, this platform raises many scientific and technical challenges on several topics such as parallelization, interpretation and rendering of dynamic partial results during a potentially very long verification, or visualization of meaningful information on high-dimensional sets of parameters.

These research axes are and will be the topics for multidisciplinary research projects (including Ph.D proposals) that will last several years and require support fundings. They are in line with the future priority research plan Océans of the CNRS that is at the moment being refined in the context of the corresponding Task Force (to which I participated). prendre en compte, en plus du non-déterminisme, les choix probabilistes. Contrairement aux travaux existants, dans lesquels des choix probabilistes pouvaient uniquement être introduits au travers des affectations de variables, notre contribution permet d'introduire ces choix probabilistes à la place de n'importe quel choix non-déterministe présent dans le système (par exemple, en plus des affectations probabilistes, dans le choix entre les évènements ou encore dans la définition des paramètres). De plus, notre théorie autorise l'introduction de ces choix probabilistes à n'importe quelle étape du raffinement du système et s'intègre avec la théorie de preuves associée au B évènementiel. Nous introduisons ainsi la notion de raffinement probabiliste et proposons une étude de la notion de convergence presque certaine dans ce cadre. La troisième contribution, d'orientation plus pratique, concerne la vérification de modèles exprimés dans le formalisme des chaînes de Markov à paramétrées (pMC). Pour cette contribution, nous partons de l'observation que la plupart des méthodes de vérification pour les pMC sont complexes et ne passent pas à l'échelle. C'est aussi le cas pour des formalismes de modélisation moins abstraits, tels que les chaînes de Markov (MC), mais il existe dans ce cas une technique statistique formelle permettant d'estimer la probabilité avec laquelle une MC satisfait une propriété donnée de manière très efficace. Cette technique, la vérification de modèle statistique, était jusqu'ici limitée aux modèles purement probabilistes et ne pouvait donc s'appliquer dans le cadre des pMC. Dans notre contribution, nous étendons la vérification de modèle statistique au cadre des pMC en s'inspirant de techniques statistiques telles que l'échantillonnage préférentiel. Nous présentons de plus un prototype d'outil, développé pour l'occasion, qui fait la preuve de concept que cette nouvelle technique est bien applicable en pratique et possède des propriétés intéressantes en termes de performances. Afin de montrer que notre technique passe à l'échelle, nous l'appliquons à un cas d'étude industriel : l'analyse du plan de vol d'un drone civil. Finalement, la dernière contribution, d'orientation plutôt pratique elle aussi, s'attaque au problème de la modélisation automatisée. En effet le processus d'écriture d'un modèle est souvent fastidieux et requiert de l'expertise humaine. Dans cette contribution, notre but est le développement de techniques automatiques permettant de créer un modèle reproduisant les caractéristiques statistiques d'un ensemble de données. Nous présentons ainsi des techniques d'apprentissage reposant sur un formalisme de modélisation à temps continu : les "Recursive Timescale Graphical Event Models" (RTGEM). Dans ce contexte, nous proposons un algorithme permettant de générer automatiquement un modèle représentant le plus fidèlement possible les données en entrée tout en satisfaisant un ensemble de propriétés de sécurité fixées. La comparaison de ce modèle avec le modèle le plus représentatif des données permet ainsi d'évaluer, de manière automatique, le degré de sécurité associé aux données en entrée.

Markov chains and associated notions

Definition 2 (Markov Chain -see p.27). A Discrete Time Markov Chain (MC for short) is a tuple M = (S, s 0 , p, L), where S is a finite set of states containing the initial state s 0 , L : S → 2 A is a labelling function, and p : S × S → [0, 1] is a probabilistic transition function such that for all s ∈ S, the function p(s, .) is a distribution (i.e., s → p(s, s ) ∈ Dist(S)).

• Given a state s ∈ S, L(s) is called the label of s.

• MC is the set containing all discrete time Markov chains.

• A run of M is a sequence w = s 1 , s 2 , . . . with p(s i , s i+1 ) > 0, for all i.

Definition 3 (Probabilistic Bisimulation -see p.28). Let M = (S, s 0 , p, L) be a Markov chain. A probabilistic bisimulation on M is an equivalence relation R on S such that for all pair of states (s 1 , s 2 ) ∈ R , we have -L(s 1 ) = L(s 2 ), and p(s 1 , T ) = p(s 2 , T ) for each equivalence class T ∈ S/R .

M 1 and M 2 are bisimilar if and only if there exists a probabilistic bisimulation R on M 1 ∪ M 2 with (s 1 0 , s 2 0 ) ∈ R .

• Given a Markov chain M , the probability of a finite sequence of states w = s 0 , s 1 , . . . , s n , written P M (w) is the product of the probabilities of the transitions transitions involved in this sequence, i.e., P M (w) = p(s 0 , s 1 ) • p(s 1 , s 2 ) • . . . • p(s n-1 , s n ).

• reach s 0 (s) = { w ∈ S * | w = s 0 , . . . s n with s n = s and s i = s ∀0 ≤ i < n } is the set of runs starting in s 0 and ending in s (without going through s before the end).

• P M (♦s) = ∑ w∈reach s 0 (s) P M (w) is the overall probability of reaching state s in M (if s = s 0 ). This notation is extended to the probability of reaching a given state label or a given atomic proposition (see p.29).

Markov chain abstractions

Definition 4 (Abstraction Model -see p.29). A Markov chain abstraction model (an abstraction model for short) is a pair (L, |=) where L is a nonempty set of models and |= is a relation between MC and L. Let P be in L and M be in MC. We say that M implements P if and only if (M , P ) belongs to |= (i.e., M |= P ).

Parametric Markov chains

Definition 5 (Parametric Markov Chain -see p.30). A Parametric Markov Chain (pMC for short) is a tuple I = (S, s 0 , P, L, X) where S, s 0 , and L are defined as for MCs, X is a set of variables (parameters) ranging over [0, 1], and P : S × S → Q X associates with each potential transition a parameterized probability.

¢

Cut here ¢ • pMC is the set containing all parametric Markov chains.

• The satisfaction relation |= p between MC and pMC is defined by M |= p I if and only if S = S , s 0 = s 0 , L = L , and there exists a valuation v of X such that p(s, s ) equals v(P(s, s )) for all s, s in S.

Interval Markov chains

Definition 6 (Interval Markov Chain -see p.31). An Interval Markov Chain (IMC for short) is a tuple I = (S, s 0 , P, L), where S, s 0 , and L are defined as for MCs, and P : S × S → I associates with each potential transition an interval of probabilities.

Definition 7 (Once-and-for-all Semantics -see p.31). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the once-and-for-all semantics, written M |= o I I , if and only if (T,t 0 , L M ) = (S, s 0 , L I ) and for all reachable states s in M and all state s ∈ S, p(s, s ) ∈ P(s, s ).

Definition 8 (IMDP Semantics -see p.31). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the IMDP semantics, written M |= d I I , if and only if there exists a mapping π from T to S such that

(1) π(t 0 ) = s 0 ,

(2) L I (π(t)) = L M (t) for all states t ∈ T , and

(3) p(t,t ) ∈ P(π(t), π(t )) for all pairs of states t,t in T , where t is reachable in M .

Definition 9 (At-every-step Semantics -see p.32). A MC M = (T,t 0 , p, L M ) satisfies an IMC I = (S, s 0 , P, L I ) with the at-every-step semantics, written M |= a I I if and only if there exists a relation R ⊆ T × S such that (t 0 , s 0 ) ∈ R , and whenever (t, s) ∈ R , we have (1) the labels of s and t correspond: L M (t) = L I (s),

(2) there exists a correspondence function δ (t,s) : T → (S → [0, 1]) such that (a) ∀t ∈ T if p(t,t ) > 0 then δ (t,s) (t ) is a distribution on S (b) For all s ∈ S, (Σ t ∈T p(t,t ) • δ (t,s) (t )(s )) ∈ P(s, s ), and (c) For all (t , s ) ∈ T × S, if δ (t,s) (t )(s ) > 0, then (t , s ) ∈ R .

Parametric interval Markov chains

• X is a finite set of parameters and v is a valuation over X.

• I(Q X ) is the set of all parametrized intervals over [0, 1].

• For all J = [ f 1 , f 2 ] ∈ I(Q X ), v(J) denotes the interval [v( f 1 ), v( f 2 )] if 0 ≤ v( f 1 ) ≤ v( f 2 ) ≤ 1 and the empty set otherwise.

¢

Cut here

¢

Properties Existential consistency

• A pIMC P = (S, s 0 , P, L, X) is existentially consistent if and only if there exists a MC M satisfying P .

• C ∃c (P ) is a CSP encoding for verifying the existential consistency of pIMC P , where the variables are as follows:

-one variable π p with domain [0, 1] per parameter p in X,

-one variable θ s s with domain [0, 1] per transition (s, s ) in { { s } × Succ(s) | s ∈ S }, and -one Boolean variable ρ s per state s in S.

Qualitative reachability

• A state label Γ is existentially reachable in pIMC P = (S, s 0 , P, L, X) if and only if there exists an implementation M of P where Γ is reachable (i.e., P M (♦Γ) > 0).

• Γ is universally reachable in P if and only if Γ is reachable in any implementation M of P .

• C ∃r (P ) is a CSP encoding, that extends C ∃c (P ), for verifying these properties. In addition to the variables of C ∃c , C ∃r contains:

-one integer variable ω s with domain [0, |S|] per state s in S.

Quantitative reachability

• C ∃r (P , Γ) is a CSP that extends C ∃r (P , Γ) with the folowing variables:

-one Boolean variable λ s per state s in S, and -one integer variable α s with domain [0, |S|] per state s in S.

• C ∃r (P , Γ) is a CSP that extends C ∃r (P , Γ) with the following variables:

-one variable γ s with domain [0, 1] per state s in S.

¢

Cut here ¢ I ( v) is the invariant, PEvts is a set of probabilistic events and Init is the initialisation event.

Semantics

Let M=( v, I( v ), PEvts, Init ) be a fully probabilistic Event-B model and σ be a valuation of its variables.

• Notation σ refers to variable valuations.

• Notation θ refers to parameter valuations when parameters are involved.

• Given a variable x ∈ v, we write [σ]x for the value of x in σ.

• Given an expression E( v), we write [σ]E( v) for the evaluation of E( v) in the context of σ.

• Given an expression E( t, v) over variables and parameters, we write [σ, θ]E( t, v) for the evaluation of E( t, v) under parameter valuation θ and variable valuation σ.

• Given a probabilistic event e i with a set of parameters t and a valuation σ of the variables, we write T e i σ for the set of parameter valuations θ such that the guard of e i evaluates to true in the context of σ and θ.

• P T e i

σ is the uniform distribution on the set T e i σ .

• E e i (x) is the set of all expressions that can be assigned to the variable x in the context of an enumerated probabilistic assignment (x := E 1 ( t, v)@ p 1 ⊕ . . . ⊕ E m ( t, v)@ p m (m ≥ 1)).

• The probability of choosing an expression E i in the above set is written P e i x (E i ) = p i . • Given a probabilistic event e i , we write Var(e i ) for the set of variables in v that are modified by the action of e i .

• If a variable x is modified by an enumerated probabilistic assignment of e i , then we write E e i (x)| σ σ,θ for the set of expressions in E e i (x) such that their evaluation in the context of σ and θ returns the value of x in the valuation σ .

• If a variable x is modified by a predicate probabilistic assignment (x : ⊕Q x ( t, v, x )), then we write V e i θ,σ (x) for the set of values x that make the predicate Q x ( t, v, x ) true when evaluated in σ and θ.

• P e i σ,θ (x, σ ) is the probability that a variable x is assigned the new value [σ ]x when executing probabilistic event e i from the valuation σ and with parameter valuation θ.

Definition 13 (Fully Probabilistic Event-B operational Semantics -see p.78). The operational semantics of a fully probabilistic Event-B model M=( v, I( v ), PEvts, Init ) is a PLTS [[M]] = (S, Acts, P, s 0 , AP, L) where S, Acts, s 0 , AP, and L are defined as in the standard LTS semantics of Event-B models (see Section 2.2.2), and P : S × Acts × S → [0, 1] is the transition probability function such that for a given state s, for all e i , s ∈ Acts × S, we have P(s, e i , s ) = 0 if e i / ∈ Acts(s) or ∃x ∈ X\{Var(e i )} st [s]x = [s ]x and otherwise • Given a pMC M and a parameter valuation v ∈ R X , we write P v for the transition probability function obtained under valuation v and M v for the resulting structure.

• If M v is a Markov chain (i.e., if P v satisfies the conditions of a Markov chain transition function), then we say that v is a valid parameter valuation w.r.t. M .

• A run ω of M is a sequence of states s 0 , s 1 , . . . such that for all i ≥ 0, P(s i , s i+1 ) = 0 (i.e., the probability is either a strictly positive real constant or a function of the parameters).

• Γ M (l) is the set of all finite runs of length l and Γ M is the set of all finite runs of a pMC.

Monte Carlo analysis

• Given a pMC M , we use the notation r : Γ M → R for the reward function that we want to estimate using Monte Carlo.

• We use the notation f for the normalization function used in importance sampling : f :

S × S → [0, 1].
We says that f is valid w.r.t. pMC M if for all states s, ∑ s ∈S f (s, s ) = 1.

• Given a pMC M and a valid normalization function f , we write M f for the MC obtained from M by replacing P by f .

• u M is the uniform normalization function, that yield uniform distributions on successor states following the original structure of M .

• Pa : Γ M → Poly(X) is a parametric reward function defined inductively as follows:

Pa(s 0 ) = 1 and for all run ω, s, s ∈ Γ M , Pa(ω, s, s ) = Pa(ω, s)P(s, s ).

• r is the parametric reward function that we use for computing the expectation of r and is defined as follows: Given any valid normalization function f and any run ω ∈ Γ M , we define r (ω) = Pa(ω)

P M f (ω) r(ω).
• Y is a random variable, evaluated on runs of M f , such that Y = r (ω) = Pa(ω) P M f (ω) r(ω). • Given an experiment consisting of a set of runs {ω 1 , . . . , ω n }, we write {Y 1 , . . .Y n } for the set of corresponding random variables.

• γ and σ 2 are the parametric functions giving, respectively, the mean value and variance of the variable Y i .

• Given an unknown value/function such as γ, we write γ for its estimation using our statistical methods.

• ϕ ϕ ϕ(z) = z -∞ exp(-x 2 /2)dx/( √ 2π) is the cumulative distribution function of the standard normal distribution N (0, 1).

• The sufficient statistics are formally written as: n t * ,l, j (x t * ) = • The "add" operator adds a non-existing edge to a model and its corresponding timescale T = (0, c], with c a constant (also called horizon).

• The "split" operator splits one interval (a, b] in the timescale of a chosen edge into two intervals (a, a+b 2 ], ( a+b 2 , b]. • The extend operator extends the horizon of a chosen edge by adding the interval (t h , 2t h ], with t h being the previous horizon.

• The symmetric backward operators are O -1 F = {reverse_add, reverse_split, reverse_extend}. • A finite consistent RTGEM learning procedure [START_REF] Gunawardana | Universal models of multivariate temporal point processes[END_REF] is to do a forward greedy search (for model construction) followed by a backward greedy search (for model refinement), both based on model selection, and using data that is faithful.

• We choose to use the Bayesian Information Criterion (BIC) [S + 78], written as follows for a model M: S t * (M) = log(p(x t * | t * ; M, λ t * ,l, j (x t * ))) -∑ l∈L

| C l | •log(t * )
where λ are the estimates of the λ parameters.

• Given a model M, the maximum likelihood estimate λ of the λ parameters for M (knowing its structure) is as follows [GM16] : λ t * ,l, j (x t * ) = n t * ,l, j (x t * ) d t * l, j (x t * ) .

Distance between RTGEM

• A timescale in an RTGEM can be represented by a vector v = [0, a, b, c, . . .] containing the successive endpoints values. Abstract: The goal of science in general is to study the functioning of complex systems in order to better understand them and potentially to be able to predict their future behaviour. An important part of this study is the realization of a model, an abstract object that faithfully represents existing knowledge about the system, and which can then be studied and simulated instead of the system itself. Unfortunately, existing knowledge about systems is often incomplete or subject to uncertainties. It is then important to include these uncertainties in the models and to develop techniques to automate the analy-sis of these models. This document presents four contributions to the field of modeling and analysis of such systems. Two of these contributions are somewhat theoreticaly-oriented, proposing modeling languages to include uncertainties in models and associated verification techniques. The other two contributions are more practical, one proposing a statistical verification technique and proving its efficiency on an industrial case study, and the other proposing a technique for automatied model construction from a data set and an analysis of its security level.

  (1) the study of a compositional specification theory for probabilistic systems called Constraint Markov Chains [CDL + 11]; and (2) the development and use of a statistical verification technique -called statistical model checking -on concrete complex systems [BBB + 10a]. 2010-2013: Postdoctoral positions During my years as a postdoctoral student, I have pursued these two lines of work. In the context of specification theories, I have enhanced the theory of Constraint Markov Chains to take into account non-determinism and modalities, yielding the formalism called Abstract Probabilistic Automata (APA) [DKL + 11a, DKL + 13]. This formalism has then been studied under all angles, yielding algebraic compositional results [DKL + 11b, DLL13, DFLL13a, DLL14, DFLL14a] and a prototype tool for their verification [DLL + 11a]. In the same line, I have also investigated other abstract formalisms for describing non-deterministic [BDF + 13] and timed [DFH + 12, DFLL13b, Introduction DFLL14b, FLDL18] systems. On the other hand, I have also pursued my work on statistical model checking, using this technique to verify different types of systems [BBB + 12] and contributing to the development of a tool [BBD + 12, NBB + 15].

  Figure 1.1 -MC M 1

Definition 4 (

 4 Markov chain Abstraction Model). A Markov chain abstraction model (an abstraction model for short) is a pair (L, |=) where L is a nonempty set of models and |= is a relation between MC and L. Let P be in L and M be in MC. We say that M implements P if and only if (M , P ) belongs to |= (i.e., M |= P ). When the context is clear, we do not mention the satisfaction relation |= and only use L to refer to the abstraction model (L, |=).

Example 6 (

 6 Example 5 continued). Consider the IMC I from Figure 1.3, the MC M 1 from Figure 1.1, the MC M 2 from Figure 1.4, and the MC M 3 from Figure 1.5. We have that M 1 satisfies I w.r.t. |= o I and we say that M 1 has the same structure as I . Trivially, we also have that M 1 satisfies I w.r.t. |= d I and |= a I .

Figure 1 Figure 1

 11 Figure 1.4 -MC M 2 satisfying the IMC I from Figure 1.3 w.r.t. |= d I

Figure 1

 1 Figure 1.8 -IMC I with three pMCs P 1 , P 2 , and P n .

Figure 1

 1 Figure 1.9 -MCs satisfying the IMC and/or one of the pMCs of Figure 1.8.

Figure 1 .Figure 1

 11 Figure 1.10 -Variables in the CSP produced by C ∃c for the pIMC P from Fig. 1.7

C

  Figure 1.18 -Comparing encoding time for the existential consistency problem

C

  Figure 1.19 -Comparing solving time for the existential consistency problem

Chapter 1 -

 1 Parametric interval Markov chains

Figure 2

 2 Figure 2.1 -How to introduce probabilities within Event-B

  e a = any t where G a (t, v) then S a (t , v) end e c = refines e a any ū where G c ( ū, w) with Wi c (t , ū, v, v ', w, w') then S c ( ū, w) end

  Figure 2.2 -Initial Event-B model of the simple P2P protocol

  Figure 2.3 -First refinement of the simple P2P protocol: step-by-step download

Figure 2

 2 Figure 2.4 -Extract of the transition system of the first refinement of the simple P2P protocol, with N=2 and K=2

  Figure 2.5 -New event introduced in the second refinement step

Figure 2 . 7 -

 27 Figure 2.7 -Probabilistic version of the simple P2P protocol

Figure 2 . 8 -

 28 Figure 2.8 -Extract of the Detailed construction of the MC of the simple P2P protocol, with N=2 and K=2

Figure 2 . 9 -

 29 Figure 2.9 -Extract of the MC of the simple P2P protocol, with N=2 and K=2

Figure 2

 2 Figure 2.10 -Probabilistic Event-B model and MC semantics illustrating the necessity of bounding event weights to ensure almost-certain convergence

Figure 2

 2 Figure 2.11 -Probabilistic Event-B model and MC semantics illustrating the necessity of bounding event parameter values to ensure almost-certain convergence

Theorem 3 .

 3 Let M=( v, I( v ), V( v), PEvts, Init ) be a probabilistic Event-B model and PEvts c ⊆ PEvts a set of convergent events. If M satisfies the above POs (1-5), then the set PEvts c almost-certainly converges.

  Figure 3.1 -pMC M 1

  2) is obtained by definition of the expected value and of the distribution of Y ; (3.3) is obtained by definition of the random variable Y ; (3.4) is direct because we only consider runs ω such that

  Figure 3.3 -Results on pMC M 1

Figure 3 . 6 -

 36 Figure 3.6 -Results for the zeroconf model obtained with our implementation (left) and taken from the PARAM website (right)

Figure 3 . 7 -

 37 Figure 3.7 -Results on the Crowds protocol model for our implementation (left) and from the PARAM website (right)

3. 5 .

 5 Application to UAV flight plan analysis will allow us (1) to analyze the drone system and detect the most important parameters, and

Figure 3

 3 Figure 3.8 -Safety zonesIn the context of software, considerations in airborne systems and equipment certification (named DO-178C) defined five levels of safety zones. We therefore chose to adapt this definition to our context and also partitioned the spatial environment of our drones into five zones, the most secure being Zone 1 and the most dangerous being Zone 5. These zones are characterized by their distance from the intended flight plan, as shown in Figure3.8.

Figure 3

 3 Figure 3.9 -Attitude coordinates Figure 3.10 -Flight Control overview

Figure 3 .

 3 Figure 3.11 -Issue on drone location and misleading positions

  [BAD + 19]), but the frequency and the accuracy of the position estimation together with details of the given flight plan (position and schedule of intermediate points) allows us to compute the effective trajectory of the drone. These computations have been taken into consideration in a global model which we now describe.

Figure 3 .

 3 Figure 3.12 -Global behaviour of the FCS

  Figure 3.14 -A simple pMC.

Figure 3 .

 3 Figure 3.15 -Impact of the choice of the normalization function f on the size of confidence intervals

Figure 4 . 1 -

 41 Figure 4.1 -Example showing the decomposition of a timeline into marked point processes

Figure 4 . 2 -

 42 Figure 4.2 -Example of a 4 labels GEM

  the graph of figure 4.2 for instance, one can identify that label B has two parents A and D, hence λ B (t | h) depends on the past history of A and D. However, label D does not have any parents and subsequently has a constant rate λ D (t | h) = λ D of occurrence.

A

  λ A,c l (h,t) B λ B,c l (h,t) C λ C,c l (h,t)

Figure 4 . 3 -

 43 Figure 4.3 -Example of a 4 labels TGEM

Figure 4 . 5 -

 45 Figure 4.5 -Example from (a) to (d) showing how an RTGEM (we omitted the timescales for a better representation) can be transformed into ordered SCCs

  2 ) and an edge e ∈ E 1 ∩ E 2 , we write v e 1 and v e 2 for the values of the respective timescales of e in G 1 and G 2 . We then write v e id = |v e 1 ∩ v e 2 |, for the number of identical endpoints in the two vectors; and v e nid = |v e 1 \ v e 2 | + |v e 2 \ v e 1 |, for the number of endpoints that are not identical in the two vectors. This allows to define an elementary distance between the two respective timescales of the edge e as follows:d(v e 1 , v e 2 ) =v e nid v e nid + v e id (4.3) endpoints and v e nid = |V e nid | for the number of unmatched enpoints. The elementary distance between v e 1 and v e 2 can now be written d * (v e 1 ,

Algorithm 1 M

 1 Proposed Strategy input: D,ϕ output: M * , ∆ 1: M o = argmax M∈RTGEM P(D | M) 2: If P(ϕ | M o ) * = find{M ∈ N (M o ), P(ϕ | M) > c} 6: ∆ = Distance(M o , M * ) 7: Return(M * , ∆) main steps. The first step (line 1) consists in learning a model M o that best represents the data D. This model is called the fittest model. The second step (line 2) consists in verifying whether the obtained model M o satisfies the security property. If this is the case, then the algorithm ends.

Chapter 4 -Figure 4 . 7 -

 447 Figure 4.7 -Mean normalized distances between learned RTGEMs with different complexities and their corresponding references, in function of the sample size t *

λFigure 4

 4 Figure 4.13 -A representation of the learned model M o OK from the set of "secure behavior" data

λ

  Figure 4.15 -A representation of the learned model M o KO 1 from the corresponding set of "not secure" data

λFigure 4

 4 Figure 4.20 -A representation of M * KO 2 1 which satisfies the security property in the neighborhood of M o KO 2

Figure 4 .( 2 (

 42 Figure 4.21 -A representation M * KO 2 2 which satisfies the security property in the neighborhood of M o KO 2

  o

6.

  The Tara Ocean GO-SEE (Global Oceans Systems Ecology & Evolution) research federation regroups 22 french and international research teams. In particular, the following organisms are involved in this federation: CNRS, CEA, Tara Ocean Foundation, Sorbonne Université, PSL, Inserm, ENS Paris, IRD, EPHE, Université d'Évry-Val-d'Essonne, Université Paris-Saclay, UPVD, AMU, Université de Toulon, École centrale de Nantes, Université de Nantes, UGA, EMBL and the Physics and Mathematics faculty of the University of Chile.

  i = l)1(c l (h i ,t i ) = j) d t * ,l, j (x t * ) = l (h i , τ) = j)dτwhere t |x t * |+1 represents the duration d between the last event occurence and the final time t * in the data. Definition 19 (Recursive timescale graphical event model -see p.133). The family of RTGEMs is defined recursively to be the finite closure of the empty model M 0 = ((L, {}), {}) under a set of allowed operators O F = {add, split, extend}.

  IMC semantics |= 1 is more general than another IMC semantics |= 2 if and only if for all IMCs I and for all MCs M if M |= 2 I then M |= 1 I . Also, |= 1 is strictly more general than |= 2 if and only if |= 1 is more general than |= 2 and |= 2 is not more general than |= 1 .

  IMC, pMC, and pIMC are three Markov chain abstraction models. In order to compare their expressiveness and compactness, we introduce the comparison operators and ≡. Let (L 1 , |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models containing respectively L 1 and L 2 . We say that L 1 is entailed by L 2 , written L 1 L 2 , if and only if all MCs satisfying L 1 satisfy L 2 modulo bisimilarity. (i.e., ∀M |= 1 L 1 , ∃M |= 2 L 2 such that M is bisimilar to M ). We say that L 1 is (semantically) equivalent to L 2 , written L 1 ≡ L 2 , if and only if L 1 L 2 and L 2 L 1 . |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models.L 1 is at least as succinct as L 2 , written L 1 ≤ L 2 , if and only if there exists a polynomial p such that for every L 2 ∈ L 2 , there exists L 1 ∈ L 1 such that L 1 ≡ L 2 and |L 1 | ≤ p(|L 2 |). Moreover, L 1 is strictly more succinct than L 2 , written L 1 < L 2 , if and only if L 1 ≤ L 2 and L 2 ≤ L 1 .

	Definition 11 introduces succinctness based on the sizes of the abstractions.
	Definition 11 (Succinctness). Let (L 1 , We start with a comparison of the succinctness of the pMC and IMC abstractions. Since pMCs
	allow the expression of dependencies between the probabilities assigned to distinct transitions
	while IMCs allow all transitions to be independent, it is clear that there are pMCs without
	any equivalent IMCs (regardless of the IMC semantics used), therefore (IMC, |= o I ) ≤ (pMC, |= p ),
	(IMC, |= d I ) ≤ (pMC, |= p ), and (IMC, |= a I ) ≤ (pMC, |= p

  that states that the |= a I and the |= o I pIMC semantics are equivalent w.r.t. existential 1.4. Qualitative properties (and universal) reachability. As for the consistency problem, we get by Proposition 1 that the three IMC semantics are equivalent w.r.t. existential (and universal) reachability. Note first that in our C ∃c encoding each ρ s variable indicates if the state s appears in the constructed Markov chain. However, the ρ s variable does not indicate if the state s is reachable from the initial state, but only if it is reachable from at least one other state (i.e., possibly different from s 0 ). Indeed, if the graph representation of the constructed Markov chain has strongly connected components (SCCs for short), then all ρ s variables in one SCC may be set to true while this SCC may be unreachable from the initial state. This is not an issue in the case of the consistency problem. Indeed, if a Markov chain containing an unreachable SCC is proved consistent, then it is also consistent without this unreachable SCC. However, in the case of the reachability problem, these SCCs are an issue. The following encoding therefore takes into account these isolated SCCs such that ρ s variables are set to true if and only if they are all reachable from the initial state. This encoding will solve the qualitative reachability problems (i.e., checking qualitative reachability from the initial state). We propose a new CSP encoding, written C ∃r , that extends C ∃c , for verifying these properties. Formally, CSP C ∃r

  contains one integer variable ω s with domain [0, |S|] per state s in S, D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

  in C ∃r ensure that ω s = k if and only if there exists in M a run from the initial state to s of length k -1 with non zero probability; and state s is not reachable in M from the initial state s 0 if and only if ω s equals to 0. Finally, constraint (8) enforces the Boolean reachability indicator variable ρ s to be set to true if and only if there exists a run with non zero probability in M from the initial state s 0 to s (i.e., ω s = 0). Let S Γ be the set of states from P labeled with Γ. Recall that C ∃r (P ) produces a Markov chain satisfying P where reachable states s are such that ρ s = true. As a consequence, Γ is existentially reachable in P if and only if C ∃r (P ) admits a solution such that s∈S Γ ρ s ; and Γ is universally reachable in P if and only if C ∃r (P ) admits no solution such that s∈S Γ ¬ρ s . This is formalised in the following proposition, whose proof is provided in [BDF + 18]. , D r ,C r ) be the CSP C ∃r (P ).-CSP (X r , D r ,C r ∪ s∈S Γ ρ s ) is satisfiable if and only if Γ is existentially reachable in P -CSP (X r , D r ,C r ∪ s∈S Γ ¬ρ s ) is unsatisfiable if and only if Γ is universally reachable in P As for the existential consistency problem, we have an exponential gain in terms of size of the encoding compared to the one we proposed in[START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]: the number of constraints and variables in C ∃r is linear in terms of the size of the encoded pIMC.Remark. In C ∃r constraints (2) inherited from C ∃c are entailed by constraints (8) and (10) added to C ∃r . Thus, in a practical approach one may ignore constraints (2) from C ∃c if they do not improve the solver performance.

	, Recall first that CSP C ∃c (P) constructs a Markov chain M satisfying P w.r.t. |= o (only for s = s 0 ) I . Informally, for each state s in M the constraints (6), (7), (9), and (10) Proposition 4. Let P = (S, s 0 , P, L, X) be a pIMC, Γ ⊆ A be a state label, S Γ = { s | L(s) = Γ }, and (X r

  and MC M 2 such that property P =1 ( X (Γ ∧ P =1 (♦Λ))) holds for M 2 , does not hold for M 1 , holds for I w.r.t. |= a I , and does not hold for I w.r.t. |= o |= a I I or M |= d I I and for all state labels Γ, there exist MCs M ≤ and M ≥ such that M ≤ |= o I I , M ≥ |= o I I , and

	I .
	best of our knowledge, there are no works addressing the same problem with the at-every-step
	semantics or showing that addressing the problem in the once-and-for-all and IMDP setting is
	sufficiently general. The following theorem fills this theoretical gap by proving that the three
	IMC semantics are equivalent w.r.t quantitative reachability. In other words, for all MCs M
	such that M

  There exists a MC M |= a pI P if and only if there exists a valuation v solution of the CSP C ∃r (P , Γ) such that for each state s ∈ S: v(λ s ) is equal to true if and only if P s M (♦Γ) = 0 and s is reachable from s 0 .Finally, we encode the equation system from[START_REF] Baier | Principles of model checking[END_REF] in a last CSP encoding that extends C ∃r .Let C ∃r (P , Γ) = (X r ∪ X , D r ∪ D ,C r ∪ C ) be such that (X r , D r ,C r ) = C ∃r (P , Γ), X contains one variable γ s with domain [0, 1] per state s in S, D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

  ⊆ R be a probability bound. Adding the constraint γ s 0 ∼ p to the previous C ∃r encoding allows one to determine if there exists a MC M |= a pI P such that P M (♦Γ) ∼ p.

	There exists a MC
	M |= a pI P if and only if there exists a valuation v solution of the CSP C ∃r (P , Γ) such that v(γ s )
	is equal to P s M (♦Γ) if s is reachable from the initial state s 0 in M , and is equal to 0 otherwise.
	Let p ∈ [0, 1]

Table 1

 1

					1.6. Prototype implementation and experiments
	Benchmarks	#nodes	#edges	#paramInBounds	#parameters
	NAND	K=1; N=2	104	147		82	4
	NAND	K=1; N=3	252	364		200	5
	NAND	K=1; N=5	930	1 371		726	7
	NAND	K=1; N=10	7 392	11 207		5 698	12
	Table 1.15 -Benchmarks composed of 4 pIMCs used for verifying quantitative properties
		Size of the Boolean Integer Real-number Boolean Linear Quadratic
	Encoding produced CSPs var.	var.	var.	constr. constr. constr.
	SotA	exponential	no	no	yes	yes	yes	no
	C ∃c	linear	yes	no	yes	yes	yes	no
	C ∃r	linear	yes	yes	yes	yes	yes	no
	C ∃r	linear	yes	yes	yes	yes	yes	yes

.16 -Characteristics of the four CSP encodings SotA, C ∃c , C ∃r , and C ∃r . 7 772 intervals.

Table 1 .

 1 21 -Comparison of solving times between qualitative and quantitative encodings.

	1.6. Prototype implementation and experiments

C ∃c encoding in SMT (time in sec) Figure 1.20 -Comparing solving time between SMT and MILP formulations pIMC C ∃c C ∃r C ∃r Benchmark #states #trans. #par. #var. #cstr. time #var. #cstr. time #var. #cstr. time

  The invariant I ( v) is a conjunction of predicates over the variables of the system specifying properties that must always hold. i is the name of the event, t = {t 1 . . . t n } represents the (optional) set of parameters of the event, G i ( t , v) is the (optional) guard of the event and S i ( t , v) is the action of the event. An event is enabled in a given valuation of the variables (also called a configuration) if and only if there exists a parameter valuation such that its guard G i ( t , v) is satisfied in this context. The action S i ( t , v) of an event may contain several assignments that are executed in parallel. An assignment can be expressed in one of the following forms:

	e i =
	any t where
	G i (t , v)

is an (optional) variant used for proving the convergence of the model, Evts is a set of events and Init ∈ Evts is the initialisation 62 2.2. Background event. then S i ( t , v) end Events An event has the following structure (see Figure on the side), where e -Deterministic assignment: x:= E( t, v) means that the expression E( t , v) is assigned to the variable x.

  The model P2P 1 given in Fig 2.2 presents an abstract Event-B specification of the protocol. It represents a general abstraction of the behaviour of the protocol with no details included. At this level, the state of the protocol is described by means of one variable DB. This variable contains a matrix which indicates for each client c∈1..N and each block b∈1..

K whether the client has already downloaded the block (DB(c →b)=finished) or not (DB(c →bc=empty). Initially, no block has been downloaded by any client.

  After reduction (removal of the intermediate states), we obtain the MC given Fig.2.9.

						2.5. Introducing probabilities in event-B models through refinement
	-the global probability of returning to (s 1 ) is 4 5 × 1 2 × 4 10 = 4 25 ,
	-the global probability of reaching (s 3 ) is 4 5 × 1 2 × 4 10 = 4 25 ,
	-the global probability of looping on (s 2 ) is 4 5 × ( 1 2 × 6 10	+ 1 2 × 6 10	) = 12 25 , where (2, 1) and
									(2,1)	(1,2)
	(1, 2) are the parameter valuations leading to these probabilistic choices.
								P2PP	1 4 |Start1DL
			1 5 |FailureDL(1,2)	1 4 |Start1DL	DB =	. . . .	, DBin =	. . . .	DB =	. . . .	, DBin =	. • . .
									1 4 |Start1DL
								1 4 |Start1DL
	DB =	. . . .	, DBin =	• . . .	DB =	. . . .	, DBin =	. • . .	DB =	. . . .	, DBin =	. . • .
	DB =	. . . .	, DBin =	. • . .	DB =	. . . .	, DBin =	. • . •	DB =	. . . .	, DBin =	. • • .

Table 3 .

 3 13 -Results of the experiments

		Model #p	10k		20k		50k
			V1	V2	V1	V2	V1	V2
	Running time M 1	28s		51-54s	142-143s
	Scenario 1	M 1 5	4.99%	5.09%	4.74%	5.10%	4.91%	4.98%
	Conf. interv.	M 1 5 ±0.85%	±0.82%	±0.55%	±0.56%	±0.36%	±0.37%
	Scenario 2	M 1 5 10.38%	10.04%	9.82%	10.05%	9.95%	9.81%
	Conf. interv.	M 1 5 ±1.15%	±1.12%	±0.79%	±0.80%	±0.51%	±0.51%
	Running time M 2	28s		53-54s	149-155s
	Scenario 1	M 2 5	5.44%	5.31%	5.61%	5.21%	5.59%	5.47%
	Conf. interv.	M 2 5 ±0.98%	±0.86%	±0.69%	±0.64%	±0.42%	±0.43%
	Scenario 2	M 2 5	10.8%	10.9%	10.8%	10.8%	10.9%	10.7%
	Conf. interv.	M 2 5 ±1.35%	±1.32%	±0.91%	±0.92%	±0.57%	±0.57%
	Running time M 3	185-190s	311-314s	612-621s
	Scenario 1	M 3 9	4.95%	5.97%	5.28%	6.62%	4.16%	5.61%
	Conf. interv.	M 3 9 ±5.22%	±5.71%	±4.71%	±6.25%	±1.86%	±4.38%
	Scenario 2	M 3 9	9.55%	9.87%	10.3%	11.3%	9.57%	10.7%
	Conf. interv.	M 3 9 ±8.40%	±7.86%	±7.04%	±7.89%	±5.29%	±3.99%

Table 4 .

 4 6 -Variations in learning time depending on the size of the sample and the complexity of the model

	80

Table 4 .

 4 9 -Results for quantitative statistical model checking of query ϕ 1 on M 1 with Monte-Carlo simulations

	δ	ε	N	Duration	p(ϕ 1 |M 1 )
	Test 1 0.02 0.01 6622	5 min 54 sec 0.9177 (> 0.80 property satisfied)
	Test 2 0.008 0.005 46808 52 min 18 sec 0.92179 (> 0.80 property satisfied)
	Test 3 0.005 0.003 130044 2 hr 8 min 0.92182 (> 0.80 property satisfied)

Table 4 .

 4 10 shows the results for this qualitative test.

	δ	α	β	m Duration P(ϕ 1 | M 1 ) > 80%
	Test 1 0.02 0.01 0.01 167 8.18 sec Property satisfied
	Test 2 0.008 0.005 0.005 275 16.21 sec Property satisfied
	Test 3 0.005 0.003 0.003 344 41.29 sec Property satisfied

Table 4 .

 4 10 -Different results for qualitative statistical model checking of query ϕ 1 on model M 1 with Sequential Probability Ratio Test.

Table 4 .

 4 An illustration of the experimental protocol that is conducted in this section erated from RTGEMs having different structures and parameters describing different behaviors on a prepaid card online service. The security query that we use in this experiment is: ϕ = 1000 (Transfer Money ⇒ Recharge (0,20] ∨ Check account (0,5] ), i.e. each time a transfer is performed, the user must have either recharged his/her account in the previous 20 time units, or checked his/her balance in the previous 5 time units. We want the model to satisfy the security property:P(ϕ | M) > 0.8.For the experiments, the statistical model checking techniques that we use are the qualitative Sequential Probability Ratio test and the quantitative estimation based on Monte-Carlo simulations. These techniques are used in parallel as previously explained. The SMC algorithms are calibrated as shown in Table4.12. Remark that the maximum number of simulations that we do (for each model) is determined by the number of simulations corresponding to the quantitative technique, in the case where SPRT does not converge faster to a result. 12 -Calibration of the SMC techniques for the following experiments First case: (OK). In the first test, we run our algorithm on the corresponding set of data to learn the model M o OK . The obtained model is shown in Figure

			Learn Model M o xxx
			from data
			Model M o xxx
			Check ϕ on M o xxx
	Stop	YES	Is ϕ satisfied?	SHD * (M * xxx , M o xxx )
			NO
			Model space exploration
			via Random Walk
		NO
			Model M * xxx
			found?
				YES
				Model M * xxx
	Figure 4.11 -

Table 4 .

 4 Table 4.14, we show the obtained results after the completion of the first test. The formal verification with SMC is repeated twenty times in order to establish mean values for the required 147 Chapter 4 -Learning and Verification of Graphical Event Models duration and number of simulations before acquiring a result. We write t for the mean duration and m for the mean number of simulations over twenty executions for all remaining tests. 14 -Results of the first test on M o OK The learned model M o OK , that best represents the data, verifies the security property in this case. Hence, no further exploration is made and the test stops. Second test: (KO 1 ). For the second test, we run our algorithm on the corresponding set of data to learn the model M o KO 1 . The obtained model is shown in Figure 4.15. The model M o KO 1

	P(ϕ | M) > 0.8	t	m
	M o OK Property satisfied 5 sec (±0.43 sec) 277 simulations (±49 simulations) M * OK ---

Table 4 .

 4 17 -Performances of the algorithm applied on the set of data KO 1The last results for this test are shown in Table4.18 after the completion of the experiment. The formal verification with SMC is repeated twenty times for both M o KO 1 and M * KO 1 in order to establish mean values for the required duration and number of simulations before converging on a result. Note that the SMC was executed only one time for all the intermediate models obtained at each step of the exploration process while constructing M * KO 1 . Third case: (KO 2 ). For the third and final test, we run our algorithm on the corresponding set of data to learn the model M o KO 2 . The obtained model is shown in Figure 4.19. The model M o KO 2does not satisfy the security property. Hence, the space exploration technique is repeated twenty separate times in order to check whether every time we find a model that satisfies the property. Again, we check whether it is the same proximal secure model that is found every time or if there are different possible ones.

	M * KO 1 which satisfies security property in the neighbor-hood of M o KO 1

Table 4 .

 4 22 -Performances of the algorithm applied on the set of data KO 2The last results for this test are shown in Table4.23. The formal verification with SMC is repeated twenty times for M o KO 2 and both M * number of simulations before acquiring a result. Note that also in this test, the SMC was executed only one time for all the intermediate models obtained at each step of the exploration process.

	τ	n Success Rate SHD * (M * KO 2 , M o KO 2 )

  P(s, e i , s ) = [s]W i ( v) ∑ e j ∈Acts(s) [s]W j ( v)

	× ∑ θ∈T s e i	x∈Var(e i ) s (θ) × ∏ P T e i	P e i s,θ (x, s )

  Modélisation et Vérification de Systèmes Incertains Mot clés : Modélisation, Vérification formelle, Systèmes probabilistes, Systèmes paramétrés Résumé : Le but des sciences en général est d'étudier le fonctionnement de systèmes complexes afin de mieux les comprendre et, éventuellement, de pouvoir prédire leur comportement futur. Une part importante de cette étude est la réalisation d'un modèle, objet abstrait qui représente fidèlement les connaissances existantes à propos du système, et qui pourra alors être étudié et simulé à la place du système lui-même. Malheureusement, il est fréquent que les connaissances existantes à propos des systèmes soient incomplètes ou sujettes à incertitudes. Il est alors important d'inclure ces incertitudes à l'intérieur des modèles et de développer des techniques afin d'automatiser l'analyse de ces modèles.Ce document présente quatre contributions au domaine de la modélisation et de l'analyse de tels systèmes. Deux de ces contributions sont à vocation plutôt théorique et proposent des langages de modélisation permettant d'inclure les incertitudes à l'intérieur des modèles ainsi que des techniques de vérification associées. Les deux autres contributions sont plus pratiques, l'une proposant une technique de vérification statistique et prouvant son efficacité sur un cas d'étude industriel, et l'autre proposant une technique de construction automatique de modèle à partir d'un ensemble de données ainsi qu'une analyse de son niveau de sécurité.

Titre : Title: Modeling and Verification of Systems with Uncertainties Keywords: Modeling, Formal verification, Probabilistic systems, Parametric Systems

The list is not exhaustive, the survey of all existing modeling formalisms is beyond the scope of this document.

https://anr.fr/Projet-ANR-14-CE28-0002

Indeed, when 0 ≤ v( f 1 ) ≤ v( f 2 ) ≤ 1 is not respected, the interval is inconsistent and therefore empty.

1.5. Quantitative properties

All resources, benchmarks, and source code are available online as a Python library at https://github. com/anicet-bart/pimc_pylib

http://www.prismmodelchecker.org/

See the category discrete-time Markov chains on the PRISM website

Note that obtaining integer integrity constraints over real-numbers can be costly.

Event-B[START_REF] Abrial | Modeling in Event-B: system and software engineering[END_REF] is a proof-based formal method used for discrete systems modeling. It is equipped with Rodin [ABH + 10], an open toolset for modeling and proving systems. This toolset can easily be extended, which makes of Event-B a good candidate for introducing probabilistic reasoning in a proof-based modeling formalism. The development process in Event-B is based on stepwise refinement: systems are typically developed progressively using an ordered sequence of models, where each model contains more details than its predecessor. Refinement allows a step-by-step description of the behaviour of systems, providing an efficient way to give a detailed description of their behaviour. So far, several research works have focused on the extension of Event-B to allow the expression of probabilistic information in Event-B models. In[START_REF] Morgan | The challenge of probabilistic event b-extended abstract[END_REF], Abrial et al. have summarised the difficulties of embedding probabilities into Event-B. This paper suggests that probabilities need to be introduced as a refinement of non-determinism. In Event-B, non-determinism occurs in several places such as the choice between enabled events in a given state, the choice of the parameter values in a given event, and the choice of the value given to a variable through some non-deterministic assignments. To the best of our knowledge, the existing works on extending Event-B with probabilities have mostly focused on refining non-deterministic assignments into probabilistic ones while other sources of non-determinism have been left untouched. In [HH07], Hallerstede et al. focus on a qualitative aspect of probability. They refine non-deterministic assignments into qualitative probabilistic assignments where the actual probability values are not specified, and adapt the Event-B semantics and proof obligations to this new setting. In[START_REF] Yilmaz | Tool support for qualitative reasoning in Event-B[END_REF], the same authors study the refinement of qualitative probabilistic Event-B models and propose a tool support inside Rodin. Other works [TTL09, TTL15, TTL10] have extended this approach by refining non-deterministic assignments into quantitative probabilistic assignments where, unlike in[START_REF] Hallerstede | Qualitative probabilistic modelling in Event-B[END_REF], the actual probability values are specified. This new proposition is then exploited in order to assess several system properties such as reliability and responsiveness.

Some aspects such as the probabilisation of standard Event-B models or the study of mixed Event-B models are deliberately left out in order to keep this document to a manageable size. Nevertheless, these topics have also been studied and are detailed in[START_REF] Amine Aouadhi | Introducing probabilistic reasoning within Event[END_REF].

In this version, we only consider constant values for the probabilities in enumerated probabilistic assignments, but this could be extended to expressions with a few adaptations.

(event/param/pWD)

sign/pWD3) for predicate probabilistic assignments and their non-emptyness is ensured by the standard feasibility POs.

This line of work has been dormant since the end of Mohamed Amine Aouadhi's thesis. In particular, our implementation of the probabilistic Even-B plugin for Rodin has been left at a

https://www.pixiel-group.com/

See https://depend.cs.uni-saarland.de/tools/param/casestudies/

See https://depend.cs.uni-saarland.de/tools/param/casestudies/ for details

http://www.gfi.world

In the following, ϕ is called a query and the term security property refers to P(ϕ | M) > c.

http://www.pilgrim.univ-nantes.fr/

Details on the random generation process are left out of this document but can be found in[START_REF] Antakly | [END_REF].

of the best (or worst) probability value, they do not provide a complete analysis of the effect of non-deterministic choice (or parametric transitions) on the satisfaction of the given property. This is something which is easily provided using our technique. Section 3.2 presents our extension of statistical model checking to parametric systems. It starts with a reminder and an adaptation of some definitions from previous chapters, and an introduction of basic definitions concerning statistical model checking. Our main theoretical contribution, the adaptation of Monte-Carlo simulation to the parametric setting, is then presented in Section 3.3. In Section 3.4, we report on our prototype implementation and its use on automata-based models, and provide some comparison to other parametric model checking tools. Our main case-study is then introduced in Section 3.5, as well as its resolution using the version of our prototype tool for Python models. In Section 3.6, we give some leads for improvement of the techniques implemented in our prototype tool. Finally, Section 3.7 concludes the chapter and Section 3.8 presents perspectives for future work.

Background

We start by recalling some basic definitions, most of them already given in Chapter 1.

Basic definitions

As usual, the set of real numbers and the set of natural numbers are respectively written R and N. In this chapter, we consider Markov Chains with no labels on states (although all our results could be extended to this setting in a straightforward manner). Here, a Markov chain is a tuple M = (S, s 0 , P) where S is a denumerable set of states, s 0 ∈ S is the initial state and P : S × S → [0, 1] is the transition probability function such that for all state s ∈ S, ∑ s ∈S P(s, s ) = 1.

Recall from Chapter 1 that a run of a Markov chain is a sequence of states s 0 , s 1 . . . such that for all i, P(s i , s i+1 ) > 0. Given a finite run ω = s 0 , s 1 , . . . s n , its length, written |ω| represents the number of transitions it goes through (including repetitions). Here |ω| = n. We write Γ M (l) (or simply Γ(l) when M is clear from the context) for the set of all finite runs of length l, and Γ M for the set of all finite runs i.e., Γ M = ∪ l∈N Γ M (l). As in Chapter 1, we use the probability measure on runs based on the sigma-algebra of cylinders (see e.g. [START_REF] Baier | Principles of model checking[END_REF]). This probability CHAPTER 4

LEARNING AND VERIFICATION OF

GRAPHICAL EVENT MODELS

The work we report in this chapter was achieved during Dimitri Antakly's thesis, in collaboration with GFI informatique group 1 . GFI is an international group that mainly proposes computer science solutions and expertise, touching a wide variety of domains (banking, road safety, supervision, cloud, AI solutions and many more). The main expectations of GFI about this thesis were to explore behavior analytics including how it can be applied in security or supervision, and to use innovative techniques in order to be distinguished from the available market solutions.

In this chapter, we tackle uncertain systems in a different form and context than those previously considered. Instead of analyzing a given model of the system of interest, we here assume that the only information we have on the system is a precomputed set of data and a security property. We therefore focus on a formalism dedicated to learning, along with its learning algorithm, and develop techniques for verifying the learned model and analyzing its security w.r.t. the given property.

Using these notations, we define the Structural Hamming Distance (SHD) between G 1 and G 2 as follows:

Each edge that is present in one RTGEM and not in the other therefore adds 1 to the SHD, while each edge that is present in both adds the elementary distance between its respective timescales. One can show that the SHD is indeed a distance metric (i.e., satisfies the distance axioms). A proof is given in [START_REF] Antakly | [END_REF].

While this distance allows to compare both the structure and the timescales of two RTGEMs using the same labels, it suffers from a "scaling" disadvantage when it comes to comparing quantitative information. For instance, consider three vectors v e 1 = [0, 1, 2], v e 2 = [0, 1.01, 1.98] and v e 3 = [0, 5, 10], representing the endpoints of three different timescales to compare using the distance measure. By computing the elementary distances between (v e 1 , v e 2 ) and (v e 1 , v e 3 ) with respect to equation 4.3, we notice that they are equal although v e 1 and v e 2 cover almost the same timescale compared to v e 3 that covers a different timescale. We have therefore proposed an extension to this SHD where the relative quantitative differences inside the timescales are taken into account in order to give a fairer distance measure.

The idea is to find matches (if existing) between the endpoints of the two timescales based on the mutual minimal absolute difference. In other words, consider two vectors v e 1 with size l and v e 2 with size k. A match is a pair (v e 1 i , v e 2 j ) such that the closest element from v e 1 i ∈ v e 1 in v e 2 is v e 2 j , and the closest element from v e 2 j ∈ v e 2 in v e 1 is the same v e 1 i . Formally, the closest element from v e

Using this notion, the set of matches for a given edge e ∈ E 1 ∩ E 2 is written V e id and defined as follows:

V

We write V e id,1 for the projection of V e id on v e 1 , i.e., V e id,1 = {v e 1,i , ∃ j, (v e 1,i , v e 2, j ) ∈ V e id }. Similarly, we write V e id,2 for the projection of V e id on v e 2 . The set of unmatched endpoints can now be defined as:

V e nid = (v e 1 \V e id,1 ) ∪ (v e 2 \V e id,2 ).

Using these sets, we enlarge the set of endpoints that are considered identical in the elementary distance between two timescales. We now write v e id = |V e id | for the number of proximally matched

Testing the proposed strategy

We now propose our main experiment, that applies our whole strategy (Algorithm 1 to three different sets of synthetic data. We start with a description of the experimental context and then present the results for each attempt.

Experimental protocol

The aim of this experiment is to perform a proof-of-concept for our proposed strategy (given in Algorithm 1). We will therefore consider three cases: (1) when the learned model satisfies the given security property (case OK), (2) when the learned model does not satisfy the security property, but is close to satisfying it (case KO 1 ), and (3) when the learned model is far from satisfying the security property (case KO 2 ).

Our pipeline is shown in Figure 4.11, where xxx ∈ {OK, KO 1 , KO 2 } and one can observe that we follow the strategy proposed in Algorithm 1. In the OK case, we use data generated from an RTGEM that satisfies the security property, while in KO 1 and KO 2 , we use data generated from models that do not satisfy the security property. The details of the models will be given further on.

Remark. A part of our strategy relies on the exploration of the neighborhood of the fittest model M o . For this experiment, we have chosen a notion of neighborhood based on the number of forward operators applied on an RTGEM. Moreover, we have chosen to perform a Random Walk on this neighborhood. In other words, we apply a limited number of random forward operators on the fittest model M o in order to produce candidate models M * . The maximum number of forward operators varies for each case (irrelevant for (OK), 25 for (KO 1 ) and 100 for (KO 2 )) and have been set according to preliminary experiments.

In order to optimize the exploration procedure, we do not choose randomly from the whole set of possible operators, but rather from a subset of those operators that might have a more significant impact on the satisfaction of the property. For the sake of completeness, if no model is found when the whole subset of operators has been sampled, we extend our search to the remaining ones. This strategy is discussed in further details in [START_REF] Antakly | [END_REF].

Experimental results for the proposed strategy

The synthetic data sets used for learning the models contain event streams with labels "Login", "Logout", "Check account", "Recharge" and "Transfer money". These sets were gen-

RÉSUMÉ EN FRANÇAIS

MAIN DEFINITIONS AND NOTATIONS

This appendix gathers the main notations and definitions introduced in the chapters of this document. It is aimed to be detached from the manuscript so the reader can use it while reading the main part of the document.

Due to the different topics treated in each chapter, remark that the definitions of some objects (Markov chains for instance) may slightly differ from one chapter to another.

Chapter 1: Parametric Interval Markov Chains

Variables, numbers and associated notions

} is a set of variables and D x , D X are their associated domains/set of domains.

• A rational function f over X is a division of two (multivariate) polynomials g 1 and g 2 over X with rational coefficients, i.e., f = g 1 /g 2 .

• Q is the set of rational numbers and Q X is the set of rational functions over X.

Valuations are extended to polynomials and rational functions.

• If S is a finite (or discrete) set, we write Dist(S) for the set of probability distributions over S, i.e., the set of functions µ : S → [0, 1] such that ∑ s∈S µ(s) = 1.

• I is the set containing all interval subsets of [0, 1].

• A is a set of atomic propositions (state labels), written in Greek alphabet.

Constraints

• An atomic constraint over X is a Boolean expression of the form f (X) g(X), with ∈ { ≤, ≥, <, >, = } and f and g two functions over variables in X. A constraint over X is a Boolean combination of atomic constraints over X.

Definition 1 (Constraint Satisfaction Problem -see p.27). A Constraint Satisfaction

Problem (CSP) is a tuple Ω = (X, D,C) where X is a finite set of variables, D = D X is the set of all domains associated to the variables from X, and C is a set of constraints over X.

Definition 10 (Parametric Interval Markov Chain -see p.35). A Parametric Interval

Markov Chain (pIMC for short) is a tuple P = (S, s 0 , P, L, X), where S, s 0 , L, and X are defined as for pMCs, and P : S × S → I(Q X ) associates with each potential transition a (parametric) interval.

• pIMC is the set containing all parametric interval Markov chains.

• Given a pIMC P = (S, s 0 , P, L, X) and a valuation v, we write v(P ) for the IMC (S, s 0 , P v , L) obtained by replacing the transition function P from P with the function P v : S × S → I defined by P v (s, s ) = v(P(s, s )) for all s, s ∈ S.

• The IMC v(P ) is called an instance of pIMC P . • The size of a pMC, IMC, or pIMC L, written |L|, corresponds to its number of states plus its number of transitions not reduced to 0, [0, 0] or ∅.

• Pred(s) = { s ∈ S | P(s , s) / ∈ { ∅, [0, 0] } } are the predecessors of a given state s (extends to sets of states).

• Succ(s) = { s ∈ S | P(s, s ) / ∈ { ∅, [0, 0] } } are the successors of a given state s (extends to sets of states).

Model comparisons

Let (L 1 , |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models containing respectively L 1 and L 2 . Chapter 2: A Probabilistic Extension for Event-B Transition Systems

• Dist(S) denotes the set of distributions over a given set S.

• A labelled transition system (LTS for short) is a tuple M =(S, Acts, s 0 , →, AP, L)

where S is a set of states, s 0 ⊆ S is the initial state, Acts is a set of actions, → ⊆ S × Acts × S is a transition relation AP is a set of atomic propositions, and L : S → 2 AP is a labelling function.

• A Probabilistic Labelled Transition System (PLTS for short) is a tuple M =(S, s 0 , Acts, P, AP, L) where S is a set of states, s 0 ∈ S is the initial state, Acts is a set of actions, AP is a set of atomic propositions, L: S → AP is a labelling function, and P: S × Acts × S → [0, 1] is the transition probability function.

• A PLTS where, for each state s ∈ S we have ∑ s ∈S,a∈Acts P(s, a, s ) = 1 is a Discrete Time Markov Chain.

Event-B

• An Event-B model is a tuple M=( v,I( v ), V( v),Evts, Init ) where v= {v 1 . . . v n } is a set of variables, I ( v) is an invariant, V( v) is an (optional) variant used for proving the convergence of the model, Evts is a set of events and Init ∈ Evts is the initialisation event.

• An event is written e i , with a set of parameters t = {t 1 . . . t n } (optional), a guard G i ( t , v) and an action S i ( t , v).

• Assignment can be written

• The action S j ( t , v) of a given event e j may contain several assignments that are executed in parallel.

Probabilistic Event-B

• Probabilistic events are equipped with weights W i ( v).

• Predicate non-deterministic assignments are replaced with predicate probabilistic assignments written x:⊕Q x ( t, v,x').

• Enumerated non-deterministic assignments are replaced with enumerated probabilistic assignments written • A Markov chain is a tuple M = (S, s 0 , P) where S is a denumerable set of states, s 0 ∈ S is the initial state and P : S × S → [0, 1] is the transition probability function such that for all state s ∈ S, ∑ s ∈S P(s, s ) = 1.

• A run of a Markov chain is a sequence of states s 0 , s 1 . . . such that for all i, P(s i , s i+1 ) > 0.

• The length of a run, written |ω| represents the number of transitions it goes through (including repetitions).

• Γ M (l) is the set of all finite runs of length l in M , and Γ M for the set of all finite runs in M .

• The probability measure on the finite runs of M is written P M .

• The probability of a finite run ω = s 0 , s 1 , . . . s n is P M (ω) = Π n i=1 P(s i-1 , s i ). • Notation r is used for reward functions r : Γ(l) → R.

• We use the notation ϕ for properties on runs (ϕ ⊆ Γ(l)).

• Reachability properties are written P M (♦ ≤l s).

• Safety properties are written Safety P M ( =l E).

• Expected reward properties are written E l M (r) = ∑ ω∈Γ(l) P M (ω)r(ω).

Parametric Markov chains

• X is a set of parameters (variables).

• Poly(X) is the set of all real (multivariate) polynomials on X.

• v ∈ R X stands for parameter valuations, and can be applied to polynomials by extension.

• A Parametric Markov chain is a tuple M = (S, s 0 , P, X) such that S is a finite set of states, s 0 ∈ S is the initial state, X is a finite set of parameters, and P : S×S → Poly(X) is a parametric transition probability function.

¢

Cut here

¢

Chapter 4: Learning and Verification of Graphical Event

Models

Statistical model checking

• We write S for a generic stochastic system.

• We use the notation ϕ for linear properties on the traces of a system.

• The random variable corresponding to the outcome of each simulation is written B i .

• The expected value of a given variable X is written E[X].

• The probability that this random variable is 1 (i.e., that S satisfies the given property

• In qualitative SMC, the parameters are θ (threshold), α (Type-I error), β (Type-II error), and [p 0 ; p 1 ], the indifference region. Usually, p 0 = pδ and p 1 = p + δ, and 2δ is the size of the symmetrical indifference region.

• In hypothesis testing, the hypothesis are H 0 : p ≤ θ and H 1 : p > θ.

• In quantitative SMC, the parameters are δ (precision) and ε (error rate).

• In all cases, the estimated value of an unknown variable p is written p.

Conditional intensity models

• We use the notation t ∈ R + for timestamps and l ∈ L for labels.

• Given a totally ordered set S and a set of marks M, a marked point process (m.p.p) can be expressed as a sequence {(s i , m i ) : i = 1, ..., n} where for all i, s i ≤ s i+1 .

• An event is a pair (t, l) consisting of a timestamp and a label.

• An event stream is an m.p.p. where S = {t i } is a set of strictly increasing timestamps and M = {l i } is a set of event labels. It is written

• t 0 = 0 and t * = t n+1 are used as conventions for a stream consisting of n events

• The length of an event stream is written |x t * | and corresponds to the number of events (n in the above item).

• The history of the event streal at time t is the set of all the events that occurred before t: h i denotes the ith history, with h i = h t i = (t 1 , l 1 ), ..., (t i-1 , l i-1 ).

• When clear from the context, we will use x to denote data (event streams), h for a given history, and t for a givem timestamp.

Definition 14 (Conditional intensity function -see p.127). Given an event stream x, a conditional intensity function λ is written as:

where N(T ) denotes the number of points (events) occurring in a time interval T in x.

Definition 15 (Conditional intensity model [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] -see p.128). Formally, a conditional intensity model (CIM) θ is a set of indexed conditional intensity functions {λ l (t | x)} l∈L . The data likelihood function is the joint density function of all the points in the m.p.p. that can be factorized into a product of all conditional intensity functions and can be written as:

where Λ l (t | h; θ) = t -∞ λ l (τ | x; θ)dτ for the data x and the indicator function 1 l (l ) is one if l = l and zero otherwise.

• In Piecewise-Constant Conditional Intensity Models (PCIMs), intensity functions are constant in given time intervals (finite number of states that depend in the time t and the observed data x).

• For each label l ∈ L, we consider a set of discrete states Σ l .

• For each label l and each state s ∈ Σ l we consider a parameter λ ls .

• We consider a mapping σ l : T × X -→ Σ l (with T the set of all possible times and X the set of all possible data) that determines the active state as a function of the chosen label and available data.

Definition 16 (Piecewise-constant conditional intensity models -see p.129). A PCIM is defined by local structures S l = (Σ l , σ l (t, x)) and local piecewise constant parameters λ ls . Let S = {S l } l∈L be the set of all discrete states and θ = {λ ls } l∈L,s∈Σ l be the set of associated parameters. The PCIM data likelihood knowing the structure and the parameters in this case can be written as:

where M ls (x) is the number of occurrence of events of type l while s is active in the event sequence x, and T ls (x) is the total duration while s was active for event type l.

¢

Cut here

¢

Graphical event models

Definition 17 (Graphical event model -see p.129). A Graphical Event Model (GEM) is a tupple G = ((L, E), θ) where (L, E) is a directed graph and θ a set of parameters. The likelihood of the data g knowing the graph and its parameters is written as:

with h |x t * |+1 the entire history of the event stream, including x t * .

• The conditional intensity functions λ l (t | h) in GEMs are not piecewise-constant.

• The Markov property in the context of m.p.p.s. implies the following property on conditional intensity functions: For all l, h and t, λ l (t | h) = λ l (t | [h] Pa(l) ), where Pa(l) are the parents of l in G.

• A timescale is a set T of half-open intervals (a, b] (with a ≥ 0 and b > a) that form a partition of some interval (c,t h ] (with c ≥ 0), where t h is the highest value of T .

• When an edge e is equiped with a timescale T , t h is called the horizon of e.

Definition 18 (Timescale graphical event model -see p.131). A TGEM M = (G, T ) consists of a GEM G = ((L, E), θ) and a set of timescales T = T e (e∈E) corresponding to the edges E of the graph of G.

• c l (h,t) is the parent count vector of bounded counts (of occurrences) over the intervals in the timescales of the parents of l at time t.

• C l is the set of all possible parent count vectors of label l.

• A given parent count vector yields a unique conditional intensity function for the corresponding parameter. The global conditional intensity functions are therefore piecewise-constant:

• In the case of TGEMs, the likelihood of the data is written as follows:

where the sufficient statistics n t * ,l, j (x t * ) and d t * l, j (x t * ) are the count of l-events and the durations, respectively, when the parent count vector was equal to j (a certain combination of parents).

¢

Cut here

¢

• Given two RTGEMs with the same set of labels L, G 1 = ((L, E 1 ), θ 1 ) and G 2 = ((L, E 2 ), θ 2 ) and an edge e ∈ E 1 ∩ E 2 , we write v e 1 and v e 2 for the values of the respective timescales of e in G 1 and G 2 .

• Consider two vectors v e 1 with size l and v e 2 with size k. The closest element from v e 1 i ∈ v e 1 in v e 2 is written cl(v e 1 i , v e 2 ) = argmin v e 2 p

(|v e 1 iv e 2 p |). • The set of matches for a given edge e ∈ E 1 ∩ E 2 is written V e id and defined as follows:

V e id = {(v e 1 i , v e 2 j ) ∈ v e 1 × v e 2 : cl(v e 1 i , v e 2 ) = v e 2 j ∧ cl(v e 2 j , v e 1 ) = v e 1 i }

• We write V e id,1 for the projection of V e id on v e 1 , i.e., V e id,1 = {v e 1,i , ∃ j, (v e 1,i , v e 2, j ) ∈ V e id }. • Similarly, we write V e id,2 for the projection of V e id on v e 2 . • The set of unmatched endpoints is V e nid = (v e 1 \V e id,1 ) ∪ (v e 2 \V e id,2 ). • The proximal distance measure between G 1 and G 2 is defined as follows:

Experimentations

• The qualitative security property is of the form P(ϕ | M) > c where ϕ is a linear property (called a query) and c is a constant.

• Our algorithm starts with data D, learns a first model M o , then (if M o does not satisfy the qualitative security property) a proximal model M * satisfying the qualitative security property.

• The empty model is written ⊥.