# Contributions to formal modelling and analysis of stochastic models
## Mémoire d'habilitation à diriger des recherches

Paolo Ballarini

► To cite this version:

Paolo Ballarini. Contributions to formal modelling and analysis of stochastic models Mémoire d'habilitation à diriger des recherches. Computer Science [cs]. CentraleSupélec, Université Paris-Saclay, 2023. tel-04335487

HAL Id: tel-04335487

https://hal.science/tel-04335487

Submitted on 11 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Contributions to formal modelling and analysis of stochastic models

Mémoire d'habilitation à diriger des recherches

Paolo Ballarini

Defended on June 6th 2023, before the committee composed of:

Eugenio Cinquemani, INRIA Grenoble Alpes reviewer

Alain Denise, Université Paris Saclay president

Giuliana Franceschinis, Università Piemonte Orientale examiner

Stefan Haar, École Normale Supérieure Saclay examiner

Dave Parker, University of Oxford reviewer

Verena Wolf, Universität des Saarlandes reviewer
Contributions to formal modelling and analysis of stochastic
models

Mémoire d'habilitation à diriger des recherches

Paolo Ballarini

Abstract

This manuscript describes a number of contributions, I had the opportunity to work on, in the field of formal methods for stochastic models. Research in formal methods is concerned with defining mathematically rigorous formalisms that allows for specification, development and verification of the correctness of complex systems. In the case of stochastic systems modelling formalisms need to account for stochasticity and verification algorithms needs to address also *quantitative* properties other than *qualitative* ones.

The first contribution included in the manuscrit overviews a statistical model checking framework, named Hybrid Automata Specification Language (HASL), which targets a general class of stochastic processes and that, by employing hybrid automata as *machineries* to specify the behavior to be analysed, leads to a powerful property specification language which surpasses the expressiveness of classical stochastic temporal logic languages such as CSL.

The second contribution is concerned with the *parametrised* extension of the stochastic model checking problem, namely the realisation of a procedure to estimate the (probabilistic) dependency of the satisfaction of a property φ w.r.t. the parameters of a stochastic model. By introducing a *satisfiability distance* measure that quantifies how distant a property φ is from being satisfied by a model’s instance and relying on the expressive of HASL property language, we adapted classical Approximated Bayesian Computation (ABC) schemes so to assess the *satisfaction probabilistic distribution* of φ against the parameter space Θ of the considered model.

The third contribution is about the application of model checking approaches to the analysis of stochastic oscillators, i.e. a prominent class of systems in systems biology. Classical temporal logic languages have proven limited success in detecting/assessing oscillations, whereas effective procedures can be obtained through automata-based model checking. The contribution we present is twofold: on one hand we show how sophisticated oscillation-related indicators can be assessed through HASL based model checking and on the other hand we show how to assess the dependency of such indicators (period and amplitude) against the parameters of a stochastic oscillator.

As a final contribution we illustrate a framework for building stochastic models that adequately reproduce periodic phenomenon from a dataset, with data about visits at a hospital Emergency Department (ED) as the reference case study. We considered different
families of stochastic processes, Markov renewal processes, Markov arrival processes (MAP)
and hidden Markov models (HMM) and through fit them so that they match meaningful
statistical indicators computed on the dataset.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Expressive statistical model checking</td><td>5</td></tr><tr><td>2.1</td><td>Preliminaries</td><td>10</td></tr><tr><td>2.1.1</td><td>Discrete event stochastic processes</td><td>10</td></tr><tr><td>2.1.2</td><td>DESP as Non Markovian Generalised Stochastic Petri Nets</td><td>11</td></tr><tr><td>2.2</td><td>Hybrid Automata Stochastic Language</td><td>13</td></tr><tr><td>2.2.1</td><td>Synchronised Linear Hybrid Automata</td><td>14</td></tr><tr><td>2.2.2</td><td>HASL expressions</td><td>25</td></tr><tr><td>2.3</td><td>Expressiveness of HASL</td><td>27</td></tr><tr><td>2.3.1</td><td>Comparison with logics for numerical stochastic model checking</td><td>27</td></tr><tr><td>2.3.2</td><td>A couple of meaningful applications</td><td>31</td></tr><tr><td>2.4</td><td>Software support: the COSMOS tool</td><td>36</td></tr><tr><td>2.4.1</td><td>A model driven code generation architecture</td><td>36</td></tr><tr><td>2.4.2</td><td>COSMOS statistical engines</td><td>37</td></tr><tr><td>2.5</td><td>Perspectives</td><td>39</td></tr><tr><td>3</td><td>A Bayesian approach to parametric model checking</td><td>41</td></tr><tr><td>3.1</td><td>Preliminaries</td><td>43</td></tr><tr><td>3.1.1</td><td>Markov Population Models</td><td>43</td></tr><tr><td>3.1.2</td><td>Linear-time temporal properties</td><td>47</td></tr><tr><td>3.1.3</td><td>Bayesian inference: the ABC method</td><td>49</td></tr><tr><td>3.1.4</td><td>Kernel density estimation</td><td>52</td></tr><tr><td>3.2</td><td>Automaton-ABC</td><td>53</td></tr><tr><td>3.2.1</td><td>Satisfiability distances for reachability problems</td><td>53</td></tr><tr><td>3.2.2</td><td>HASL specifications for satisfiability distances</td><td>57</td></tr><tr><td>3.2.3</td><td>Automaton-ABC: ABC with satisfiability distance</td><td>60</td></tr><tr><td>3.2.4</td><td>Estimation of the satisfaction probability function</td><td>63</td></tr><tr><td>3.3</td><td>Experiments</td><td>63</td></tr><tr><td>3.3.1</td><td>Enzymatic reaction system</td><td>63</td></tr></table>
<table><tbody><tr><td>3.4</td><td>Perspectives</td><td>67</td></tr><tr><td>4</td><td>Formal analysis of stochastic oscillators</td><td>69</td></tr><tr><td>4.1</td><td>Introduction</td><td>69</td></tr><tr><td>4.2</td><td>Formal analysis of stochastic oscillators</td><td>70</td></tr><tr><td>4.2.1</td><td>Temporal logic characterisation of oscillations</td><td>71</td></tr><tr><td>4.2.2</td><td>Partition based noisy periodicity</td><td>74</td></tr><tr><td>4.2.3</td><td>Peaks-detection based noisy periodicity</td><td>78</td></tr><tr><td>4.2.4</td><td>Case study: circadian clock</td><td>84</td></tr><tr><td>4.3</td><td>Tuning stochastic oscillators</td><td>89</td></tr><tr><td>4.3.1</td><td>An automaton for the distance from a target period</td><td>90</td></tr><tr><td>4.3.2</td><td>Case studies</td><td>90</td></tr><tr><td>4.4</td><td>Perspectives</td><td>94</td></tr><tr><td>5</td><td>Data-driven predictive modeling of periodic phenomena</td><td>97</td></tr><tr><td>5.1</td><td>Statistical analysis of an emergency department dataset</td><td>98</td></tr><tr><td>5.2</td><td>Deriving adequate models of patient arrivals</td><td>101</td></tr><tr><td>5.2.1</td><td>Markovian renewal process approximation</td><td>102</td></tr><tr><td>5.2.2</td><td>Markovian arrival processes</td><td>104</td></tr><tr><td>5.2.3</td><td>Hidden Markov processes</td><td>106</td></tr><tr><td>5.3</td><td>Empirical validation of derived models</td><td>108</td></tr><tr><td>5.3.1</td><td>Petri Net encoding of HMM arrival models</td><td>109</td></tr><tr><td>5.3.2</td><td>KPI assessment</td><td>113</td></tr><tr><td>5.4</td><td>Perspectives</td><td>116</td></tr><tr><td>6</td><td>Other contributions</td><td>117</td></tr><tr><td>7</td><td>Conclusion and perspectives</td><td>124</td></tr></tbody></table>
# Chapter 1

# Introduction

Modern engineering is a wide, pervasive and multi-disciplinary field which revolves around the notion of *system*, where a system may be thought of as aggregation of interacting parts (i.e. objects, things, components) which are combined by nature or humans to form a *whole* that performs a certain *function*. Examples of systems can be found in just about any possible field. They may be distinguished according with their kind, as in, for example, *natural systems* (e.g. the system of biochemicals reactions which controls the Eukaryotic cell's reproduction cycle), or *human-made systems*, (e.g. embedded systems consisting of interacting electronic devices that allow an aircraft to be safely flown, or smart devices, as in modern *Internet of Things* context allowing remote control for home utilities), or *abstract systems* (e.g. an ecosystem of beings interacting with each other for preserving their species). The development of model-based methodologies to support the design, validation, performance analysis, testing and maintenance is fundamental in system engineering. A model is a mathematical, abstract representation of a system, that, depending on its nature, may serve diverse purposes the main of which is to reliably mimic the system's dynamics, while possibly allowing for performance analysis, and/or verification of requirements. As our society relies more and more on complex software controlled critical systems the need to employ formal methods to design, validate and certify them becomes paramount.

## Model-based methods

My work is mainly focused on model-based methods, namely, the area of research whose goal is the development of formal modelling frameworks which allows one to unambiguously build a model of a real-life system and, based on rigorous mathematical foundations, assess the model's properties. The overall goal is to establish whether the modelled system functions correctly, which, depending on the nature of the model, may entail the verification of some *qualitative* requirement and/or the evaluation of some *quantitative* performance indicator.
**Model-checking.** Amongst model-based methods *model checking*, quickly became a prominent one as it gives the modeller the possibility to automatically verify a system's requirements, formally expressed as temporal logic properties, against a system's model. The seminal work of Pnueli [Pnu77] on the verification of sequential and concurrent programs set the ground for the formalisation of algorithms for model-checking of non-probabilistic systems introduced by Emerson and Clarke [EC80] and also, independently, by Quelles and Sifakis [QS82] (resulting in the 2007 Turing award awarded to the authors) and eventually leading to many extensions targeting different classes of models including probabilistic ones [ASSB96, ASSB00, BHHK03a]. The attractiveness of model checking is mainly due to its practical simplicity as, while keeping the conception of the model separated from that of the requirements (i.e. *separation of concerns*), it provides the user with a *push-button* machinery for certifying that the system is correct (i.e. if the model *satisfies* the requirements).

**Classes of models.** The variety of real-life systems calls for models to account for different features including timing constraints, probabilities and unknown parameters.

Timing constraints are crucial for safety-critical systems, examples of which can be found in domains such as aircraft controlling, nuclear power plants and the automotive industry, to mention a few. Timed automata [AD94] and timed Petri, for instance, are formalisms suitable to account for time dependent behavior and where timing requirements are expressed through dense-time variables that evolve all at the same speed and that can be used to capture either the delay since the last occurrence of an event or how long a given transition has been enabled for.

Other kind of systems instead are characterised by a *probabilistic behaviour*, that is, a behaviour which can not be *certainly* predicted in advance. We may distinguish between systems that are *probabilistic by design* as, for example, WiFi protocols which employ probabilistic algorithms to reduce *packet collisions*, as opposed to systems that are *probabilistic by nature*, as the dynamics of the real-life phenomenons they depend upon can only be expressed through some probability law, like for example for failure rates of physical components.

Finally most often models depend on a number of parameters whose value is not necessarily known a priori. For example, for a model of queuing system the maximum number of clients standing in a queue for a given service might not necessarily be known, or, for a model of a manufacturing system the probability that a given workstation breaks down in the next time unit may not be known. Whenever some parameters are unknown the model becomes *parametric* and any analysis technique also needs to be adapted to account for the effect that parameters have on the behavior. In model checking terms this entails adapting verification algorithms so to assess the effect that parameters have on the satisfaction of the model's requirements.

**Trading off expressiveness with tractability.** The realm of model-based techniques is
affected by two opposing forces i) the *expressiveness* of the modelling framework and of the
corresponding analysis techniques and ii) the tractability of the models. These two aspects
are *in opposition* as, generally speaking, the more expressive the modelling framework, the
more complex and computationally costly are the mathematical tools necessary to analyse
its behaviour.

In model checking terms, for example, such an antagonism entails considering to what extent the relevant characteristic of a system’s behaviour can be accounted for by the mod-eller through a formally expressed property. This depends both on the kind of model (e.g. transition-system, timed-automaton, Markov-chain) and of the mathematical means on top of which an algorithm for verifying the formulae of the corresponding property language can be defined.

Contents of this thesis

Over the last decade I have been focusing on research in the area of complex systems taking
into account several amongst the above described aspects. My contributions are twofold:
*methodological/theoretical* as well as *applicative*. From a methodological perspective the
core of my work is in the area of probabilistic systems and revolves around the probabilistic
model checking problem, that is, looking at how to improve on one hand the scalability
and on the other the expressiveness of approaches for the verification of probabilistic sys-
tems. Although the probabilistic model checking has been at the core of my interests,
progressively my activities diversified towards other areas, including that of parameter in-
ference as well as that of introducing novel timed-automata based formalisms. From an
applicative perspective the focus of my work has been mainly in demonstrating the benefit
that the methodological novelties I have introduced brought in different application areas
including that of biological systems, that of manufacturing systems and that of wireless
communication networks.

**Chapter 2** gives an overview of an important contribution in the field of probabilistic model checking. Specifically we introduce a *statistical model checking* (SMC) framework which relies on a novel automata-based formalism, namely the Hybrid Automata Stochastic Logic (HASL), for expressing properties of a probabilistic model. This set of works started through a collaboration with Serge Haddad, Nihal Pekergin and Marie Duflot which was then further extended to Benoit Barbotand resulted with the introduction of the HASL formalism [BDD$^{+11}$c, BDD$^{+11}$a, BBD$^{+15}$a] as well as with the development of the COSMOS software tool [BDD$^{+11}$a, BB22, COS].

In **Chapter 3** we report on a newly introduced approach for parametric verification of temporal properties against stochastic models, i.e., a framework which addresses the problem of estimating the satisfaction probability of a given temporal logic property in function of a model’s parameters. This is achieved through adaptation of the Approximation
Bayesian Computation (ABC) method, a technique for calibrating a model's parameters according to observations of the actual system under study. The adaptation we introduced is based on the notion of satisfiability distance<fnref n="1" /> used to quantify “how far” (the model instance corresponding to) a parameter set is from satisfying the considered property. In practical terms the framework uses stochastic simulation for sampling paths of a model's instance and an HASL monitor automaton (of Chapter 2) for measuring the satisfiability distance of the sampled trajectories. This contributions have been developed as part of the PhD thesis work of Mahmoud Bentriou [Ben21] which I supervised with Paul-Henry Cournède and have led to the publications of papers [BBC19, BBC21].

In **Chapter 4** we present contributions addressing the problem of formal analysis of stochastic oscillators. Although established mathematical approaches, such as Fourier transform and autocorrelation analysis are commonly employed for the analysis of the periodic signals here we face the problem from a different perspective: we look for a *logical* characterisation of periodicity which can be used *i)* to assess oscillation related indicators against a discrete event stochastic model and *ii)* to search the parameter space of a stochastic oscillator so that the resulting model meets given periodicity constraints (e.g. the average duration of the period). Inspired by the PhD thesis work of David Spieler [Spi13] we extend the characterisation of *noisy periodicity* by introducing rich hybrid automaton specifications for assessing both the period and the amplitude of oscillations. Based on those automata we then further adapt the automata-based adaptation of ABC parameter inference outline in Chapter 3 to the case of oscillatory models. The contributions presented in this chapter are based on the following papers [BMM09, BMR12, BD15, Bal15a].

In **Chapter 5** we describe a number of contributions in the field of performance analysis of stochastic models in different domains, including that of biological systems, of telecommunication networks and of manufacturing systems. The overall goal of this chapter is to demonstrate, by taking into account diverse type of systems, to what extent one can take advantage of formal model-based approached to address relevant performance analysis questions.

Although conclusive remarks concerning research directions on each topic are given at the end of the respective chapters in **Chapter 7** a few general perspectives I wish to consider as directions to guide my future research work are discussed.

Finally in **Chapter 6** we give a brief overview of other contributions I have worked on and that are not described in this manuscript. Those contributions span rather heterogeneous domains, including that of *compositional approaches to stochastic model analysis*, that of *modelling formalisms for timed, stochastic, concurrent systems*, that of *improvements to the approximate Bayesian computation approach*, that of *formal performance modelling of manufacturing systems* and also of *wireless network protocols*. A brief description of each of these other contributions is given together pointers with the corresponding publications.

<footnote n="1">¹The dual of *space-time robustness* used in approaches for assessing the robust satisfaction of MITL formulae over non-probabilistic models [DM10].</footnote>
# Chapter 2

## Expressive statistical model checking

Since its introduction in the late-70s/early-80s model checking became a prominent technique for the automatic verification of a system's properties. Given a property formally stated through a temporal logic *property language*, algorithms are provided that allow one for establishing whether the property is *satisfied* by a given model of the considered system. The popularity of model checking is mainly due to three factors: (1) the ability to express specific properties by formulas of an appropriate logic, (2) the firm mathematical foundations based on automata theory and (3) the simplicity of the verification algorithms which has led to the development of numerous tools.

With its seminal work [Pnu77] Pnueli set the ground for model checking by introducing the verification of *linear-time temporal logic (LTL)* properties against discrete-state, untimed models of computer programs. LTL semantics gave rise to so-called *linear-time reasoning* as properties are evaluated w.r.t. *linear sequences of states* (i.e. *paths*) that represent possible executions of the system under analysis. The original work of Pnueli was further developed by both Emerson and Clarke [EC80] and Quelles and Sifakis [QS82], with the introduction of the *computation tree logic (CTL)*, a temporal logic for *branching time* properties, i.e. properties that allows for taking into account possible *branching* along a model's execution, and hence are evaluated against *trees* (rather than *paths*) issued by a model. Both LTL and CTL properties are inherently *qualitative* as they evaluate to either *true* or *false* against a *transition system model*. Example of properties that can be verified through LTL or CTL model checking are 1) "absence of concurrent access to a mutually exclusive shared resource" (safety), or 2) "possibility to access a shared resource infinitely often" (liveness). Such kind of properties are called *qualitative* as, by definition, they cannot include any form of quantification. For example, one can state that *"every message that is sent will eventually be received"* but cannot quantify within which delay it will be received.
Notice that LTL and CTL are distinct property languages (some LTL properties have no CTL equivalent) with a non-empty intersection (some LTL properties have a CTL equivalent and vice versa) and that they are both subsumed by the CTL* [EH86] temporal logic.

**Probabilistic model checking.** As of the mid 90s researchers started to get interested in extending the model checking approach to the verification of probabilistic models. The seminal work of Hansson and Jonsson [HJ90] extended CTL model checking to the verification of discrete-time Markov chains (DTMC), through the introduction of the PCTL temporal logic, namely the probabilistic version of CTL. PCTL allows modellers to express properties which refer to the probability of (discretely) time-bounded events (e.g. *once a message is sent there is at least a 95% probability it will be delivered within 5 time units*). Later on Aziz [ASSB00] first, and Baier *et al.* [BHHK03a] then, further extended CTL model checking to continuous-time Markov chains (CTMCs) by introducing the Continuous Stochastic Logic (CSL) hence allowing for reasoning about the probability of continuously time-bounded events (e.g. *when the queue contains 5 jobs there is at most a 25% chance of a further job arrival within the next 2.75 seconds*). As CTL is a *branching time logic*, so are PCTL and CSL, which means that formulae from both PCTL and CSL are interpreted over the tree consisting of all trajectories unfolded from the considered Markov chain model (either a DTMC for PCTL or a CTMC for CSL). Verifying that a PCTL or CSL specification $\phi$ holds against a given Markov chain requires assessing the probability measures of the sub-tree that satisfies $\phi$, which is well-defined based on the characterisation of the probability space for the trajectories of Markov chains [Kul16]. Specifically, formulae of a probabilistic temporal logic, require comparing the probability that a condition $\phi$ holds against a probability bound $p$, denoted $P_{\sim p}[\phi]$, with $\sim \in \{<, \le, \ge, >\}$ and $p \in [0, 1]$. The *qualitative* fragment of PCTL/CSL is that resulting when the considered probability bound $p$ is either 0 or 1, conversely the *quantitative* fragment is that in which $p \in ]0, 1[$, they are called CSL model checking has then been further extended by Kwiatkowska, Norman and Parker [KNP07a] giving the modeller the possibility to attach (state and/or transition) *rewards* to Markov chain models and to assess *reward based properties* all of which is supported by PRISM [PRI] arguably the most popular probabilistic model checking tool.

**Numerical probabilistic model checking of Markovian models.** Probabilistic model checking mainly targeted Markov chains models only and consisted of *exhaustive* procedures (i.e., algorithms requiring the complete storage of the Markov chain's *transition matrix*), which is in line with the actual notion of *verification* that inherently entails the truth of a formula must be established by considering a model in its entirety. The definition of *exhaustive* probabilistic model checking algorithms for Markov chain models is relatively straightforward thanks to two main factors: 1) the *memoryless property* of exponential distribution which characterises Markovian models and 2) the (somewhat) limited expres-
siveness of the considered temporal logics (e.g. PCTL and CSL). From an algorithmic
point of view the verification of a temporal logic formula against a Markov chain model
is obtained by a combination of mathematical methods that depend on the kind of for-
mula. For example for DTMC models the verification of *untimed* PCTL properties boils
down to solving of a system of linear equations whereas for *time-bounded* PCTL properties
it amounts to matrix-vector multiplications. On the other hand for CTMC models the
verification of *untimed* CSL formulae reduces to a PCTL problem on the corresponding
*embedded DTMC*, whereas the verification of *time-bounded* CSL properties boils down to
a *transient-state analysis* problem for CTMC which is numerically approximated through
*uniformisation*.

**Beyond Markovian models.** Whenever the timing of stochastic events that charac-
terise the considered system's cannot be reasonably assumed to follow an exponential
distribution then the underlying stochastic process is no longer Markovian and hence the
definition of probabilistic model checking algorithms becomes cumbersome. In this con-
text the *transient stochastic state classes* method, by Horváth, Paolieri, Ridi and Vicario
(i.e., a method that allows for keeping track of the evolution in time of the *transient-state*
distribution for non-Markovian models) has been showed an effective means for adapting
algorithms for the verification of time-bounded properties against Generalised Semi Markov
Process (GSMP) models [HPRV11].

**The curse of scalability.** Independently of the targeted class of models and of the expressiveness of the property language the main limitation of *exhaustive approaches* to model checking is their scalability: as the space complexity of exhaustive algorithms depends on the size of the model they can only be applied to models of “reasonable” size, thus‘ ruling out many real-life applications.

**Statistical model checking.** In order to overcome the scalability issue of *exhaustive*
algorithms, as of the mid 2000s, a novel kind of *non-exhaustive* approach, named *statistical*
*model checking* (SMC), has been introduced by the pivotal work of Younes [YS06]. The
main idea with SMC is to employ a finite number of model’s *simulations* in order to obtain
an approximated answer to the question of whether a probabilistic temporal logic formula
$P_{\sim p}[\psi]$ (where $p \in [0, 1]$ and $\sim \in \{<, \leq, \geq, >\}$) holds against a stochastic model $\mathcal{M}$ which
boils down at deciding whether the probability $Pr(\mathcal{M}, \psi) \sim p$ of those paths satisfying
$\psi$ comply with $\sim p$. In this respect there exist two families of SMC approaches: those
which aim at *deciding* whether $Pr(\mathcal{M}, \psi) \sim p$ without actually estimating $Pr(\mathcal{M}, \psi)$ (by
application of *hypothesis testing*), and those which aim at computing an estimate $p'$ of
$Pr(\mathcal{M}, \psi)$ (by application of *confidence interval* methods) and hence establish whether
$Pr(\mathcal{M}, \psi) \sim p$ by comparing the estimate $p' \sim p$.

The clear advantage of SMC is that one does not even need to build (hence store) a
system’s model, it suffices that a syntactical representation, in terms of some modelling
language (equipped with an operational semantics), is available so that a stochastic simula-
tion algorithm can be devised for sampling (finite) trajectories from the model’s state-space.
Each trajectory sampled from M either satisfies (i.e. is accepted) or does not satisfy (i.e.
is rejected) ψ and therefore it corresponds to a realisation of a Bernoulli experiment. As a
consequence assessing the probability of ψ against M in the context of SMC corresponds
to estimating the distribution of a Binomial random variable which, in turns, boils down
to computing the ratio of sampled trajectories that satisfy ψ.

**SMC approaches.** Ever since its introduction a number of SMC formalisms, algorithms and tools have been proposed and their effectiveness demonstrated through analysis of large-scale systems in a variety of fields, including computer networks, manufacturing systems, and biological processes [AP18]. Younes opened the way by introducing an *hypothesis testing* based approach to decide, with arbitrary approximation, whether the probability that a CSL property *φ* is satisfied by a model *M* comply with a threshold *p*. An alternative solution to the same kind of problem has then been proposed by Peyrronet *et al*. [HLMP04] by introducing a *randomized approximation scheme* to establish whether the probability that an LTL property *φ* is satisfied by a DTMC match a given constraint *p* (which has then been integrated amongst the statistical engines of PRISM model checker). Alternatively confidence interval approaches have been proposed to solve the SMC *estimation* problem (i.e. the problem of estimating what is the probability that *φ* is satisfied by *M*) with dedicated efficient sequential schemes discussed by Jegourel *et al*. [JSD18] which have been adopted in many SMC tools (PRISM included).

**Expressiveness: Numerical versus Statistical model checking.** The use of sta-
tistical methods instead of numerical ones allows for relieving the limitations inherent to
stochastic model checking procedures based on numerical methods, both in terms of the
family of stochastic models as well as of the type of properties that can be considered. If
numerical model checking (NMC) is mainly targeting Markovian models (because of the
nice mathematical properties of the Exponential distribution), statistical model checking
(SMC) permit one to trespass the Markovian perimiter, hence to use a very wide range of
distributions, and to *synchronise* such a model with an automaton that includes linearly
evolving variables, complex updates and constraints. In terms of the properties that one
can target SMC allows for a much increased expressiveness. In fact since with NMC the
verification of temporal logic properties boils down to either a *transient-state* or *steady-
state* distribution computation, the kind of properties that one can formulate are limited:
NMC algorithms can only cope with the evaluation of the probability of paths which are
characterised by the temporal evolution of state conditions (something which goes under
the name of *state-based reachability problems*). With SMC, on the other hand, one is com-
pletely freed from the *state-based* reachability perimeter which constrains NMC, therefore
one can easily combine state-conditions with other-conditions, e.g. constraints on given
indicators evaluated along the path (e.g. such as waiting time, number of clients in a
system, production cost of an item), for characterising the target random variable whose
mean value is the target of the SMC verification process. As a result SMC is naturally
suited for easily conceive and compute sophisticated performability parameters.

**In quest for performance oriented property languages.** Despite their success the vast majority of SMC frameworks are based on temporal logic languages (such as PCTL or CSL) that are limited to *state-based* properties, that is, properties that use *state-conditions* as criteria to characterise the behaviour (i.e. the executions) of interest. For example, in a systems consisting of 2 queues where different kind of jobs arrive and compete to access a shared resource an example of property based on *state-conditions* is “*what is the probability that the shared resource is busy when the number of jobs in the first queue does not exceed N?*” In standard performance evaluation problems, however, it makes sense to allow for a wider, *performance-oriented*, class of properties to be considered, i.e., properties that may use *quantitative statistics* (evaluated along a path) as criteria for selecting the executions of interest. For instance, the average length of a waiting queue during a busy period or the mean waiting time of a client are typical examples of criteria that cannot be expressed by the quantitative logics that use *state-conditions* as the only kind of criteria to characterise the executions whose probability shall be assessed.

**Outline of contribution.** This chapter reports on our contribution in the field of expressive SMC. We introduced a novel, hybrid-automata-based, property specification language named *Hybrid Automata Specification Language* (HASL) [BDD$^+$10, BBD$^+$15b] for expressing temporal-dependent performance indicators referred to a discrete-event stochastic model and we introduced the statistical model checking procedure to assess them together with an associated software tool named COSMOS [BDD$^+$11a, COS] to use it. HASL is the evolution of the time-automata (TA) extension of CSL, namely the CSL$^{TA}$ [DHS07] logic (by Donatelli, Haddad and Sproston), in that it replaces TA which, within CSL$^{TA}$ are used as a machinery for specifying the behaviour of interest that one want to study on a system, with linear hybrid automata (LHA). Employing LHA in place of TA opens up to a much increased expressiveness, as more sophisticated criteria for specifying the behavior of interest can be formulated through an LHA. As it turns out the family of properties that can be taken into account through HASL go beyond those of most popular temporal logic languages rendering HASL a specification language more suited for performance analysis as demonstrated by several applications presented in Chapter 4 and Chapter 5.
## 2.1 Preliminaries

### 2.1.1 Discrete event stochastic processes

We target a very broad class of discrete-state, continuous-time stochastic processes which we generically refer to as Discrete Event Stochastic Processes (DESP), but which is often referred to as generalised semi-Markov processes (GSMPs) [Gly83, ACD91]. Differently from the most common probabilistic model checking frameworks such as those supported by the PRISM [KNP11] and Storm [HJK+22] tools, with DESPs we get rid of any restriction about the type of probability distribution allowed for modelling events’ delay, therefore we are free to leave the realm of memorylessness and consider also processes with memory. Syntactically speaking this comes to the price of a more cumbersome formalism as we need to 1) specify the probability distribution associated to each event of the process, 2) account for (time) memory in the specification of an event’s occurrence and 3) disambiguate between concurrently occurring events, which, within DESPs, may have a non null probability to occur.

**Definition 2.1.** A Discrete Event Stochastic Process (DESP) is a tuple $\mathcal{M} = (S, \pi_0, E, \text{Ind}, \text{enabled}, \text{delay}, \text{choice}, \text{target})$ where $S$ is a countable set of states, $\pi_0 \in \text{dist}(S)$ is the initial state distribution, $E$ is a finite set of events, $\text{Ind} : S \to \mathbb{R}$ is a set of state indicator functions, $\text{enabled} : S \to 2^E$ denotes the events enabled in each state, $\text{delay} : S \times E \to \text{dist}(\mathbb{R}^+)$ characterises the probability distribution of the delay of occurrence of an event $e$ in state $s$, $\text{choice} : S \times 2^E \times \mathbb{R}^+ \to \text{dist}(E)$ characterises the probability distribution over the set of concurrently occurring events in a given state and $\text{target} : S \times E \times \mathbb{R}^+ \to S$ describes which state is entered when from state $s$ an enabled event $e \in \text{enabled}(s)$ occurs with a delay $d$.

Given a state $s \in S$ of a DESP $\mathcal{M}$, $\text{enabled}(s)$ is the set of events enabled in $s$. For an event $e \in \text{enabled}(s)$, $\text{delay}(s, e)$ is the distribution of the delay between the enabling of $e$ and its possible occurrence. If $d \in \mathbb{R}^+$ is the earliest delay in some configuration of the process with state $s$, and $E' \subseteq \text{enabled}(s)$ the set of events with minimal delay, $\text{choice}(s, E', d)$ describes how the conflict is randomly resolved: for all $e' \in E'$, $\text{choice}(s, E', d)(e')$ is the probability that $e'$ will be selected among $E'$ after waiting for the delay $d$. The function $\text{target}(s, e, d)$ denotes the target state reached from $s$ on occurrence of $e$ after waiting for $d$ time units.

**DESP semantics.** As the dynamics of a DESP is not necessarily *memoryless* we need to carry on extra information about a system’s state along the unfolding of an execution. Therefore a DESP path amount to a sequences *configurations*, where a configuration consists of a triple $(s, \tau, \text{sched})$ with $s \in S$ the current state, $\tau \in \mathbb{R}^+$ the current time and $\text{sched} : E \to \mathbb{R}^+ \cup \{+\infty\}$ the function that describes the occurrence time of each enabled event $(+\infty$ if the event is not enabled). From a given configuration $(s, \tau, \text{sched})$ a path
is unfolded by iteratively selecting the next event $e \in enabled(s)$ to occur together with
its occurrence delay $d$ leading to the next configuration $(s, \tau, sched) \xrightarrow{d} (s', \tau', sched')$.
The next transition to fire is determined through the following steps. First function sched
determines $E^{min} = \{e \in enabled(s) | \forall e' \in enabled(s), sched(e) \le sched(e')\}$, i.e., the set
of enabled events with minimal delay $d_m = sched(e) - \tau, \forall e \in E'$. Then, if $E^m$ contains
more than 1 element, the probability distribution choice$(s, E^{min}, d_m)$ is used to randomly
pick the actual next event $e_m \in E_m$ to fire. The next state $s'$ is then determined by fir-
ing the previously selected next event to fire $e_m$, i.e. $s' = target(s, e_m, d_m)$ and also by
updating the configuration's time accordingly, i.e., $\tau' = \tau + d_m$. Finally the schedules of
currently enabled events is updated by application of function sched as follows: first for all
events $e' \ne e_m$ that are still enabled, sched($e'$) is maintained (i.e. sched'($e'$) = sched($e$),
$\forall e' \in enabled(s') \cap enabled(s)$, such that $e' \ne e_m$) while for all other "newly" enabled event
$e' \in enabled(s')$, a new delay $d'$ is sampled according to the distribution delay($s'$, $e'$) and
sched($e'$) is set to $\tau + d + d'$. Finally for all disabled events $e' \notin enabled(s')$ the sched'($e'$)
is set to $+\infty$.

**Operational semantics of a DESP.** The characterisation of the operational semantics of a DESP reduces to showing that a sound definition of probability for a generic path $\sigma : (s_0, 0, sched_0) \xrightarrow{d_1}_{e_1} (s_1, \tau_1, sched_1) \xrightarrow{d_2}_{e_2} \dots$ of a DESP exists. This, in turns, boils down to formally characterising a number of relevant families of random variables, namely: the family $e_1, \dots, e_n, \dots$ (events associated to $\sigma$) with support the set of events $E$, the family $s_0, \dots, s_n, \dots$ (the states of $\sigma$) with support the set $S$ of $M$, the family $\tau_0 \le \dots \le \tau_n \le \dots$ (the occurrence time of the events on $\sigma$) with support $\mathbb{R}^+$, the family $\{sched(e)_0, \dots sched(e)_n \dots\}$ with support $\mathbb{R}^+ \cup \{+\infty\}$ (the events' queue along $\sigma$) and the family $\{E_0^{min}, \dots E_n^{min} \dots\}$ (the set of events with minimal delay along $\sigma$). The characterisation of such families of random variables is inductive w.r.t. the length $n$ of the $\sigma$ and is formally given in [BBD$^+$15b].

2.1.2 DESP as Non Markovian Generalised Stochastic Petri Nets

For the sake of practicality within HASL we adopt the non-Markovian extension of Gen-
eralised Stochastic Petri Net (NMGSPNs) as high level modelling language to represent
DESP models. This has two main motivations: 1) to prevent the user from the tedious
and error prone need to specify a DESP by explicit listing of all of its elements (as of Def-
inition 2.1) and 2) to take advantage of the Petri nets semantics which is particularly well
suited to yield highly efficient stochastic simulation algorithms. In essence a NMGSPNs is
a generalised stochastic Petri nets [AMBC$^+$95], that is a Stochastic Petri net formalism
which combines stochastic timed transitions with probabilistic immediate ones, extended
with the possibility to associate generic probability distributions (non necessarily Expo-
nential ones) with stochastic timed transitions, and therefore, according to Definition 2.1,
must be equipped with necessary means to probabilistically disambiguate between events
that might occur at the same time.

**NMGSPN models and token game.** A NMGSPN model is a bipartite graph consist-
ing of *places* (circles) and *transitions* (bars) nodes. Places may contain *tokens* (representing
the state of the modelled system) while transitions indicate how tokens “flow” within the
net (encoding the model dynamics). The state of a NMGSPN (hence of the DESP it rep-
resents) consists of a *marking* indicating the distribution of tokens among the places (i.e.
how many tokens each place contains). A transition *t* is enabled whenever every *input
place* of *t* contains a number of tokens greater than or equal to the multiplicity of the
corresponding (input) arc. An enabled transition may *fire*, consuming tokens (in a number
indicated by the multiplicity of the corresponding input arcs) from its input places, and
producing tokens (in a number indicated by the multiplicity of the corresponding output
arcs) in its *output places*. Transitions can be either *timed* (denoted by thick bars) or *imme-
diate* (denoted by thin filled-in bars, see Figure 2.1). Transitions are characterised by: (1)
a *distribution* which randomly determines the delay before firing it; (2) a *priority* which
deterministically selects among the transitions scheduled the soonest, the one to be fired;
(3) a *weight*, that is used in the random choice between transitions scheduled the soonest
with the same highest priority. With the original GSPN formalism [AMBC$^+$95] the delay
of timed transitions is assumed *exponentially* distributed, whereas with our NMGSPN it
can be given by any distribution. Thus a NMGSPN timed-transition is characterised by a
tuple: $t = (\mathit{type}, \mathit{par}, \mathit{pri}, w)$, where $\mathit{type}$ indicates the type of distribution (e.g. uniform),
$\mathit{par}$ indicates the parameters of the distribution (e.g. $[\alpha, \beta]$), $\mathit{pri} \in \mathbb{R}^+$ is a priority assigned
to the transition and $w \in \mathbb{R}^+$ is used to probabilistically choose between transitions occur-
ring with equal delay and equal priority. The information associated with a transition (i.e.
$\mathit{type}, \mathit{par}, \mathit{pri}, w$) is exploited in different manners depending on the type of transition. For
example for a transition with a continuous distribution the priority ($\mathit{pri}$) and weight ($w$)
records are superfluous (hence ignored) since the probability that the schedule of the corre-
sponding event is equal to the schedule of the event corresponding to another transition is
null. Similarly, for an immediate transition (denoted by a filled-in bar) the specification of
the distribution type (i.e. $\mathit{type}$) and associated parameters ($\mathit{par}$) is irrelevant (hence also
ignored). Therefore these unnecessary informations are omitted in Figure 2.1.

**Example 1 (Mutual exclusive shared resource).** The NMGSPN model of Figure 2.3 (top) describes the behaviour of an open system where two classes of clients (processes) (namely 1 and 2) compete to access a shared resource (e.g. the memory). Class *i*-clients ($i \in \{1, 2\}$) enter the system according to a Poisson process with parameter $\lambda_i$ (corresponding to the exponentially distributed timed transition $Arr_i$ with rate $\lambda_i$, represented by thick empty bars). On arrival, clients cumulate in places $Req_i$ waiting for the resource to be free (a token in place $Free$ witnessing that the resource is available). The exclusive access to the shared resource is regulated either deterministically or probabilistically by the *priority*
<figure><img src="image_1.png" /><figcaption>Figure 2.1: The NMGSPN description of a shared resource system.</figcaption></figure>

(priᵢ) and the weight (wᵢ) of immediate transitions *Start*₁ and *Start*₂. Thus in presence of a competition (i.e. one or more tokens in both *Req*₁ and *Req*₂) a class *i* client wins the competition with a class *j* ≠ *i* client with probability 1 if priᵢ > priⱼ, and with probability wᵢ/(wᵢ+wⱼ) if priᵢ = priⱼ. The occupation time of the resource by a class *i* client is assumed to be uniformly distributed within the interval [αᵢ, βᵢ] (corresponding to the thick filled-in transitions *Serv*ᵢ). Thus on firing transition *Serv*ᵢ the resource is released and a class *i* client leaves the system. Notice that graphically Exponential, respectively non-Markovian, timed transitions are represented by thick empty, respectively filled-in, bars).

## 2.2 Hybrid Automata Stochastic Language

Alur and Henzinger's hybrid automata [ACHH92, ACHH93] are a popular formalism for modelling linear hybrid systems, i.e. systems consisting of interacting digital and analog components. Although reachability problems for hybrid automata are, generally speaking, undecidable, Henzinger et al. proved that under certain conditions reachability decidability can be established [HKPV98] leading to the development of dedicated model checkers such as the HyTech tool [HHWT97]. Here we consider hybrid automata from a different perspective, that is, as a language to specify properties of a system rather than as a language to model a system. To this aim we introduce the Hybrid Automata Stochastic Language (HASL) [BDD$^+$11b, BBD$^+$15a], the basic idea it being to employ a tailored linear hybrid automaton (LHA) $\mathcal{A}$ as a *monitor* for selecting relevant paths of a DESP model $\mathcal{M}$, based on sophisticated selection criteria. This allows one to conceive very sophisticated performance indicators of a model $\mathcal{M}$, by means of a dedicated automaton $\mathcal{A}$ and to assess them through *synchronisation*, that is, through a process by which paths $\sigma \in \text{Path}(\mathcal{M})$ are processed *on-the-fly* and are either accepted or rejected by $\mathcal{A}$, while relevant statistics are maintained in the automaton's variables and are eventually used to produce estimates
<figure><img src="image_2.png" /><figcaption>Figure 2.2: HASL-SMC schema: sampled paths are filtered by a LHA and the accepted ones are used for a confidence interval estimate of the target measure Z.</figcaption></figure>

of the performance indicators of interest.

In practical terms a HASL formula $\phi$ consists of two components $\phi = (\mathcal{A}, Z)$, where $\mathcal{A}$ is the *monitor* LHA and $Z$ is an expression that specifies the quantity to be evaluated (based on the $\mathcal{A}$'s data variables).

**HASL model checking scheme.** The synchronisation between a DESP model $\mathcal{M}$ and a LHA $\mathcal{A}$ naturally leads to definition of a statistical model checking procedure for the HASL formalism which is outlined in Figure 2.2. A HASL model checker takes a DESP model $\mathcal{M}$ and a HASL formula $(\mathcal{A}, Z)$ as inputs together with two auxiliaries parameters, i.e., the confidence level $(\epsilon)$ and the width $(\delta)$ of the confidence interval estimation process. The LHA part of a HASL formula is used to select, via synchronisation, the paths sampled on the DESP model and while so doing the (real-valued) variables of the LHA are updated with relevant statistics of the synchronising path. Such statistics are then employed for producing a confidence-interval estimation of the mean value of the considered expression $Z$ (i.e. the second part of the HASL formula). $Z$ is the quantity the modeller wants to assess w.r.t. the paths accepted by LHA.

### 2.2.1 Synchronised Linear Hybrid Automata

In the context of HASL model checking we introduce an adaptation of Alur and Henzinger's hybrid automata [ACHH92] which we refer to as synchronised LHA the main difference being that HASL automata are designed to interact with a DESP model, hence they must fulfil certain syntactical as well as semantical constraints. From a syntactical perspective there are two elements that impact the syntax of a HASL hybrid automaton definition: first, an automaton edges are now labeled either with events of the DESP model, indicating that the edge traversal takes place only through synchronised occurrence of a corresponding DESP event or, alternatively, through a special label $\sharp$, indicating that the edge traversal happens autonomously without synchronisation with any DESP event, and second, an automaton's variables can evolve at a rate which depends on the current state of the
DESP model (i.e. the rate of a variable is given by a DESP *indicator* function). From a
semantical perspective we need to establish a number of *determinism-related* constraints
so to guarantee the unicity of synchronisation, that is, the fact that there exists always a
single execution of the synchronised product process (*M* × *A*). All of these aspects are
covered by the following Definition.

**Definition 2.2.** A synchronised linear hybrid automaton, referred to a DESP model $\mathcal{M}$ with states $S$, events $E$ and indicator functions $\text{Ind} : S \to \mathbb{R}$, is defined by a tuple $\mathcal{A} = \langle E, L, \Lambda, \text{Init}, \text{Final}, X, \text{flow}, \to \rangle$ where: $E$ is a finite alphabet of events (the events of $\mathcal{M}$); $L$ is a finite set of locations; $\Lambda : L \to \text{Prop}$ denotes the locations' invariants<fnref n="1" />; $\text{Init}, \text{Final} \subset L$ denotes the initial, resp. final, locations; $X = (x_1, \dots, x_n)$ is a $n$-tuple of data variables, and:

*   *flow :* $L \mapsto \text{Ind}^n$ associates with each location one indicator per data variable representing the evolution rate of the variable in the location.

*   $\rightarrow \subseteq L \times ((2^E \times \text{Const}) \uplus (\{\#\} \times \text{ICons})) \times \text{Up} \times L$ is a set of edges, where *Const*, (resp. *ICons*), denotes the set of *linear constraints*, (resp. *left-closed constraints*) consisting of a boolean combination of inequalities of the form $\sum_{1 \le i \le n} \alpha_i x_i + c < 0$ where $\alpha_i, c \in \text{Ind}$ are DESP indicator functions and $< \in \{=, <, >, \le, \ge\}$, $\uplus$ denotes the disjoint union and $\text{Up}$ denotes the set of $n$-tuples of update functions with the following form $x_k = \sum_{1 \le i \le n} \alpha_i x_i + c$ where again $\alpha_i, c \in \text{Ind}$.

Furthermore $\mathcal{A}$ fulfils the following conditions.

*   **Initial determinism:** $\forall l \neq l' \in \text{Init}, \Lambda(l) \wedge \Lambda(l') \Leftrightarrow \text{false}$. This must hold whatever the interpretation of the indicators occurring in $\Lambda(l)$ and $\Lambda(l')$.

*   **Determinism on events:** $\forall E_1, E_2 \subseteq E \text{ s.t. } E_1 \cap E_2 \neq \emptyset, \forall l, l', l'' \in L, \text{ if } l'' \xrightarrow{E_1,\gamma,U} l \text{ and } l'' \xrightarrow{E_2,\gamma',U'} l' \text{ are two distinct transitions, then either } \Lambda(l) \wedge \Lambda(l') \Leftrightarrow \text{false} \text{ or } \gamma \wedge \gamma' \Leftrightarrow \text{false.} \text{ Again this equivalence must hold whatever the interpretation of the indicators occurring in } \Lambda(l), \Lambda(l'), \gamma \text{ and } \gamma'$.

*   **Determinism on $\sharp$:** $\forall l, l', l'' \in L, \text{ if } l'' \xrightarrow{\sharp, \gamma, U} l \text{ and } l'' \xrightarrow{\sharp, \gamma', U'} l' \text{ are two distinct transitions, then either } \Lambda(l) \wedge \Lambda(l') \Leftrightarrow \text{false} \text{ or } \gamma \wedge \gamma' \Leftrightarrow \text{false.}$

*   **No $\sharp$-labelled loops:** For all sequences $l_0 \xrightarrow{E_0, \gamma, U_0} l_1 \xrightarrow{E_1, \gamma, U_1} \dots \xrightarrow{E_{n-1}, \gamma_{n-1}, U_{n-1}} l_n$ such that $l_0 = l_n$, there exists $i \le n$ such that $E_i \ne \sharp$.

In the remainder we write $l \xrightarrow{E',\gamma,U} l'$ to denote an edge $(l, E', \gamma, U, l') \in \rightarrow$ of an LHA,
furthermore we refer to an edge such as $l \xrightarrow{\gamma,E',U} l'$, where $E' \subseteq E$ is a subset of events,

<footnote n="1">¹Where *Prop* ⊂ *Ind* is the subset of DESP indicator functions that evaluates to {0, 1}</footnote>
as a synchronising edge as opposed to an edge such as $l \xrightarrow{\gamma, \# , U} l'$ which we refer to as an autonomous edge.

**Example 2.** The two automata $\mathcal{A}_1$ and $\mathcal{A}_2$ in Figure 2.3 are designed to be synchronised with the NMGSPN model of mutual exclusion of Example 1. Both $\mathcal{A}_1$ and $\mathcal{A}_2$ consists of

<figure><img src="image_3.png" /><figcaption>Figure 2.3: Examples of two LHA referred to mutual exclusion NMGSPN model</figcaption></figure>

2 locations $L=\{l_0, l_1\}$, an initial one $Init=\{l_0\}$ and a final one $Final=\{l_1\}$. They both refer to the same event set $E = \{Arr_1, Start_1, Serv_1, Arr_2, Start_2, Serv_2\}$ i.e. that of the transitions of the NMGSPN model of mutual exclusion and both are equipped with 3 data variables, though $\mathcal{A}_1$ has a clock variable ($t$) whereas $\mathcal{A}_2$ is not equipped with any clock. The flow of variables is conventionally null in the final location ($Flow(l_1) = (0, 0, 0)$) as once in the final location the synchronisation ends$^2$, whereas is non-null (for certain variables)

<footnote>$^2$in fact the flow of variables in any final location could equivalently be set to an arbitrary value as it cannot affect the synchronisation process given the synchronisation stop precisely on entering a final location</footnote>
in the initial location $l_0$: variable $x_1$ (for both $\mathcal{A}_1$ and $\mathcal{A}_2$) evolves at a rate given by the marking of place $Req_1$ of the NMGSPN model, whereas the clock $t$, for $\mathcal{A}_1$, naturally has constant rate of 1. As for the edges $\mathcal{A}_1$ has 2 synchronising self-loop edges on $l_0$, $l_0 \xrightarrow{\{Serv_1\},t<T,\{x_2++\}} l_0$, which synchronises with the occurrences of the $Serv_1$ transition of the NMGSPN model (and in so doing it stores the occurrences of $Serv_1$ in variable $x_2$, and $l_0 \xrightarrow{E\{Serv_1\},t<T,\emptyset}$ $l_0$ which synchronises with any other transition (without registering any information in the variables of $\mathcal{A}_1$) plus an autonomous edge $l_0 \xrightarrow{\sharp,t=T,\{x_1/=T\}} l_1$ which establishes the end of the synchronisation as soon as the clock is $t=T$ (where $T$ is a parameter of $\mathcal{A}_1$), therefore paths accepted by $\mathcal{A}_1$ are *time-bounded*. Similarly $\mathcal{A}_2$ has 3 synchronising self-loop edges on $l_0$, one for counting the occurrences of $Serv_1$ (stored in $x_2$), one for counting the occurrences of $Start_1$ (stored in $x_3$) and the remaining for ignoring any other transition occurring during the synchronisation. $\mathcal{A}_2$ autonomous edge establishes the ending of the synchronisation as soon as $k$ (where $k$ is a parameter of $cal\mathcal{A}_2$) occurrences of the $Serv_1$ transition have been observed, therefore paths accepted by $\mathcal{A}_2$ are *event-bounded*.

## LHA semantics Intuitively

**Real-valued variables with linear flow.** An LHA is equipped with a *n*-tuple of real-valued variables $x_1, \dots, x_n \in \mathbb{R}$ (called data variables) whose value evolves with a linear rate that depends both on the location of the automaton and on the current state of the DESP<fnref n="3" />. More precisely, a *flow function* associates with each location *l* of the automaton an *n*-tuple of real-valued *indicator functions* (one for each variable) which define how the value of each variable evolves, in function of the current state of the DESP, while the LHA spends time in *l*. For example $flow_i(l)(s)$ (where $flow_i(l)$ is the *i*-th component of $flow(l)$) denotes the flow of variable $x_i$ while the LHA is in *l* and the DESP in *s*.

**Autonomous versus synchronised transitions.** An edge of a LHA is characterised by 3 components: (*eventSet*, *constraint*, *update*), where *eventSet* denotes either a subset of *E*, the event set of the DESP, or the special symbol $\sharp$ (denoting a pseudo-event that is not included in the event set *E*), *constraint* is a boolean proposition composed of linear inequalities on the LHA variables $x_i$ and *update* denotes a set of equations describing how the value of variables $x_i$ is updated when the edge is traversed. Edges are enabled (depending on the LHA variables and possibly also on the state of the DESP) as soon as the *constraint* condition is satisfied. Depending on *eventSet* we distinguish between two type of edges: *autonomous* and *synchronised*. When enabled an autonomous ($\sharp$) edge can be traversed *autonomously* (i.e. without synchronisation with the occurrence of a corresponding DESP event). On the contrary a synchronised edge is associated to a subset of

<footnote n="3">³the value of $x_i$ evolves in a location of the LHA according to a rate defined by a real-valued linear function of some model's *indicator*</footnote>
events *eventSet* ⊆ *E*, meaning that, when enabled, it can be traversed only through synchronisation with the occurrence of a DESP event *e* ∈ *eventSet*. Notice that autonomous edges have priority on synchronised ones: if both are enabled the autonomous edge will be traversed before the synchronous. Because of the co-existence of autonomous and synchronised edges an LHA is thus capable of taking into account, on one hand the system behaviour, through synchronised transitions, but also to autonomously move the automaton to another location, through autonomous transitions, as soon as specific conditions captured by the data variables *x*ᵢ take place.

A *constraint* of an LHA edge is a boolean combination of inequalities of the form $\sum_{1 \le i \le n} \alpha_i x_i + c < 0$ where $\alpha_i, c \in \mathbb{Ind}$ are indicators, and $< \in \{=, <, >, \le, \ge\}$. The set of constraints is denoted by *Const*. Given a location and a state, an expression of the form $\sum_{1 \le i \le n} \alpha_i x_i + c$ linearly evolves with time. An inequality thus gives an interval of time during which the constraint is satisfied. We say that a constraint is left closed if, whatever the current state *s* of the DESP (defining the values of indicators), the time at which the constraint is satisfied is a union of left closed intervals. These special constraints are used for the “autonomous” edges, to ensure that the first time instant at which they are satisfied exists. We denote by $\mathbb{IConst}$ the set of left closed constraints.

An *update* is more general than the reset of timed automata. Here each data variable can be set to a linear function of the variables' values. An update *U* is then an *n*-tuple of functions *u*₁, ..., *u*ₙ where each *u*ₖ is of the form *x*ₖ = Σ₁≤i≤n αᵢ*x*ᵢ + *c* where the αᵢ and *c* are indicators. The set of updates is denoted by *Up*.

Finally each location of an LHA is also associated with an *invariant* (through a labelling function Λ), i.e. a boolean valued proposition built on top of DESP indicators. A location invariant plays also a role in the enabling of LHA edges (autonomous and synchronised): given in the current location *l* the constraint *c* of an edge *l* $\xrightarrow{-.c,u}$ *l*' is satisfied, the edge *l* $\xrightarrow{-.c,u}$ *l*' is actually enabled only if the invariant of Λ(*l*') will also be satisfied after traversing *l* $\xrightarrow{-.c,u}$ *l*' (hence after applying the update *u* to the LHA variables).

**Remark 1** (Unique synchronisation). The four constraints outlined in Definition 2.2 let the automata we consider deterministic in the following (non usual) sense. Given a path σ of a DESP, there is at most one synchronisation with the linear hybrid automaton. The three first constraints ensure that the synchronised system is still a stochastic process. The fourth condition disables “divergence” of the synchronised product, i.e. the possibility of an infinity of consecutive autonomous events without synchronisation.

**Remark 2** (Linearity). We stress that the linearity restriction that concerns different elements of an LHA (i.e. the guards and the updates of an arcs as well as the flow of the variables) can be relaxed, as long as those elements do not concern autonomous transitions. Polynomial evolution of constraints could easily be allowed for synchronised edges for which we would just need to evaluate the expression at a given time instant. Since the best
algorithms solving polynomial equations operate in PSPACE [Can88], such an extension
for autonomous transitions cannot be considered for obvious efficiency reasons.

**Example 3.** Figure 2.4 depicts two slight variants of an LHA for measures related to the time the shared resource of the mutual exclusion model is occupied by class 1 clients or by class 2 clients. Both automata employ 2 variables: $x_0$ that counts the number of completed occupation by either class of clients and $x_1$ dedicated to measuring the difference between the time the resource is occupied by clients of type 1 minus the time it is occupied by clients of type 2, and to this aim variable $x_1$ has flow 1 (resp. -1) when the resource is used by class 1 (resp. 2) clients, and 0 when the resource is not used. The termination condition, given in the constraint of arcs leading to the final location $l_3$, is defined so that synchronisation ends on observing the $k$-th occurrence of either $Serv_1$ or $Serv_2$ (corresponding to the completion of occupation of the shared resource by either type of client), i.e. $k$ being the sum of completed occupation of the resource. The most peculiar characteristic of the automata in Figure 2.4 is that they both use location invariants, (i.e. state-dependent proposition, corresponding to function $\Lambda$ of Definition 2.2) associated with the non-final locations: proposition **Acc1** $> 0$ (there is at least a token in place **Acc1**) associated to $l_1$, **Acc2** $> 0$ (there is at least a token in place **Acc2**) associated to $l_2$ and **Free** $> 0$ (there is no token neither in **Acc1** nor in **Acc2**) associated with the initial location $l_0$. Hence, starting from location $l_0$ (assuming initially the resource is free), no matter which precise event occurs in the system, the automaton will switch from $l_0$ to $l_1$ and $l_2$ depending on which class of clients has access to the resource. The fact that for example three different transitions labelled with $E$ without any constraint are available in location $l_0$ does not induce non determinism as only one of these transitions is possible at a time thanks to the location invariants. The synchronisation terminates in location $l_3$ on observation of the $k$-th departure of (any) client from the system and on condition that type 1 clients have used the resource longer than type 2's (constraint $x_1 \ge 0$ for $A_{3a}$) or viceversa on condition that type 2 clients have used the resource longer than type 1's (constraint $x_1 < 0$ for $A_{3b}$). This allows to assess, through a corresponding HASL expression $Z \equiv PROB()$ (see Section 2.2.2) the probability that type 1 (resp. type 2) clients use the resource longer than type 2 (resp. type 1).

**Operational Semantics: synchronised process $M \times A$.** The formal characterisation of such synchronisation process boils down to showing that the synchronised process is itself a DESP (equipped with absorbing accept/reject states) and by providing its formal definition (Definition 2.3) in function of the synchronising components, that is the DESP model $M$ and the LHA $A$.

**Notations.** Given a DESP model $M = \langle S, \pi_0, E, Ind, enabled, delay, choice, target \rangle$ and an LHA automaton $A = \langle E, L, \Lambda, Init, Final, X, flow, \rightarrow \rangle$ in the remainder we adopt the
<figure><img src="image_4.png" /><figcaption>Figure 2.4: An LHA to assess the probability that class 1 (resp. class 2) clients use the shared resource longer than class 2’s. (resp. class 1’s).</figcaption></figure>

following notations:

* LHA variables valuations: we denote $\nu: X \to \mathbb{R}$ a valuation of the LHA variables and $Val$ the set of all possible valuations. For $\nu \in Val$ and $x_i \in X$ we denote $\nu(x_i)$ the value of $x_i$ in $\nu$. We denote $\nu = 0$ the valuation that assigns 0 to each variables $x_i \in X$.

* LHA constraints valuations: let $\nu \in Val$ be a valuation, $s \in S$ a state of $M$ and $exp = \sum_{1 \le i \le n} \alpha_i x_i + c$ an expression related to variables $x_i \in X$ and indicators $\alpha_i, c \in Ind$ we denote $exp(s, \nu)$ the interpretation of $exp$ w.r.t. $\nu$ and $s$ where $exp(s, \nu)$ is defined by $exp(s, \nu) = \sum_{1 \le i \le n} \alpha_i(s)\nu(x_i) + c(s)$.

* LHA constraints satisfaction: let $\nu \in Val$ be a valuation, $s \in S$ a state of $M$ and $\gamma \equiv exp \prec 0 \in (Const \cup |Const)$ a constraint built on top of an expression $exp = \sum_{1 \le i \le n} \alpha_i x_i + c$, we write $(s, \nu) \models \gamma$ (i.e. $(s, \nu)$ satisfies $\gamma$) if $exp(s, \nu) \prec 0$.
* LHA updates valuation: let $\nu \in Val$ be a valuation, $s \in S$ a state of $M$ and $U = (u_1, \dots, u_n) \in Up$ be an update (appearing on some edge) of the LHA we denote by $U(s, \nu)$ the valuation defined by $U(s, \nu)(x_k) = u_k(s, \nu)$ for $1 \le k \le n$.

* state proposition satisfaction: Let $s \in S$ be a state of $M$ and $\varphi \equiv \sum_{1 \le i \le n} \alpha_i + c < 0$ be a state proposition (i.e. a constraint built on an expression that involves only indicators of $M$ and no variables of the LHA). We write $s \models \varphi$ (i.e. $\varphi$ is satisfied in $s$) if $\varphi(s) = true$.

* Time elapsing: Given a state $s \in S$ of $M$ a non final location $l \in L$ and a valuation $\nu$ of $A$, we define the effect of time elapsing by: $Elapse(s, l, \nu, \delta) = \nu'$ where, for every variable $x_k$, $\nu'(x_k) = \nu(x_k) + flow_k(l)(s) \times \delta$.

* Earliest autonomous delay: Given a state $s \in S$ of $M$ a non final location $l \in L$ and a valuation $\nu$ of $A$ we denote Autdel(s, l, v) the minimum delay which enables an outgoing autonomous edge originating in l,

$$ Autdel(s, l, \nu) = \min(\delta | \exists l \xrightarrow{\sharp, \gamma, U} l' \wedge s \models \Lambda(l') \wedge (s, Elapse(s, l, \nu, \delta)) \models \gamma) $$

**Remark:** Whenever Autdel(s, l, v) is finite, we know that there is at least one executable transition with minimal delay and, thanks to the “determinism on $\sharp$” (Definition 2.2), we know that this transition is unique. We call the transition with such earliest delay the earliest and we denote Next(s, l, v) its target location and Umin(s, l, v) its update.

**Definition 2.3 (M×A synchronised process).** Let $M = \langle S, \pi_0, E, Ind, enabled, delay, choice, target \rangle$ be a DESP and $A = \langle E, L, \Lambda, Init, Final, X, flow, \rightarrow \rangle$ a corresponding LHA we define $M' = \langle S', \pi'_0, E', Ind', enabled', target', delay', choice' \rangle$ as the DESP corresponding to the synchronisation of $M$ with $A$ as follows:

* $S' = (S \times L \times Val) \uplus \{\perp\}$ among which $(S \times Final \times Val) \uplus \{\perp\}$ are the absorbing states.

* $\pi'_0(s, l, v) = \begin{cases} \pi_0(s) & \text{if } (l \in Init \wedge s \models \Lambda(l) \wedge v = 0) \\ 0 & \text{otherwise} \end{cases}$ and $\pi'_0(\perp) = 1 - \sum_{s \in S, l \in L, v \in Val} \pi'_0(s, l, v)$.

Note that this definition gives a distribution since, due to “initial determinism” (Definition 2.2), for every $s \in S$, there is at most one $l \in Init$ such that $s \models \Lambda(l)$.

* $E' = E \uplus \{\sharp\}$

* $Ind' = \emptyset$. In fact $Ind'$ is useless since there is no more synchronisation to make.
• $enabled'(s, l, v) = \begin{cases} enabled(s) \cup \{\sharp\}, & \text{if } Autdel(s, l, v) \neq \infty \\ enabled(s), & \text{otherwise} \end{cases}$

• $delay'((s, l, v), e) = \begin{cases} delay(s, e), & \text{if } e \in enabled(s) \land \sharp \notin enabled'(s, l, v) \\ Dirac(Autdel(s, l, v)) & \text{if } \sharp \in enabled'(s, l, v) \end{cases}$

Notice that the Dirac distribution centred on the earliest autonomous delays indicates the priority of autonomous transitions over synchronised ones.

• $choice'((s, l, v), E', d)(e) = \begin{cases} 1, & \text{if } (\sharp \in E' \land e = \sharp) \\ 0, & \text{if } (\sharp \in E' \land e \neq \sharp) \\ 0, & \text{if } e \notin E' \\ choice(s, E', d)(e) & \text{otherwise} \end{cases}$

Again this is coherent since, as soon as $\sharp \notin E'$, then $E'$ is a subset of $enabled(s)$ on which choice is thus defined.

• target for occurrence of a synchronising event $e \in E$

$$target'((s, l, v), e, d) = \begin{cases} (target(s, e, d), l', v'), & \text{if } (\exists l \xrightarrow{E', \gamma, U} l' \in \mapsto) \land (e \in E' \cap enabled(s)) \\ \land Elapse(s, l, v, d) \models \gamma \\ \land target(s, e, d) \models \Lambda(l') \\ \land v' = U(Elapse(s, l, v, d)) \\ \bot & \text{otherwise} \end{cases}$$

Notice that due to the determinism constraints there is at most one such transition.

• target for occurrence of an autonomous event $\sharp$

$$target'((s, l, v), \sharp) = \begin{cases} (s, l', v'), & \text{if } (\exists l \xrightarrow{\sharp, \gamma, U} l' \in \mapsto) \land (\sharp \in enabled'(s, l, v)) \\ \land l' = Next(s, l, v) \\ \land v' = Umin(s, l, v)(Elapse(s, l, v, Autdel(s, l, v))) \\ \bot & \text{otherwise} \end{cases}$$

• In order to get a DESP, for any absorbing state one adds a special event **tick** only enabled in absorbing states with Dirac distribution on value 1 (so that time diverges in these states).

**Summary remarks.** Roughly speaking the formal characterisation of synchronised ($M \times A$) process given above (Definition 2.3) tells us that as long as the automaton is not in a final state, the synchronisation of $M$ with $A$ waits for the first (synchronised or autonomous) transition to occur. If an autonomous transition occurs then only the location
of the automaton $\mathcal{A}$ and the valuation $\nu$ of its variables change, whereas, if a synchronising event of $\mathcal{M}$ occurs then, either the LHA can take a corresponding transition and the system goes on with the next transition or the system goes to a dedicated rejecting state $\perp$ implying the immediate end of the synchronisation. In case of a conflict between two transitions, an autonomous and a synchronised one, the autonomous transition is taken first. Furthermore note that, by initial determinism, for every $s \in S$ there is at most one $l \in \text{Init}$ such that $s$ satisfies $\Lambda(l)$. In case there is no such $l$ the synchronisation starts and immediately ends up in the additional state $\perp$. Determinism on events (resp. on $\#)$ ensures that there is always at most one synchronised (resp. autonomous) transition fireable at a given instant.

<figure><img src="image_5.png" /><figcaption>Figure 2.5: Example of synchronisation of a NGGSPN model with a LHA.</figcaption></figure>
**Example 4 (LHA and DESP synchronisation).** Figure 2.5 depicts an example of synchronisation between the NMGSPN model $\mathcal{M}$ (top left) of the mutual exclusion system (Example 1) with the LHA $\mathcal{A}_1$ of Example 2 (top right). The example shows a path $\sigma$ (sampled from $\mathcal{M}$ and whose projection $\sigma_{Req_1}$ over the $Req_1$ place is depicted as the blue step-like plot in the bottom picture) and the corresponding path $\sigma_{Req_1} \times \mathcal{A}_1$ induced by $\sigma$ on the synchronised product process $\mathcal{M} \times \mathcal{A}_1$. The bottom picture includes also a plot (i.e. red linear plot) of the corresponding time evolution (induced by synchronisation with $\sigma$) of the $x_1$ variable of $\mathcal{A}_1$ that is the variable storing the value of the integral of the number of class 1 processes waiting for the shared resource.

We quickly recall that $\mathcal{A}_1$ uses a clock variable $t$ to measure time elapsing and a variable $x_1$ which, while spending time in $l_0$, register the integral of class 1 processes waiting to get access to the shared resource ($Flow(l_0)(x_1) = Req_1$). The initial state of path $\sigma$ is assumed to be ($Req_1: 1$, $Req_2: 1$, $Free: 1$, $Acc_1: 0$, $Acc_2: 0$)$^4$, i.e. we assume that initially 1 process of each class is waiting to access the shared resource and that the resource is free hence the initial corresponding configuration of the $\mathcal{M}_{Req_1} \times \mathcal{A}_1$ product process is $(1, l_0, [t=0, x_1=0, x_2=0])$ denoting that the marking of $Req_1$ is 1, the location of $\mathcal{A}_1$ is $l_0$ and that in the initial valuation all $\mathcal{A}_1$ variables are null $[t=0, x_1=0, x_2=0]$. The initial state of $\sigma$ does not change until the first event ($Arr_1$) occurs after 0.5 time units yielding a change to state $(2, 1, 1, 0, 0)$, i.e. an increment in the tokens in $Req_1$ (the number of class 1 processes waiting increases by 1). While $\sigma$ is sojourning in the initial state the automaton spends time in $l_0$ hence $x_1$ increased is value with rate $Req_1=1$ hence at $t=0.5$, its value is $x_1=0.5$, which corresponds to the transition to the second configuration of $\mathcal{M}_{Req_1} \times \mathcal{A}_1$, that is $(2, l_0, [t=0.5, x_1=0.5, x_2=0])$. Notice that as long as $t < 4$ the autonomous transition $l_0 \xrightarrow{\sharp, t=4, x_1/4} l_1$ of $\mathcal{A}_1$ is not enabled hence the synchronisation cannot stop. The next transition in $\sigma$ (again $Arr_1$) happens after further 1.5 time units, at that moment the clock is $t=2$, while $x_1$ continued to grow but at a greater rate $\dot{x}_1 = 2$ as for $t \in [0.5, 2]$ the marking of place is $Req_1 = 2$. Therefore the occurrence of the second $Arr_1$ events yields the third configuration of $\mathcal{M}_{Req_1} \times \mathcal{A}_1$, that is $(3, l_0, [t=2, x_1=3.5, x_2=0])$. The following event in $\sigma$ is the occurrence of the immediate transition $Start_2$ representing the fact that a class 2 process starts to occupy the shared resource: notice that the occurrence of $Start_2$ although yields a synchronisation with the automaton does not actually change the configuration of the product process as neither the location nor the variables of $\mathcal{A}_1$ are affected by the occurrence of $Start_2$. The synchronisation goes on in a similar fashion and ends with an autonomous transition through which the final location $l_1$ of $\mathcal{A}_1$ is reached from configuration $(2, l_0, [3.5, 7.5, 1])$ after a delay of 0.5 units, which corresponds with the clock reaching the value $t=4$ hence fulfilling the constraint for the $l_0 \xrightarrow{\sharp, t=4, x_1/4} l_1$ of $\mathcal{A}_1$ edge. By traversing such edge $x_1$ is updated to $x_1 = 8/4 = 2$ which indeed corresponds with the mean number of class 1 process that have been waiting to access the resource in

<footnote>$^4$we represent states of $\mathcal{M}$ as 5-tuples with the following correspondence with the places name ($Req_1$, $Req_2$, $Free$, $Serv_1$, $Serv_2$)</footnote>
the time interval $[0, 4]$.

**Remark 3** (Almost-sure termination of the synchronisation process). A necessary condition for establishing a statistical model checking framework based on the HASL formalism is the guarantee that a stochastic simulation algorithm that (based on the operational semantics outlined in Definition 2.3) samples paths from the synchronised $M \times A$ product process always terminates. Generally speaking HASL model checking consists of an iterative procedure for samping infinite path of $M \times A$ however the quantity $Z$ of interest (see Section 2.2.2) is evaluated on the finite prefix that ends when the path reaches an absorbing state. In the general case, it would be possible that an infinite path does not admit such a prefix. Here we assume that given a DESP $M$ and a HASL formula $(A, Z)$, with probability 1, the synchronising path generated by a random execution path of $M$ reaches an absorbing configuration. This semantical assumption can be ensured by structural properties of $A$ and/or $M$. For instance, the time bounded Until of CSL guarantees this property. As a second example, the time unbounded Until of CSL also guarantees this property when applied on finite CTMCs where all terminal strongly connected components of the chain include a state that either fulfils the target sub-formula of the Until operator or falsifies the continuing sub formula. This (still open) issue is also addressed in [SVA05a, HJB+10].

### 2.2.2 HASL expressions

The second component of an HASL specification is an expression, denoted $Z$, given by grammar (2.1). $Z$ is associated to a LHA $A$ and expresses the target measure whose confidence interval should be estimated based on the paths accepted by $A$.

$$ 
\begin{align} 
Z &::= \mathbb{E}(Y) \mid Z + Z \mid Z \times Z \mid Pdist \nonumber \\ 
Pdist &::= PDF(Y, step, start, stop) \mid CDF(Y, step, start, stop) \mid PROB() \tag{2.1} \\ 
Y &::= c \mid Y + Y \mid Y \times Y \mid Y/Y \mid last(y) \mid min(y) \mid max(y) \nonumber \\ 
y &::= c \mid x \mid y + y \mid y \times y \mid y/y \nonumber 
\end{align} 
$$

There are two main types of expressions $Z$: $\mathbb{E}(Y)$ (where $\mathbb{E}^5$ indicates mean value of, i.e. the first moment of) and $Pdist$ (indicating a probability distribution or probability value expression). $Y$ represent a random variable built on top of algebraic combination of some path operators applied to an algebraic expression composed of LHA variables $y$, i.e. $last(y)$ (i.e. the last value that $y$ has at the end of an accepted path, $min(y)(resp.max(y))$, the min (resp. max) value of $y$ along an accepted path. Conversely $Z$ expressions of $Pdist$ type include $PDF(Y, step, start, stop)$, which allows for estimating the PDF of random variable $Y$

<footnote n="5">⁵Notice that within the COSMOS tool $\mathbb{E}$ is actually denoted $AVG$, therefore HASL expressions may equivalently be formulated using either $\mathbb{E}(Y)$ or $AVG(Y)$.</footnote>
computed by discretisation of the support set [start, stop] ⊂ R≥0 in (stop−start)/step sub-intervals of size step and similarly CDF(Y, step, start, stop). Finally expression PROB() allows for estimating the probability that a path is accepted, otherwise said PROB() is used to estimating the probability of the paths event set represented by the considered automaton A. Notice that one can then define sophisticated probability measures by considering specific conditions as constraints of the edge(s) leading to an accepting location of A.

Notice that since Z can be a composed expression obtained by algebraic operators (e.g. Z + Z, Z × Z) built on top of the first moment of Y (AVG(Y)) it is possible to take into account diverse significant characteristics of Y (apart from its expectation) as the quantity to be estimated, including, for example, the variance Var[Y] = E[Y²] − E[Y]² and the covariance Covar[Y₁, Y₂] = E[Y₁Y₂] − E[Y₁]E[Y₂]. Furthermore we point out that for efficiency reasons, in the implementation of the COSMOS software tool (Section 2.4), we have considered a restricted version of grammar (2.1), where products and quotients of data variables (e.g. x₁ × x₂ and x₁/x₂) are allowed only within the scope of the LAST operator (i.e. not with MIN, MAX, INT). Indeed, allowing products and quotients as arguments of path operators such as MAX or MIN requires the solution of a linear programming problem during the generation of a synchronised M×A path which, although feasible, would considerably affect the computation time.

**Example 5 (HASL expressions).** Referring to the automata $\mathcal{A}_1$ and $\mathcal{A}_2$ of Example 2.3 and $\mathcal{A}_{3a}$ and $\mathcal{A}_{3b}$ of Example 3 let us consider a few examples of complete HASL specifications together with their informal description.

• $\phi_{1a} \equiv (\mathcal{A}_1, \mathbb{E}(\mathrm{last}(x_1)))$: the mean number of class 1 processes waiting to access the
shared resource within the time interval $[0, 4]$

• $\phi_{1b} \equiv (\mathcal{A}_1, \mathbb{E}(\mathrm{max}(x_2)))$: the max number of class 1 processes that got access the
shared resource within the time interval $[0, 4]$

• $\phi_{2a} \equiv (\mathcal{A}_2, \mathbb{E}[\mathrm{LAST}(x_1)/k])$: the upper bound on the waiting time of the first $k$ instances of class 1 processes arrived in the system.

• $\phi_{2b} \equiv (\mathcal{A}_2, \mathbb{E}[\mathrm{LAST}(x_1)/\mathrm{LAST}(x_3)])$: the upper bound on the waiting time of the first $k$ instances of class 1 processes arrived in the system.

• $\phi_{3a} \equiv (\mathcal{A}_{3a}, \mathbb{E}[\mathrm{PROB}()])$: the probability that type 1 clients have used the shared resource longer than type 2's, given that $k$ clients (of any kind) have used the shared resource.

• $\phi_{3b} \equiv (\mathcal{A}_{3b}, \mathbb{E}[\mathrm{PROB}()])$: the probability that type 2 clients have used the shared resource longer than type 1's, given that $k$ clients (of any kind) have used the shared resource.
**Remark 4** (Existence of expectation $\mathbb{E}(Y)$). We emphasise that the (conditional) expectation of a path random variable is not always defined. There are two obvious necessary conditions on the synchronised product $(\mathcal{M} \times \mathcal{A})$: (1) almost surely the random execution ends either in a final state of the LHA or in the rejecting state, and (2) with positive probability the random execution ends in a final state of the LHA. However these conditions are not sufficient. Different restrictions on the path formula ensure the existence of expectations. For instance, when the formula only includes bounded data variables and the operator *INT* and the division are excluded, the expectation exists. Divisions may be allowed when the path expression is lower bounded by a positive value: for example w.r.t. formula $\phi_{2b}$, $LAST(x_3)$ is lower bounded by $LAST(x_2)$ which is lower bounded by *k*. The existence of regeneration points of the synchronised product may also entail the existence of such expectations. We do not detail the numerous possible sufficient conditions for the expectation of a path expression *Z* to exist but for all applications discussed is this manuscrit the existence of such expectation can be easily shown.

## 2.3 Expressiveness of HASL

As arguably the most relevant feature of the HASL formalism is its expressiveness it is important to address how it compares with respect to related formalisms. In the remainder we briefly overview temporal logic based frameworks for the verification of stochastic models, referring to [BBD$^{+}$15a] for a more in depth discussion.

### 2.3.1 Comparison with logics for numerical stochastic model checking

**CSL and variants.** A number of temporal logic formalisms have been introduced for model checking of CTMCs. These include the seminal Continuous Stochastic Logic (CSL), first introduced by Aziz *et al.* [ASSB00] and further extended by Baier *et al.* [BHHK03a]. CSL extended CTL *branching-time* state-based reachability reasoning to the realm of CTMCs, by 1) introducing densely time-bounded temporal operators, hence allowing for time-bounded reachability verification, and 2) by allowing for steady-state reasoning through the introduction of a dedicated steady-state modality. In an effort to extending model checking with performance analysis capabilities several variants of CSL have been introduced that operate on *reward enriched Markov models*, such as Baier *et al.* ’s CSRL [BHHK00], that extends the CSL approach to Markov reward models, i.e. CTMCs with a single reward on states or Cloth *et al.* ’s [CKKP05] that further extended that approach with impulse rewards on actions. Kwiatkowska *et al.* [KNP07b] then combined those advances by equipping the PRISM model checker [PRI] with complete state and transition reward-based analysis facilities which is achieved by giving the possibility to the modeler to enrich CTMCs with multiple state and/or transition rewards and to analyse them through dedicated *reward based operators* which allow for assessing their expected value either at given instant of time, at steady-state or cumulated over finite paths.
Another significant evolution of CSL is Baier et al. asCSL logic [BCH+07] through which the time-bounded Until of CSL is replaced by a regular expression with a time interval constraint. These path formulas can express elaborated functional requirements as in CTL* but the timing requirements are still limited to a single interval globally constraining the path execution.

**Remark 5 (CSL expressiveness).** Despite the considerable effort to evolve CSL logic the expressiveness of all CSL-like formalisms is limited to state-based, singly time-bounded reachability reasoning, that is: with CSL we can only consider state conditions and at most a single time interval as constraining factors to characterise the paths whose probability measure we wish to account for.

**Timed automata based CTMC model checking.** In order to increase the expressiveness of CSL-like temporal reasoning timed-automata based model checking of CTMCs has been proposed. Haddad et al. introduced the overall idea with the logic CSL<sup>TA</sup> [DHS09], whereby path formulae are defined by a single-clock deterministic timed automaton (DTA). This logic has been shown to be strictly more expressive than CSL and also more expres- sive than asCSL when restricted to path formulas. A generalisation of CSL<sup>TA</sup> has been introduced by Chen et al. [CHKM09], deterministic timed automata with multiple clocks are considered and the probability for random paths of a CTMC to satisfy a formula is shown to be the least solution of a system of integral equations. The cost of this more expressive model is both a jump in the complexity as it requires to solve a system of partial differential equations, and a loss in guaranty on the error bound.

**Remark 6 (Timed-automata temporal logic expressiveness).** The employment of (single or multiple-clock) timed-automata as paths specifiers increases the expressiveness of temporal reasoning as one can now not only combine state and event conditions but also consider multiple time-intervals as criteria to specify the paths of interest, therefore effectively moving beyond CSL expressiveness. However the proposed timed automata based model checking approaches do not provide any support for reward-based analysis as clocks being the only kind of variable supported by timed automata they do not allow to neither measure any reward nor consider cumulated rewards as criteria to characterise the paths of interest.

**M(I)TL.** Several logics based on linear temporal logic (LTL) have been introduced to consider timed properties, including Metric (Interval) Temporal logic in which the Until operator is equipped with a time interval. Chen et al. [CDM11] have designed procedures to approximately compute desired probabilities for time bounded verification, but with complexity issues. The question of stochastic model checking on (a sublogic of) M(I)TL properties, has also been tackled see e.g. [ZC13].
**Remark 7** (Limited expressiveness for numerical model checking). All of the above mentioned CSL and timed-automata logics have been designed so that numerical methods can be employed to assess the probability measure of a formula. This very constraint is at the basis of their limited expressive scope which has two aspects: first the targeted stochastic models are necessarily CTMCs; second the expressiveness of formulas is constrained by decidability/complexity issues. Furthermore the evolution of stochastic logics based on CTL seems to have followed two directions: one targeting *temporal reasoning* capability (in that respect the evolutionary pattern is: CSL $\rightarrow$ asCSL $\rightarrow$ CSL<sup>TA</sup> $\rightarrow$ multi-clock DTA), the other targeting *performance evaluation* capability (evolutionary path: CSL $\rightarrow$ CSRL $\rightarrow$ CSRL+impulse rewards). A unifying approach is currently not available, thus, for example, one can calculate the probability for a CTMC to satisfy a sophisticated temporal condition expressed with a DTA, but cannot, assess performance evaluation queries at the same time (i.e. with the same formalism).

**Hybrid automata based model checking.** As HASL is inherently based on stochastic simulation, it naturally allows for releasing the constraints imposed by logics that rely on numerical solution of stochastic models and therefore broaden the class of target stochastic models to DESP which includes, but is not limited to, CTMCs. From an expressiveness point of view, the use of LHA, i.e. a generalisation of DTA, allows for generic variables, which include, but are not limited to, clock variables. That opens up to a much wider expressiveness as by combining real-valued, boolean, integer-valued variables one can now take into account sophisticated reward measure that can not only be estimated *per se* but, most importantly, can be used as selection criteria of the paths of interest. This means that sophisticated temporal conditions as well as elaborate performance measures of a model can be accounted for in a single HASL formula, rendering HASL a unified framework suitable for both model-checking and performance and dependability studies. Note that the nature of the (real-valued) expression *Z* (2.1) generalises the common approach of stochastic model checking where the outcome of verification is (an approximation of) the mean value of a certain measure (with CSL, asCSL, CSL<sup>TA</sup> and DTA a measure of probability).

More specifically it is also worth noting that the use of data variables and extended updates in the LHA enables to compute costs/rewards naturally. The rewards can be both on locations and on actions. First using an appropriate flow in each location of the LHA, possibly depending on the current state of the DESP we get “state rewards”. Then, by considering the *update expressions* on the edges of the LHA, we can model sophisticated “action rewards” that can either be a constant, depend on the state of the DESP and/or depend on the values of the variables. Thus HASL extends the possibilities of CSRL (and extensions [LKKP05]). The extension does not only consist of the possibility to define multiple rewards (that can be handled, for example, through the reward-enriched version of CSL supported by the PRISM tool) but rather of their use. First several rewards can be used in the same formula, and last but not least these rewards have a more active role, as they can not only be evaluated at the end of the path, but they can also take an
important part in the selection of enabled transitions, hence of accepted paths. It is for
example possible and easy to characterise the set of paths along which a reward reaches a
given value and after that never goes below another value, a typical example that neither
PRISM-CSL nor CSRL can handle.

**Limitations of HASL.** Finally, we briefly discuss two features that are available in the above mentioned stochastic logics but not in HASL. First HASL does not allow to properly model nesting of probabilistic operators. The key reason is that this nesting is meaningful only when an identification can be made between a state of the probabilistic system and a configuration (comprising the current time and the next scheduled events). While this identification was natural for Markov chains, it is not possible with DESP and general distributions that have no memoryless property, and therefore this operation has not been considered in HASL. Furthermore, even for Markovian systems, the complexity of the statistical method on formulas with nesting is quite high [YS06] as the verification time per state along a path is no longer constant.

A similar problem arises for the steady state operator. The existence of a steady state distribution raises theoretical problems, except for finite Markov chains. With HASL we allow for not only infinite state systems but also non Markovian behaviours. However, when the DESP has a regeneration point, various steady state properties can be computed considering the sub-execution between regeneration points.

In conclusion it is worth noting that these limitations are rather due to the verification
method (statistical in our case) and to the expressiveness of the model (allowing non
Markovian systems) than to a particular tool. All this information given, particularly
concerning nested and steady state formulae, we can now state and prove our claim about
the respective expressiveness of HASL, CSRL and CSL<sup>TA</sup>:

**Proposition 1.** Given a non nested transient CSRL formula $P_{\succdot q}\phi$ and a system described as a Markov Reward Model, it is possible to build an *LHA* to estimate the probability $p$ for $\phi$ to hold, and then decide whether it fulfils the bound required (i.e. $p \succdot q$ with $\succdot \in \{\langle, \rangle, \leq, \geq\}$).

Proof. See [BBD<sup>+</sup>15a]

**Proposition 2.** Given a non nested transient CSL<sup>TA</sup> formula $P_{\succdot q}A(\phi_1, ..., \phi_n)$ and a system described as a Continuous Time Markov Chain, it is possible to build an LHA to estimate the probability $p$ for an execution to be accepted by the DTA $A(\phi_1, ..., \phi_n)$, and then decide whether it fulfils the bound required (i.e. $p \succdot q$ with $\succdot \in \{\langle, \rangle, \leq, \geq\}$).

Proof. See [BBD<sup>+</sup>15a]
### 2.3.2 A couple of meaningful applications

The expressive power of the HASL formalism has been demonstrated through a number of applications in different domains, including that of biological modelling [BGGM14, BD15], that of performance analysis of manufacturing systems [BDD$^+$11d, BH21a, BH21b], that of verification of wireless network protocols [BBV19, MMB$^+$17], in medicine [BKMP15] and in autonomous vehicle safety analysis [BBDH17] (a list of publications related to HASL is available at [COS]). To illustrate HASL expressiveness below we briefly overview a couple of meaningful applications.

#### Analysis of gene expression models with stochastic delayed dynamics

Gene expression is a fundamental biological intra-cellular process by which proteins are synthesised from a *gene*, i.e., a sequence in the DNA. It consists of two main phases: *transcription*, i.e., the copying of a sequence in the DNA strand by an RNA polymerase (RNAp) into an RNA molecule, followed by *translation*, i.e., the process by which proteins are synthesised from the (transcribed) RNA sequence. The rate of expression of a gene is usually regulated at the stage of transcription, by repressor molecules that can bind to the operator sites (generally located at the promoter region of the gene) and then inhibit transcription initiation. Evidence suggests that this is a highly stochastic process since usually, the number of molecules involved, e.g. transcription factors and promoter regions, is very small, ranging from one to a few at a given moment, which motivates the application of stochastic modelling to the analysis of gene expression dynamics. Ribeiro *et al.* [RZK06] argued that the dynamics of chemical reactions involved in the gene expression process exhibits a *stochastic delayed* nature, meaning that the appearance of a reaction's product may be subject to a further (non Markovian) stochastic delay after the actual occurrence of the reaction which, in turns, is selected stochastically in a standard manner, i.e., through Gillespie's SSA [Gil77].

In [BMR12] we formally analysed the dynamics of a single gene reaction network taken from [RZK06] whose stochastic Petri net encoding is depicted in Figure 2.6. The six reactions are modelled by small subnet blocks consisting of few timed Exponential timed transitions (depicted as empty rectangles) and non-Markovian timed transitions (grey filled rectangles). Place *Pro* represent the single promoter area of the gene and a token in it indicates the promoter is not occupied. Expression begins with the initialisation reaction $R_1$, whereby an RNAp binds to a promoter (*Pro*), which remains unavailable for more reactions until reaction *transc* occurs. Following reaction *transc* both the promoter and the ribosome binding site *RBS* become available although the production of an RNA is not terminated until a further delay given by (the Gamma distributed) reaction *termin* (notice that the number of RNA molecules produced corresponds with the marking of place *RNA*). As soon as the *RBS* site of the RNA is completed (through reaction *transc*) the *translation process* begins by binding (*reaction transl*) of a ribosome molecule *Rib* to the *RBS*. This process terminates only after further delays: a Gamma distributed timed
<figure><img src="image_6.png" /><figcaption>Figure 2.6: Stochastic Petri net model of the single-gene expression network</figcaption></figure>

transition representing the release of the *Rib* molecule from the *RBS*, another Gamma distributed timed transition modelling the termination of translation hence the synthesis of the product proteins *P* and finally a Dirac distribution for freeing the *RBS*. The model is completed by few further Markovian transitions: one modelling the RBS decay (*R*<sub>4</sub>) and a pair of transitions (*R*<sub>5</sub> and *R*<sub>6</sub>) modeling the periodic occupation of the promoter by a repressor molecule which inhibits the whole expression process.

<figure><img src="image_7.png" /><figcaption>Figure 2.7: Sustainment of translation during repression intervals.</figcaption></figure>

In order to analyse the dynamics of gene expression process we conceived a few prop-
erties which depend on non trivial random variables such as:

* $\phi_1$: *"duration of translation process given that translation ends within a repression interval"*.

* $\phi_2$: "probability of the number of translation events occurring in between two consecutive transcription events given at least *N* translations occur between two consecutive
transcriptions”.

Random variable $\phi_1$ quantifies the inertia (sustainment) of the translation phase in terms of how long translation keep going since expression is turned off (through repression). The overall idea is illustrated in Figure 2.7 (left) where a *repression interval* is the time interval between the occurrence of a repression event ($R_5$) and the corresponding un-repression event ($R_6$) and a “positive” repression interval is one in which the translation process actually completely terminated (i.e. by releasing a *P* molecule and detaching the *Rib* molecule from the *RBS*) whereas a “negative” interval is one in which translation did not terminate. Figure 2.7 (right) shows the average value of $\phi_1$ measured over a time-bounded horizon (through a dedicated HASL formula on COSMOS) in function of RBS decay rate.

<figure><img src="image_8.png" /><figcaption>Figure 2.8: Completed translations between consecutive transcriptions.</figcaption></figure>

On the other hand random variable $\phi_2$ is meant as an indicator of the relative frequency of the transcription and translation processes. Figure 2.8 (left) illustrates intuitively the meaning of measures related to $\phi_2$, where a “positive” interval is one in which at least $N = 2$ transcriptions event occurred as opposed to a “negative” interval is one in which less than $N$ transcriptions occurred. Figure 2.7 (right) shows the value of $\phi_2$ measured (through a dedicated HASL formula on COSMOS) in function of $N$ and w.r.t. to unrepressed and differently repressed setting of the single-gene model.

We stress that $\phi_1$ and $\phi_2$ are meaningful witnesses of the expressive power of the HASL formalism: an equivalent formulation in terms of popular temporal logic formalisms such as, e.g., the reward-extension of CSL supported by the PRISM model checker, could hardly be obtained.

## Analysis and computation of passage-time measures

In performance modelling the *passage time distribution* is a specific type of performance index particularly useful when reasoning about properties related with Service Level Agreements (SLA) or safety requirements. Generally speaking a passage time measure is con-
cerned with tracking the duration for a single instance of a process to complete, for example the time it takes for an order to be treated (archived) since it is placed on a selling platform (see model in Figure 2.9). In practice passage time measures can be formulated as the distribution of the time required to reach a *goal* state from any *entry* state possibly without hitting any *forbidden* or equivalently by matching the occurrence of a *start* with a corresponding *end* event. Notice that classical performance measures based on mean values, like the average response time, are not sufficient to estimate a passage time distribution, whose calculation requires dedicated methods and tools [Kul95, DHK04].

<figure><img src="image_9.png" /><figcaption>Figure 2.9: A NMGSPN model for an order-handling process.</figcaption></figure>

From a modelling standpoint measuring of passage time distribution induces some difficulties, as it requires the ability to track the start and end of a specific process instance. In Petri nets terms this corresponds to the ability of tracking a specific token (e.g. the one corresponding to a specific order) as it flows through the Petri net up until it reaches the place that correspond to the end of the monitored process instance<fnref n="6" />. That lead to the extension of standard Petri net formalisms as is the case with the Tagged extension of GSPN [BDPF09] (TGSPN) where tokens can be *tagged* in order to follow their traversal of the net. The evaluation of passage time measures on TGSPN models boils down to a *transient analysis* problem of the underlying CTMC.

Alternatively passage times can be accounted for by combining a “standard” model with an automata-based *monitor* where the latter is used as a machinery to track the beginning and end of the activity of interest. Examples of this approach are the Extended Stochastic Probes (XSP [CG08], operating on PEPA models [Hil05]), Path Automata (PA,

<footnote n="6">⁶Notice that tracking tokens in a Petri net is not trivial as soon as several tokens, representing parallel instances of the monitored process, are permitted.</footnote>
operating on Stochastic Activity Networks [OS99]) or Probe Automata (PrA [ABDF11]
operating on GSPN models, where PrA are a kind of deterministic, untimed automata).

In [ABB+13] we introduced the use of the HASL formalism for the specification of passage time measures over NMGSPN models and showed that HASL allows for specifying and assessing more sophisticated passage time measures than those supported by Probe Automata, more specifically conditional passage times measures obtained by a combination of time, state and events constraints.

For example, referring to the order-handling model of Figure 2.9, we introduced
the following HASL-expressed conditional passage time measures which we assessed (see
Figure 2.10) through a number of experiments run on the COSMOS tool (see Section 2.4):

* **Measure w1:** *the CDF of the passage-time for an ordered good to be delivered.*

* **Measure w2:** *the CDF of the passage-time for an ordered good to be delivered given that it was out-of-stock.*

* **Measure w3:** *the CDF of the passage-time for an ordered good to be delivered given that it is not out-of-stock and that the total delay for checking its availability and shipping it does not exceed K.*

* **Measure w4:** *the CDF of the passage-time for an ordered good to be delivered given that the total delay for reordering and updating the stock does not exceed K.*

<figure><img src="image_10.png" /><figcaption>Figure 2.10: Passage times for the workflow properties *w*₁ to *w*₄ assessed over time interval [0, 100].</figcaption></figure>

It is worth stressing that apart from the increased expressiveness in [ABB+13] we also
demonstrated the effectiveness of a general purpose performance analysis framework, such
as HASL model checking, in treating of a type of problem that otherwise is solved through customised formalisms. In fact assessment of passage time measures on a GSPN model through the Probe Automata approach is obtained through a quite cumbersome procedure which involves first constructing the tangible reachable graph of the GSPN, than composing it by considering the product with the Probe Automaton and finally solving the transient state distribution of the product model.

## 2.4 Software support: the COSMOS tool

In order to provide software support to the HASL formalism we developed COSMOS<fnref n="7" /> [BDD+11a, COS, BB22] a software platform for statistical model checking of HASL formulae. COSMOS is based on the SMC scheme outlined in Figure 2.2 hence it implements the stochastic simulation algorithm for the synchronised product process ($M \times A$) and supports different statistical methods for estimating the target measure $Z$ of a HASL formula $\phi \equiv (A, Z)$, including confidence interval estimation as well as hypothesis testing.

### 2.4.1 A model driven code generation architecture

COSMOS is implemented in C++ and relies on the Boost libraries for random number generation functionalities. In an effort to obtain a performant tool COSMOS is designed according to a *model driven code generation* scheme (Figure 2.11): the inputs $M$ (NMGSPN) and $A$ (LHA) are parsed in order to generate an efficient customised C++ implementation of the simulation engine of the synchronized product $M \times A$. The generated code is linked

<figure><img src="image_11.png" /><figcaption>Figure 2.11: COSMOS’s model-driven code generation scheme</figcaption></figure>

with a library containing parts of the simulator that are independent of $M$ and $A$. This library contains the main function that determines the next event to occur by means of an event heap and the generated code. In addition a pseudo random generator computes delays for the new events that are put into the heap.

<footnote n="7">⁷COSMOS is an acronym of the french sentence “Concept et Outils Statistiques pour le Modèles Stochastiques” whose english translation would sound like: “Tools and Concepts for Statistical analysis of Stochastic MOdels”.</footnote>
**Parallel threads of execution.** When an experiment is launched on a given model-automaton pair ($M, A$) COSMOS generate the C++ code for the corresponding SMC engine, compiles it and launches several threads of the resulting executable code in parallel that repeatedly generate trajectories and send back the evaluation of the formulas on these trajectories. COSMOS aggregates these evaluations and stops the simulation depending on the selected statistical method (see below). The code generation and compilation time is generally negligible compared to the simulation runtime.

### 2.4.2 COSMOS statistical engines

COSMOS is equipped with a number of statistical engines that are applied depending on the nature of the quantity to be estimated. Specifically *confidence interval* techniques, are used for *quantitative experiments*, i.e. experiments whose goal is to estimate the mean value of a given (real-valued) random variable, as opposed to *hypothesis testing* methods, which are applied for *qualitative experiments*, i.e. experiments whose goal is to establish how the mean value of random variable compares w.r.t. a given threshold.

**Confidence interval methods.** Confidence Interval (CI) methods works by collecting a number of samples $Z_1, \dots, Z_n$ of an unknown random variable, in order to produce an estimate of an interval which is likely to contain the exact mean value $\mu_Z$ of the considered random variable ($Z$ in the case of HASL). The accuracy of the estimate is controlled w.r.t. two complementary aspects: i) the *confidence level* ($\alpha \in (0, 1)$) which expresses how reliable the produced estimate is and ii) the *admitted error bound* (i.e. expressed as the width $\delta$ of the resulting interval). This means that if we repeatedly estimate the interval for a given $\theta$ we are guaranteed that the (possibly different) resulting intervals will contain the actual value of $Z$ in a proportion corresponding to $(1 - \alpha)$. Given an $n$-sized sample $Z_1, \dots, Z_n$ the general form of the $100(1 - \alpha)\%$ confidence-interval for the expected value $\mu_Z$ of $Z$, denoted $CI_{\mu_Z}^{\alpha}$, is:

$$CI_{\mu_Z}^{\alpha} = (PE_{\mu_Z}) \pm EB_Z$$

where $PE_{\mu_Z}$ is a Point Estimate of $\mu_Z$ and $EB_\alpha$ is the Error Bound, which corresponds to the semi-width of the CI interval, i.e. $EB_\alpha = \delta/2$. The sample mean $\bar{Z} = \frac{\sum_{i=1}^n Z_i}{n}$ is used as (unbiased) Point Estimator of $\mu_Z$. The execution of CI experiment depends on three parameters: the sample size $n$, the confidence level $\alpha$ and the error bound $\delta$. Generally speaking one can fix two of them and get the corresponding value of the remaining one in function of the previous two, however the actual choice of the two parameters to be fixed depend on the nature of variable to be estimated: in some case the sample size must be fixed a priori.

**Static sample size estimation.** With respect to CI estimation COSMOS distinguish between i) *indicator variables* (Bernoulli variables) used for evaluation of unknown prob-
abilities, ii) *bounded variables* used for evaluation of proportions, ratios, mean number of clients in a system with finite capacity and iii) *general variables* without any additional knowledge. For *indicator variables* and *bounded variables* COSMOS allows the user to choose 2 out of 3 of the CI parameters and automatically establishes the value of the remaining one in function of the fixed two parameters. For *general variables* however the sample size must be provided leaving free choice of one among the two remaining parameters. In general when the user fixes the sample size *n* and the confidence level *α* COSMOS compute the corresponding error bound based on the following statistical methods:

* Indicator variables: the error bound is obtained through the Clopper-Pearson method [CP34].

* Bounded variables the error bound is obtained through the Chernoff-Hoeffding method [Hoe63].

* General variables: an asymptotically correct error bound is obtained as approximation of a normal distribution with unknown mean and variance.

**Dynamic sample size estimation.** On the other hand, when estimating a general variable, the user may choose a *priori* the confidence *α* and the error bound *δ* and in this case cases COSMOS applies *sequential CI scheme* through which the number of samples (i.e. trajectories) is established at runtime depending on a stopping condition (normally the convergence of the sample standard deviation). More specifically COSMOS supports the following CI sequential schemes: Chow and Robbins [CR65] Clopper-Pearson [CP34] and Chernoff-Hoeffding [Hoe63].

**Hypothesis testing.** In case the user is interested in a qualitative formula, i.e. in establishing whether the average value of a given quantity *Z* fulfils a given constraint then COSMOS adopt a hypothesis testing approach. Specifically if the considered variable is a Bernoulli of unknown mean *p*, and given two probabilities *p*₀ < *p*₁ (the maximal probability of false positive and true negative results), COSMOS applies the Sequential Probability Ratio Test (SPRT) [Wal45] method (i.e. an optimal sequential test method) for deciding whether *p* ≥ *p*₀ or *p* ≤ *p*₁ holds.

## Related Tools

Numerous tools that support statistical model checking have been developed over the years, some of which also support numerical model checking engines. These include COSMOS [BDD^+11c], PLASMA [JLS12], PRISM [KNP11], UPPAAL [BDL^+12], MARCIE [HRS13], STORM [HJK^+22], APMC [HLP06], YMER [You05], MRMC [KZH^+09] and VESTA [SVA05b]. An in depth comparison of these tools would certainly be an interesting task in which one would need to take into account several aspects including the
supported family of models, the solution engines, the expressive power of the property
language but also how they perform w.r.t. some benchmark case studies. In an attempt
to tackle such a goal in [BBD$^+$15c] we provided a tentative comparative analysis limited
to the YMER, PRISM, UPPAAL, PLASMA, MARCIE and COSMOS tools.

## 2.5 Perspectives

The HASL model checking framework effectively widen the ability of modellers to analyse the systems under study. By employing a rich automata-based formalism as a means to monitor a model’s executions one may conceive sophisticated properties (by freely combining state-conditions, event-conditions and reward-based conditions) that allow one to capture complex aspects of a system’s dynamics that could not be accounted for with “classical” temporal logic formalisms. A further positive consequence of such expressive power is that it yields a “full” *separation of concerns* meaning that modelling and properties conception are independent efforts in the sense that a model does not need to be enriched with complementary elements (rewards and/or dedicated monitor variables) which are otherwise necessary to overcome the deficits of less expressive formalisms and whose goal would be to register relevant statistics that one wants to account for during model checking.

If the expressive power is undoubtedly the strongest point in favour of the HASL formalism there are a number of factors worth considering for improving things. First the unsurprising cost of such an expressive power is the *unfriendliness* of the formalism: indeed the design of a hybrid automaton to capture a given temporal dynamics is a delicate and error prone task requiring expertise. In this respect it would be worth to come up with a workaround, in terms of a user friendly (syntax based) language to express a target property of the system under study (e.g. “what is the mean period of oscillation of species A?”) which is suitable to automatic translation into a corresponding (HASL) hybrid automaton monitor. The conversion of temporal logic formulae into a corresponding (timed) automata formalism has been discussed by Donatelli, Sproston and Haddad [DHS07] as well as by ourselves [BBD$^+$15a]. It is known that commonly used not-nested, (bounded/unbounded) reachability properties based on the *Until* temporal operator (and its variants) can straightforwardly be encoded in HASL terms while encoding of formulae involving nested temporal operators is more cumbersome. Generally speaking an in depth treatment of this matter is missing. It would therefore be worth figuring out to what extent one can outline an approach to automatically generate hybrid automaton monitors that represents relevant properties expressed through a more user friendly syntax.

Another aspect worth considering in order to evolve the HASL framework concerns the speed up of the statistical model checking engine. Although, differently from the memory constrained numerical model checking scheme, SMC allows for taking into account very large (even infinite state) models, the runtime of a SMC experiment can be huge given that sampling of a model’s path through *exact* stochastic simulation can be very computational
intensive (depending on the stochasticity as well as of the size of the model). In recent
time I got interested in investigating whether, within a SMC framework, one could take
advantage of accelerated, *approximated* stochastic simulation algorithms, e.g., those based
on the τ-leaping scheme [Gil01]. Roughly speaking that entails investigating what is the
effect of the error induced by *approximated* paths on the (confidence-interval) estimation of
the mean value of the property φ which is targeted by the SMC experiments. A preliminary
prototype implementation developed as part of a student’s research project seems to provide
promising results in this respect.
# Chapter 3

## A Bayesian approach to parametric model checking

The automated verification of formally expressed properties by model checking have proved a successful achievement in theoretical computer science as demonstrated by numerous applications in heterogenous domains. One limitation of model checking though is that it does not come with natural means to analyse how the satisfaction of the considered property $\phi$ is affected by the model's parameters. Given that a system's model $\mathcal{M}$ normally depends on a set of parameters<fnref n="1" /> $\theta = (\theta_1, ..., \theta_n) \in \mathbb{R}^n$ what model checking does is to establish whether a given model's instance $\mathcal{M}_\theta$ (i.e. the instance of $\mathcal{M}$ obtained by replacing the parameters with the values $\theta \in \mathbb{R}^n$) complies with $\phi$, denoted $\mathcal{M}_\theta \models \phi$. A trivial approach to study how the satisfaction of $\phi$ is affected by parameters $\theta \in \mathbb{R}^n$ would be to systematically run a model checking experiment for every possible instance $\mathcal{M}_\theta$. Given the parameter space is normally dense such a *brute force* exhaustive approach would hardly be feasible other than being very inefficient given the computational cost for iterating single model checking experiments $\mathcal{M}_\theta \models \phi, \forall \theta \in \mathbb{R}^n$.

Here we are interested in *parametric model checking* of continuous-time probabilistic systems (particularly Markov chains) and therefore parameters correspond with the rates of the probability distributions that characterise the random delays of state transitions. The relevant question we wish to answer is: given a stochastic model with parametric rates for which sets of rate values a given temporal logic specification is satisfied? There are two families of approaches to answer such question: those that rely on *numerical methods* as opposed to those that rely on *statistical methods*.

**Numerical approaches.** Numerical approaches for model checking parametric CTMCs

<footnote n="1">¹Example of which may be, the rate at which clients arrive in a given queue of a queueing system, or the maximum size of a mailbox in a message receiving server, or the probability at which a production machine in a manufacturing system breaks down in the next time unit, or yet the rate at which a given protein is synthesised when its corresponding gene is not inhibited.</footnote>
rely on basic findings in CTMC model checking [BHHK03a], that is: the verification of time-
bounded CSL specifications against a CTMC instance can be reduced to a transient-state
distribution problem whose solution, in turn, is approximated through uniformisation. To
deal with parametric CTMC a parametrised version of uniformisation has been proposed
leading to algorithms for parameter synthesis. Initially limited to qualitative non-nested
CSL specifications [HKM08], (i.e. boolean valued formulae of the form $P_{\sim p}[\phi]$ where $\phi$ is a
non-nested path formulae and $\sim \in \{<, \le, \ge, >\}, p \in [0, 1]$) they have then been extended to
the entire CSL spectrum including quantitative formulae (i.e. formulae of the form $P_{=?}[\phi]$.
whose verification requires assessing the probability that condition $\phi$ is met) [CDP+17].
The positive aspects of parametrised-uniformisation synthesis schemes is that they im-
plement an exhaustive search of the parameter space yielding regions in which either a
qualitative CSL formula is guaranteed to be satisfied or, alternatively, regions in which
the probability of a quantitative CSL formula is ensured to fall within a bounded interval
(of arbitrary width). The downside however is their computational cost<fnref n="2" /> and hence their
inability to scale up other than their limitation to the Markov chains realm.

**Statistical approaches.** Statistical approaches to parametric verification constitute a computationally cheaper alternative to numerical ones, the basic idea being to obtain an approximation of the *satisfaction probability* of *φ* through non-exhaustive exploration of the parameter space. The prominent contribution in this still relatively young area of research is that of smoothed model checking [BMS16] in which Bortolussi and Sanguinetti demonstrated that the satisfaction function of a Metric Interval Temporal Logic (MITL) formula against a parametric CTMC is (under mild conditions over the nature of the parameters) a smooth function of the CTMC parameters. This leads to algorithms where an arbitrarily precise approximation of the satisfaction function is derived via a Gaussian Process based on its smoothness.

**Outline of contribution.** In this chapter we present a statistical approach for approxi-
mating the *satisfaction probability* of temporal properties in function of a stochastic model's
parameters. Differently from smoothed model checking our approach relies on Bayesian
statistics, specifically on the Approximated Bayesian Computation (ABC) method, i.e.,
a well established method that given a set of *observations* of the modelled system and
an initial *prior* distribution over the model's parameter space allows to approximate the
*posterior* distribution by selecting parameters based on their *their distance* from the sys-
tem's *observations*. Here we adapt ABC schemes by replacing the *observations* with a
formal *machinery* (consisting in a linear hybrid automata of the HASL formalism) through
which we assess *how far a set of parameters is from satisfying the considered temporal
property*, hence letting ABC algorithms select only those parameters that match the de-
sired behaviour. Since we use HASL linear hybrid automata as machinery to assess the
distance of parameters we name Automaton-ABCthe parametric model checking adapta-

<footnote n="2">²in the general case parametrised uniformisation boils down to solving a nonlinear programming problem.</footnote>
tion of ABC. Differently from standard ABC schemes, the automaton-ABC procedure is *observations free*, in that it does not require any dataset in order to drive the parameter selection process, however it relies on the existence of a notion of *distance* between the paths of a model instance $M_\theta$ and the property $\phi$ which is the target of the parameter search. We demonstrate the effectiveness of the automaton-ABC scheme by introducing a formal definition of distance for a classical time-bounded reachability problems and we give the corresponding distance-monitor LHA automata. This contribution has been at the core of Mahmoud Bentriou PhD which I co-directed with Paul-Henry Cournède and whose results have been published in [BBC19, BBC21].

**Structure of the chapter.** The chapter is organised as follows: Section 3.1 gives an overview of necessary background material including the class of parametric Markov models which we consider, the property language we refer to, the classical ABC algorithms which we aim to extend and finally the Kernel estimation method which we resort to for reconstructing the satisfaction probability of a property $\varphi$ from the posterior distribution issued by the Automaton-ABC algorithm. In Section 3.2 we introduce the notion of satisfiability distance for reachability problems and introduce the novel Automaton-ABC framework based on such distance measure is developed while in Section 3.3 the effectiveness of the framework is demonstrated through a couple of case studies.

# 3.1 Preliminaries

## 3.1.1 Markov Population Models

Before looking at the basics of how the ABC method works we start by overviewing the class of models we are going to consider, namely that of (parametric) Markov population models (pMPM), a (parametrised) form of continuous-time Markov chain (CTMC) [Kul16, BC19] suitable for modelling of *population processes*, i.e. systems whose states represent the number of individuals of different *species* and whose transitions correspond to adding/re-moval of individuals. Since we target modelling of biological systems we also overview the Chemical Reaction Networks (CRNs) formalism as a means for expressing population models.

**Definition 3.1** (Markov Population Model). A Markov population model (MPM) for $n \in \mathbb{N}$ population species is a triple $\mathcal{M} = (S, Q, \pi_0)$ such that:

* $S \subseteq \mathbb{N}^n$ is a countable set of states

* $\mathbf{Q}: S \times S \rightarrow \mathbb{R}$: is the infinitesimal generator matrix (with $\mathbf{Q}(s_i, s_j) = -\sum_{j \neq i} \mathbf{Q}(s_i, s_j)$)

* $\pi_0: S \rightarrow [0, 1]$ is the initial state probability distribution<fnref n="3" /> (i.e. $\sum_{s \in S} \pi_0(s) = 1$).

<footnote n="3">³Whenever $\exists s_o \in S : \pi_0(s_o) = 1$ we use $\mathcal{M} = (S, Q, s_o)$ to denote an MPM.</footnote>
Because of the memoryless property Markov models allows for a compact analytical expression of the probability for timed state transitions: given that the system is in state *s* at time *t*, the probability of observing a transition to state *s'* within time *t'* is $Pr((s, t) \to (s', t')) = P(s, s') \cdot (1 - e^{-Q(s, s') \cdot (t'-t)})$ where $E(s) = \sum_{s' \neq s} Q(s, s')$ is the *exit rate* of state *s*, and $P(s, s') = Q(s, s')/E(s)$ is the (time-independent) probability of jumping from *s* to *s'*. Such basic property induces naturally a probability measure over the space of paths of a MPM model.

Whenever matrix *Q* is dependent on a *d*-dimensional vector of parameters $\theta \in \Theta \subseteq \mathbb{R}^d$
we talk of *parametric* MPM (pMPM).

**Paths of a MPM.** Since our final goal is to be able to estimate $Pr(\varphi|\mathcal{M}_{\theta})$, i.e. the probability that a temporal property $\varphi$ is satisfied in function of the parameters of a parametric MPM $\mathcal{M}_{\theta}$, we first recall that $Pr(\varphi|\mathcal{M}_{\theta})$ corresponds with the probability of those paths of $\mathcal{M}_{\theta}$ that satisfy $\varphi$ and then revise how, thanks to the memoryless property of MPMs, such probability enjoys a nice analytical expression leading to a compact a measure of probability<fnref n="4" /> on the set of paths [BHHK03b].

A path of an MPM model $\mathcal{M}_\theta$ is a (possibly infinite) sequence $\sigma = s_0 \xrightarrow{t_0} s_1 \xrightarrow{t_1} \dots \xrightarrow{t_{n-1}} s_n \dots$, with $t_i \in \mathbb{R}_{>0}$ being the sojourn-time in state $s_i \in S$. We denote $\text{Path}_{\mathcal{M}_\theta}$ the set of all possible paths of $\mathcal{M}_\theta$ originating in the initial state $s_0$. For $\sigma \in \text{Path}_{\mathcal{M}_\theta}$, $i \in \mathbb{N}$ and $t \in \mathbb{R}_{>0}$, we denote $\sigma[i] = s_i$ the $i$-th state of $\sigma$, $\delta(\sigma, i) = t_i$ the sojourn-time of $\sigma$ in the $i$-th state, $\sigma_i$ the suffix of $\sigma$ starting at state $\sigma[i]$, $T_k = \sum_{i=0}^k \delta(\sigma, i)$ the sum of the sojourn times up to state $k$, $\sigma@t$ the state of $\sigma$ at time $t$ and $\sigma[t]$ the $t$-shifted suffix of $\sigma$, i.e. the suffix of $\sigma$ that starts at time $t$. Formally $\sigma[t] = \sigma[k+1] \xrightarrow{t-T_{k+1}} \sigma_k$ where $k$ is the greatest index such that $T_k \le t$. For example, for $\sigma = s_0 \xrightarrow{0.25} s_1 \xrightarrow{0.5} s_2 \xrightarrow{0.15} s_3 \xrightarrow{1} \dots$ we have $\sigma[1] = s_1$, $\delta(\sigma, 2) = 0.15$, $T_1 = 0.75$, $T_2 = 0.9$ and $\sigma[0.8] = s_2 \xrightarrow{0.1} s_3 \xrightarrow{1} \dots$ and $\sigma[1.5] = s_3 \xrightarrow{0.4} \dots$. Notice that trajectories of an MPM are *càdlàg* (i.e. step) functions of time. In order to indicate the reaction event that yielded a transition of the path of a CRN model, we sometimes adopt the following notation $\sigma = s_0 \xrightarrow{0.25} s_1 \xrightarrow{0.5} s_2 \xrightarrow{0.15} s_3 \xrightarrow{1} \dots$, where $R_{ij}$ indicates that reaction $R_j$ occurred on the $i$-th transition of the path.

**Paths probability space.** For $s_0, s_1, \dots, s_k$ a sequence of states of $\mathcal{M}_\theta$ such that $\mathcal{Q}(s_i, s_{i+1}) > 0 \ (0 \le i < k)$ and $I_0, \dots, I_{k-1}$ a sequence of non-empty time intervals in $\mathbb{R}_{\ge 0}$, we let $C(s_0, I_0, s_1, \dots, I_{k-1}, s_k)$ be the *cylinder set* consisting of all paths $\sigma \in \text{Path}_{\mathcal{M}}(\bar{\mathcal{M}})$ such that $\sigma[i] = s_i \ (i \le k)$ and $\delta(\sigma, i) \in I_i \ (i < k)$. Furthermore we let $\mathcal{F}(\text{Path}_{\mathcal{M}})$ denote the smallest $\sigma$ algebra containing all sets $C(s_0, I_0, s_1, \dots, I_{k-1}, s_k)$. The probability measure on

<footnote n="4">⁴Notice that the notion of probability measure for paths of an MPM naturally extends to parametric MPMs which we target in the remainder.</footnote>
$\mathcal{F}(Path_M)$ is inductively defined by:

$$ \Pr_M(C(s_0, I_0, \dots, I_{k-1}, s_k)) = \begin{cases} 1, & \text{if } k=0 \\ \mathbf{P}(s_{k-1}, s_k) \cdot (e^{-E(s_{k-1}) \cdot t} - e^{-E(s_{k-1}) \cdot t'}), & \text{Pr}_M(C(s_0, I_0, \dots, I_{k-2}, s_{k-1})), \text{ otherwise} \end{cases} \quad (3.1) $$

where $t = inf(I_{k-1})$ and $t' = sup(I_{k-1})$. In essence $\Pr_M$ states that for a MPM $M$ the probability of a path is given by the product of the probability to observe each constituent transitions $s_i \to s_{i+1}$ with a delay that falls in the corresponding binding interval $I_i$. In terms of vocabulary a subset of trajectories of $Path_M(M)$ may be referred to as an event of $M$ and its probability is given by $\Pr_M$. In the realm of probabilistic model checking, temporal logic languages provide the modeller with a powerful language for characterising relevant events of an MPM model in terms of *formulae* (i.e. formal properties). The probability that a temporal logic property $\varphi$ is satisfied by an MPM model $M$ is defined in terms of the probability measure $\Pr_M$ (see Definition 3.5).

In the remainder we consider Chemical Reaction Networks (CRNs) as a formalism for expressing population models.

**Definition 3.2 (Chemical Reaction Network).** A (parametric) chemical reaction network (pCRN) with $n$ species and $m$ reaction channels is a triple $N_\theta = (\mathcal{X}_n, \mathcal{R}_m, \theta, \mathbf{X}^0)$ defined as follows:

* $\mathcal{X}_n = \{X_1, \dots, X_n\}$ is a set of species

* $\mathcal{R}_m = \{R_1, \dots, R_m\}$ is a set of reaction channels where each $R_j$ ($j \in \{1, \dots, m\}$) is characterised by an equation with the following form:

$$ R_j : \sum_{i=1}^{n} a_{ij}^{-} X_i \xrightarrow{k_j} \sum_{i=1}^{n} a_{ij}^{+} X_i $$

where $a_{ij}^{-}, a_{ij}^{+} \in \mathbb{N}$ are the stoichiometric coefficients of the reaction's *reactants*, respectively, *products*. Furthermore $R_j$ is characterised by a pair $R_j : (\nu_j, \eta_j)$ with

$$ - \nu_j = [\nu_{1j}, \dots, \nu_{nj}] \text{ the change vector,} $$

$$ - \eta_j : \mathbb{N}^n \times \Theta \to \mathbb{R}_{\ge 0} \text{ is the propensity function of } R_j. $$

* $\theta = [\theta_1, \dots, \theta_d]$ is a $d$-dimensional vector of parameters affecting the kinetic rate of the reaction channels, with $\theta \in \Theta \subset \mathbb{R}^d$.

* $\mathbf{X}^0 \in \mathbb{N}^n$ is the initial state
Although in the literature CRNs are often inherently mapped on their continuous-deterministic semantics (i.e., a system of differential equations) here we focus on their discrete-stochastic semantics, hence we assume the dynamics of a CRN to correspond with a (càdlàg) step-function governed by its reactions channels $R_j$. Assuming the system is in state $\mathbf{X} \in \mathbb{N}^n$ at time $t \in \mathbb{R}_{\ge 0}$ reaction $R_j : (\nu_j, \eta_j)$ may occur, moving the system to state $\mathbf{X}' = \mathbf{X} + \nu_j$, at time $t' > t$, with the delay $t' - t$ which is stochastically dependent on both the current state $\mathbf{X}$ and the actual value of the parameters $\theta$.

**Remark.** For the sake of simplicity in our framework we assume reactions to obey the mass-action law. That means that the propensity functions are proportional to the product of the non-null stoichiometric coefficients of a reaction's reactants.

**Definition 3.3** (pMPM semantics of a CRN). A pMPM model $\mathcal{M}_{\theta} = (S, Q_{\theta}, \pi_0)$ of a CRN $\mathcal{N}_{\theta} = (\mathcal{X}_n, \mathcal{R}_m, \theta)$ is defined as follows:

* $S \subseteq \mathbb{N}^N$ is a countable set of states whose elements are vectors $\mathbf{X} = [X_1, \dots, X_n] \in S$ where $X_i$ is the population of the $i$-th species.

* $\mathbf{Q}_{\theta}: S \times S \to \mathbb{R}$: is the infinitesimal generator matrix whose entries are defined as:

$$ \mathbf{Q}_{\theta}(\mathbf{X}, \mathbf{Y}) = \begin{cases} \sum_{\{R_j | \mathbf{X} + \nu_j = \mathbf{Y}\}} \eta_j(\mathbf{X}, \theta), & \text{if } \mathbf{X} \neq \mathbf{Y} \\ -\sum_{\mathbf{Z} \neq \mathbf{X}} \mathbf{Q}(\mathbf{X}, \mathbf{Z}), & \text{otherwise} \end{cases} $$

* $\pi_0: S \to [0, 1]$ is defined as $\pi_0(\mathbf{X}^0) = 1$

Notice that by construction, the non-diagonal entries of $\mathbf{Q}_{\theta}$ are given by the sum of the propensities of those reactions whose occurrence leads the CRN to move from state $\mathbf{X}$ to state $\mathbf{Y}$. This is in line with the semantics of Markovian events, according to which the distribution of the minimum between a set of concurrent exponentially distributed reactions is itself exponentially distributed with rate given by the sum of the rates of the racing events.

**Example 6** (CRN of infection spreading). As a first example of CRN let us consider the SIR compartmental model [KM27], which describes the spread of infectious disease among a constant population. The CRN for the SIR is defined as $\mathcal{N}_{\text{SIR}} = (\{S, I, R\}, \{R_1, R_2\}, \{k_i, k_r\})$ where species S represents the susceptible individuals, I the infected and R the recovered ones. The system's dynamics is given by two reactions channels encoded by chemical equations (3.2).

$$ R_1: S + I \xrightarrow{k_i} 2I \qquad R_2: I \xrightarrow{k_r} R \qquad (3.2) $$

Reaction $R_1$ describes the infection step: a susceptible person meets an infected person and gets infected. Reaction $R_2$ models the recovering step: infected may become immune from the disease. The parameter vector of the model is $\theta = (k_i, k_r)$. The CRN of the
SIR yields a finite-state MPM with the following kinds of state-dependent transitions. For $\mathbf{X} = (\mathbf{X}_S, \mathbf{X}_I, \mathbf{X}_R)$ a state such that $\mathbf{X}_S > 0 \wedge \mathbf{X}_I > 0$ two kinds of transitions are possible, i.e. $\mathbf{Q}_\theta((\mathbf{X}_S, \mathbf{X}_I, \mathbf{X}_R), (\mathbf{X}_S-1, \mathbf{X}_I+1, \mathbf{X}_R)) = \mathbf{X}_S \cdot \mathbf{X}_I \cdot k_i$ and $\mathbf{Q}_\theta((\mathbf{X}_S, \mathbf{X}_I, \mathbf{X}_R), (\mathbf{X}_S, \mathbf{X}_I-1, \mathbf{X}_R+1)) = \mathbf{X}_I \cdot k_r$. For states such that $\mathbf{X}_S=0 \wedge \mathbf{X}_I > 0$ only one transition is possible i.e. $\mathbf{Q}_\theta((\mathbf{X}_S, \mathbf{X}_I, \mathbf{X}_R), (\mathbf{X}_S, \mathbf{X}_I-1, \mathbf{X}_R+1)) = \mathbf{X}_I \cdot k_r$, whereas any state such that $\mathbf{X}_I=0$ is absorbing.

**Regions of an MPM.** In the remainder, we refer to the notion of region associated with an MPM. A region, respectively a time-bounded region, of an n-dimensional MPM is any subset of $\mathbb{N}^n$, respectively $\mathbb{N}^n \times \mathbb{R}_{\ge 0}$, characterised by a collection of hyper-rectangles of dimension no larger than n. A region is elementary if it is characterised by a single hyper-rectangle. For example for a bi-dimensional MPM with state space $\mathbf{X} = \{X_1, X_2\}$, $R_1 \equiv [[1, 2]] \times \mathbb{N}$ is an elementary region ($X_1$ in $[[1, 2]]$ while $X_2$ is unbounded), $R_2 \equiv [[0, 3]] \cup [[5, 8]] \times [[5, \infty[ is a non-elementary region ($X_1$ is either in $[[0, 3]]$ or $[[5, 8]]$, $X_2$ is larger than 5), whereas $TR_1 \equiv ([1, 2] \times \mathbb{N}) \times [0.2, 1.41]$ is an elementary time-bounded region (similar to $R_1$, but with the supplemental condition that the time is in $[0.2, 1.41]$).

### 3.1.2 Linear-time temporal properties

To express properties of a MPM we refer to the Metric Interval Temporal Logic (MITL [MN04]), i.e. a linear-time temporal logic which allows for stating time-bounded reachability problems by combining state-conditions (expressed through inequalities on state variables) through classic time-bounded temporal operators. We choose MITL in order to formally set up reachability problems, however since we eventually employ hybrid automata as *meters* to measure the *satisfiability distance* of a MPM instance from a MITL formula, in practice this gives us the possibility to trespass the MITL expressiveness (e.g. by analysing oscillations properties).

**MITL syntax.** MITL formulae are terms of the following grammar:

$$ \varphi ::= \top | \mu | \neg \varphi | \varphi_1 \wedge \varphi_2 | \varphi_1 \mathbf{U}^{[t_1, t_2]} \varphi_2 $$

where $\top$ stands for the true formula, $\mu$ denotes an atomic proposition (i.e. an inequality built on top of model's state-variables), $\neg$ and $\wedge$ are the basic negation and conjunction connectives of propositional logic and $\mathbf{U}^{[t_1, t_2]}$ is the time-bounded *until* temporal operator with $[t_1, t_2] \subseteq \mathbb{R}_{\ge 0}$ being the bounding interval.

We consider the truth of an MITL formula to be established w.r.t. to a path $\sigma$ of an MPM model through a so-called *satisfaction relationship*, denoted $|=$.

**MITL semantics for temporal formulae.** For $\sigma \in \text{Path}_{M_\theta}$ a path of MPM model $M_\theta$, $t \in \mathbb{R}_{\ge 0}$ a time instant the satisfaction relation $|=$ of MITL temporal formula is defined as follows:
$$
\begin{align*}
(\sigma, t) &\models \top &\Leftrightarrow& \text{true} \\
(\sigma, t) &\models \mu &\Leftrightarrow& \sigma@t \models \mu \\
(\sigma, t) &\models \neg\varphi &\Leftrightarrow& (\sigma, t) \not\models \varphi \\
(\sigma, t) &\models \varphi_1 \wedge \varphi_2 &\Leftrightarrow& (\sigma, t) \models \varphi_1 \text{ and } (\sigma[t, t) \models \varphi_2 \\
(\sigma, t) &\models \varphi_1 \mathbf{U}^{[t_1, t_2]}\varphi_2 &\Leftrightarrow& \exists t' \in [t+t_1, t+t_2] : \\
&&& (\sigma, t') \models \varphi_2 \wedge \forall t'' \in [t, t'], \\
&&& (\sigma, t'') \models \varphi_1
\end{align*}
$$

Intuitively MITL semantics states that an atomic proposition $\mu$ is satisfied by a path $\sigma$
as of time $t$ if the state condition $\mu$ is satisfied in the state in which $\sigma$ is at time $t$ ($\sigma@t \models \mu$).
On the other hand, a time-bounded until formula $\varphi_1 \mathbf{U}^{[t_1,t_2]} \varphi_2$ is satisfied by $\sigma$ as of time
$t$ if and only if $\varphi_2$ is satisfied by $\sigma$ as of a future time instant $t'$ which is no further than the
time-bounding interval, (i.e. $t' \in [t+t_1, t+t_2]$) while $\varphi_1$ is sustainedly satisfied beforehand
(i.e. $\forall t'' \in [t, t']$).

As usual, we consider two derivations of the time-bounded until operator: the time-
bounded *eventuality* $\mathbf{F}^{[t_1,t_2]}\varphi \equiv \top \ \mathbf{U}^{[t_1,t_2]}\varphi$, which stands for “at some point within $[t_1,t_2]$ $\varphi$ is satisfied” and the time-bounded *globally* $\mathbf{G}^{[t_1,t_2]}\varphi \equiv \neg\mathbf{F}^{[t_1,t_2]}\neg\varphi$ which stands for “$\varphi$ is always satisfied within $[t_1,t_2]$”. In the remainder we assume that a path $\sigma \in \text{Path}_{\mathcal{M}_{\theta}}$ satisfies an MITL formula $\varphi$, denoted $\sigma \models \varphi$, if it does so starting from $t=0$, i.e. $\sigma \models \varphi \iff \sigma[0] \models \varphi$. Furthermore, we restrict our focus to the non-nested fragment of MITL, i.e., we consider only formulae such that the operands of a temporal modality are boolean combinations of atomic propositions $\mu$. While bearing a definite limitation in terms of expressiveness, this constraint still allows us to treat the most common reachability problems.

**MITL formulae and MPM regions.** In the remainder we refer to the notion of untimed/time-bounded reachability regions associated with a MITL propositional formula. Specifically we point out that MITL propositional formulae induce *untimed* regions over the state space of an MPM. For example, for a bi-dimensional MPM the formula $\mu_1 \equiv x_1 \ge 1 \wedge x_1 \le 2$ induces the region $R_1 \equiv [[1,2]] \times \mathbb{N}$ while $\mu_2 \equiv [(x_1 \le 3) \lor (x_1 \ge 5 \land x_1 \le 8)] \land x_2 > 4$ induces the region $R_2 \equiv [[0,3]] \cup [5,8] \times [[5,\infty[$. By a slight abuse of vocabulary, we say that two formulae $\mu_1, \mu_2$ are *disjoint* if the corresponding regions are. In the remainder, we assume regions are characterised by MITL propositional formulae in disjunctive normal form (DNF).

On the other hand we may associate a time-bounding interval to a MITL propositional
formula and in this case we talk of time-bounded reachability region.

**Definition 3.4 (Time-bounded reachability region).** Given an *n*-dimensional MPM population model $\mathcal{M}_\theta$ with state-space $S \subseteq \mathbb{N}^n$, a propositional formula $\mu$ and a time-interval $I \subseteq \mathbb{R}_{\ge 0}$ we define $S_\mu^I \subset S \times \mathbb{R}_{\ge 0}$ the time-bounded reachability as:

$$
S_{\mu}^{I} = \{(s, t) \in S \times \mathbb{R}_{\ge 0} \mid s \models \mu \wedge t \in I\}
$$

If $I = [0, \infty)$ then we omit the time interval and use $S_{\mu} = \bigcup_{i}^{n} \{s \in S \mid s \models \mu\}$
Note that depending on $\mu$, we distinguish between *elementary reachability regions*, i.e. regions that consist of a single hyperrectangle (e.g. $\mu \equiv x_1 > 2 \land x_1 \le 4$), as opposed to non-elementary regions, i.e. regions that consist of the union of several hyperrectangles (e.g. $\mu \equiv x_1 < 2 \lor x_1 > 4$).

**Definition 3.5 (Satisfaction probability of an MPM).** For $M_\theta$ an MPM and $\varphi$ an MITL formula, the satisfaction probability of $\varphi$ w.r.t. $M_\theta$ is defined as:

$$Pr(\varphi|M_\theta) = \Pr_{M_\theta}(\{(\sigma, 0) \models \varphi, \sigma \in Path_{M_\theta}\})$$

where $\Pr_{M_\theta}$ is the probability measure (3.1) induced by $M_\theta$ over $Path_{M_\theta}$.

Notice that assessing $Pr(\varphi|M_\theta)$ is indeed the target of probabilistic model checking and that its value can be obtained either exactly through *numerical model checkers* [KNP02, DJKV17] or being approximated through *statistical model checkers* [You05, LST16, BDD$^{+}$11a, SVA05a]. However since here we are interested in estimating how the satisfaction probability of $\varphi$ changes with the model's parameters we introduce the notion of satisfaction probability function.

**Definition 3.6 (Satisfaction probability function).** Let $(M_\theta)_{\theta \in \Theta}$ be a parametric MPM and $\varphi$ a MITL formula, the function:

$$f_{\varphi}: \Theta \rightarrow [0, 1] \\ \theta \rightarrow Pr(\varphi|M_{\theta})$$

is called the satisfaction probability function.

We will see how through an adaptation of the ABC method we manage to obtain an approximation of $f_{\varphi}$ for a parametric MPM model $M_\theta$.

### 3.1.3 Bayesian inference: the ABC method

Given a parametric model $M_\theta$ and an MITL formula $\varphi$, our goal is to estimate the satisfaction probability function $f_\theta$, i.e. the function that characterises how the probability that $\varphi$ is satisfied by $M_\theta$ varies w.r.t. the parameter $\theta \in \Theta$.

Given a set of observations $y_{obs}$ Bayesian inference is concerned with learning probability distributions over the parameters space $\Theta$ of a parametrised model $M_\theta$. Relying on the Bayesian interpretation of probability, initial beliefs on the parameters expressed via a *prior distribution* $\pi$ over $\Theta$ are progressively updated, via $y_{obs}$, based on the model's dynamics encoded by a *likelihood function* $p(y_{obs}|\theta)$ which expresses how probable $y_{obs}$ is to be observed given the model's parameters $\theta \in \Theta$. The output of the inference procedure
consists of the *posterior distribution* $\pi(\theta|y_{obs})$ computed over $\Theta$ based on the observed data $y_{obs}$ and which is defined by

$$ \pi(\theta|y_{obs}) = \frac{p(y_{obs}|\theta)\pi(\theta)}{\int_{\theta'} p(y_{obs}|\theta')\pi(\theta') d\theta'} $$

For realistic models the likelihood function $p(y_{obs}|\theta)$, necessary to derive the posterior distribution, is too expensive to compute or even intractable. For this reason *likelihood-free* methods, such as Approximate Bayesian Computation [MPRR12, SFB18], have been introduced to obtain an approximation, denoted $\pi_{ABC,\epsilon}(\theta|y_{obs})$, of the posterior distribution $\pi(\theta|y_{obs})$. With ABC schemes the likelihood is approximated by matching trajectories $y$ sampled, through simulation, from the likelihood $y \sim \pi(y|\theta)$ with the observations $y_{obs}$ via a distance metric $\rho(y, y_{obs})$. This lead first to the simple ABC *rejection sampling* algorithm which has then been improved yielding the faster converging ABC-SMC algorithm.

**ABC rejection sampling.** The simplest form of ABC, known as *rejection sampling*, operates by iteratively sampling parameters $\theta'$ from the prior distribution $\theta' \sim \pi(\cdot)$. For each $\theta'$ a trajectory $y' \sim p(\cdot|\theta')$ is simulated from the corresponding model instance $M_{\theta'}$ and $\theta'$ is accepted if $\rho(y', y_{obs}) \le \epsilon$ (i.e. if it is sufficiently matching observations $y_{obs}$) or rejected if $\rho(y', y_{obs}) > \epsilon$, where $\epsilon \in \mathbb{R}_{\ge 0}$ represents a chose tolerance. The accepted parameters $\theta_i$ together with the corresponding traces $y'_i$ give samples $(\theta_i, y_i)$ drawn from the joint distribution: $\pi_{ABC,\epsilon}(\theta, y|y_{obs}) \propto 1_{A_{\epsilon}(y, y_{obs})}p(y|\theta)\pi(\theta)$ where $A_{\epsilon}(y, y_{obs}) = \{y \in \mathcal{Y} | \rho(y, y_{obs}) \le \epsilon\}$ (with $1_{A_{\epsilon}(y, y_{obs})}$ denoting the indicator function representing the set of traces whose distance from $y_{obs}$ is within the tolerance $\epsilon$). $\pi_{ABC,\epsilon}$ approximates the posterior distribution: the smaller the $\epsilon$, the closer the simulations $y_i$ are to the observations $y_{obs}$, the better the approximation.

**Algorithm 1 ABC rejection sampling**

**Require:** $N$: number of particles, $y_{obs}$ observations, tolerance $\epsilon$, distance $\rho$, summary statistics $\eta$

**Ensure:** $(\theta_i)_{1 \le i \le N}$ drawn from $\pi_{ABC,\epsilon}$
**for** $i = 1:N$ **do**

**repeat**

$$ \theta' \sim \pi(\cdot) $$

$$ y' \sim p(\cdot|\theta') $$

**until** $\rho(\eta(y'), \eta(y_{obs})) \le \epsilon$

$$ \theta_i \leftarrow \theta' $$

$$ y_i \leftarrow y' $$

**end for**

Notice that in Algorithm 1, $\eta: \mathcal{Y} \to S \subset \mathbb{R}^{k_1}$ represents summary statistics$^5$ computed on

<footnote>$^5$The choice of summary statistics is a crucial point in ABC (see for example [ANJB10]).</footnote>
the observations $y_{obs}$ and on the simulated trace $y'$, while $\rho: S \times S \to \mathbb{R}^+$ is a distance in the space of summary statistics.

**ABC Sequential Monte Carlo method.** The chosen value of $\epsilon$ is crucial for the performance of the simple ABC algorithm: a small $\epsilon$ is needed to achieve a good approximation. However, this may result in a high rejection rate leading to cumbersome computations. To overcome this issue, the more elaborate algorithm known as ABC Sequential Monte Carlo (ABC-SMC), has been proposed [BCMR09]. It is an SMC based approach [DMDJ06] through which a population of $N$ particles is iteratively sampled with increasing accuracy until the targeted level of accuracy $\epsilon_M$ is obtained. At the first iteration, the particles are initialised through the simple ABC rejection sampling algorithm using a large enough $\epsilon_1$ to limit the computation cost. $\epsilon_1$ possibly equals infinity, which is equivalent to only sampling from the prior distribution. Then, at each step $i$, $i = 2, \dots, M$, the particles are moved by a transition kernel $K(\cdot|\cdot)$ (for example, a Gaussian one [DMDJ06]) until they match the tighter, next level, approximation constraint $\epsilon_i$. At iteration $M$, we finally get $N$ particles that fulfil the desired approximation $\epsilon_M$. Some ad-hoc strategies are proposed to find a proper sequence $(\epsilon_i)_{1\le i\le M}$ ensuring an efficient convergence towards the posterior distribution.

<caption>**Algorithm 2 ABC Sequential Monte Carlo**</caption>

**Require:** $N$: number of particles, $y_{obs}$, $(\epsilon_i)_{1\le i\le M}$, $\rho$, $\eta$
**Ensure:** $(\omega_j)_{1\le j\le N}$, $(\theta_j)_{1\le j\le N}$ weighted samples drawn from $\pi_{ABC,\epsilon_M}$
Iteration $i = 1$: find $(\theta_j^{(i)})_{1\le j\le N}$ with ABC rejection sampling with tolerance $\epsilon_1$
$\omega_j \leftarrow \frac{1}{N}$
**for** $i = 2: M$ **do**
&nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1: N$ **do**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**repeat**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Take $\theta'_j$ from $(\theta_j^{(i-1)})_{1\le j\le N}$ with probabilities $(\omega_j)_{1\le j\le N}$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\theta_j^{(i)} \sim K(\cdot|\theta'_j)$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y' \sim p(\cdot|\theta_j^{(i)})$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**until** $\rho(\eta(y'), \eta(y_{obs})) \le \epsilon_i$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\omega_j \leftarrow \frac{\pi(\theta_j^{(i)})}{\sum_{j'=1}^{N} \omega_{j'}^{(i-1)} K(\theta_j^{(i)}|\theta_{j'}^{(i-1)})}$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Normalize $(\omega_j)_j$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**

**Remark:** ABC is a *distance related class of methods*. A basic point about ABC algorithms is that, by definition, they all rely on a *distance metric* (between the simulations $y'$ and observations $y_{obs}$). It is precisely based on this characteristic that we extend ABC
statistical inference to the model checking problem. In Section 3.2, we first introduce the notion of satisfiability distance, which quantifies how far a model instance $M_\theta$ is from satisfying a temporal logic property $\varphi$ and then we exploit such a novel distance for adapting ABC schemes to the estimation of $f_\theta$, i.e. the satisfaction probability function for property $\varphi$. In practice we do so by plugging a distance automaton in the ABC procedure and use it as a machinery to assess how far trajectories issued by an MPM model with parameter $\theta$ are from satisfying an MITL formula $\varphi$. Furthermore, as we will demonstrate, distance automata yield a null distance for any simulation that satisfy the considered formula $\varphi$ therefore with the HASL-based extension of ABC, we actually estimate the ABC-posterior using a zero tolerance, i.e. with $\epsilon = 0$.

### 3.1.4 Kernel density estimation

Given that ABC methods only deliver samples $(\theta_i)_i$ drawn from the ABC posterior distribution $(\pi_{ABC,\epsilon}(\theta|y_{obs}))$ and that our end goal, instead, is to obtain an approximation of a continuous probability density function, namely the probability satisfaction $f_\varphi$, we need to resort to methods for reconstructing a density function given a set of samples. Specifically we resort to kernel density estimation (KDE), which allow one to derive an approximation $\hat{\pi}$ of an unknown probability density function $\pi$ given a finite number of samples of a random variable. The approximation $\hat{\pi}$ is obtained as the sum of the application of a kernel function $K$ to the samples, i.e. a continuous function which quantify the contribution, in terms of probability mass, brought by a sample to the density $\pi$ to be estimated (i.e. essentially, a kernel is a manner to weight data samples).

**Definition 3.7 (Kernel function).** A function $K: \mathbb{R} \to \mathbb{R}_{\ge 0}$ is called a kernel if:

$$1. \int_{\mathbb{R}} K(u)du = 1$$

$$2. \forall u \in \mathbb{R}, K(u) = K(-u)$$

Based on kernel functions, one can define a kernel density estimator [Sil86].

**Definition 3.8 (Kernel density estimator).** Let $\theta^{(1)}, \dots, \theta^{(N)}$ be $N$ i.i.d samples from an unknown density $\pi$ on $\Theta$. The kernel density estimator $\hat{\pi}$ associated with a kernel function $K$ on $\mathcal{X}$ is:

$$\forall \theta \in \Theta, \hat{\pi}(\theta) = \frac{1}{Nh} \sum_{i=1}^{N} K\left(\frac{\theta - \theta^{(i)}}{h}\right)$$

$K$ is a kernel function, $h$ is a smoothing parameter called the bandwidth$^6$.

<footnote>$^6$Where $h$ is either a scalar, a vector or a matrix, depending on the dimension of $\Theta$ and the choice of $K$.</footnote>
The rationale behind estimator $\hat{\pi}(\theta)$ is that each value $\theta$ contributes to the probability mass over the whole set $\Theta$, with the idea that the further $\theta$ is from the observation $\theta^{(i)}$, the lower is its contribution to the probability mass. The selection of the kernel $K$ and the calibration of the bandwidth $h$ are crucial points in KDE and several possibilities exist [Sil86], [CD18]. In Section 3.2.4 we adapt the KDE approach to the Automaton-ABC method more specifically by employing the minimisation of the Least Squares Cross-Validation as criterion for selecting the bandwidth $h$ jointly with Gaussian and Beta kernels [Che99].

## 3.2 Automaton-ABC

We introduce an adaptation of the ABC methodology, named automaton-ABC, aimed at approximating the satisfaction probability function $f_{\varphi}$ (Definition 3.6) of a MITL formula $\varphi$. It relies on the following four aspects: i) the formalisation of the notion of satisfiability distance expressing how far a model's path is from satisfying the considered formula $\varphi$ ii) the introduction of the corresponding HASL meter automaton to measure such distance, iii) the definition of a novel ABC algorithm whose convergence is driven by assessment of the satisfiability distance and, finally, iv) the derivation of the normalisation constant through which the KDE approximation of the probability satisfaction function $f_{\varphi}$ is obtained from the ABC posterior density. In the remainder we present the automaton-ABC approach taking into account simple, non-nested, time-bounded reachability problems only. Its extension to the full MITL spectrum could intuitively be straightforwardly obtained as long as the notion of satisfiability distance can be reasonably extended to the entire MITL spectrum.

### 3.2.1 Satisfiability distances for reachability problems

The basic idea behind the notion of satisfiability distance is that a time-bounded temporal logic specification $\varphi$, which depends on atomic propositions $\mu_i$ and time intervals $I_j \subseteq \mathbb{R}_{\ge 0}$, identifies a set of spatio-temporal regions $S_{\mu_i}^{I_j} \subseteq S \times \mathbb{R}_{\ge 0}$ (see Definition 3.4) that a path $\sigma \in \text{Path}_{M_{\theta}}$ should (partially or totally) traverse in order to satisfy $\varphi$. Therefore one can sensibly introduce a *satisfiability distance*, denoted $d: \text{Path}_{M_{\theta}} \times \Phi \rightarrow \mathbb{R}_{\ge 0}$, which quantifies how far a path $\sigma$ is from satisfying $\varphi$ by taking into account some form of Euclidean distance of $\sigma$ from the relevant regions $S_{\mu_i}^{I_j}$. To characterise $d(\sigma, \varphi)$ we consider the following guidelines: i) $d(\sigma, \varphi)$ should evaluate to zero whenever $\sigma \models \varphi$ and ii) $d(\sigma, \varphi)$ should favour the convergence of the ABC algorithm.

Before formally introducing it (Definition 3.9), we illustrate, in Figure 3.1, the notion of distance for the three kind of (non-nested) temporal MITL formulae, namely $\mathbf{F}^I\mu$ (eventual reachability), $\mathbf{G}^I\mu$ (global reachability) and $\mu_1\mathbf{U}^I\mu_2$ (conditional reachability), assuming, for simplicity, that $\mu, \mu_1$ and $\mu_2$ corresponds each to a single elementary region, hence that
both $\mathbf{F}^I\mu$ and $\mathbf{G}^I\mu$ induce a simple satisfiability region $S_\mu^I$, whereas $\mu_1\mathbf{U}^{[t_1,t_2]}\mu_2$ induce a compound satisfiability region given by $S_{\mu_1}^{[0,t_1]} \cup S_{\mu_2}^{[t_1,t_2]} \cup S_{\mu_1}^{[t_1,t_2]}$.

<figure><img src="image_12.png" /><figcaption>Figure 3.1: Examples of trajectories with zero-distance (left) and positive distance (right) from an F, a G and a U region (positive distances are depicted in red).</figcaption></figure>

Plots in the top row of Figure 3.1 depict cases in which a path $\sigma$ has to have a null satisfiability distance from a given formula. For $\mathbf{F}^I\mu$ distance has to be 0 if and only if $\sigma$ has at least one point traversing region $S_\mu^{[t_1,t_2]}$ (Figure 3.1a), for $\mathbf{G}^{[t_1,t_2]}\mu$ if and only if within $t \in [t_1, t_2]$ all points $\sigma@t$ fall in $S_\mu^{[t_1,t_2]}$ (Figure 3.1b) while for $\mu_1\mathbf{U}^{[t_1,t_2]}\mu_2$ the distance is 0 if and only if there exists $t' \in [t_1, t_2]$ such that $\sigma@t'$ is in $S_{\mu_2}^{[t_1,t_2]}$ while is consistently in $S_{\mu_1}^{[0,t']}$ beforehand (Figure 3.1c)<fnref n="7" />.

The bottom row of Figure 3.1 instead depicts cases in which a path $\sigma$ has to have a positive distance, illustrating the rationale we adopted in the characterisation of the metric $d(\sigma, \varphi)$ given in Definition 3.9. Specifically for $\mathbf{F}^{[t_1,t_2]}\mu$ (Figure 3.1d), assuming $\sigma$ does not enter $S_\mu^{[t_1,t_2]}$, the distance corresponds with the minimal Euclidian distance between $\sigma$ and $S_\mu^{[t_1,t_2]}$ (Equation (3.3)), whereas for $\mathbf{G}^{[t_1,t_2]}\mu$ (Figure 3.1e), assuming $\sigma$ leaves $S_\mu^{[t_1,t_2]}$, the distance corresponds with the volume of the hyperrectangles delimited by the segments of $\sigma$ that lie outside $S_\mu^{[t_1,t_2]}$ (Equation (3.4)). Finally, for $\mu_1\mathbf{U}^{[t_1,t_2]}\mu_2$ (Figure 3.1f) the distance of $\sigma$ (Equation (3.5)) bears 3 components: $d(\sigma, \mathbf{G}^{[0,t_1][\mu_1]})$ which accounts for the fact that a path satisfying $(\mu_1\mathbf{U}^{[t_1,t_2]}\mu_2)$ must never leave region $S_{\mu_1}$ before $t_1$, $d(\sigma, \mathbf{F}^{[t_1,t_2]}\mu_2)$ that accounts for the fact that $\sigma$ must enter region $S_{\mu_2}$ within $[t_1, t_2]$ and $d(\sigma, \mathbf{G}^{[t_1,t']}(\mu_1 \vee \mu_2))$

<footnote n="7">⁷Note that for some $t'' > t' \sigma@t'' \notin S_{\mu_1}$ but this is uninfluential, as the satisfaction of $\mu_1\mathbf{U}^{[t_1,t_2]}\mu_2$ by $\sigma$ is established as of what happens before $t'$ not what happens after.</footnote>
that accounts for the fact that there must be a time $t'$ where $\sigma$ switch from region $\mu_1$ to region $\mu_2$ directly, i.e. without spending time in any intermediate $S_{-(\mu_1 \lor \mu_2)}$ region: if that is not the case (i.e. if $\sigma$ within $[t_1, t_2]$ has points in the complementary region $S_{-(\mu_1 \lor \mu_2)}^{[t_1, t_2]}$) then Equation (3.5) yields a positive value which accounts for the sum of the minimal distances of each such point from either regions $S_{\mu_1}$ or $S_{\mu_2}$.

**Definition 3.9 (Satisfiability metrics).** Given a path $\sigma \in Path_M$ of a $n$-dimensional MPM $M$ a closed time-bounding interval $[t_1, t_2] \subset \mathbb{R}_{\ge 0}$, an elementary propositional formulae $\mu, \mu_1, \mu_2$, we define the distance $d(\sigma, \varphi)$ from the satisfiability region for the following kinds of temporal formulae:

$$1) \varphi \equiv \mathbf{F}^{[t_1, t_2]} \mu$$

$$d(\sigma, \mathbf{F}^{[t_1, t_2]} \mu) = \begin{cases} d_e((\sigma@t_l(t_2), t_l(t_2)), S_{\mu}^{[t_1, t_2]}) & \text{if } t_l(t_2) < t_1 \land \sigma@t_l(t_2) \not\models \mu \\ \min(d_e(\sigma@t_l(t_1), S_{\mu}^{[t_1, t_2]}), \min_{t \in [t_1, t_2]} d_e((\sigma@t, t), S_{\mu}^{[t_1, t_2]})) & \text{otherwise} \end{cases} \quad (3.3)$$

where $t_l(t^*) = \min\{t \in [0, t^*] : \forall t' \in [t, t^*], \sigma@t' = \sigma@t\}$ is the time instant of the last jump occurred on $\sigma$ before $t^*$ and $d_e((s, t), S_1 \times T_1) = \min_{t' \in T_1, s' \in S_1} \sqrt{(t - t')^2 + \sum_{i=1}^d (s[i] - s'[i])^2}$ denotes the euclidean distance of a point $(s, t) \in S \times \mathbb{R}_{\ge 0}$ from the closest point of a time-bounded region $S_1 \times T_1$.

For non-elementary propositional formulae $\mu \equiv \bigvee \mu_i$, we define the distance:

$$d(\sigma, \mathbf{F}^{[t_1, t_2]} \bigvee \mu_i) = \min_i d(\sigma, \mathbf{F}^{[t_1, t_2]} \mu_i)$$

where $\mu_i$ are elementary formulae.

$$2) \varphi \equiv \mathbf{G}^{[t_1, t_2]} \mu$$

$$d(\sigma, \mathbf{G}^{[t_1, t_2]} \mu) = \int_{t_1}^{t_2} d_e(\sigma@t, S_{\mu})dt \quad (3.4)$$

where $d_e(s, S_1) = \min_{s' \in S_1} \sqrt{\sum_{i=1}^d (s[i] - s'[i])^2}$ denotes the euclidean distance of a point $s \in S$ from the closest state of a state space subset $S_1 \subseteq S \subseteq \mathbb{N}^d$. $d(\sigma, \mathbf{G}^{[t_1, t_2]} \mu)$ is the integral of the Euclidean distance from $S_\mu$ of any point of $\sigma$ that occurs within $[t_1, t_2]$,
Similarly, for non-elementary propositional formulae, we define the distance

$$d(\sigma, \mathbf{G}^{[t_1, t_2]} \lor \mu_i) = \min_i d(\sigma, \mathbf{G}^{[t_1, t_2]} \mu_i)$$

where $\mu_i$ are elementary formulae.

$$3) \varphi \equiv \mu_1 \mathbf{U}^{[t_1, t_2]} \mu_2$$

$$d(\sigma, \mu_1 \mathbf{U}^{[t_1, t_2]} \mu_2) = d(\sigma, \mathbf{G}^{[0, t_1]} \mu_1) + d(\sigma, \mathbf{F}^{[t_1, t_2]} \mu_2) + d(\sigma, \mathbf{G}^{[t_1, t_{min}]} (\mu_1 \lor \mu_2)) \quad (3.5)$$

where $t_{min} = \min(\arg\min_{t \in [t_1, t_2]} d_e(\sigma@t, S_{\mu_2}))$ is the earliest time corresponding to the closest point between $\sigma$ and region $\mu_2$.

**Satisfiability metric and ABC convergence aspects.** In order to be employed in ABC frameworks, satisfiability distances shall account for the convergence of the ABC algorithms, in the first place by ensuring that each path is ranked with an as large a value of distance as the path is further from satisfying the considered formula, which, indeed (3.3), (3.4) and (3.5) do. Experimental evidence showed that the convergence of ABC algorithms for eventual formulae $\mathbf{F}^{[t_1, t_2]}\mu$ is affected by a peculiar aspect of a model's path $\sigma$, that is, the presence/lack of jumps within the bounding interval $[t_1, t_2]$. Specifically if, within $[t_1, t_2]$, $\sigma$ does not contain any jump (and lies outside region $S_{\mu}^{[t_1, t_2]}$) then it is more convenient (from a convergence standpoint) that the distance $d(\sigma, \mathbf{F}^{[t_1, t_2]}\mu)$, i.e. (3.3), is set to the Euclidian distance between region $S_{\mu}^{[t_1, t_2]}$ and the point entered by $\sigma$ at the last jump occurred before entering $[t_1, t_2]$ even if within $[t_1, t_2]$ the path is actually closer to $S_{\mu}^{[t_1, t_2]}$ (e.g. Figure 3.2a). Such an aspect ensures that the ABC-driven parameter search is not mislead by anomalous situations such as, e.g. parameters $\theta \in \Theta$ that yield a model $\mathcal{M}_{\theta}$ for which there is a non-null probability of reaching an absorbing state in $S_{-\mu}$ before $t_1$. On the other hand, if $\sigma$ contains jumps in $[t_1, t_2]$, then it is more convenient that $d(\sigma, \mathbf{F}^{[t_1, t_2]}\mu)$ is set to the Euclidian distance of the closest point amongst the point corresponding to the last jump before $t_1$ and those corresponding to jumps occurring in $[t_1, t_2]$ (e.g. Figure 3.2c).

An essential aspect of the satisfiability metrics introduced in Definition 3.9 is their soundness, i.e., to prove that they yield a null distance $d(\sigma, \varphi) = 0$ if and only if $\sigma \models \varphi$.

**Proposition 3 (Soundness of satisfiability metric for elementary F, G and U).** For $\sigma \in \mathcal{P}ath_M$ a path of an MPM $M$ and $\mu, \mu_1, \mu_2$ are MITL propositional formulae in DNF then

$$\sigma \models \mathbf{F}^{[t_1, t_2]} \mu \iff d(\sigma, \mathbf{F}^{[t_1, t_2]} \mu) = 0$$

$$\sigma \models \mathbf{G}^{[t_1, t_2]} \mu \iff d(\sigma, \mathbf{G}^{[t_1, t_2]} \mu) = 0$$

$$\sigma \models \mu_1 \mathbf{U}^{[t_1, t_2]} \mu_2 \iff d(\sigma, \mu_1 \mathbf{U}^{[t_1, t_2]} \mu_2) = 0$$

where $d$ is as per Definition 3.9.

*Proof.* see [BBC21] □
<figure><img src="image_13.png" /><figcaption>Figure 3.2: Automaton $\mathcal{A}_F$ (top) for measuring the distance of a path $\sigma$ for an eventual property concerning observed species $X_O$ and examples (bottom) of measured distance $d$: positive distance (a), null distance (b), selection of the minimum distance (c) in case of presence of jumps in $[t_1, t_2]$ and evolution of the computed distance $d$ along a path (d).</figcaption></figure>

### 3.2.2 HASL specifications for satisfiability distances

Having established the notion of satisfiability metric for non-nested MITL reachability properties we introduce hybrid automata *meters* for assessing them. For simplicity in the remainder we present automata for mono-dimensional reachability properties (i.e. temporal formulae built on top of elementary conditions such as $\mu \equiv x_1 \le x_O \le x_2$, where $x_O$ denotes the population of an observable quantity *O* and $x_1 < x_2 \in \mathbb{N}$), bearing in mind that distance automata for formulae based on *n*-dimensional regions are straightforward adaptations of those discussed below.

Meter automaton $\mathcal{A}_F$. Automaton $\mathcal{A}_F$ (Figure 3.2) is designed to measure the distance (3.3) of a path $\sigma$ from the region associated with $F^{[t_1, t_2]}(x_1 \le x_O \le x_2)$, i.e. the region corresponding to the observed species $x_O \in [x_1, x_2]$ within time $t \in [t_1, t_2]$. It uses 4 variables: *d* (computed distance), *t* (current time along the path), *n* (population of the observed species *O* after the most recent occurrence of a reaction) and *n'* (population of *O* before the most recent occurrence of a reaction). The synchronisation of $\sigma$ with $\mathcal{A}_F$ is managed through a number of mutually exclusive *autonomous* transitions (from $l_1$ to $l_3$), plus a single *synchronised* transition (from $l_3$ to $l_1$), which results in the automaton looping between $l_1$ and $l_3$ up until a termination condition is fulfilled. It is straightforward to show that $\mathcal{A}_F$ complies with the HASL determinism constraints (Chapter 2) and
therefore the synchronisation of an arbitrary path $\sigma$ yields a unique path in the product process $M \times A_F$. Specifically, synchronisation of $\sigma$ with $A_F$ works as follows. At the start ($l_0 \to l_1$) the distance is initialised to $d := \infty$ and the initial value of the observed species are stored in $n:=x_O$. Once in $l_1$, the analysis of $\sigma$ begins and is driven by seven mutually exclusive autonomous transitions. If initially $\sigma$ is inside the region (and this include even initially with $t=0$ in case $t_1=0$ too), then transition $l_1 \xrightarrow{\sharp, t_1 \le t \ge t_2 \wedge (x_1 \le n \le x_2), \{d:=0\}} l_2$ occurs immediately and the synchronisation stops with distance $d=0$. On the other hand if, while in $l_1$, the path has not entered $[x_1, x_2]$, distance $d$ must be computed depending on different conditions (that correspond to 4 mutually exclusive autonomous transitions linking $l_1 \to l_3$). Specifically: in case $t < t_1$ (i.e. $\sigma$ has not yet temporally reached the time interval $[t_1, t_2]$) then either $\sigma$ has entered $[x_1, x_2]$, in which case $d$ is correctly set to $d := 0^8$ through firing of $l_1 \xrightarrow{\sharp, (t \le t_1) \wedge (x_1 > n \wedge n > x_2) \wedge \{d:=0, n' := n\}} l_3$ or, $\sigma$ has not entered $[x_1, x_2]$ and then $d$ is set to the Euclidian distance of the current point of the path from the nearest corner of the region (either $(t_1, x_1)$ or $(t_1, x_2)$) through firing of $l_1 \xrightarrow{\sharp, (t \le t_1) \wedge (x_1 > n \wedge n > x_2) \wedge \{d:=\min(\sqrt{(t-t_1)^2+(n-x_2)^2}, \sqrt{(t-t_1)^2+(n-x_1)^2})\}} l_3$.

On the other hand if $t \ge t_1$, in accordance with (3.3), the distance of the current point of the path is either: i) left unchanged (by firing of $l_1 \xrightarrow{\sharp, (t \ge t_1) \wedge (n=n') \wedge (x_1 > n \wedge n > x_2) \wedge \emptyset} l_3$), if the last occurred reaction had not produced a jump w.r.t. the observed species (i.e. $n' = n$), or, conversely, ii) to the minimum between the previous value of $d$ and the distance of the current point from $[x_1, x_2]$ (by firing of $l_1 \xrightarrow{\sharp, t \ge t_1 \wedge (n \ne n') \wedge (x_1 > n \wedge n > x_2) \wedge \{d:=\min(d, \min(|n-x_1|, |n-x_2|))\}} l_3$) if the last occurred reaction did produce a jump w.r.t. O.

**Meter automaton** $\mathcal{A}_G$. Automaton $\mathcal{A}_G$ (Figure 3.3) is designed to measure the distance of a path $\sigma$ w.r.t. to a formula $G^{[t_1, t_2]}(x_1 \le x_O \le x_2)$, based on (3.4). It uses the same variables as $A_F$ (hence $d$ stores the measured distance corresponding with the integral of the segments that, within $[t_1, t_2]$, fall outside the region) plus an extra timer $t'$, to measure the duration of a segment falling outside the region within $[t_1, t_2]$, and a boolean flag $in$, which is set to true if the last segment of the path originates in $[t_1, t_2]$ outside of the region $[x_1, x_2]$. $in$ is used to distinguish cases where the path is out of the region $[x_1, x_2]$ with $t' < t_1$ and a new event occurs after a time $t'' > t_1$, in order to add $(t''-t_1) * \min(|n-x_1|, |n-x_2|)$ instead of $(t''-t') * \min(|n-x_1|, |n-x_2|)$. After the initialisation of variables ($l_0 \to l_1$), analysis begins in $l_1$: for events occurring before $t_1$, we distinguish two cases. If $\sigma@t \in [x_1, x_2]$, the distance is set to zero ($l_1 \to l_3$ top arc). Otherwise, $d$ is the distance of $\sigma@t$ from $[x_1, x_2]$ otherwise ($l_1 \to l_3$ midway arc). Indeed, if, for example, the next jump of $\sigma$ happens at $t > t_2$, then the final distance is given by $d \cdot (t_2 - t_1)$ ($l_1 \to l_2$ bottom arc). For events occurring at $t \in [t_1, t_2]$, if $\sigma@t \notin [x_1, x_2]$ (sequence $l_1 \to l_4 \to l_1$), the distance is incremented

<footnote n="8">⁸this is because MPM paths are *càdlàg* functions of time, if the next reaction occurs at time $t \ge t_1$ then it is certain that the current path has at least one point within the considered region hence the path will be then accepted ($l_1 \xrightarrow{\sharp, (d=0) \wedge (t \ge t_1)} l_2$) and the finally measured distance will be $d=0$).</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 3.3: Automaton $\mathcal{A}_G$ for global property.</figcaption></figure>

by the surface defined by the path segment (of duration $t'$) laying outside $[x_1, x_2]$ and the closest border of $[x_1, x_2]$. The distance is left unchanged if $\sigma@t \in [x_1, x_2]$ (sequence $l_1 \to l_3 \to l_1$).

Having introduced the meter automata for formulae $\varphi \equiv \mathbf{F}$ and $\varphi \equiv \mathbf{G}$ we prove their soundness. i.e. we demonstrate that the value stored in the automaton's variable $d$ at the end of the synchronisation with a path $\sigma$ corresponds with the distance $d(\sigma, \varphi)$ given in Definition 3.9.

**Proposition 4.** Let $X_O$ be a species of an MPM model $\mathcal{M}$, $\mathcal{A}_F$ be the distance LHA corresponding to the MITL reachability formula $\varphi \equiv \mathbf{F}^{[t_1, t_2]} x_1 \le x_O \le x_2$ and $\sigma \in \text{Path}_\mathcal{M}$ be a path of $\mathcal{M}$, then:

$$ \text{last}(\mathcal{A}_F(\sigma).d) = d(\sigma, \mathbf{F}^{[t_1, t_2]} x_1 \le x_O \le x_2) $$

$$ \text{last}(\mathcal{A}_G(\sigma).d) = d(\sigma, \mathbf{G}^{[t_1, t_2]} x_1 \le x_O \le x_2) $$

where $\text{last}(\mathcal{A}_F(\sigma).d)$ (resp. $\text{last}(\mathcal{A}_G(\sigma).d)$) is the value stored in variable $d$ of automata $\mathcal{A}_F$ (resp. $\mathcal{A}_G$) at the end of the synchronisation between $\mathcal{A}_F$ (resp. $\text{last}(\mathcal{A}_G(\sigma).d)$) and $\sigma$ while $d(\sigma, \mathbf{F}^{[t_1, t_2]} x_1 \le x_O \le x_2)$ (resp. $d(\sigma, \mathbf{G}^{[t_1, t_2]} x_1 \le x_O \le x_2)$) is the distance of $\sigma$ from $\mathbf{F}^{[t_1, t_2]} x_1 \le x_O \le x_2$ ($\mathbf{G}^{[t_1, t_2]} x_1 \le x_O \le x_2$) as per the metric in Definition 3.9.

*Proof.* See [BBC21]. □

**Meter automaton for other formulae.** We stress that automata for other kind of temporal formulae such as $\mathbf{U}$ and the conjunction of a $\mathbf{G} \wedge \mathbf{F}$ have been obtained [Ben21] following a similar approach to that outlined for $\mathbf{F}$ and $\mathbf{G}$ formulae. Specifically the automaton for MITL formula: $\mathbf{G}^{[t_1, t_2]}(x_1 \le x_O \le x_2) \wedge \mathbf{F}^{[t_3, t_4]}(x_3 \le x_{O'} \le x_4)$ with $t_2 \le t_3$ (i.e. the $\mathbf{G}$ region precedes the $\mathbf{F}$ region) while $x_1, x_2, x_3, x_4 \in \mathbb{N}$ and $x_O$, resp. $X_{O'}$, denotes the population of species $\mathbf{O}$, resp. $\mathbf{O'}$ is shown to correspond to a combination of automata $\mathcal{A}_G$ and $\mathcal{A}_F$ [Ben21].
### 3.2.3 Automaton-ABC: ABC with satisfiability distance

Based on the notion of satisfiability distances and on the corresponding automata meters to measure them, we introduce the Automaton-ABC algorithms, i.e. the adaptation of the ABC-algorithms to estimate the satisfaction probability of a MITL formula $\varphi$ by a parametric MPM $(M_\theta)_\theta$. The basic idea behind Automaton-ABC schemes is to inject a model-automaton synchronisation procedure through which trajectories of the product process $M_\theta \times A_\varphi$ are sampled and rejected if the distance measured by $A_\varphi$ through synchronisation is non-null. The accepted parameter samples $\theta_i$ constitute a sample to be used for estimating Automaton-ABC posterior distribution which we denote $\pi_{\varphi-ABC}$. We point out the estimation of the Automaton-ABC posterior distribution, differently from classical ABC schemes, is no longer computed as a limit approximation (i.e. $\lim_{\epsilon \to 0} \hat{\pi}_{ABC,\epsilon}(\cdot|y_{obs})$), but rather as an estimation of the exact distribution, since paths are accepted exclusively if their distance to the satisfiability regions for the considered formula $\varphi$ is zero.

**Automaton-ABC rejection sampling.** Algorithm 3 introduces a modified version of the rejection sampling ABC Algorithm 1 adapted to satisfiability distances. The algorithm takes as inputs a parametric MPM $(M_\theta)_{\theta \in \Theta}$, a prior distribution $\pi$ over $\Theta$ and a satisfiability distance automaton $A_\varphi$ corresponding to a MITL formula $\varphi$. The workflow is as for classical ABC rejection sampling hence at each iteration, a parameter $\theta'$ is drawn from the prior $\pi(\cdot)$, however now a path $\sigma'$ is sampled from the product process $M_{\theta'} \times A_\varphi$ and the corresponding satisfiability distance $d(\sigma', \varphi)$ is stored in one variable of $A_\varphi$; $\theta'$ is accepted only if the distance from $\varphi$ is $d(\sigma', \varphi)=0$. (i.e. if $\sigma' \models \varphi$ by Proposition 3).

**Algorithm 3** Automaton-ABC rejection sampling

**Require:** $(M)_{\theta \in \Theta}$ a pMPM, $\pi(\cdot)$ prior, $A_\varphi$ satisfiability distance automaton for MITL formula $\varphi$, $N$ number of particles

**Ensure:** $(\theta_i)_{1 \le i \le N}$ drawn from $\pi_{\varphi-ABC}$
for $i = 1:N$ do

repeat

$$ \theta' \sim \pi(\cdot) $$

$$ d(\sigma', \varphi) \sim (M_{\theta'} \times A_\varphi) $$

until $d(\sigma', \varphi) = 0$

$$ \theta_i \leftarrow \theta' $$

end for

**Automaton-ABC sequential montecarlo.** The Automaton-ABC rejection sampling (Algorithm 3) suffers from slow convergence in an even stronger fashion than its standard counterpart given the strict acceptance condition (i.e. only parameters that yield a null satisfiability distance are accepted). To improve convergence we introduce the sequential Monte Carlo version of Automaton-ABC (Algorithm 4), through which we insert
a decreasing (distance) tolerance $\epsilon \in \mathbb{R}_{>0}$ for accepting parameters. In essence we rank
parameters according to the satisfiability distance of the corresponding paths and accept
only parameters whose corresponding satisfiability distance is within tolerance $\epsilon$ (even if
the corresponding paths do not necessarily satisfy $\varphi$). This can lead to a faster convergence
of the parameter inference process hence resulting in an algorithm with a reduced *runtime*
(w.r.t. Algorithm 3) but that still yields samples from the posterior satisfiability density
function $\pi_{\varphi-ABC}$ (as per Proposition 1).

<caption>**Algorithm 4 Automaton-ABC Sequential Monte Carlo**</caption>

**Require:** $(\mathcal{M})_{\theta \in \Theta}$ a pMPM, $\pi(\cdot)$ prior, $\mathcal{A}_{\varphi}$ satisfiability distance automaton for MITL formula $\varphi$, $N$ number of particles, $\alpha \in (0, 1)$, $K$ kernel distribution

**Ensure:** $(\omega_j, \theta_j)_{1 \le j \le N}$ weighted samples drawn from $\pi_{\varphi-ABC}$

$$ (\theta_j^{(1)})_{1 \le j \le N} \sim \pi(\cdot) $$

$$ \forall j \in 1, \dots, N, d_j \sim (\mathcal{M}_{\theta_j^{(1)}} \times \mathcal{A}_{\varphi}) $$

$$ \epsilon \leftarrow \text{quantile}(\alpha, (d_j)_{1 \le j \le N}) $$

$$ (\omega_j)_{1 \le j \le N}^{(1)} \leftarrow \frac{1}{N} $$

$$ i \leftarrow 2 $$

**while** $\epsilon > 0$ **do**

**for** $j = 1:N$ **do**

**repeat**

Take $\theta'_j$ from $(\theta_j^{(i-1)})_{1 \le j \le N}$ with probabilities $(\omega_j^{(i-1)})_{1 \le j \le N}$

$$ \theta'_j \sim K(\cdot | \theta'_j) $$

$$ d_j \sim (\mathcal{M}_{\theta_j^{(i)}} \times \mathcal{A}_{\varphi}) $$

**until** $d_j \le \epsilon$

$$ \omega_j \leftarrow \frac{\pi(\theta_j^{(i)})}{\sum_{j'=1}^{N} \omega_{j'}^{(i-1)} K(\theta_j^{(i)} | \theta_{j'}^{(i-1)})} $$

**end for**

Normalise $(\omega_j^{(i)})_j$

$$ \epsilon \leftarrow \text{quantile}(\alpha, (d_j)_{1 \le j \le N}) $$

$$ i \leftarrow i + 1 $$

**end while**

The Automaton-ABC SMC algorithm (Algorithm 4), takes the same inputs as the Automaton-ABC rejection sampling (Algorithm 3), i.e. a parametric MPM $\mathcal{M}_{\theta}$, a satisfiability distance automaton $\mathcal{A}_{\varphi}$, the size $N \in \mathbb{N}$ of parameters' sample to be drawn at each iteration and a prior distribution $\pi(\cdot)$, plus two additional ones: a kernel distribution $K$ and a hyperparameter $\alpha \in ]0, 1[$ used to control how fast the tolerance $\epsilon$ decreases along with the iterations. The algorithm operates iteratively following the ABC-SMC schemes of Algorithm 2 yet with a few peculiarities. Initially, $N$ parameters $(\theta_j^{(1)})_{1 \le j \le N}$ are drawn from the prior $\pi(\cdot)$, their corresponding satisfiability distances $(d_j)_{1 \le j \le N}^{(1)}$ are computed (through synchronised simulations $\mathcal{M}_{\theta_j^{(1)}} \times \mathcal{A}_{\varphi}$). Parameters $\theta_j^{(1)}$ of this initial sample
are all accepted (regardless of their distance $(d_j)^{(1)}$) and coupled with weights purposely uniformly set to $(\omega_j)_{1\le j\le N}^{(1)} \leftarrow \frac{1}{N}$ to form the first draw from the posterior distribution $\pi_{\varphi-ABC}$. Furthermore $\theta_j^{(1)}$ are the used to establish the first tolerance level $\epsilon$ which is set to the $\alpha$-quantile of the distances $(d_j)_{1\le j\le N}^{(1)}$ (i.e. $\epsilon \leftarrow \text{quantile}(\alpha, (d_j)_{1\le j\le N}^{(1)})^9$). Notice that the first tolerance level $\epsilon$ may actually result in a quite large value of distance $d_j^{(1)}$ as no control is imposed on the selection of the initial sample $\theta_j^{(1)}$. The convergence of the parameter search process is then driven through the iterative loop. At each iteration $i \ge 2$ a new sample of parameters $\theta_j^{(i)}$ is obtained by moving parameters $\theta_j^{(i-1)}$ accepted at the previous iteration through the kernel distribution $K$. Each newly sampled parameters $\theta_j^{(i)}$ is accepted only if its corresponding distance $d_j^{(i)}$ (assessed through synchronised simulation $M_{\theta_j^{(i)}} \times A_{\varphi}$) is below the tolerance level $\epsilon$ established at previous iteration. The accepted $N$ parameters $\theta_j^{(i)}$ are then used to establish the next tolerance level $\epsilon$ as the $\alpha$-quantile of the distances $(d_j^{(i)})_{1\le j\le N}$. The iterations end as soon as the last tolerance level $\epsilon = 0$ is reached. The introduction of several steps with positive decreasing tolerances leads to an efficient exploration of the parameter space driven by $A_{\varphi}$.

**Theorem 1.** For $(M_\theta)_{\theta \in \Theta}$ a parametric MPM, $\varphi$ an MITL formula and $\pi$ a prior distribution over the parameter set $\Theta$, the $(\theta_i)_{1\le i\le N}$ issued by either Algorithm 3 or Algorithm 4 are drawn from a density function $\pi_{\varphi-ABC}$:

$$ \pi_{\varphi-ABC}(\theta_i) = f_{\varphi}(\theta_i) \cdot \frac{\pi(\theta_i)}{K} $$

where $f_{\varphi}: \Theta \to [0, 1]$ (i.e. $f_{\varphi}(\theta) = \text{Pr}(\varphi | M_\theta)$) is the probability satisfaction function of $\varphi$ (Definition 3.6) and $K \in \mathbb{R}_{\ge 0}$ is a positive constant.

*Proof.* See [BBC21].

<figure><img src="image_15.png" /></figure>

Theorem 1 links the Automaton-ABC algorithms (Algorithm 3 and 4) with the satisfaction probability function $f_{\varphi}$ (Definition 3.6). It establishes that with Automaton-ABC the regression of a smooth function, i.e. the likelihood function targeted by classical ABC schemes, is transformed into the regression of a probability density function, i.e. the probability satisfaction function $f_{\varphi}$. In this respect we observe that each parameter $(\theta_i)_{1\le i\le N}$ sampled (through the Automaton-ABC method) from the $\pi_{\varphi-ABC}$ is informative as it corresponds to a path that verifies $\varphi$. Furthermore, as $\pi_{\varphi-ABC}$ is a probability density function, the relative position of accepted parameters is highly informative w.r.t. the estimation of the satisfaction probability function: the denser a subset of accepted parameters, the higher the satisfaction probability function over the subset.

<footnote n="9">⁹Therefore $\epsilon$ is set to the distance $\epsilon = d_k$, with $1 \le k \le N$ such that the empirical cumulative distribution computed over $(d_j)_{1\le j\le N}^{(1)}$ is $\ge \alpha$.</footnote>
### 3.2.4 Estimation of the satisfaction probability function

Following Theorem 1 we have that:

$$f_{\varphi}(\theta_i) = K \cdot \frac{\pi_{\varphi-ABC}(\theta_i)}{\pi(\theta_i)}$$

This implies that if from the parameter's samples $(\theta^{(i)})_{1\le i\le M}$ we obtain by application of either Automaton-ABC algorithm we can derive an estimate of the ABC posterior density $\pi_{\varphi-ABC}(.)$ and if we can then further estimate the constant $K$ then we inherently obtain an estimate of the satisfaction probability function $f_{\varphi}(.)$ as the product $K \cdot \pi_{\varphi-ABC}(.)$ divided by the (known) prior $\pi(.)$ (the one which has been employed in the application of the Automaton-ABC algorithm).

**Estimation of the $\varphi - ABC$ posterior distribution.** To obtain an estimate of the Automaton-ABC posterior distribution $\pi_{\varphi-ABC}(.)$ we apply kernel density estimation (Section 3.1.4) to the samples $(\theta_i)_{1\le i\le N}$ issued by an Automaton-ABC algorithm. Specifically we considered two different kinds of kernels: Gaussian and Beta [Che99]. Beta kernels are useful when we have to estimate densities over bounded supports with positive probabilities on the boundaries, but are more computationally expensive for the calibration of the bandwidth. The optimal bandwidth is obtained by Least Squares Cross-Validation minimisation [Sil86].

**Estimation of $K$.** An estimate of $K$ can be obtained directly as the ratio between a single-point estimation of $Pr(\varphi|M_{\theta^*})$ and $\pi_{\varphi-ABC}(\theta^*)$ where $\theta^* \in \Theta$ is a specifically chosen value in the parameter space. Since we propose to estimate $Pr(\varphi|M_{\theta^*})$ through (standard) statistical model checking, $\theta^*$ should be chosen so that the runtime to statistically model check $M_{\theta^*} \models \varphi$ (hence to estimate $Pr(\varphi|M_{\theta^*})$) is not excessively high, i.e. which is the case as long as $Pr(\varphi|M_{\theta^*}) \gg 0$. Therefore, based on the estimate of the posterior distribution $\pi_{\varphi-ABC}$ one should chose $\theta^*$ wisely in a high probability region of $\pi_{\varphi-ABC}$. Alternatively, in order to improve kernel density estimation stability, one could select several $\theta^*$, estimate the corresponding constants $K^*$ and compute their mean value.

## 3.3 Experiments

We demonstrate the effectiveness of parametric model checking via the automaton-ABC method through a few examples of biological systems.

### 3.3.1 Enzymatic reaction system

We consider a simple model of enzymatic reaction (ER) system in which a *substrate* species $S$ is converted into a *product* $P$ through mediation of an *enzyme* $E$. The dynamics is given
by chemical equations (3.6) which depend on the parameters $\theta=\{k_1, k_2, k_3\}$, i.e. the kinetic rate constants of reactions $R_1, R_2, R_3$.

$$R_1 : E + S \xrightarrow{k_1} ES \quad R_2 : ES \xrightarrow{k_2} E + S \quad R_3 : ES \xrightarrow{k_3} E + P \quad (3.6)$$

The dynamics of the ER system is such that the totality of the substrate (initially $S_0 = 100$) is converted into the product at a speed dependent on parameters $\theta$. For example with $\theta=(1, 1, 1)$, the totality of $S$ is converted to $P$ before time = 5 (Figure 3.4 left), whereas with a tenfold speed reduction in the formation of the $ES$ complex and synthesis of $P$ (i.e. $\theta=(0.1, 1, 0.1)$), only about 30% of $S$ has been converted at time = 5 (Figure 3.4 right).

<figure><img src="image_16.png" /><figcaption>Figure 3.4: Trajectories of the ER system with $\theta_{\text{left}}=(1, 1, 1)$, $\theta_{\text{right}}=(0.1, 1, 0.1)$.</figcaption></figure>

**Satisfaction probability function estimation.** As a case study we consider the estimation of the satisfaction probability function $f_{\varphi_i}$ for two sets of time-bounded reachability formulae $\varphi_1, \varphi_2, \varphi_3$ (Figure 3.6) that we use for mono-dimensional experiments and $\varphi_4, \varphi_5, \varphi_6$ (Figure 3.7) which we use for bi-dimensional experiments. While $\varphi_1, \varphi_2, \varphi_3, \varphi_4$ and $\varphi_5$ refer to “simple” time-bounded reachability conditions, i.e. they are concerned with the probability that the population of a single species (i.e. species $P$ for $\varphi_1, \varphi_2, \varphi_3, \varphi_4$ and $E$ for $\varphi_5$) enter a specific elementary spatio-temporal regions (i.e. $R_1, R_2, R_3, R_4$ and $R_5$), $\varphi_6$ consists of a conjunction of time-bounded conditions involving two different species, namely $E$ and $P$, thus corresponds with a combination of 2 elementary regions one constraining the $E$-projection the other concerning $P$-projection of paths issued by the ER model. Before discussing the outcome of $f_{\varphi_i}$ estimation we point out the effect that parameter $k_3$ has on the dynamics of the $P$ species by showing few paths (Figure 3.5 left) sampled from the ER model by varying parameter $k_3 \in \{10, 20, 50\}$ while keeping $k_1 = k_2 = 1$ (i.e. equally coloured plots correspond to the same parameter set). Intuitively one can point out that parameter set $\theta_1:(1, 1, 50)$ induces a high probability that trajectories enter region $R_1$ (which corresponds to formula $\varphi_1$), while $\theta_2=(1, 1, 20)$ induces a high probability on region $R_2$ and $\theta_3=(1, 1, 10)$ on region $R_3$.
<figure><img src="image_17.png" /><figcaption>Figure 3.5: Projection of trajectories of ER system along P-dimension, and three spatial-temporal regions with corresponding MITL formulae encoding.</figcaption></figure>

**Mono-dimensional experiments.** Figure 3.6 illustrates the estimated posterior distribution $\pi_{\varphi_i-ABC}$ and corresponding satisfaction probability function $f_{\varphi_i}$ in function of parameter $k_3$ (1D experiments) obtained by application of the automaton-ABC method (Algorithm 3 and 4) to formulae $\varphi_1$, $\varphi_2$ and $\varphi_3$.

<figure><img src="image_18.png" /><figcaption>Figure 3.6: Weighted histograms of posterior distribution $\pi_{\varphi_i-ABC}$ for ER system obtained with Automaton-ABC-SMC using constant $k_1 = k_2 = 1$, prior distributions $\pi_{k_3}(\cdot) \sim U(0, 100)$ and using 1000 particles. The red plot depicts the estimated satisfaction probability function obtained through kernel density estimation method, while the blue plot depicts the “true” satisfaction probability function estimated on a selection of points through statistical model checking (at 95% confidence semi-interval width 0.01). In red: the satisfaction function estimated through kernel density estimation method.</figcaption></figure>

The marginal for $\pi_{\varphi_1-ABC}$ exhibits a rather uniform profile with the 95% credibility interval that $\varphi_1$ is satisfied for $k_3 \in [20, 100]$ (approximately) whereas the marginals for $\pi_{\varphi_2-ABC}$ and $\pi_{\varphi_3-ABC}$, instead, result in narrower 95% credibility intervals with $k_3 \in [15, 50]$ ($\varphi_2$) resp. $k_3 \in [5, 25]$ ($\varphi_3$). The estimation of the satisfaction probability function (red plot in Figure 3.6) from the posterior $\pi_{\varphi_i-ABC}$ has been obtained using Beta kernels for $\varphi_1$ (to handle boundary constraints) and Gaussian kernels for $\varphi_2, \varphi_3$. **Bi-dimensional**
**experiments.** Figure 3.7 depicts the results of the 2D experiments on the ER system for formulae $\varphi_4, \varphi_5, \varphi_6$. The triangular profile of the joint posterior in experiments for $\varphi_4$ and $\varphi_5$ (computed with $k_3=1$ and $\pi_{k_1}(\cdot), \pi_{k_2}(\cdot) \sim U(0,100)$) indicates that only very low values of $k_1$ ($k_1 \le 0.015$ for $\varphi_4$, $k_1 \le 1$ for $\varphi_5$) combined with rather high-values of $k_2$ (i.e. $k_2 \in [40, 100]$ for $\varphi_4$, $k_2 \in [50, 100]$ for $\varphi_5$) result in trajectories entering R4, resp. never leaving R5, which means that the algorithm managed to catch the correlation between the parameters. This is intuitively correct in both cases, in fact R4 corresponds to a very low synthesis of P which is not compatible with fast creation of the ES complex (i.e. only very small $k_1$ are not ruled out) and even the compensation effect obtained by fast decomplexation (i.e. large $k_2$) will not suffice for trajectories to stay in R4. Similarly, R5 limits the speed of the initial decrease of E (which initially is $E_0=100$), to 50 within $t \le 0.8$, which again is compatible only with slow ES complexation and cannot be compensated by fast decomplexation. The R6 experiment caught an even more important correlation between parameters, in fact the posterior for $\varphi_6$ is contained in that for $\varphi_5$ (which is expected because if a trajectory verifies $\varphi_6$ then it also verifies $\varphi_5$).

**Algorithm selection remarks.** The selection of either version of the Automaton-ABC algorithm, i.e. the Automaton-ABC rejection sampling or Automaton-ABC-SMC algorithm, depends very much on the computational complexity of the estimation problem, which in turns depends on the proportion between the support of the posterior distribution $\pi_{\varphi_i-ABC}(\cdot)$ to be estimated and the support for the considered prior distribution $\pi(\cdot)$, i.e. $|\pi_{\varphi_i-ABC}(\cdot)|/|\pi(\cdot)|$, where we denote $|\pi_{\varphi_i-ABC}(\cdot)|$ (resp. $|\pi(\cdot)|$) the size of the support of $\pi_{\varphi_i-ABC}(\cdot)$ (resp. $\pi(\cdot)$). Since there is no way to predict $|\pi_{\varphi_i-ABC}|$ the only way to proceed is by empirically monitoring the runtime of the considered experiment settings. For the mono-dimensional experiments, i.e. those concerning formulae $\varphi_1, \varphi_2$ and $\varphi_3$, we noticed no notable differences of performance between Automaton-ABC rejection sampling (Algorithm 3) and Automaton-ABC-SMC (Algorithm 4): both can be applied and results in comparable runtime. This is because of the fairly large ratio $|\pi_{\varphi_i-ABC}(\cdot)|/|\pi(\cdot)|$ between the support of the obtained posterior distributions, and the size of the sampling support given the prior for the single parameter is $k_3 \sim U(0,100)$. Conversely for the bi-dimensional experiments ($\varphi_4, \varphi_5$ and $\varphi_6$), we point out that Automaton-ABC rejection sampling algorithm is not worth using because of its extremely large runtime. In fact for those experiments the support for the prior distribution is $[0, 100]\times[0, 100]$ while the region with non-null probability of satisfying the considered formula corresponds with, in the worst case (i.e. the support of $\pi_{\varphi_4-ABC}$ in Figure 3.7 top left), the area of a triangle with sides roughly 90 and 0.03. Therefore the probability of sampling a bi-dimensional parameter ($k_1, k_2$) in the resulting posterior-distribution is roughly $|\pi_{\varphi_4-ABC}(\cdot)|/|\pi(\cdot)| \approx \frac{90 \cdot 0.03}{2} \cdot \frac{1}{100 \cdot 100} \approx 10^{-4}$. This results in an infinitesimal probability of drawing from the prior $N=1000$ particles that fall in such a narrow distribution let alone the fact that even a parameter sampled in the obtained distribution could produce paths that don't satisfy $\varphi$. By adding several transitional steps with Algorithm 4, the problem becomes treatable. The results for $\varphi_4$ required about $2 \times 10^5$ simulations of the model (i.e. the highest number of simulations
for experiments on the ER model).

**Further case studies.** We tested the Automaton-ABC method on other systems including the popular SIR model of infection spreading and a model of intracellular viral infection (see [BBC21]).

## 3.4 Perspectives

We introduced a novel approach for parametric verification of temporal logic specifications against probabilistic models based on an adaptation of the ABC parameter inference scheme. As the original ABC scheme is based on a notion of *distance* (between the system's observations and an instance of the model's parameters) the parametric verification extension relies on the adaptation of the distance notion to the *satisfiability distance*, which quantifies how far a model's instance is from satisfying a given temporal logic condition. We have then shown that dedicated hybrid automaton monitor can be defined and plugged in within an adapted ABC scheme in order to assess on-the-fly the satisfiability distance of a temporal property, yielding an approximation of the satisfaction probability distribution for the considered property. The presented approach suffers from lack of generality as, in its current formulation, the satisfiability distance is custom-defined and restrained to non-nested time-bounded temporal specifications, rather than covering the entire syntactical spectrum of a temporal logic formalism as done by Donzé and Maler in similar approaches that target the robustness of MITL formulae against deterministic (ODEs) models [DM10]. On the other hand the fact we employ hybrid automata to monitor the paths issued by a model's instance allows us to widen the realm of parametric verification problems by taking into account much more sophisticated behaviours such as e.g. that of tuning of stochastic oscillators, as discussed in Chapter 4.

One aspect worth investigating in order to evolve this contribution is how to generalise the satisfiability distance adaptation of the ABC scheme to the entire MITL spectrum. Intuitively that would entail working out a reformulation of the satisfiability distance definition based on the recursive syntax of MITL temporal formulae, similar to Donzé and Maler's approach to tackle robustness of MITL properties [DM10]. Given a recursive definition of satisfiability distance can be obtained then it'd also necessary to figure out how it can be conveniently encoded in hybrid automata terms.