Maria Climent- Pommeret

Nicolas Dandrimont

Jonathan Lajus

Loïc Hélouet

Leïla Kloul

Émilie Baudet

Gaëlle Carrette

Romain Gatelet

Mamadou Gaye

Agathe Herrou

Arnaud Martinet

Sophie Pires Doliget

Théo Marchal

Jacqui Samlone

Vincent Lejeune

François Verbeck

Antoine Roche

Louise Xie

Emmanuel Arrighi

Peggy Boulay

Hugues Moretto-Viry

Francis Hulin-Hubard

Linda Zhao

Jean Van

MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS

This PhD thesis could not have been achieved without the help of

Introduction iii I Context

INTRODUCTION

Context

The population growth in most countries, alongside with the increased distance from home to activities, have led to an ever increasing demand in transportation systems. While train and steam buses exist from the industrial age, the invention of personal cars has drastically changed the behaviour of the travelers. More precisely, with the introduction of the Ford Model T in 1908, automobiles became affordable for the middle-class, introducing a better option than public transportation. Unfortunately, the increasing car traffic has triggered some issues: driving induces strain, the increasing number of accidents is a major preoccupation for public safety, and the circulation conditions are degrading. In order to address these issues, a range of advanced driver-assistance systems have been introduced thanks to breakthrough in computation power and electronics. Today, smartphones and mobile networks enable communication while traveling. However, drivers have to focus on the road circulation. So the next step is to make personal cars capable of performing the entire process of of going to some destination: the so-called autonomous vehicles (or self-driving cars). It could even be a part of public transportation, such as the planned night shuttle between Massy and Polytechnique [START_REF] Chevallier | Essonne : bientôt des navettes sans chauffeur entre Massy et le plateau de Saclay[END_REF]. Waymo and Tesla are known to be pioneers in fully driver-less vehicles ( [START_REF] Hawkins | Waymo is first to put fully self-driving cars on US roads without a safety driver[END_REF], [START_REF] Reuters | Tesla's Autopilot to get 'full self-driving feature[END_REF]), but it is now a concern of all car manufacturers. Obviously, these companies want to sell these new vehicles to the general public. But this requires ensuring that such vehicles do not increase the risk of accidents, nor pose a threat to the life of their owners. Since the question of safety has an high societal impact and is difficult to solve, most companies want to let their vehicles drive enough kilometers, in order to have sufficient data and design improvements. However, this comes with a huge legislation problem: they must have the right to drive those possibly dangerous cars. In France, the autonomous vehicles are limited to small portions of motorways [START_REF] Rfi | Emmanuel Macron veut faciliter les tests de véhicules autonomes[END_REF], which can only be a benchmark for Trafic Jam Controllers. In the United States, the companies started to test their iii INTRODUCTION autonomous vehicles in California and Texas, but this had led to several accidents and the public opinion became uncomfortable with the matter.

Similarly to the crash tests for the safety of drivers and passengers during a crash (which is now entirely done, in Europe, via EURO NCAP test benches [START_REF]EUROpean New Car Assessment Programme[END_REF]), this requires a shift to more abstract verification methods. Such approach has two advantages: first, it does not take as much time nor money to get the target precision on the crash rate (10 -9 crashes every 100km driven), and most importantly it does not involve taking any vital risk, and formal verification offers better guarantees than testing [START_REF] Berard | Systems and Software Verification: Model-Checking Techniques and Tools[END_REF].

Related works. Since an autonomous vehicle is very complex by nature, any direct verification task is extremely difficult. Hence, researchers have defined and tackled many subproblems, like: the effects of weather perturbation on sensors [START_REF] Rasshofer | Influences of weather phenomena on automotive laser radar systems[END_REF], the detection of road signs [START_REF] Saadna | An overview of traffic sign detection and classification methods[END_REF] and other entities [START_REF] Sun | On-road vehicle detection: a review[END_REF], inter-vehicle communication [START_REF] Sichitiu | Inter-vehicle communication systems: a survey[END_REF], human-vehicle interaction [START_REF] Matthews | Intent communication between autonomous vehicles and pedestrians[END_REF], design of controllers for road parts [START_REF] Bérard | Collision avoidance in Intelligent Transport Systems: towards an application of control theory[END_REF] and for the single vehicle.

In this thesis, we focus on the verification of such controllers. Several related works will be presented in Chapter 2. They can be classified according to two criteria: • Geometrical versus logical. In the former case, the behaviour of the vehicles in a small time interval is modeled by the set of all positions they occupy, using differential equations, and a check is performed on the interserction of these sets; In the latter case, the expected properties on vehicles and roads are formalised as logical constraints, which are then solved by existing tools. • Synthesis versus verification. In the former case, a controller satisfying the expected properties is automatically generated. In the latter case, the expected properties are checked on an already existing controller.

Objectives

This PhD thesis takes place in the IRT SystemX project Simulation pour la Sécurité du Véhicule Autonome [START_REF]Simulation of Autonomous Vehicle Safety -IRT SystemX[END_REF], whose aim is to provide simulation methods ensuring the safety of autonomous vehicles. This project includes as partners various automobile manufacturers such as PSA Peugeot Citroën and Renault, but also suppliers such as Continental and Valeo. It is also in partnership with software developers, such as ANSYS, Oktal, Apsys-Airbus. In this thesis, we propose a statistical modelchecking based approach which could be used by automobile manufacturers who design controllers of autonomous vehicles. More precisely, our objectives are as follows.

Part

Formal Specifications of Systems

The first step of the process shown in Figure 1.1 consists in specifying a formal model of the system expressed with some formalism. So here, we present some of them. Transition Systems which is the basic one is not used for design. However it is provided as a semantic for more concise formalisms and for defining the truth value of a formula. Still at the semantical level, Discrete Event Stochastic Processes may be viewed as an extension of Transition Systems by integrating probabilities and time. At a syntactical level, we introduce Petri Nets and their extensions which are widely used for concurrent and distributed systems. Finally, we present Hybrid Systems which allow to describe a continuous dynamics and thus to express cyberphysical systems.

Transition Systems

As discussed above, transition systems is the standard low-level formalism. It mainly consists of a set S of states (which may be uncountable) and a transition relation →⊆ S × S describing state changes. We first present the particular case of Kripke structures, which have a finite number of states equipped with labels representing atomic properties that hold in the state.

Definition 1 (Kripke structure)

Let AP be a set of atomic propositions. A Kripke structure is a 4-uple K = (S, I, R, Σ, L) where:

• S is a finite set of states;

• I ⊆ S is the subset of initial states;

• R ⊆ S × S is a transition relation satisfying:

∀s ∈ S : ∃s ∈ S : (s, s ) ∈ R

• Σ = 2 AP is the set of subsets of atomic propositions;

• L : S → Σ is a labelling function

Remark. The restriction on the transition relation means that for every state there always exists a transition starting from that state. In fact this restriction can be easily overcomed by adding a special idle state with a loop and such that for all deadlock state (i.e. a state without outgoing transitions) one adds a transition to the idle state. Moreover a special property, say deadlock, would only hold in the idle state in order to reason about it in formulas. Often, one also includes in the definition of a transition system a subset of initial states S 0 which is sometimes a singleton {s 0 } . Example 1. A Kripke structure modeling a (very simple) coffee machine is shown in Figure 1.2; the labels on arcs are only provided here for readability of the example, but are not a part of the Kripke structure. The atomic property p means that the machine is busy and the atomic property c means that the coffee is filling the cup. If there is no fairness condition, this machine may prepare coffee indefinitely. An execution or run of a Kripke structure K is an infinite sequence of states σ = s 0 , . . . , s n , . . . such that, for every i ∈ N : (s i , s i+1 ) ∈ R. This run generates an (infinite) word w = L(s 0 )

• • • L(s n ) • • • ∈ Σ ω , called the trace of σ.
The language of K, denoted by L(K), is the set of words generated by runs of K.

Remark. When S 0 is specified one restricts the runs to those starting in S 0 .

Example 2.

∅ ∅ ∅ {p, q}

Discrete Event Stochastic Processes

We extend transitions systems by adding probabilities and time to the formalism.

Notation 3. We denote by P (e) the probability of an event e, and by P (e|e ) the conditional probability of e given e has happened. We denote by E(X) the expectation of the (numerical) random variable X. Given a set of basic events E equipped with a σ-algebra, we denote by dist(E) the set of distributions over E.

Remark. Observe that the word event is overloaded. Within a σ-algebra it denotes a measurable subset while in a transition system it denotes what triggers the state change. In the sequel, the context should clarify which notion is used.

Behaviors for stochastic discrete event systems can be modeled via two families of random variables: • (S n ) n∈N taking values in the system state space. The random variable S 0 represents the initial state of the system, and S n the state after the n th event. • (T n ) n∈N taking values in R + . The random variable T 0 corresponds to the time elapsed before the first event, and T n to the time between the n th and the (n+1) th event. Note that there may be zero delay between two events, like in the case of a sequence of instantaneous instructions. Note that this definition implies that events are instantaneous, and that there might be a null delay between two events. Example 3. state time q 3 q 4 q 6 q 7 q 12 τ τ 0 τ 0 + τ 1 The figure above illustrates a discrete event system where the first realizations of its random variables are defined as follows:

• S 0 = q 4 , S 1 = q 4 , S 2 = q 6 , S 3 = q 3 , S 4 = q 12 et S 5 = q 7 ; • and T 0 = τ 0 , T 1 = τ 1 , T 2 = T 3 = 0.

According to physical constraints, the only restriction for these families of random variables is to forbid Zeno phenomenon, where an infinite number of actions can occur in a finite duration. This is mathematically equivalent that almost surely ∞ n=0 T n = +∞. Example 4. In this section, we use the Tandem Queue Network (TQN) illustrated in Figure 1.5 as a running example. This system consists of two queues: the event in represents the arrival of a client in the first queue; move indicates the move of a client from the first to the second queue; finally, out represents a client leaving the system. While the first event could happen at any moment, the other events require the presence of at least one client in their corresponding queue. have not yet specified its probabilistic behavior, this high level model leads to the transition system shown in Figure 1.6 which a countable number of states, each state being specified by the numbers of clients in the two queues.

0,0 0,1 0,3 0,3

. . . For modeling purposes, it is useful to have a syntactical formalism for which the semantics are described by the two random variable families (S n ) n∈N and (T n ) n∈N previously defined.

Definition 4 (Discrete Event Stochastic Process)

A Discrete Event Stochastic Process (DESP) is a tuple D = (S, π 0 , E, Ind, enabled, delay, choice, target) where: • S is a (possibly infinite) set of states; • π 0 ∈ dist(S) is the initial distribution over states; • E is a finite set of events; • Ind is a set of state indicators, i.e. functions S → R that may be constants; • enabled : S → P(E) indicates the set of enabled events for each state, satisfying: for all s ∈ S, enabled(s) = ∅; • delay : S × E → dist(R + ) is a mapping defined for pairs (s, e) such that s ∈ S and e ∈ enabled(s); • choice : S × P(E) × R + → dist(E) is a mapping defined for triples (s, E , d) such that E ⊆ enabled(s) and the possible outcomes of the distribution are restricted to e ∈ E ; • target : S × E × R + → S is a mapping indicating the state changes triggered by events, defined for triples (s, e, d) with e ∈ enabled(s).

The set enabled(s) is the set of events enabled from state s. The delay function, when applied to a pair (s, e) gives the distribution of the delay between the activation of e in state s and its possible occurrence. The function choice gives an additional probabilistic distribution in the case of several events having the same due date, to randomly solve conflicts. The function target denotes the (random) destination state given the source state, the occurrence of an event and the sojourn time in the source state. The set of functions Ind are related to the observed behavior of the model and will be used for specifying formulas.

Example 5. In the case of the tandem queues (see Figure 1.5), the set of states is N × N, π 0 = Dirac ((0, 0)) meaning that the initial distribution is concentrated on state (0, 0), E = {in, move, out}. As discussed above,

enabled(s) =          {in}
when s = (0, 0), {in, out} when ∃n ∈ N + such that s = (0, n), {in, move} when ∃n ∈ N + such that s = (n, 0), {in, move, out} otherwise.

Let event in follow an exponential law of parameter λ denoted by exp(λ) with density x → λe -λx and events move and out follow exponential laws of respective parameters µ 0 and µ 1 , we have, for s ∈ S and e ∈ enabled(S):

delay(s, e) =      exp(λ)
if e = in, exp(µ 0 ) if e = move, exp(µ 1 ) if e = out.

10CHAPTER 1. MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS

The function choice is irrelevant as we use continuous laws: the simultaneous occurrence of two events has null probability. We might for example use choice(s, E , t) = U(E ) (a uniform distribution). Finally, as can be seen in if e = out.

Observe that here the target function does not depend on the sojourn time.

Definition 5. A configuration of a DESP is a triple (s, τ, sched) where s is the current state, τ ∈ R + the current time and sched : E → R + ∪ {∞} the function describing the times at which each enabled event will occur (+∞ if the event cannot be activated).

Informal semantics of a DESP. Starting from a given configuration (s, τ, sched) of a DESP, its dynamics corresponds to an infinite loop, where each iteration consists of the following steps: 1. Using sched, one get the set E = {e ∈ enabled(s) | ∀e ∈ enabled(s) : sched(e) ≤ sched(e )} the set of enabled events with minimum delay. The corresponding delay d is then sched(e)-τ for all e ∈ E . If E contains several elements, we then use a random variable of distribution choice(s, E , d). The corresponding sampling provides an event e, which will be activated. 2. The next state s is then defined by s = target(s, e, d) and sched is updated: for all event e = e still enabled, the value of sched(e ) is kept. For any other event e that is enabled in s , we sample a new delay d using delay(s , e ) and then sched(e ) = τ + d + d . For all disabled event e in s , sched(e ) = +∞. The initial configuration is obtained by sampling s using the distribution π 0 , τ is the initial time (τ = 0) and sched is obtained by sampling every event delay of enabled(s).

A particular case of DESP corresponds to the classical formalism of Discrete Time Markov Chain (DTMC). Definition 6. A square matrix P ∈ R S×S is stochastic (or a transition matrix) if all its coefficients are non-negative and the sum of coefficients over each line is equal to 1, i.e.: ∀i, j ∈ S, P(i, j) ≥ 0 and j∈S P(i, j) = 1

Definition 7 (Discrete-Time Markov Chain)

A sequence (X k ) k≥0 of random variables with values in S is a Markov chain of transition matrix P and initial law π 0 if:

• P (X 0 = i) = π 0 [i] for every i ∈ S;

• for all k ∈ N and i ∈ S:

P (X k+1 = i | X k = i k , . . . X 0 = i 0 ) = P (X k+1 = i | X k = i k ) = P(i n , i).
When the initial distribution is a Dirac distribution, the initial state is usually denoted s 0 . The sequence (X k ) k∈N is a stochastic process. This definition highlights the characteristic property of Markov chains: the future evolution of the process is independent from its past evolution; it only depends on its current state.

Definition 8 (Graph of a DTMC )

A DTMC is represented by a weighted directed graph:

• The set of vertices is S.

• Ther is an edge i → j with weight P(i, j) if P(i, j) > 0 .

In addition to the interest of a graphical representation (see Figure 1.7) which may be more readable for small DTMC, the analysis of such graphs leads to characterizations of properties of the states like transience, recurrence, periodicity, etc. Implicitly in DTMC, an event occurs every time unit which implies that w.r.t. time elapsing DTMC are not memoryless. A contrario Continuous Time Markow Chains (CTMC) are "fully" memoryless since the time between two events follows an exponential distribution whose rate only depends on the current state.

Definition 9 (Continuous Time Markov Chain)

Two sequences (X k ) k≥0 and (T k ) k≥0 represent a CTMC of transition matrix P, initial distribution π 0 and rate {λ i } i∈S if: • The sequence (X k ) k∈N is a DTMC with transition matrix P and initial distribution π 0 ; • for all k ∈ N and τ ∈ R + :

P (T k ≤ τ |X k = i k , . . . , X 0 = i 0 , T 0 ≤ τ 0 . . . T k-1 ≤ τ k-1 ) = P (T k ≤ τ |X k = i k ) = 1 -e -λ i k τ 12CHAPTER 1.

MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS

The DTMC associated with the CTMC is called the embedded DTMC. The graph of a CTMC is the one of embedded DTMC where the states have for labels their rate (see Figure 1.8).

Petri Nets

Carl Adam Petri, a German mathematician, has defined during his PhD [START_REF] Adam | Kommunikation mit Automaten[END_REF] a formalism for describing relations between conditions and events, thus modeling the behavior of distributed discrete event systems. In this section, we give a short presentation of this model, which is now called Petri net, and of some of its extensions.

Definition 10 (Petri net)

A Petri net is a 5-tuple N = (P, T, W -, W + , m 0 ) where P ∩ T = ∅ such that: • P is a non-empty finite set of places; • T is a non-empty finite set of transitions; • W -: P × T → N is the forward incidence function;

• W + : P × T → N is the backward incidence function; • m 0 ∈ N P is the initial marking of the net.

A Petri net is represented by a bipartite graph where: • vertices are places depicted as circles and transitions depicted as rectangles; • there is an edge from place t to transition t labelled by

W -[p, t] if W -[p, t] > 0;
• there is an edge from transition t to place p labelled by

W + [p, t] if W + [p, t] > 0;
• Inside place p, one writes either the number m 0 (p) of tokens or represent this number (when enough small) by bullets. When the label of an edge is missing, it means that this label is 1. The net shown in Figure 1.9 models the tandem queue where the marking of p 1 (respectively p 2 ) indicates the number of clients in the first (respectively second) queue of the tandem.

Example 6.

A marking of a Petri net N is a mapping m ∈ N P that associates an integer with each place of the net.

A transition

t ∈ T is firable in a marking m if ∀p ∈ P : m(p) ≥ W -(p, t).
The firing of such a transition from m leads to a new marking m defined by: ∀p ∈ P : Example 7. In Figure 1.9, the initial marking is m 0 = (0, 2) and a possible sequence of transitions would be:

m (p) = m(p) -W -(p, t) + W + (p, t).
(0, 2) in -→ (1, 2) out --→ (1, 1) move ---→ (0, 2)
Transitions in and out are firable from initial marking m 0 while move is not.

It is possible to add (quantitative) time in Petri nets, with two main approaches: • either Time Petri nets [START_REF] Berthomieu | Modeling and Verification of Time Dependent Systems using Time Petri Nets[END_REF] in which transitions are equipped with intervals constraining their firing; • or Timed Petri nets [START_REF] Parosh | Timed Petri Nets and BQOs[END_REF] where clocks are added to tokens, and edges are equipped with guards.

The first approach can be extended to obtain stochastic Petri nets, by defining a probabilistic distribution over the time interval of the transition delay. In addition, one has to solve conflicts between simultaneous occurrences of transitions and this is done with the help of priorities and weights.

Definition 11 (Stochastic Petri net)

A stochastic Petri net is a 8-tuple N = (P, T, W -, W + , m 0 , Λ, Π, w) where:

• (P, T, W -, W + , m 0 ) is a Petri net; • Λ : T → dist(R)
is a function that associates a distribution with each transition; • w : T → R a function that associates a weight with each transition; • Π : T → N a function that associates a priority with each transition.

14CHAPTER 1. MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS

Remark. This definition could be extended by considering Λ as a mapping from markings to rates.

Stochastic nets can be seen as a particular case of DESP with states N P , the initial distribution over states π 0 = Dirac(m 0 ), no indicator, and: d) is a distribution defined by the following procedure: select the subset T max = arg max{Π(t) | t ∈ T } and then sample some t in T max with distribution

• enabled(m) = {t ∈ T | ∀p ∈ P : W -(p, t) ≤ m(p)}; • delay(m, t) = Λ(t); • choice(m, T ,
w(t) t ∈Tmax w(t ) ; • and target(m, t, d) = m with m (p) = m(p) -W -(p, t) + W + (t, p) for all p in P .
We now present some high level Petri nets, in which tokens are labelled with information and transition can be fired in multiple ways Let us first recall the a natural extension of sets: multisets. A multiset may contain multiple occurrences of the same element.

Definition 12.

A multiset is a pair (S, µ) where S is a set and µ : S → N is a mapping from the elements of the set to natural numbers. Notation 13. We denote by Bag(S) set of multisets using the support set S. A multiset µ is usually written as a symbolic sum s∈S µ(s) • s with term omitted when µ(s) = 0 and factor omitted when µ(s) = 1. In order to define the incidence functions W -and W + , it is necessary to introduce a typed variable per color class occurences of the cartesian product of Domain(t). They are used to build expressions with values in Bag(Domain(P )), For instance, with Domain(T ) = C 1 × C 1 × C 2 and Domain(P ) = C 1 × C 2 such a mapping could be:

C 1 × C 1 × C 2 → Bag(C 1 × C 2 ) < x 1,1 , x 1,2 , x 2,1 > → 2• < x 1,1 , x 2,1 > + < x 1,2 , x 2,1 > whose type of variable x i,j is C i .

16CHAPTER 1. MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS

The incidence functions of this example are:

W -(Get, P 1 ) = Processes × Ressources → Bag(Processes) < x, y > → < x > W -(Get, P 2 ) = Processes × Ressources → Bag(Processes) < x, y > → < y > W + (Get, P 3 ) = Processes × Ressources → Bag(Processes × Ressources) < x, y > → < x, y >
Stochastic Petri nets and high-level Petri nets can be combined to define stochastic high-level Petri nets which is indeed the input formalism of the tool Cosmos used in our case studies.

Hybrid systems

Hybrid systems [START_REF]Hybrid systems[END_REF] combine continuous dynamics, i.e. evolution of variables according to flow functions (possibly described by differential equations) in control locations, and discrete jumps between these locations, equipped with guards and variable updates.

Example 10. A classical example is the thermostat depicted in Figure 1.11, controlling a heater by turning it on or off to adjust the room temperature represented by variable x, maintaining it between x min and x max . When the heater is off, the temperature decreases according to differential equation ẋ = -Kx and when the heater is on, it increases according to ẋ = K(x h -x), where x h and K are parameters for the heater and the room.

ON ẋ = K(x h -x) x ≤ x max OFF ẋ = -Kx x ≥ x min x = x max stop x = x min start Figure 1.11: Room temperature controller
Formally, hybrid automata over R can be defined as follows:

Definition 15 (Hybrid automaton)

A hybrid automaton of dimension n is a tuple H = (S, Σ, E, Dyn, Inv, G, R),

where:

• S is a finite set of states;

• Σ is a finite alphabet of events;

• E ⊆ S × Σ × S is a finite set of edges;

• Dyn is the dynamics assigning to each state a set of differential equations over R n ; • Inv assigns to each state a subset of R n called invariant; • G assigns to each edge a subset of R n called guard;

• R assigns to each edge a subset of {1, . . . , n} called reset.

In the example above

, n = 1, S = {OFF, ON}, Σ = {stop, start}, E = {(OFF, start, ON), (ON, stop, OFF)}, Dyn(OFF) = { ẋ = -Kx}, Dyn(ON) = { ẋ = K(x h -x)}, Inv(OFF) = {x ∈ R | x ≥ x min }, Inv(ON) = {x ∈ R | x ≤ x max }, G(OFF, start, ON) = {x min }, G(ON, stop, OFF) = {x max } and R(OFF) = R(ON) = ∅.

Definition 16

The semantics of hybrid automaton H is defined as the transition system

T H = (Q, →) where: • the set of configurations is Q = S × R n ;
• the transition relation contains two types of steps: For this expressive formalism, where the associated transition system has an uncountable state space, most verification questions are undecidable and in particular the reachability problem. Analysis results were obtained for subclasses, usually by building a finite abstraction based on some bisimulation equivalence, preserving a specific class of properties, like reachability or those expressed by temporal logic formulas.

-discrete step: (s 1 , y 1 ) → (s 2 , y 2 ) iff there exists e = (s 1 , a, s 2 ) in E with y 1 ∈ G(e) and y 2 [i] = 0 for i ∈ R(e) and y 2 [i] = y 1 [i] otherwise. -continuous step: (s 1 , y 1 ) → (s 2 , y 2 ) iff s 1 = s 2 and
Observations. While this definition seems rather general, it suffers from some limitations. For instance, the specification of the differential equations only depends on the current state and not on the whole history. In Chapter 3, a more general formalism based on Simulink will be detailled.

There are several ways to introduce probabilities in such systems: either by associating weights with discrete transitions or associating random delays with states. In Chapter 4, we propose another way to do so.

Formal Specifications of Properties

Generally properties are specified using some logic. In the case of dynamical systems temporal logics that express properties related to runs of sytems are the most appropriate ones. Following our presentation of formalisms, we begin by temporal logics for discrete event systems and then we introduce some of their extensions introduced to capture features of stochastic and/or hybrid systems.

Linear Temporal Logic

The aim of temporal logics is to express properties over system executions. Two approaches are possible: linear logics like LTL (Linear Temporal Logic) which consider sequences independently from each other, and tree logics such as CTL (Computation Tree Logic) which consider the execution tree and thus can be used to reason over branchings.

We first focus on LTL, for which a formula is interpreted over a system execution trace.

Definition 17 (Syntax of LTL). Given a set AP of atomic propositions, LTL formulas are defined by the following grammar:

ϕ ::= a | ϕ 1 ∨ ϕ 2 | ¬ϕ | X ϕ | ϕ 1 U ϕ 2 , a ∈ AP
This logic is an extension of propositional logic; we first remark that by using De Morgan rules, it is possible to redefine the ∧ operator.

(true) can be defined as a ∨ ¬a and ⊥ (false) as a ∧ ¬a. The X and U modalities have the following interpretations: • X ϕ means that at the next moment, ϕ holds; • ϕ 1 U ϕ 2 means that ϕ 2 will be true at some point in the future, and ϕ 1 holds until then. More formally:

Definition 18 (Semantics of LTL)

The satisfaction semantics of an LTL formula ϕ over an infinite word w = w 0 . . . w n . . . ∈ (2 AP ) ω and a position i of this word, denoted w, i |= ϕ, is inductively defined by:

w, i |= a ∈ AP if a ∈ w i w, i |= ϕ 1 ∨ ϕ 2 if (w, i |= ϕ 1 ) or (w, i |= ϕ 2 ) w, i |= ¬ϕ if w, i |= ϕ w, i |= X ϕ if w, i + 1 |= ϕ w, i |= ϕ 1 U ϕ 2 if ∃j ≥ i : (w, j |= ϕ 2 ) ∧ (∀i ≤ k < j : w, k |= ϕ 1 )
The classical operators ♦ and are defined by: • ♦ ϕ ≡ U ϕ: the formula ϕ will be satisfied eventually; • ϕ ≡ ¬ ♦ ¬ϕ: at every moment the formula is satisfied (there is no time in the future where ϕ is not satisfied) Moreover, for an LTL formula ϕ and a word w, we say that w satisfies ϕ if w, 0 |= ϕ, written w |= ϕ. We denote by L ϕ the set of words w satisfying ϕ. Definition 19. Let K be a Kripke structure, and ϕ an LTL formula. We say that K |= ϕ if, and only if, for any execution σ of K, the trace w of σ satisfies ϕ.

Remark. Said otherwise, K |= ϕ if and only if

L(K) ⊆ L ϕ .
Similarly to regular languages (and finite automata), it is possible to characterise the language L ϕ with a class of automata called Büchi automata.

Definition 20 (Büchi Automaton)

A Büchi automaton is a 6-uple B = (S, T, I, F, Σ, L) with: • S a set of states, and T ⊆ S × S a set of transitions; • I ⊆ S the set of initial states, F ⊆ S the set of recurring states; • Σ an alphabet and L : S → Σ a labelling function.

Definition 21.

An execution σ of a Büchi automaton is an infinite sequence of states s 0 , . . . , s n , . . . such that, for each i ∈ N : (s i , s i+1 ) ∈ T . The trace of σ is the word w = L(s 0 )...L(s n ).... An execution is accepting if the set {i ∈ N | s i ∈ F } is infinite (it goes infinitely often in a recurring state). A word w is recognized by B if there exists an accepting execution σ with trace w. We denote by L(B) the set of words recognized by B.

A usual question on languages and automaton is the emptiness problem.

Proposition 1. Let B be a Büchi automaton. It is possible to check in linear time if the language recognized by B is empty. Furthermore, the emptiness problem is NLOGSPACE-complete.

There exist several algorithms answering this question; a comparison between some of them has been published by Andreas Gaiser and Stefan Schwoon [START_REF] Gaiser | Comparison of algorithms for checking emptiness on buechi automata[END_REF].

Let AP be a set of atomic propositions. We extend Definition 20 with a labelling function using the alphabet P(P(AP )), but where words are still considered using P(AP ).

Definition 22 (Extension of definition 21)

. Let Σ = P(AP ) and let B be a Büchi automaton over alphabet P(Σ). Let w ∈ Σ ω and let σ be an accepting execution of B. Then w is recognized by σ if and only if: ∀i ∈ N :

w i ∈ L(s i ).
Using this extension, we can now state the following theorem, with L ϕ the set of words w satisfying ϕ recognized by a Büchi automaton:

Theorem 2

Let ϕ be an LTL formula, over a set of atomic propositions AP . Then, we can build in exponential time a Büchi automaton B ϕ over the alphabet 2 2 AP which exactly recognises the set of words satisfying ϕ:

L(B ϕ ) = L ϕ .
Several algorithms perform the construction from ϕ to B ϕ , notably Paul Gastin and Denis Oddoux's one [START_REF] Gastin | Fast LTL to Büchi Automata Translation[END_REF].

Remark.

In order to shorten state labelling notations and allow for the use of the transformation of [START_REF] Gastin | Fast LTL to Büchi Automata Translation[END_REF], we use boolean formulas over AP : we recursively define a valuation function Val : Bool(AP ) → 2 Σ = 2 2 AP as follows:

Val( ) = 2 AP Val(p) = Σ p = {A ∈ 2 AP | p ∈ A} Val(ϕ ∧ φ) = Val(ϕ) ∩ Val(φ) Val(ϕ ∨ φ) = Val(ϕ) ∪ Val(φ) Val(¬ϕ) = Σ \ Val(ϕ)
For instance, with AP = {p, q, r}, we have Val(p ∧ ¬r) = {{p}, {p, q}}. We use this notation to define sets of P(AP ).

Example 11. Let p and q be atomic propositions. The Büchi automaton shown in Figure 1.12(a) recognizes ♦ p (at each time, there exists a time in the future where p is satisfied), and the Büchi automaton shown in Figure 1.12(b) recognizes p U q.

Computational Tree Logic

We now focus on CTL which considers the system execution tree, and is used to reason over branchings.

Definition 23 (Syntax of CTL). Given a set AP of atomic propositions, CTL formulas are defined by the following grammar: In this case, the temporal operators X and U are in the scope of a path quantifier, E being interpreted as for some path while A means for all paths.

ϕ ::= a | ϕ 1 ∨ ϕ 2 | ¬ϕ | E X ϕ | A X ϕ | E ϕ 1 U ϕ 2 | A ϕ 1 U ϕ 2 , a ∈ AP p ¬p (a) p q (b)
A CTL formula is interpreted over a state of the execution tree of a Kripke structure K = (S, I, R, Σ, L). We write Exec(s) for the execution tree starting from state s ∈ S, and σ = s 0 s 1 . . . with s 0 = s for a path in this tree, with σ(i) = s i for the state at position i.

Example 12.

s 0

s 1 s 0 s 2 s 1 s 0 s 2 . . .

Definition 24 (Semantics of CTL)

The satisfaction of a CTL formula over a state is inductively defined by:

s |= a ∈ AP if a ∈ L(s) s |= ϕ 1 ∨ ϕ 2 if (s |= ϕ 1 ) or (s |= ϕ 2 ) s |= ¬ϕ if s |= ϕ s |= E X ϕ if ∃σ ∈ Exec(s) : σ(1) |= ϕ s |= A X ϕ if ∀σ ∈ Exec(s) : σ(1) |= ϕ s |= E ϕ 1 U ϕ 2 if ∃σ ∈ Exec(s), ∃j ≥ 0 : (σ(j) |= ϕ 2 ) ∧ (∀0 ≤ k < j : σ(k) |= ϕ 1 ) s |= A ϕ 1 U ϕ 2 if ∀σ ∈ Exec(s), ∃j ≥ 0 : (σ(j) |= ϕ 2 ) ∧ (∀0 ≤ k < j : σ(k) |= ϕ 1 )
Like for LTL, the standard operators ♦ and are defined by: • E ♦ ϕ ≡ E U ϕ. A state s statisfies this formula if there exists a path starting from s where ϕ will be satisfied eventually ; We continue the example on the Kripke structure K shown in Figure 1.2. We have:

• K |= A ♦ c: there exists a loop between states s 0 and s 1 , which avoids the only state where c is satisfied, s 2 ; • K |= A X E X c: the single successor of s 0 is s 1 which has for possible successor s 2 satisfying c; • K |= E X A X c: state s 1 has also for possible successor s 0 not satisfying c.

Hybrid Automata Stochastic Logic

Logics LTL and CTL are appropriate for discrete event systems but cannot express usual quantitative timed requirements like a bounded response time for requests. In addition, they cannot express any satisfaction probability which is mandatory for risk evaluation. So we now present stochastic extensions to LTL and CTL that have mainly be introduced for numerical model checking (following the analyses of [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF]) and we end this section by the HASL logic suited for statistical model checking of probabilistic hybrid systems and integrated in the tool Cosmos.

Continuous Stochastic Logic.

In [START_REF] Aziz | Modelchecking continuous-time markov chains[END_REF], Continuous Stochastic Logic (CSL) has been introduced and the decidability of the verification problem over a finite continuous-time Markov chain (CTMC) has been established. CSL extends the branching time reasoning of CTL to CTMC models by replacing the discrete CTL path-quantifiers All and Exists with a continuous path-quantifier that expresses that the probability of CTMC paths satisfying a given condition fulfils a given bound.

Several variants have been introduced and studied such as CSRL [START_REF] Baier | On the logical characterisation of performability properties[END_REF], that extends CSL to take into account Markov reward models, i.e. CTMCs with a single reward on states or possibly on actions [START_REF] Cloth | Model checking Markov reward models with impulse rewards[END_REF], that is used to specify an interval (on the accumulated reward) on the Until or Next operator. asCSL, introduced in [START_REF] Baier | Model checking action-and state-labelled Markov chains[END_REF] replaces the interval time constrained Until of CSL 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. In the logic CSLTA [START_REF] Donatelli | Model checking timed and stochastic properties with CSL T A[END_REF], path formulas are defined by a single-clock deterministic timed automaton. This logic has been shown strictly more expressive than CSL and also more expressive than asCSL when restricted to path formulas.

DTA.

In [START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF], 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.

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. [START_REF] Chen | Time-bounded verification of ctmcs against real-time specifications[END_REF] 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. [START_REF] Zuliani | Bayesian statistical model checking with application to stateflow/simulink verification[END_REF].

Observe that all of the above mentioned logics have been designed so that numerical methods can be employed to decide about 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 → asCSL → CSLTA → DTA), the other targeting performance evaluation capability (evolutionary path: CSL → CSRL → 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 Stochastic Logic. The logic HASL (Hybrid Automata Stochastic Logic) describes properties via two elements:

• a Linear Hybrid Automaton (LHA) which contains numerical variables and is synchronised with the observations of the analyzed stochastic process; • and an expression specifying the quantity to be evaluated.

Example 14.

0 ẋ1 :1 ẋ2 :0 Size<10

1 ẋ1 :1 ẋ2 :1 Size≥10 2 E; x 1 ≤100 E; x 1 ≤100 E; x 1 ≤100 E; x 1 ≤100
; x 1 = 1 0 0 ; r : = x 2 1 0 0 (n, m) → n + m. This formula contains an automaton with three states ( 0 , 1 , 2 ) and an expression E(LAST (r)) which describes the proportion of time during which too many people are in the queues (over 10 simultaneous clients), using a overall duration of 100 time units. This automaton handles three variables:

; x 1 = 1 0 0 ; r : = x 2 1 0 0 E(LAST(r))
• x 1 representing the elapsed time; • x 2 representing the time during which the queues were too full; • and r representing the time proportion over which the queues were too full.

The initial state 0 corresponds to the system running without overfilling. In this state, only x 1 evolves over time. In state 1 , the variable x 2 also evolves, measuring the time during which the queues are overfilled. The measure ends in state 2 , with the maximum simulation time.

There are two types of transitions: those who are fired by a process event, shown in blue, and the autonomous transitions fired by a condition over variables (in this case, the execution time) represented in red.

In the expression E(LAST (r)), the operator LAST is used to evaluate r at the end of an execution, and E the expectation (over successful runs).

We now describe more formally the syntax and semantics of HASL. In the LHA, the transitions are equipped by constraints that form guards, and variable updates.

Definition 25. A constraint is a boolean combination of comparisons of the form

1≤i≤n α i x i + c 0 where x i are variables, α i , c are indicators and ∈ {=, <, > , ≤, ≥}. We denote by Const the set of constraints. Definition 26. An update is a tuple of functions u 1 , . . . , u n such that every u k is of the form u k = 1≤i≤n α i x i + c where α i and c are indicators. We denote by Up the set of updates.

In order to synchronize with the model of the system, we define Linear Hybrid Automata: Definition 27

A Linear Hybrid Automaton (LHA) is a tuple A = (E, L, Λ, Init, Final, X, flow, →) where: • E is a finite alphabet of events; • L is a finite set of locations; • Λ : L → Prop is a location labelling function;

• Init is a subset of L containing the initial locations; • Final is a subset of L containing the final locations; • X = (x 1 , . . . , x n ) is a n-uple of data variables; • flow : L → Ind n is a function which associates with each location one indicator per data variable representing the evolution rate of the variable in this location. The projection of flow on its i th component is denoted by

flow i ; • the transition relation → is a subset of L× (2 E × Const) ({#} × Const) ×
Up × L, where denotes the disjoint union. A transition ( , E , γ, U, ) ∈→ is written

E ,γ,U ----→ .
Furthermore, the following conditions are required:

• Initial determinism : ∀ 1 = 2 ∈ Init, Λ( 1 ) ∧ Λ( 2 ) ⇔ ⊥; • Determinism on events : ∀E 1 , E 2 ⊆ E st E 1 ∩ E 2 = ∅ : ∀ , 1 , 2 , si E 1 ,γ 1 ,U 1 -----→ 1 et E 2 ,γ 2 ,U 2 ) -----→ 2 are two distinct transitions, then either Λ( 1 ) ∧ Λ( 2 ) ⇔ ⊥ or γ 1 ∧ γ 2 ⇔ ⊥; • Determinism on # : ∀ , 1 , 2 ∈ L si #,γ 1 ,U 1 ----→ 1 et #,γ 2 ,U 2 ----→ 2 are two distinct transitions, then either Λ( 1 ) ∧ Λ( 2 ) ⇔ ⊥ or γ 1 ∧ γ 2 ⇔ ⊥; • No #-labelled loops: for any sequence 0 E 0 ,γ 0 ,U 0 -----→ 1 E 1 ,γ 1 ,U 1 -----→ . . . E n-1 ,γ n-1 ,U n-1 ---------→ n such that 0 = n , ithere exists i ≤ n such that E i = #.
There are two kinds of transitions in this definition: • transitions with label #, which are spontaneously fired with respect to the automaton's variable values; • transitions with a label containing a subset of events, which are fired upon synchronisation with the model to evaluated: the occurence of a transition of the model "sends" the corresponding event to the automaton.

A configuration of the hybrid automaton consists of a state and a valuation of its variables.

The synchronised product of the stochastic process of the model and the hybrid automaton is also a stochastic process. The state of the synchronised product is a pair, composed of a configuration of the net and a configuration of the hybrid automaton.

A change of state in this synchronised product is produced either by the firing of a transition of the model (synchronised with an edge of the automaton), or by the firing of an autonomous transition from the current automaton state. As previously, the event to be considered is the earliest one. If a transition firing occurs in the model and cannot be synchronised with the automaton, the process ends in a failure state. The only other way to end this new process is to reach a final state of the automaton.

Example 15. In the example of figure 1.14, the set of events E is {in, move, out}. The two autonomous transitions are used to stop the synchronisation after 100 time units by reaching the final state 2 , whatever the previous state. Conditions on state 0 and 1 distinguish the state in which tandem queues are in: either they are overfull (Size ≥ 10) or there is a normal load. Since ẋ1 = 1 during the whole run of the automaton, variable x 1 measures the elapsed time (and therefore stops the execution when needed). The variable x 2 is only modified in state 2 , to measure the duration of overfull load in the tandem queues.

We now focus on the evaluation of a performance index over a trajectory of the synchronised product between the LHA and the model. This index is obtained via expressions over state variables:

Definition 28

HASL expressions are defined by the following grammar:

Z ::= c | E(Y ) | Z + Z | Z -Z | Z × Z | Z/Z Y ::= c | Y + Y | Y × Y | Y /Y | LAST(y) | MIN(y) | MAX(y) | INT(y) | AVG(y) y ::= c | x | y + y | y × y | y/y
where x is a state variable (in the set X) of the automaton, and c a constant.

In this definition (from bottom to top): • the y expressions are arithmetic expressions over state variables and constants; • path operators are then applied on the previous expressions: minimum, maximum, integral with respect to time, last value. New arithmetic expressions are then built from these values. The result is a random variable Y associated with the random trajectory of the synchronised product.

• The expectation of these random variables form performance indices that we can then combine with arithmetic operations to define new indices. These expectations are actually conditional ones since the trajectories have to be accepted. A variant for these expressions is to consider directly as a performance index, the probability that the trajectory is accepted.

Example 16. To continue again the example of figure 1.14, by construction trajectories always end in the final state 2 (after 100 time units), where we consider the last value of r = x 2 100 , which is set by the transitions leading to 2 . The HASL expression associated with this LHA is used to measure the fraction of time in which the queue was overloaded, for simulations of 100 time units.

Compared to the logics introduce above, HASL is the most expressive one thanks to (1) the selection of trajectories, (2) the existence of multiple indicators, and (3) the arithmetic operators.

Verification Methods

Verifying Kripke Structures

To check whether a Kripke structure K satisfies an LTL formula ϕ, a classical method consists in trying to exhibit a counter-example: more precisely, a particular execution σ, whose trace w satisfies ¬ϕ. We have seen in section 1.2.1 that a Büchi automaton can be used to generate the whole set of words (or executions) satisfying a LTL formula. Using this result, one only checks the emptiness of the synchronous product of B ¬ϕ with the Kripke structure K to prove that there is no such counter example.

The synchronised product of a Kripke structure K and the Büchi automaton B ¬ϕ is also a Kripke structure:

Definition 29 (Kripke structure and Büchi automaton product)

Let AP be a set of atomic propositions. Let K = (S K , I K , R K , L K ) be a Kripke structure (which states are labelled by elements of P(AP )), and B = (S B , T B , I B , F B , P(P(AP )), L B ) a Büchi automaton over the alphabet P(P(AP )). We define the product Büchi automaton P(S, T, I, F, Σ, L) as follows:

• its set of states is

S = {(s, s ) ∈ S K × S B | L K (s) ∈ L B (s )}
The states (s, s ) ∈ S are labelled by L(s, s ) = L K (s).

• its initial states are I = (I K × I B ) ∩ S, and the recurring states are F = (S K × F B ) ∩ S;

• its alphabet Σ is P(AP );

• its set of transitions is defined using:

((s 1 , s 1 ), (s 2 , s 2 )) ∈ T ⇐⇒ (s 1 , s 2 ) ∈ R K ∧ (s 1 , s 2 ) ∈ T B
We denote this product by P = K ⊗ B. Proof : Let Σ = P(AP ) and w be a word recognized by the product Büchi automaton.

There exists an accepting execution σ = (s 0 , s 0 ) . . . (s n , s n ) . . . which recognizes w. We have, for every i ∈ N : ((s i , s i ), (s i+1 , s i+1 )) ∈ T and w = L(s i , s i ) = L K (s i ). Furthermore, using the definition of S, we have w ∈ L B (s i ).

• Using the definition of T , we have for every i ∈ N : (s i , s i+1 ) ∈ R K and by construction of S, s 0 , . . . , s n , . . . is an infinite sequence of states of the Kripke structure K, which means K recognizes w. • Similary, for every i ∈ N: (s i , s i+1 ) ∈ T B and s 0 , . . . s n . . . is an infinite sequence of states of the Büchi automaton. As σ is an accepting execution in the product automaton, the set {i

∈ N | (s i , s i ) ∈ F } is infinite. Moreover, as F = (S 1 × F B ) we have for all t ∈ S 1 : (t, s i ) ∈ F ⇔ s i ∈ F B .
We can deduce that {i ∈ N | s i ∈ F B is also infinite. Finally, using the definition of S, w = L K (s i ) ∈ L B (s i ) which means that w is recognized by the Büchi automaton B

We can show similarly that if the word w is recognized by both automata, then it is recognized by the product automaton.

Proposition 4

Let K be a Kripke structure, and ϕ a LTL formula. Then the problem: «K |= ϕ?» is decidable, with temporal complexity |K| × 2 |ϕ| and belongs to PSPACE.

Proof : This proposition has been proved by Moshe Y. Vardi and Pierre Wolper in 1986 [START_REF] Vardi | An Automata-Theoretic Approach to Automatic Program Verification[END_REF].

Example 17. We want to check if ϕ = (p ⇒ X ♦ q) holds on the Kripke structure shown in Figure 1.3. We first compute ¬ϕ ≡ ¬ (p

⇒ X ♦ q) ≡ ♦(¬(p ⇒ X ♦ q)) ≡ ♦(p ∧ X (¬ ♦ q)) ≡ ♦(p ∧ X ( (¬q))
). Then we build a Büchi automaton that recognizes ¬ϕ and afterwards the product automaton. This operation is described in Figure 1 It is enough to read the product automaton: as it contains an accepting loop of recurring states reachable from the initial state, the formula ϕ is not satisfied over the Kripke structure. Even better, we immediately obtain on the Kripke structure a counter-example: reaching the p, q state and then looping over the three first states of the Kripke structure.

Numerical model-checking

Numerical model-checking of stochastic models has been intensively studied. Here we only give some relevant examples of these techniques and then discuss the limitations of numerical model checking. Let us start with time-bounded reachability for DTMC which requires to reach a subset of states S + in at most u steps. Such a property can inductively defined w.r.t. u as follows.

Definition 30 (Time-Bounded Reachability for DTMC )

Given a DTMC C = (S, S 0 , P), a set of states S + ⊂ S and a positive integer u, we define the vector of probabilities µ u ∈ [0, 1] S as:

µ u (s) =          1 if s ∈ S + 0 if u = 0 ∧ s ∈ S + s ∈S P(s, s )µ u-1 (s ) if u > 0 ∧ s ∈ S +
We now switch to CTMCs, and define µ τ as the vector of probabilities to reach 30CHAPTER 1. MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS an absorbing state s + in τ time units by µ τ (s) = π τ,s (s + ). Using the uniformization method [START_REF] Jensen | Markov chains as an aid in the study of markov processes[END_REF], one assumes w.l.o.g. that the exit rate is independent from the state and equal to some λ. Thus by conditionning on the number of steps during an interval [0, τ ], the time-bounded reachability may be defined as follows.

Definition 31 (Time-Bounded Reachability for CTMC ) Given a CTMC C, a state s + ∈ S and a positive real τ , we define the vector of probabilities µ τ ∈ [0, 1] S by:

µ τ (s) = e -λτ n∈N (λτ ) n n! µ n (s)
In order to generalize it to a CTL-type formulas, one considers the formula as a tree, and evaluates those starting from the leaves. Given a sub-formula consisting of only one probabilistic operator at the root, and no nested probabilistic operators, the evaluation of U , X and formulas can be reduced to time-bounded or time unbounded reachability problems, as sub-formula are states formulas. Once subformulas have been evaluated on each state, they can be replaced by a proposition in the upper sub-formulas, until the root is reached. For CTL-type formulas, one proceeds similarly to section 1.3.1: the formula is translated into an automaton, then the synchronized product of the stochastic system with the automaton is built. This product is still a stochastic system, and the model-checking problem is thus reduced to time-unbounded reachability. However, if the formula is large, the size of the automaton and thus the size of the product may become problematic. The reader may refer to chapter 10 of Principles of Model Checking [START_REF] Baier | Principles of Model Checking (Representation and Mind Series[END_REF] for more details about the model-checking of probabilistic systems.

Most numerical model-checking algorithms rely on matrix vector multiplication, usually taking advantage of the sparsity of transition probability matrices. However, these computations are difficult to parallelize. The tool Marcie [START_REF] Heiner | Marcie -model checking and reachability analysis done efficiently[END_REF] features a parallel implementation of this method, but still a non negligible part of the computation cannot be performed in parallel.

In order to perform numerical model-checking, strong probabilistic hypotheses on the stochastic system to be analyzed are mandatory. The easiest case is to suppose that the system is Markovian and has a finite state space. Note that the system of linear equations produced during the computation can still be huge, and that the computation of solutions may be difficult and therefore require iterative methods. This implies that only approximate solutions of such systems can be obtained.

In the case of using a time automaton for the specification of the formula, there exist efficient techniques only if there is at most one clock in the automaton, that are described in [START_REF] Donatelli | Model Checking Timed and Stochastic Properties with CSL TA[END_REF] and [START_REF] Barbot | Efficient ctmc model checking of linear real-time objectives[END_REF]. More generally, numerical approaches can be adapted to deal with non determinism [START_REF] Vardi | Automatic verification of probabilistic concurrent finite state programs[END_REF] and to tackle model-checking problems on Markov Decision Processes by alternating probabilistic transitions with non-deterministic ones.

This approach has been efficiently implemented in several tools with different types of Markovian systems, such as:

• Prism [START_REF] Kwiatkowska | PRISM: Probabilistic Symbolic Model Checker[END_REF] deals with Discrete and Continuous Time Markov Chains, probabilistic automatas and Markov Decision Processes; • Uppaal [START_REF] Behrmann | Developing UPPAAL over 15 years[END_REF] deals with various types of automata;

• GreatSPN [START_REF] Baarir | The greatspn tool: Recent enhancements[END_REF] and Marcie [START_REF] Heiner | Marcie -model checking and reachability analysis done efficiently[END_REF] deal with Stochastic Petri nets.

Statistical model-checking

The alternative to numerical approach when dealing with large systems is to use statistical methods. The statistical model checking relies on a Monte Carlo algorithm to estimate the probability of interest. More precisely, given a linear temporal logic formula, one defines a random variable X which takes 1 as value when the trajectory satisfies the formula and 0 otherwise. Thus, this variable follows a Bernouilli law, and one computes its expected value E(X). The Monte Carlo algorithm simulates a large number N of trajectories, and counts the number of those trajectories satisfying the specification. An estimation of the probability is obtained as the ratio of the number of successful trajectories on the total number of trajectories. The random variable Z is defined as the mean of N independent copies of X: Z = 1 Let X and Y be two random variables of distribution UNIF(0, 1). Let Z be a random variable such that an observation z is 1 when the realisations x and y satisfy x 2 + y 2 ≤ 1 and 0 otherwise. This definition means that z is equal to 1 if and only if (x, y) belong in the quarter-circle of radius 1 (see Figure 1.16), which has an area of π 4 . Thus E(Z) = π 4 .

UNIF(0,1)

e 1 UNIF(0,1) e 2 (a) ẋ = 1 ẏ = 0 s 0 ẋ = 0 ẏ = 1 s 1 s 2 s 3 e 1 e2 ; x 2 + y 2 ≤ 1 ; r := 4 e 2 ; x 2 + y 2 > 1 ; r := 0 E(LAST(r)) (b)
Figure 1.17: Computing π using a HASL formula and a stochastic net.

A trajectory is obtained by synchronizing the LHA shown in Figure 1.17(b) with the stochastic net shown in Figure 1.17(a), producing random values for x (the firing delay of e 1 ) and y (the firing delay of e 2 ) respectively. The last value r is 4 if x 2 + y 2 ≤ 1, and 0 otherwise. In both case, the last synchronisation lead to a final state. The expression E(LAST (r)) corresponds to the expected value of r at the end of a trajectory, which is equal to π, four times the value of E(Z).

It is easy to parallelize this method: it suffices to run several "simulators" of the system on different processes (either on different processors or machines), and take the mean result of trajectories of all simulators. The main difficulty comes from choosing the right random number generator so that all generated trajectories are independent one from each other. The only sequential operation being the computation of the mean value and the confidence interval. This implies there is almost no additional computing cost in using parallel computations. These statistical methods can be naturally extended to performance evaluation: instead of computing a probability, we compute the expected values of random variables whose values depend on trajectories of the systems.

Statistical procedures

We now detail the different proccedures that are available in Cosmos for the evaluation of HASL expressions: • Sequential hypothesis testing [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF]. This procedure checks whether a probability is above a threshold. Parameters of this procedure are the probability of an error for a positive answer and a negative answer and the width of the indifference region. When the value of the probability is outside the indifference interval, the probability of an error is bounded by the parameter corresponding to the answer.

• Chernoff-Hoeffding bounds [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF]. This static method requires three related parameters, each of them can be determined by the two others. These parameters are the interval width, the confidence level and the number of samples. It outputs a confidence interval whose width satisfies the requirement and where the probabilistic guarantee is exact. It applies to estimate the expectation of a bounded random variable. • Chow-Robbins bounds [START_REF] Yuan | On the asymptotic theory of fixed-width sequential confidence intervals for the mean[END_REF]. This sequential method requires two parameters:

the interval width and the confidence level. It outputs an interval whose width satisfies the requirement and where the probabilistic guarantee is asymptotic w.r.t. the width of the interval. It applies to estimate the expectation of a random variable, when no known bound is available. • Gaussian approximation. This static method requires two parameters. The number of samples has to be given. The second parameter is either the confidence level or the interval width, one of these determining the other one. It ouputs an interval whose width satisfies the requirement and where the probabilistic guarantee is asymptotic w.r.t. the number of samples. It applies to estimate the expectation of a random variable. It is based on the central limit theorem. • Clopper-Pearson bound [START_REF] Clopper | The use of confidence or fiducial limits illustrated in the case of the binomial[END_REF] This static method computes confidence intervals for binomial distributions. It takes as input three parameters, the total number of samples, the confidence level and the number of successful samples and outputs a confidence interval for the probability of a sample to be successful.

Comparison of numerical and statistical methods.

The statistical methods have several advantages compared to the numerical methods:

• the memory required to simulate one trajectory of the system is usually very small (and the possible values of the random variables are rarely large), making the memory requirement of a statistical model-checker very low. A contrario, since models and formulas may introduce a large size increase of either the probability matrix or the model/formula product, memory is a bottleneck of numerical methods; • we have seen in section 1.3.2 that numerical methods are only usable on a very restricted set of models. A contrario, statistical methods only need an operational semantic of the stochastic model; • as already mentioned, it is straightforward to parallelize the simulation process and it comes with no meaningful additional computation cost. A contrario, parallelization is much more difficult for numerical model-checking; • complex linear time properties, expressed for instance using the HASL logic (see section 1.2.3), can be efficiently evaluated using statistical model-checkers while numerical model checkers require restricted temporal logics. However, statistical methods have several drawbacks such as:

• the computation time may be a problem in the case of tight confidence intervals, that one needs for precise results. In general, dividing by two the width of the confidence interval requires performing four time more simulations; • it is not suited for the evaluation of logics based on state formulas; • as it is possible to simulate only finite trajectories, a statistical model-checking procedure cannot evaluate the unbounded until operator φ U ψ; a workaround is described in [START_REF] He | A bounded statistical approach for model checking of unbounded until properties[END_REF]; • a crude Monte carlo estimator is unsuitable to compute very small probabilities, such as the probability of a rare event; this specific problem is widely discussed in [START_REF] Barbot | Acceleration for statistical model checking[END_REF].

Tools for statistical model-checking

In this section, which is only a slight improvement of the Related Tools and Tool Evaluation sections of [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF], we list and compare a certain number of tools that can be used for statistical model-checking. Figures of this section are taken directly from the journal article.

Related tools.

Most of these tools have been briefly mentioned with respect to their numerical model-checking features.

Cosmos1 [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF] is a statistical model-checker for the Hybrid Automata Stochastic Logic over high-level stochastic Petri nets. It has been developed first during Hilal Djafri's PhD [START_REF] Djafri | Numerical and statistical approaches for model checking of stochastic processes[END_REF] at LSV, ENS Cachan. It was improved, adding HASL support in [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF], and rare event handling in Benoît Barbot's PhD [START_REF] Barbot | Acceleration for statistical model checking[END_REF] still at LSV, ENS Cachan. As several improvements have been made during this PhD, the tool is presented more thoroughly in Chapter 4. It is currently maintained by Benoît Barbot.

Prism2 [START_REF] Kwiatkowska | PRISM: Probabilistic Symbolic Model Checker[END_REF] is a tool for performing model checking on probabilistic models, that has been used for numerous applications [START_REF]Prism -case studies[END_REF]. It can perform numerical modelchecking on discrete and continuous-time Markov chains, Markov decision processes and probabilistic timed automata. Its statistical part only deals with Markov chains. The Prism language defines these probabilistic systems with a synchronised product between reactive modules and can, in this manner, describe large systems in a compact way. The property language subsumes several well-known probabilistic temporal logics, including CSL and PCTL.

Uppaal3 [START_REF] Behrmann | Developing UPPAAL over 15 years[END_REF] is a verification tool that includes many formalisms, such as timed automata, or timed games. It supports automata-based and game-based verification techniques, and has shown its ability to analyse large scale applications. A statistical-model checker engine has been added in 2011 in the form of an extension4 (here called Uppaal-smc), enriching Uppaal with the ability to verify timed systems with stochastic semantics. The specification language, PLTL, is an adaptation of LTL with path operators substitued for quantifiers and with bounded Until ( U ) operator.

Plasma5 [START_REF] Jegourel | A PLatform for High Performance Statistical Model Checking -PLASMA[END_REF] is a platform dedicated for statistical model-checking. It accepts the Prism language for its models, extended with more general distributions and a dedicated biological language. Its property specification language is a restricted version of PLTL, with a single threshold operator. Furthermore, it is built with a plugin system giving the ability for developers to extend it. It can also be integrated, via a library, to another software.

Ymer 6 [74] is a statistical model checker for CTMCs and generalised semi-Markov processes described using the PRISM language. Its property specification language is a fragment of CSL without the steady-state operator but including the unbounded Until ( U ) operator. It seems to be no more maintained.

Marcie • the first model is a Tandem Queuing Network (TQN) taken from [START_REF] Hermanns | Multi terminal binary decision diagrams to represent and analyse continuous time Markov chains[END_REF]. The interarrival time of clients in the tandem queue follows an exponential distribution, the service of the first queue follows a Cox 2 distribution and the service of the second queue follows an exponential distribution. A state of this system is composed of the number of clients in each queue and the state of service of the first queue when non empty. • the second model is a model of dining philosophers.

Parameters.

For the experiments, the following parameters were used: the queue capacity N = 5, the arrival rate λ = 20, the service rates on the first phase in the first queue µ 1 = 0.2 and µ 1 = 1.8 (the latter corresponding to clients without a second service phase) and the service rate of the second phase in the first queue is µ 2 = 2. The service rate of the second queue is κ = 4.

The dining philosopher model is a mutual exclusion problem where N philosophers are sitting around a table. They can decide to eat by taking two forks, that are shared with their right and left neighbours. A contention problem may arise due to the sharing of resources. For the experiments, the rate of all exponential distribution is chosen equal to 10.

Performance indices.

For the tandem queue system, we consider the following time-bounded reachability measure: the probability φ T QS that the first queue of the tandem gets full within time T . For the dining philosophers, it is the probability φ DQM of reaching a deadlock state before N philosophers eat. Such a deadlock occurs when all philosophers have taken one fork. Those properties can be straightforwardly encoded in CSL and in HASL: equivalent verification experimentations can be performed by all tools.

Experiment settings. The following statistical parameters have been set: the confidence level is 0.95, and the width of the confidence interval 0.005. For the hypothesis testing, the probability of error is 0.005 and the width of the indifference region is set to 0.001. The other parameters have been set to default values. The tools must generate a large number of trajectories to be able to fulfill these parameters. Most tools can take advantage of parallelisation; however the experiments have been executed on a single core of a 2.4 GHz Intel Core 2 Duo.

Comments. Figure 1.18(a) shows the runtime for the dining philosopher models as a function of the number of philosophers. At the time of [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF], Cosmos is the fastest tool using Chernoff-Hoeffding bounds. For a hundred philosophers, Marcie is 1.4 times slower, Uppaal-smc is 1.5 times slower, Plasma is 1.9 times slower, and Prism using APMC 8 is 2.5 times slower. Among the tools using sequential procedures, the two versions of Prism have similar run times, Cosmos being up to 1.9 times faster.

Figure 1.18(b) shows the runtime comparison with different time bounds T for the tandem queue system experiment. Two kinds of behaviours for tools can be distinguished, depending on the applied statistic method. For the first one, that corresponds to the Chernoff-Hoerding method, the simulation time is increasing with the time bound T . About 295 000 trajectories are required to obtain the specified confidence interval. For the second one, corresponding to sequential confidence interval methods, the required number of sampling decreases when Accuracy comparison. Finally, to assess the accuracy of Cosmos, its output has been compared to the one produced by Prism via both its numerial engine and its statistical engine (with confidence level 99.99% and 0.01 interval width). The results indicate that Cosmos and the statistical engine of Prism are comparably accurate with the estimated intervals. Both always contains the value obtained using the numerical engine of Prism. This numerical value is also used to perform a coverage test of Cosmos: the ratio of simulation that return a confidence interval containing the real value is always close to the confidence level.

CHAPTER 2

APPROACHES FOR THE CONTROL OF AUTONOMOUS VEHICLES

Abstract. The control of autonomous vehicles triggers numerous issues: how to deal with the mechanical constraints of the vehicles? how to take into account the uncontrolled behaviour of other vehicles or pedestrians? how to specify several road situations such as crossing, or entrance ramps? what are the relevant safety indices? etc. Given the huge scope of these tasks, researchers have tackled subproblems, by identifying significant case studies. The goals related to these subproblems can be classified in two different ways:

• must the controller be automatically generated (synthesised) or manually designed and then verified? • are the involved methods based on geometry or logic? Given the diversity and the number of related publications, we only present a selective bibliography which provides relevant examples for each kind of approach. We emphasize that most of the methods do not support probabilities and, as such, cannot take uncertainties into account and quantify risk. Objective. The goal is to check whether the trajectories given by another module (the planner) are robust with respect to noise and vehicle deviations.

Assumptions. In this communication, one assumes that checking evasive maneuvers is done in two steps:

• every vehicle computes periodically its spatiotemporal trajectory in the current situation, and broadcasts it to other vehicles; • then, every controller checks whether this set of trajectories is collision-free. In this case, these trajectories are selected. In the other case, the process is iterated with additional information. The study is restricted to the analysis of this second phase.

While broadcasting full trajectories may be a time and space-consuming task, it provides a more accurate information than the one only given by the position and speed of other vehicles (see Figure 2.1, all figures from this section are extracted from the communication). Indeed, the latter case would require a probabilistic model of the behaviour of these vehicles. In Figure 2.1(b), the support of random trajectories is illustrated by the colored cones, and it results in a probability of collision. This is addressed in [START_REF] Eidehall | Statistical threat assessment for general road scenes using monte carlo sampling[END_REF] and [START_REF] Althoff | Model-based probabilistic collision detection in autonomous driving[END_REF], and induces a more costly and conservative controller design. Computing trajectories with disturbances. In order to check for robustness, the model is enlarged with perturbations, leading to a deviation from the reference trajectory, which is specified by the following equation:

ẋ ∈ Ax + Bu (2.1)
where:

• x is the state vector, containing the deviations to the reference trajectory and their derivatives; • u is the input vector, containing the curvature of the reference trajectory and the steering angle of the vehicle; • A and B are interval matrices that represent the uncertainty over the trajectory. The actual trajectory is described by:

ẋ(t) = Ax(t) + Bu(t) (2.2)
where A and B are matrices compatible with A and B, and u(t) ∈ U .

A temporal discretisation with time parameter r is used to solve this equation system. It consists in iteratively building the set of possible positions in the time interval [kr, (k + 1)r] based on the set of possible positions in [(k -1)r, kr]. Once this discretisation is done, the solution of the homogenous equation is computed then the inhomogenous part is added. Such a process provides a sequence of polytopes whose union is the occupancy set over the trajectory.

If the occupancy sets do not intersect over a time interval [0, t f ] then the manoeveurs of all vehicles are safe.

Case of wrong-way driver. The usefulness of this method is illustrated with the case study of a vehicle going the wrong way. This vehicle puts at risk two autonomous vehicles, on a three-lane road. To reduce the risk of accident, both vehicles plan a lane-changing manoeuver, as shown in Figure 2.2.

In the case of the collision detection over random rectangles, the proposed approach is faster by a factor 16.04. For verifying planned trajectories, the factor is 8.9.

Logic-based approaches 2.2.1 Control for Collision Avoidance [25]

Objective. The goal is to build a robust controller on a motorway segment.

Assumptions. Every vehicle of the segment is controlled. During the control phase, at most one vehicle may lose control and behave in a random way.

n lanes    a 2 a 0 a 1 horizontal positions

3).

A position is defined by its coordinates (x, y) with y ∈ 0, n -1 the lane and x ∈ 0, -1 the horizontal position.

Physical parameters (such as the maximal speed vmax, maximal acceleration amax, minimum entry speed vmin and minimum delay between two entering vehicles) are assumed to be related to the motorway section. Moreover, it is supposed that state and position changes are synchronous. We give below a specification for relevant information of a vehicle and the configuration of a section:

A vehicle present in the section is a tuple a = a.x, a.y, a.v, a.c where: (a.x, a.y) is the position of a; a.v ∈ 0, vmax is its current horizontal speed; and a.c ∈ B a boolean denoting whether the vehicle is controllable. A configuration of the section is a tuple s = d, A where A is a finite set of vehicles, and d is an integer array indexed by the lanes such that d[i] denotes the time elapsed since the last entrance on lane i when d[i] < dmin, and d[i] = dmin otherwise.

A possible configuration is shown in Figure 2.3. We consider a motorway segment of one kilometer, with two lanes and = 100 positions. There are at most 25 vehicles on each lane, due to the safety distance. The speed values are in 0, 4 , which corresponds to a maximum speed of 40m.s -1 . The possible acceleration values are -1, 0 and 1.

Principle. This situation can be seen as a turn-based game between two players, environment and controller. It is a zero-sum game as the environment tries to force a collision and the controller wants to avoid it. The next paragraphs describe the behaviour of the system in a formal way. Vehicle movement constraints. Let a = a.x, a.y, a.v, a.c be a vehicle, its next state is written as a.x , a.y , a.v , a.c , where:

Description of the

• a.x = a.x + a.v (and if a.x > then the vehicle is removed); • a.y ∈ [a.y -1, a.y + 1] ∩ [0, n -1] if a.v > 0 (a.y = a.y otherwise); • a.v ∈ [a.v -amax, a.v + amax] ∩ [0, vmax].
The controllability status is not modified by a vehicle move. The speed of the vehicle can change in a non-deterministic manner within the requirements of the motorway segment. If the vehicle has a positive speed, it can also change lane.

The system evolves from a state s = d, A to a state s = d , A in two steps: • the environment step: a vehicle can be marked as non-controlled (if there is no such one), some vehicles may be added to the segment, and the non-controlled vehicle is moved (if it exists); • the controller step: all controlled vehicles are moved.

Let s = d, A be a configuration. Then s 1 = d 1 , A 1 can be reached by a environment transition by the following operations:

1. If no vehicle in A is uncontrollable, then the environment may select a ∈ A and set a.c = ⊥; 2. For every lane i, the value d

1 [i] = min(d(i] + 1,
a 2 .y = a 2 .y and a 1 .y = a 2 .y) : if a 1 is behind a 2 ,
then either a 1 stays behind or the two vehicles stay in their respective lanes (only if they were in separate lanes); 3. If a 1 .x = a 2 .x then a 1 .y = a 1 .y and a 2 .y = a 2 .y : if the two vehicles start from the same horizontal position, they stay in their respective lanes. These additional conditions help to avoid risky behaviours, as shown in Figure 2.4 extracted from [START_REF] Bérard | Collision avoidance in Intelligent Transport Systems: towards an application of control theory[END_REF].

It is well known that for such games, there exist positional winning strategies. If the initial state is winning for the controller, its strategy gives us an implementation of the controller; in the other case, the environment strategy provides a failure scenario that may help to modify the settings of the problem in such a way that the controller is winning. Reachability games are solved in polynomial time by a saturation algorithm, that we now describe.

Let S be the set of all configurations and let s 0 be the initial configuration, in which the motorway section is empty. For s ∈ S we define:

• Succ e (s) the set of successors of s by a environment transition; • Succ c (s) the set of successors of s by a controller transition.

The set S fail = {s | Succ c (s) = ∅} denotes the set of configurations in which the controller can not avoid a collision and we define S * = S \ S fail .

A strategy is a mapping f :

S * → S such that f (s) ∈ Succ c (s). The correspond- ing set of states G f (s 0 ) = G c f (s 0 ) ∪ G e f (s 0 ) is defined inductively by: • s 0 ∈ G c f (s 0 ); • if s ∈ G c f (s 0 ) then ∀s ∈ Succ e (s), s ∈ G e f (s 0 ) • if s ∈ G e f (s 0 ) ∩ S * then f (s) ∈ G c f (s 0 ). It is a winning strategy if G e f (s 0 ) ∩ S fail = ∅.
Informally, the standard algorithm searching for a winning strategy is as follows.

Let Sbad c (resp. Sbad e ) be the set of configurations from which the controller will lose (resp. the environment will win) if it its turn to play.

1. These sets are initially defined by: Sbad e = S fail and Sbad c = ∅; 2. These two sets are enriched using the following rules: a Let s be a configuration such that ∀s ∈ Succ c (s), s ∈ Sbad e , then Sbad c = Sbad c ∪ {s}; b Let s be a configuration such that ∃s ∈ Succ e (s), s ∈ Sbad c , then Sbad e = Sbad e ∪ {s}. 3. The algorithm stops if either s 0 ∈ Sbad c (in which case there is no winning strategy for the controller), or when no rule can be applied anymore. This algorithm is polynomial with respect to the size of the system. However, such a system is generally very big, which requires to define adapted data structures for the implementation.

Example 19. During the first step, the controller state q 4 that can only lead to a crash is added to Sbad c . Then, the environment states q 2 and q 3 leading to q 4 ∈ Sbad c are added to Sbad e . Finally, the controller state q 1 which can only lead to Sbad e is added to Sbad c . The green states are safe states.

q 1 q 3 q 2 q 4 q crash

Conclusion.

Other works have implemented the algorithm described above. In particular, in [START_REF] Zhang | Conception Semi-automatique de Contrôleurs avec VeriJ[END_REF] the objective is to generate a controller for a system described by a Java program.

Verification of Smart Intersections [59]

Objective. The case study presented in the chapter 6 of this PhD thesis focuses on Cooperative Intersection Collision Avoidance Systems (CICAS) that would enhance intersections with controllers that communicate with "smart" vehicles.

More specifically, it focuses on the CICAS for Stop-Sign Assist implementation (CICAS-SSA), where a driver crossing the road needs to decide when it enters in the crossing, depending on the gaps between other vehicles and the crossing. The verification objective is to ensure that, if the subject vehicle 1 driver follows the system advice, there is no collision.

Assumptions. The other vehicles behave according to the expected behaviour on a motorway.

Principles. The input of the method is the strategy of the intersection controller, which is specified as follows: the intersection controller allows the subject vehicle to enter the crossing only if all the oncoming vehicles are far enough away (beyond a fixed distance) from the intersection. Note that it is an advice, and the vehicle may wait until the next one.

A hierarchical heterogeneous tree (see Figure 2.6 extracted from [START_REF] Akshay | Multi-Model Heterogeneous Verification of Cyber-Physical Systems[END_REF]) is built: the idea is to decompose the problem into subproblems. The basic subproblems (i.e. located at the leaves of the tree) are modelled with an appropriate formalism and solved by a dedicated method. The other subproblems combine subproblems already handled and are then solved by compositional methods. More specifically, the compositional methods are either conjunctive, disjunctive, or disjunctive with switching (see below, node 3j).

We only present some subproblems of the case study, one for each kind of composition.

Node 01. The root presents a case of conjunction. The first level nodes correspond to independent aspects of the system: • in Node 11, the goal is to estimate the driver's response time;

1 the vehicle at the stop-sign, served by CICAS-SSA 48CHAPTER 2. APPROACHES FOR THE CONTROL OF AUTONOMOUS VEHICLES • in Node 12, the goal is to do worst-case execution time analysis on the computational model -estimating the delay in which the response of the intersection controller is given; • in Node 14, the goal is to make sure that the error of the sensing subsystem is always bounded, within a given range of positions and velocities; • in Node 15, the goal is to find the upper bound on the communication delay between the sensing subsystem and the computation subsystem • in Node 13 (see below), the goal is to establish that the subject vehicle is not in the intersection at the same time as another vehicle, for given values on computation time, communication time, sensing errors and driver response time. The Node 01 is a general case of the Node 13, that is completed with the analysis from Node 11, 12, 14 and 15. Therefore, this conjunctive abstraction completes the verification task of Node 01 (if the verification tasks of all first level nodes have been handled).

Node 13. This node corresponds to a case of disjunction. The goal described in this node is to ensure that the subject vehicle is not in the intersection at the same time as another vehicle. The verification model M 13 (shown in Figure 2.7) is consists of two hybrid automata components: • the Major Road component models the dynamics of oncoming vehicles, where

x represents the position of the nearest2 oncoming vehicle from the intersection (called the principal other vehicle). The value of x is always negative, by the definition of the intersection (which is located at 0, and x is increasing). Note that when a vehicle reaches the intersection, another vehicle becomes this principal other vehicle and x decreases to its position. • the SV components models the subject vehicle, which has to commit to a direction: either going straight, or turning left or right. It also represents the controller, through the conditions to move from a conflict state to a clear state. The position variables are distinguished with respect to the decision of the vehicle. The collision-freedom specification S 13 can be defined by the temporal logic formula 2.3 where constants are chosen based on a typical highway intersection geometry. Using behaviour relations to ensure that there is no variable conflicts between the three submodels, the authors show that there are heterogeneous implications and that finally if all these submodels verify their respective specifications, M 13 satisfies S 13 , as depicted in the verification tree shown in Figure 2.6.

¬((x = 0∧0 < y s < 4.5)∨(x = 0∧0 < y r < 170)∨(x = 0∧0 < y l < 180) (2.
Node 3j. This node corresponds to a case of a disjonction with mode switching. The goal described by this node is to ensure that the subject vehicle and a specified car are never in the intersection at the same time (in the case of a single lane). Figure 2.8 shows an inter-model switching between: • M 41 , initially safe, where the incoming vehicle is far enough from the intersection, and the subject vehicle can cross the intersection if it chooses to; • M 42 , initially unsafe, where the incoming vehicle is too close and the subject vehicle must stop. Observe that the designer (or the safety engineer) must provide this inter-model switching. This reduces the dimension (number of vehicles involved) to one of the single-lane multi-vehicle Major Road models to one for single-lane single-principal other vehicle. The different subject vehicle actions for the two models can be analyzed independently for safety while the system is in those modes.

Conclusion.

This application chapter tackles a high-level specification (the subject vehicle must not collide with other vehicles) by splitting it into subproblems, and applying the various methods developed earlier in the thesis. This early separation of various aspects of the verification model is useful in keeping the analysis complexity at check. This proof of concept may be applied to other objectives, including the vehicle industry, which heavily uses model-based designs. However, this kind of methods is limited: it is not always possible to do such a decomposition, and the mode switching operator requires additional modeling.

Formal Approach to Vehicle Coordination [6]

Objective. The goal of this paper is to provide a formalism describing coordination protocols between autonomous vehicles such that it can be encoded as a constraint satisfaction problem and solved by standard tools, like Z3.

Assumptions. Each vehicle is equipped with an autonomous controller and can communicate with other vehicles. Figure 2.9: System Architecture Principle. The system consists of (see Figure 2.9, all figures are extracted from [START_REF] Asplund | A formal approach to autonomous vehicle coordination[END_REF]) a control module represented by an automaton, the coordination protocol, other cars, the network and of the environment. This work adopts the point of view of a single vehicle. The behaviour of the other vehicles is modelled by a set of constraints overapproximating their possible trajectories. This set of constraints, with the property, is then given to a solver (such as Z3).

Case Study: Simple Intersection. The case study consists of a four way intersection. The entity A is approaching the intersection and is assumed not to turn. There are four regions for this entity in relation to the intersection: "far away", "close", "in intersection" and "passed" (see Figure 2.10). Conclusion. The solver approach has been validated by the case study, partly solved using Z3. The entire model is composed of 825 lines of SMT-lib code (including comments), and the verification by Z3 took 14 seconds on a Dell optiplex 990 with a 3.4 GHz Intel Core i7 processor, using 109MB of memory and generating 965 000 equations. There are limitations on the expressivity of the formulas (euclidian distance has not been implemented by the authors). Moreover, some features of the approach are not sufficiently described in order to convince the reader that it may be applied to a more general context.

CHAPTER 3 AN OPERATIONAL SEMANTICS FOR SIMULINK

Abstract. We briefly mentioned Simulink ® in section 1.1 as a more complex model than hybrid automata. This tool from MathWorks implements a modelling language and simulation engine, largely used in the French automotive industry for the design of controllers for autonomous vehicles. However, the model behaviours are not defined formally. For this reason, several works proposed formal translations from (subsets of) Simulink blocks to other models like hybrid automata [START_REF] Tiwari | Formal Semantics and Analysis methods for Simulink Stateflow Models[END_REF][START_REF] Agrawal | Semantic Translation of Simulink/Stateflow Models to Hybrid Automata Using Graph Transformations[END_REF], or languages like Lustre [START_REF] Tripakis | Translating discrete-time Simulink to Lustre[END_REF]. Other works directly define exact semantics [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF] or operational semantics [START_REF] Bouissou | An operational semantics for simulink's simulation engine[END_REF] for Simulink. We follow the latter approach and propose semantics for a significant fragment of Simulink. We proceed in two steps: we first develop an exact version, and then enrich it with effective procedures with the aim to integrate it into the model-checker Cosmos.

Contents

Syntax

We first introduce a formal syntax for Simulink block diagrams, that we call here Simulink models or, shortly, SK-models. An SK-model consists of a set of blocks connected by wires supporting signals. More precisely, blocks act as operators that transform their input signals into output signals.

Types, signals and operators

Types are associated with signals and operators. We denote by Type the set of types used in SK-models.

Definition 32 (Basic types and constructors)

Basic types are a subset of Type containing subsets of the set R of real numbers:

• booleans B = {0, 1};

• (un)signed integers, either represented with varying number of bits (int8, uint8, int16, uint16, int32, uint32) or with usual sets N or Z; • floating numbers (double, single) or the entire set of real numbers R. A constructor is a function Type n → Type, where n ∈ N is the arity of the constructor.

Remark. This definition covers both exact (N, Z, R) and operational (other types) semantics.

Example 20. The constructor Tuple 3 multiplexes three signals into one, like in Tuple 3 (N, R, R).

Signals will be evaluated over a time domain, which will usually be a finite interval of R denoted by T = [t init , t end ] and called the simulation interval. Except for a finite number of discontinuities, the signals are smooth functions because they are solutions of differential equations over subintervals:

Definition 33 (Signals)

A signal of type Tp ∈ Type is a right-continuous piecewise smooth function s : T → Tp, admitting a left limit denoted by s(t -) for each t ∈ T.

We denote by type(s) the type of s and by Sig T (Tp) the set of signals of type Tp defined on T.

Discontinuities are handled like in discrete event systems by applying instantaneous operations to the left-limit s(t -) to obtain the value s(t).

Definition 34 (Operators)

An operator is a function op : Sig T (Tp 1 ) × . . . × Sig T (Tp m ) → Sig T (Tp) such that for each t ∈ T, op(s 1 , . . . , s m )(t) only depends on the restriction over [t init , t] of the signals

(s i ) 1≤i≤m .
This restriction is natural as an operator cannot predict the future values of its input signals.

Blocks

A block contains a tuple of operators generating output signals from input signals. Blocks are classified using three criteria:

(i) Whether they are continuously evaluated. If not, they are called discrete, either asynchronous or synchronous. In the latter case, a sampling delay must be provided; (ii) Whether there is a latency for the evaluation of inputs. A block without latency is called immediate and a non zero latency is either positive or infinitesimal, where the output value at time t only depends on the inputs over [t init , t[ (which is the case for integration); (iii) Whether the output value depends on threshold crossing by an input signal, called critical input and denoted by i c . In this case, the threshold values (v i ) i∈I must be specified, as a countable increasing sequence without accumulation points. Such blocks are called threshold blocks.

Definition 35 (Block type)

A For instance, the generic block type Add can be defined with (Add i ) i∈N where i is the arity.

Example 21. We illustrate the definition of block type with standard examples: The Switch block (Figure 3.1(b)) is a threshold block (over its input i c ) with exactly one discontinuity point. Its operator is defined by: op(i

1 s v 0 i (a): Integrator >0 i 1 i c i 2 (b): Switch 1 z v 0 i (c): Unit Delay
1 , i c , i 2 )(t) = if i c (t) v then i 1 (t) else i 2 (t)
where the comparison operator and the threshold value v are the additional parameters. It is a continuous and immediate block. On the figure, is > and v = 0. The Unit Delay block (Figure 3.1(c)) samples its input every δ, with a δ latency. It is a discrete block with initial value v 0 as additional parameter.

Block list

We now give a more extensive description of the set of blocks considered in this work. ) is an immediate block adding its two inputs:

op 1 (i 1 , i 2 )(t) = i 1 (t) + i 2 (t)
. It corresponds to the instance Add 2 of the generic block type Add introduced earlier.

This generic block type can be extended to (Add) sign where sign is a tuple of signs in {+, -}. For this new generic block type, the operator will be defined by op

1 ((i k ) 1≤k≤m )(t) = m k=1 sign k • i k (t). Block Product (Figure 3.3(b)
) is an immediate block that computes the product of its two inputs: op Block Delay is a synchronous discrete block, with three parameters: the sampling delay δ, the delay length n and the initial value v 0 . It samples and delays its input signal and its behaviour is represented by the operator op 1 (i 1 )(t) = i 1 (t init + t-t init δ -n δ). This block usually has a positive latency, unless n = 0 in which case it is immediate.

1 (i 1 , i 2 )(t) = i 1 (t) × i 2 (t). It can
Block Unit Delay is particular case of the Delay block, where the delay length is fixed to n = 1. Hence, the operator is defined by op

1 (i 1 )(t) = i 1 (t init + t-t init δ -1 δ).
op(i 1 )(t) = v 0 if t -t init < r i(t -r) otherwise.
Block Integrator is a continous block with infinitesimal latency. It has one parameter: the initial value v 0 . Its operator is op(i 1 )(t) = v 0 + t t init i 1 (τ )dτ .

e) Threshold blocks

>0 i 1 i c i 2 (
(i 1 , i c , i 2 )(t) = if i c (t) v then i 1 (t) else i 2 (t).
More generally, we define a generic block type Switch Tp 1 ,Tp 2 with Tp 1 the type of the input signals i 1 , i 2 and the output signal, and Tp 2 the type of the conditional input. If the output type Tp 1 is a boolean type, then the block is an asynchronous discrete block.

Block Relay (Figure 3.6(b)) is an immediate asynchronous discrete block with four parameters: v on , v off , s on , and s off . It represents a bang-bang controller which outputs v on when it is activated and v off otherwise. It has one input and one output and its current state can be viewed as an internal signal. Its unique operator is defined as follow:

op(i)(t) =          v on if i(t) > s on v off if i(t) < s off op(i)(t -) if t > t init v off otherwise.
Block Floor (Figure 3.6(c)) is an immediate asynchronous discrete block without parameter. It represents the exact floor function with the countable set Z of threshold points, with: An event operates as a block with the exception that its output is not a signal as it does not fulfill the right-continuity condition. The input and output signals of a Chart are considered as variables of the state machine, respectively denoted by X I and X O . We write X = X I X O for the set of variables and type(X) = x∈X type(x). A valuation is an element v of type(X), also written as a tuple v = (v(x)) x∈X .

op(i)(t) = n if n ≤ i(t) < n + 1. f) (Stateflow) Chart
In order to specify transitions, we need a syntax to define their guards:

Definition 38.
For the set X of variables and the set E of events, the set G(X, E) of guards is defined by the following grammar:

ag ::= F | x c | x x | τ g ::= ag | g ∧ g | g ∨ g where F ⊆ E is a subset of events, τ ∈ Q + represents a time constraint, x ∈ X, ∈ {<, ≤, =, ≥, >, =}.
Variable updates are generally specified by MATLAB code. Therefore, in order to avoid ad hoc notations, we define them in a semantical way. An update modifies the values of output variables according to the current values of all variables. Definition 39. An update is a mapping a : type(X) → type(X) such that for every x ∈ X I : a(v)(x) = v(x). The set of updates is denoted by Up(X).

We now have all the ingredients necessary to define a Chart:

Definition 40 ((Stateflow) Chart)

A Chart is a tuple S = (Q, q 0 , X, E, ∆, Pri, init) where:

• Q is a finite set of states and q 0 ∈ Q is the initial state; • X = X I X O is a finite set of variables composed of a set X I of input variables and a set X O of output variables; • E is a finite set of events; In order to evaluate guards and updates, we need the notion of environment. An environment is composed of a valuation for every signal, an activation condition and the time elapsed since the last state change: Definition 41. An environment of a Chart S is a tuple env = (v, b, d) where v is a valuation, b an activation and d a duration. A configuration of a Chart is a pair (q, env) where q is a state of S and env an environment.

• ∆ ⊆ Q × G(X, E) × Up(X) × Q is the transition function; • Pri : ∆ → N is an injective priority function; • init ∈ Up(X) is the initial action. A transition δ = (q 1 , g, a, q 2 ) ∈ ∆ is written q 1 g,a -→ q 2 .

Example 22.

Note that q and d are stored as internal signals of the Chart block.

Definition 42 (Evaluation of a guard). Given an environment env = (v, b, d) and a guard g, we define env |= g inductively by:

• env |= F if ∃e ∈ F : b e ; • env |= x c if v(x) c; • env |= x y if v(x) v(y); • env |= τ iff d ≥ τ ; • env |= g 1 ∧ g 2 iff env |= g 1 and env |= g 2 ; • env |= g 1 ∨ g 2 iff env |= g 1 or env |= g 2 .
Definition 43 (Evaluation of an update). Let a ∈ Up(X) and env = (v, b, d). We define:

a(env) = (v , b, 0) where ∀x ∈ X : v (x) = a(v)(x).
Example 23. In the Chart shown in Figure 3.9, the update function defined by the transition from q 1 to q 2 is a : type(x) × type(y) → type(x) × type(y) such that a((α, β)) = (α, α). Note that if a variable is not specified on the transition, its value is unchanged. Definition 44. A transition δ = q 1 g,a -→ q 2 ∈ ∆ is enabled in a configuration (q 1 , env) if env |= g. The configuration reached by firing δ is (q 2 , a(env)). The global transition is written (q 1 , env) δ -→ (q 2 , a(env)). We have (q, env)

δ = ⇒ (q , env ) if: • (q, env) δ - → (q , env );
• there is no δ ∈ ∆ such that Pri(δ ) < Pri(δ) enabled in (q, env). We define ⇒= ∪ δ∈∆ δ = ⇒ and ⇒ * the transitive reflexive closure of ⇒.

We now explain how, on activation, a Chart block updates its output and internal signals. We adopt here the maximal step semantics meaning that all possible transitions are fired to obtain the next configuration. We assume here that every loop between states contains either at least one time constraint or exclusive conditions on input signals.

Definition 45. We define (q, env) max = = ⇒ (q , env ) if and only if (q, env) ⇒ * (q , env ) and there is no enabled transition from (q , env ). The resulting configuration (q , env ) is denoted by SuccMax((q, env)).

Let (q , env ) = SuccMax((q, env)) with env = (v , b , d ). We write:

• q = SuccMax st (q, env); • d = SuccMax dur (q, env); • v (x) = SuccMax x (q, env) for all x ∈ X O .

In order to completely specify the semantics of a Chart, we introduce the sequence (t i ) i≥0 of activation times containing (in increasing order) the initial time and the times where at least one event is activated.

The time sequence can be determined from the specification of the input signals, as illustrated by the following example:

Example 24.
e The event signals e 1 and e 2 in Figure 3.10 are computed respectively from input signals i 1 and i 2 . The sequence of activation times, in this example, is (t init , t 1 , t 2 , t 3 ).

From the time sequence (t i ) i≥0 , the operator of the Chart provides the corresponding sequence of configurations (q i , (v i , b i , d i )) i≥0 . Given (q i , (v i , b i , d i )) at time t i , the configuration at time t i+1 is defined as follows:

• First, the configuration at t - i+1 is built without changing the inputs:

(q i , (v i , b i , d i + t i+1 -t i ))
• Then, the environment for t i+1 is obtained by copying output signals and updating the valuation of inputs signals alongside with their events:

(q i , (v i , b i+1 , d i + t i+1 -t i )) • This finally yields the configuration: (q i+1 , (v i+1 , b i+1 , d i+1 )) = SuccMax((q i , (v i , b i+1 , (d i + t i+1 -t i ))))
Note that we only have considered a subset of the capabilities of Stateflow Charts, which includes their main features. For instance, we do not consider nested Stateflow Charts.

g) Summary

We give a summary of blocks in alphabetical order: Simulink models have an underlying graph representation which has a semantic interest.

+ + 1 z -k Add Constant

Definition 47

The graph G M of the SK-model M = (B, L) is a labelled directed (multi-)graph defined by: • the set B of nodes • for every link ( B , o , i, B ) ∈ L an edge from B to B labelled by (o, i).

Example 26. The graph of the SK -model from Figure 3.11 is given in Figure 3.12: The application of block operators in order to obtain a trajectory requires temporal dependency constraints between blocks and more specifically between immediate blocks. The following definition allows us to identify incorrect dependencies:

B 1 B 2 B 3 B 4 B 5 (1, 1) (1, 1) (1, 1) (1, 1) (1, 1) (1, 2)

Definition 48

An immediate cycle of an SK-model M is a cycle in G M along immediate blocks.

An SK-model is correct if it does not contain any immediate cycle.

Example 27. We illustrate the notions above with the SK -model shown in Figure 3.13(a), which contains an immediate cycle with immediate block B 2 linked to itself. It would represent the equation x(t) = x(t) + 1 for all t ∈ T, which is impossible.

1 As shown in the next proposition, the correctness property is easily decided.

B 1 + + B 2 (

Proposition 5. SK-model correctness is decidable in linear time.

Proof : It is done with a topological sort on its graph restricted to immediate blocks.

The decision procedure could also use a depth-first search algorithm but the topological sort is used later (definition 54) in the notion of block ordering. Like in Simulink, from now on, we only consider correct models.

Exact semantics

Since an SK -model represents a deterministic hybrid system, its semantics should be a unique trajectory. It should be clear that such a semantics cannot be operational since the trajectory is evaluated over an interval in R. Moreover, even with discrete sampling points, the possible presence of an Integrator block forbids the exact computation of these values. Furthermore, due to the general form of implicit differential equations, there is no guarantee that a trajectory exists.

Principle. The trajectory of an SK -model over its simulation interval T = [t init , t end ], if it exists, can be obtained informally as follows. The interval [t init , t end ] is split into a finite sequence of contiguous subintervals. The intermediate interval bounds correspond to either threshold crossing or discrete block sampling. Within each interval, the trajectory is the solution of a differential equation system depending on the signal values obtained until then.

Differential equations of an SK -model

In this section, we first define the backward graph of a block in an SK -model. We also define the notions of mode (for threshold blocks) and history (for blocks with latency). We use these definitions to explain how to obtain the differential equations that represent the behaviour of the SK -model.

Example 28. Using Add, Product and Integrator blocks, it is easy to represent polynomial differential equations systems by SK -models. It requires exactly as many Integrator block as variables in the system. For instance the system { ẋ = x 2 + xy, ẏ = xyz, ż = y 2 } might be described as shown in Figure 3.14: The specification of differential equations also depends on the mode of threshold blocks. We denote by Th(M) the set of threshold blocks of M.

1 s 1 s 1 s x y z + + ẋ × × × × × × ×

Definition 50 (Mode)

Given a threshold block with threshold values {v i } i∈I , these discontinuities define a partition of R into contiguous intervals. A mode of M is a mapping associating an interval with each block in Th(M).

To complete the specification of the differential equation system, we still have to define which functions should be substituted for each non-immediate block. There are three types of non-immediate blocks B to be considered: • B is a block without inputs. In this case, the function associated with each output of B is substituted in the differential equation; • B is an integration block. In this case, the integration variable associated with the output of B is substituted in the differential equation; • B is a positive latency block. In this case, the substituted function should be given by an external mechanism, which requires the following definition.

Let Lat(M) be the set of positive latency blocks of M.

Definition 51 (History).

A history of M is a mapping associating a function h (B,o) with each output o of a block B of Lat(M).

We can now define the differential equation system obtained once a mode and a history are fixed.

Definition 52 (Differential equations of a model)

Let m be a mode, and h an history.

init = 0 B 1 -1 B 2 1 s init = 1 B 3 ẋ2 ẋ3 ẋ1

Example 32.

Recall the SK -model of Figure 3.11, for which the necessary backward graphs were given in Figure 3.14. For this model, the equations are ż = h (3) , the output of block B 3 , and either ẏ = 1 if mode

= {] -∞, 0[} or ẏ = z if mode = {[0, +∞[}.

Interval Partitioning

Several factors may impact the partition of the simulation interval T = [t init , t end ]: • Sampling blocks modify their output instantly at each sampling step, hence creating a discontinuity which requires splitting the simulation interval at their sampling times; • Positive latency blocks signals must be known during subinterval where differential equations are fixed. Hence the length of a subinterval cannot be greater than the minimal positive latency.

• Signals on critical input of threshold blocks should satisfy the mode specified by the differential equation to solve. A new mode corresponds to a new splitting point of the simulation interval. We introduce the static next step function next s in order to deal with the two first factors, as well as δ max , a bound on the length of the subintervals, which corresponds to the delay specified by an SK-model. Given t the lower bound of a subinterval, next s (t) gives the (maximal) upper bound of the subinterval that satisfies the first two conditions above.

Definition 53 (Static Next Step)

Let M be a SK-model, let ∆ be the set of its sampling steps and let R be the set of its positive latencies. We set δ lat = min(δ max , min(R)), and we define the static next step value next s (t) for time t < t end by:

δ samp (t) = min{pδ | δ ∈ ∆, p ∈ N, pδ > t} next s (t) = min(δ samp (t), t + δ lat , t end )
When the interval changes, we need to reevaluate the output signals of some of the sampling blocks. This requires an order for the evaluation of blocks, which is also needed to evaluate signals which are not output of an Integrator block after differential equations have been solved.

Definition 54 (Block Ordering)

The block ordering BO of a (correct) SK-model M is a total order which extends the following order on blocks:

1. blocks without input; 2. blocks with positive latency, by increasing latency; 3. blocks with infinitesimal latency; 4. immediate blocks in topological order.

Trajectories

We now give detailed explanations about trajectories, beginning with conditions under which a trajectory exists for an SK -model.

Definition 55 (Trajectory)

The vector w is a trajectory of the SK-model M over T = [t init , t end ] if there exists (i) a sequence (t i ) i∈ 0,N with t 0 = t init and t N = t end and for each i: t i+1 ≤ min τ >t i (τ | ∃j ≤ i, r ∈ R : τ = t j + r), (ii) a minimal sampling delay ε T > 0, which will be a lower bound on the length of sub-intervals, and (iii) a sequence (m i , h i ) i∈ 0,N -1 of modes and histories, such that: 0. The initial mode m 0 and history h 0 are those of the initial values of the model. Moreover, for 0 < i < N , h i on [t i , t i+1 [ agrees with w on [t 0 , t i [.

1. For all i < N , the differential equation system ˙ x(t) = F m i ,h i (t, x(t)) admits a solution on the interval [t i , t i+1 + ε T ] which coincides with w on [t i , t i+1 [ and is consistent with m i on ]t i , t i+1 [. The vector w agrees with the operations of the remaining blocks.

2. For all 0 ≤ i < N -1, w(t i+1 ) corresponds to the application to w(t - i+1 ) of the operators of active discrete blocks (i.e. such that t i+1 is a sampling time). The solution of ˙

x(t) = F m i ,h i+1 (t, x(t)) on [t i+1 , t i+1 + ε T ] is consistent with m i+1 .
In this definition, we use the sampling delay ε T to deal with mode changes. Another approach, described in [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF], uses non-standard analysis which implies handling infinitesimals.

Example 33. Returning to the example of Figure 3.11, the history is h [START_REF] Parosh | Timed Petri Nets and BQOs[END_REF][START_REF] Agrawal | Semantic Translation of Simulink/Stateflow Models to Hybrid Automata Using Graph Transformations[END_REF], and the mode m (1) 1 =]0, +∞[ must agree with the value of z at time 1 + ε T . Hence with t 2 = 2, the differential equations are ż(t) = t -1 and ẏ = z, which yields z(t) = (t -1) 2 /2 and y(t) = (t -1) 3 /6 + 1.

(3) 1 (t) = y(t -1) = t -1 on
In some cases, there is no trajectory:

Example 34.

The differential equation ẋ = 1 + x 2 with x(0) = 0, easily represented with a Simulink model, has the unique maximal solution tan(x) over ] -π 2 , π 2 [. Hence, there is no trajectory over any interval strictly containing ] -π 2 , π 2 [.

Unicity of a trajectory.

With suitable hypotheses on the operators, we prove that if a trajectory exists, it is unique. Unfortunately, as is often the case for hybrid systems, the existence of a trajectory is an undecidable problem. This semantic looks like it is non-deterministic as it seems to depend on how the interval is splitted, the value of ε T and the solutions of the succesive differential equations. However, with the hypothesis made on operators (which are piecewise smooth), we have:

Proposition 6

An SK-model has at most one trajectory.

The proof of this proposition is based on the following theorem ([26] theorem 1.8.2 (page 121)): Theorem 7. Let f : (R × R) → R be a locally Lipschitz function with respect to its second variable and let I be an interval of R. If there are two exact solutions ϕ 1 and ϕ 2 : I → R of the differential equation ẋ = f (t, x), and if they are equal at a point t 0 ∈ I, then ϕ 1 and ϕ 2 are equal into the whole interval I.

A trajectory is defined by several differential equations, for which we need the local-Lipschitz condition to hold:

Lemma 8

Let w be a trajectory with its sequences (t i ), (m i ) and (h i ), then for every i ∈ 0, N -1 : F m i ,h i is locally Lipschitz with respect to the second variable.

Proof : From definition 52 we remark that all the intermediate functions (including the history functions) are C ∞ since they are defined inductively by solutions of previous differential equations, except in the initialisation where the history functions are constants and for blocks without input.

Then, the derivative of the solutions are bounded and F m i ,h i are locally Lipschitz with respect to the second variable.

Proof of proposition 6 :

Let us consider two trajectories w and w of the SKmodel, with (t i ) i∈ 0,N and (t i ) i∈ 0,N the two sequences characterising the underlying time interval partitions and ε T , ε T the respective minimal sampling delays. We define (u i ) i∈ 0,M as the time sequence obtained by merging these sequences with u 0 = t init = t 0 = t 0 . Note that (u i ) still satisfies for each i:

u i+1 ≤ min τ >u i (τ | ∃j ≤ i, r ∈ R : τ = u j + r). We denote by (m i , h i ) i∈ 0,M -1
and (m i , h i ) i∈ 0,M -1 the modes and histories over the intervals [u i , u i+1 [ induced by the two original trajectories.

We prove by induction that, for all i ∈ 0, M -1 , m i = m i , h i = h i , and w(u i ) = w (u i ).

Initialisation. By definition, as they are defined using the initial values of the model, we have m 0 = m 0 , h 0 = h 0 and, w(u 0 ) = w (u 0 ).

Induction. Let k ∈ 1, M -1 . Using the induction hypothesis, we have m k-1 = m k-1 and h k-1 = h k-1 and w(u k-1 ) = w (u k-1 ). Thus, the differential equation ˙ x(t) = F m k-1 ,h k-1 (t, x)
is the same for both trajectories over [u k-1 , u k [, hence, using lemma 8 and theorem 7, those trajectories coincide over this interval. This implies that h k = h k and that the limits w(u - k ) = w (u - k ). The values w(u k ) and w (u k ) are defined by application of the operator of active discrete blocks, hence they are equal. Let ε = min(ε T , ε T ). Since the modes m k and m k are determined using the solution of ˙

x = F m k-1 ,h k (t, x) over [u k-1 , u k + ε],
they are equal. Therefore, we have:

h k = h k , m k = m k and w(u k ) = w(u k ).
We have shown that for every i ∈ 0, M -1 , the trajectories are equal over [u i , u i+1 [. The trajectories are thus equal over [t init , t end [. We also have w(u M ) = w (u M ) using the same argument as in the recursion part of the proof. Thus, the two trajectories w and w are equal on the whole interval [t init , t end ] and unicity holds.

Failure cases. There are several cases where a trajectory does not exist: • as seen above, if the model contains a differential equation without solution over its time interval; • if the solution violates the mode constraints related to some threshold block (see example 35); • or if there are an infinite number of discontinuities (see the model

s on = 1, s off = 0, v on = -1, v off = 1.
Seeing as B 1 always output 1 or -1 and B 4 is the absolute value of the output of B 1 multiplied by B 7 (which is always positive), computing the value of B 7 requires to study the behaviour of the mode shown Figure 3.19.

Lemma 9

Let M be the model shown Figure 3.18, and let

N ∈ N. Let t k = 1 -1 2 k for all k ∈ [0, N [. Then, a(t) = (k + 1) + 2 k+1 (t -t k ) b(t) = k + 1 c(t) = 2 k+1
Proof : We proceed by induction. Let

ε T = 1 2 N +2
Initially (at t 0 = t init = 0) we have a(t 0 ) = 1 and by direct implication b(t 0 ) = 1 and c(t 0 ) = 2 = 2 1 . Over [0, 1 2 [ a(t) = 1 + 2t, b(t) = 1 and c(t) = 1. There is no crossing of the thresholds of the block B 6 . However at t 1 = 1 2 , a(t 1 ) = 2 corresponds to some crossing of B 6 . The next step of the solution stays consistent over [ 1 2 , 1 2 + ε T ] with the mode of B 6 . Also, this means that a(t 1 ) = 2, b(t 1 ) = 2, and c(t 1 ) = 2 2 .

We assume the proposition true for

k ∈ [0, N -1[. In that case, a(t k+1 ) = k + 2, b(t k+1 ) = k + 2 and c(t k+1 ) = 2 k+2 . For any t ∈ [t k+1 , t k+2 [ we have a(t) = k + 2 + 2 k+2 (t -t k+1 ), which implies b(t) = k + 2 and c(t) = 2 k+2 .
There is no crossing over this interval of the threshold of the block B 6 . At t k+2 , a(t k+2 ) = k +3 and we have directly b(t k+2 ) = k +3 and c(t k+2 ) = 2 k+3 . Similarly, the next step of the solution stays consistent over [t k+2 , t k+2 + ε T ]. Hence, the proposition is true for k + 1.

Finally, the input of B 5 is the absolute value of the input of block B 3 . However, the input of B 3 changes signs whenever B 3 reaches 1 (with increasing values) or 0 (with decreasing values). Hence,

Proposition 10

For all k ∈ N : over [t k , t k+1 [:

x(t) = (k mod 2) + (1 -2(k mod 2))2 k+1 (t -t k )
Note that there is an infinite number of discontinuities from the threshold crossings of blocks B We used a discrete controller operating in a continuous environment to define the exact semantics of SK -models. This controller is limited by its reaction capacities, that are implicitly specified by ε T . Hence, when going from an interval of the simulation to another, the controller needs to keep a mode to decide whether it applies to the next differential equation.

It is well-known that questions on hybrid systems trajectories are usually undecidable (e.g. threshold crossing). Here, even the problem of existence of a trajectory over a finite interval is undecidable.

Proposition 11

Trajectory existence problem for SK-models is undecidable.

Proof : The proof is a reduction of the halting problem for two-counter machines (which is undecidable [START_REF] Minsky | Computation: Finite and Infinite Machines[END_REF]) to the trajectory existence problem for SK -models.

Recall that a two-counter machine consists of a finite sequence of labelled instructions, which handle two counters c and d and end with a special instruction with label Halt. The other instructions have one of the two forms below with x ∈ {c, d} representing one of the two counters:

1. : x := x + 1; goto 2. : if x = 0 then goto else x := x -1 goto Since a (long) sequence instructions of the first type can set the counters at any initial values, we assume that the counters have initial value zero. The behaviour of such a machine is described by a (possibly infinite) sequence of configurations:

< 0 , 0, 0 >< 1 , c 1 , d 1 > • • • < i , c i , d i > • • •
where c i and d i are the respective values of counters c and d, and i the label of the current instruction after the i th instruction.

The reduction is done by encoding a two-counter machine as a Stateflow Chart, which will be activated by the signal x of the model shown in Figure 3.18. Then, by adding an output to this Chart which is connected to the enabler model, we ensure that the final model has a trajectory if, and only if, the two-counter machine halts.

B 1 × × × B 2 1 s B 3 init = 0 |i| B 4 1 s B 5 init = 1 B 6 2 i B 7

Figure 3.21: The resulting model

The model built by this operation is depicted in Figure 3.21. The output of block B 6 is also connected to the Chart to ensure that only one step of the two-counter machine is computed at each activation.

The Stateflow Chart is constructed as follows:

• for each instruction label , the Chart has a corresponding state state( ); • for each counter c and d, a variable named similarly is chosen; • two additional variables are added: s indicating whether the Chart reached the Halt instruction (which is initialised to 1 and is an output of the Chart) and n a variable that will be used to prevent the Chart from doing more than one step at each activation; • two inputs, e being the event input (which is the output of block B 3 above) and m the output of the B 6 block which will be used to check the current iteration;

• for each instruction : x := x + 1 goto an arc state( ) ----------→ -------------→ state( ); • for each arc leading into the state state(Halt), we append s ← 0 to its list of instructions.

m>n,{x←x+1;n←m} -------------→ state( ); • for each instruction : if x = 0 then goto else x := x -1 goto , two arcs state( ) m>n∧x=0,n←m
Now, it has been seen (proposition 10) that the block B 3 generates a peak at each 1 -1 2 ( 2n+1) , n ∈ N. This sequence has an infinity of values in [0, 1] meaning that the Stateflow Chart will be activated as many times as possible. When the Chart reaches the state(Halt) state, the output B 2 is forced to 0 and this stops the output of B 3 from changing further. The model Figure 3.21 has a trajectory if, and only if, the Chart halts, which is if and only if the two-counter machine halts. Thus, the trajectory existence problem for SK -models is undecidable.

t i σ 0 = next s (t i ) σ 1 = next int (t i , σ 0 ) t i+1 = next z (t i , σ 1 )
Observe that this proof uses discrete asynchronous blocks.

Operational semantics

The exact semantics relies on solving differential equations and determining threshold crossings. It is well-known that such operations are not effective. We now define an approximate semantics that will be the base of a simulation engine for SKmodels. Moreover, it should be possible to give guarantees on the approximation towards exact semantics. General principle. The approximate semantics is based on the iterative construction of a subinterval partition of [t init , t end ] = N -1 i=0 [t i , t i+1 ]. In order to control errors induced by approximations, we again define ε T as the minimal stepsize and add an accuracy parameter ε V satisfying: |x| < ε V ⇒ x 0. The approximate semantics will replace complete trajectories by an array mapping, for each output o of each block B, its values W B,o [i] at each step t i , 0 ≤ i ≤ N . Iteration step. Let t 0 , . . . , t i be known, as well as the values W B,o [j] for each output o of each block B, and 0 ≤ j ≤ i. The first step is to determine t i+1 . This value is obtained by applying the next step function next s : σ 0 = next s (t i ), and then possibly lowering this value:

1. First, by applying a variable-step integration method, such as Runge-Kutta-Felhberg 1 [START_REF] Fehlberg | Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems[END_REF], between t i and σ 0 over all Integration blocks, giving σ 1 = next int (t i , σ 0 ); 2. Second, by detecting the mode change of threshold blocks, by linear approximations, producing t i+1 = next z (t i , σ 1 ). In both cases, it is necessary to compute the values corresponding to backwards graphs of concerned blocks at intermediate points generated by these methods.

Variable-step Integration in SK -models. The chosen method, ODE45, is based on the simultaneous computation of two approximations by fourth and fifth-order Runge-Kutta methods, as described in [START_REF] Fehlberg | Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems[END_REF]. If the difference between the two approximations is above ε V , then the stepsize is reduced with a minimum value of ε T and the procedure repeated.

However, in order to use Runge-Kutta methods, we need to explicitly know differential equations, while here they are implicitly specified by the SK -model. This requires to adapt these methods, evaluating backward graph blocks to the intermediate times of the Runge-Kutta method.

Example 36. In order to give a general idea, we show an application of the classical fourth-order Runge-Kutta method on a simple Simulink model. This can be adapted to any other Runge-Kutta method.

Classical fourth-order Runge-Kutta method. We consider the differential equation ẏ = f (t, y) with initial value y 0 = y(t 0 ). The classical RK4 method of stepsize h iteratively computes the values y n = y(t 0 + nh) values:

y n+1 = y n + h 6 (k 1 + 2k 2 + 2k 3 + k 4 )
where

         k 1 = f (t n , y n ) k 2 = f (t n + h 2 , y n + h 2 k 1 ) k 3 = f (t n + h 2 , y n + h 2 k 2 ) k 4 = f (t n + h, y n + hk 3 )
At each step, the output signal of an Integration block is estimated using initial value and an elapsed time multiplied by the derivative. To simplify notations, we will symbolically denote k j = f (τ j , z j ) the equation corresponding to step j. Note that this example is simplified as the block outputs never directly depends on the current time (except for integrals) and there is no delay functions nor threshold block. In a Simulink model, the function represented by f results from the computation of the backward graph of the input block of an Integrator block. We consider the Simulink model of Figure 3.23 as an example. Initial values of B 1 and B 3 are respectively 1 and 0. We will simulate the model over interval [0, 10] with a maximal stepsize 0.5. The exact semantics of output signals are cos(t) for B 1 , and cos(t) for B 2 andsin(t) for B 3 . We compute the backward graphs for blocks B 3 and B 2 (inputs of B 1 and B 3 ), that are respectively shown in figures 3.24(a) and 3.24(b). We denote k j,i the coefficient k j associated to block B i . Similar notations will be used for f and other blocks. Using those backwards graphs, we can unwrap the definition of equations k j,i = f (τ j , z j,i ). We then have k j,1 = z j,3 for block B 1 and k j,3 = z j,1 for block B 3 . We can then apply the RK4 method to the Simulink model. We then describe the calculations for t 0 = 0 and h = 1 2 . For the first step, we have

1 s B 1 init = 1 -1 B 2 1 s B 3 init = 0
1 s B 3 (a): Backward graph of B 3 , in- put of B 1 1 s B 1 -1 B 2 (b): Backward graph of B 2 , in- put of B 3
k 1,1 = f (t 0 , z 1,1 ) = W 3,0 = 0. Likewise, k 1,3 = f (t 0 , z 1,3 ) = -W 1,0 = -1.
The other steps are described in the following table: Final values are approximations of cos(0.5) andsin(0.5) exactly to 10 -3 .

1 k 1,1 = f 1 (0, y 0,3 ) = y 0,3 = 0 k 1,3 = f 3 (0, y 0,1 ) = -y 0,1 = -1 2 k 2,1 = f 1 ( 1 4 , 0 + 1 4 k 1,3 ) = 1 4 k 1,3 = -1 4 k 2,3 = f 3 ( 1 4 , 1 + 1 4 k 1,1 ) = -1 -1 4 k 1,1 = -1 3 k 3,1 = f 1 ( 1 4 , 0 + 1 4 k 2,3 ) = 0 + 1 4 k 2,3 = -1 4 k 3,3 = f 3 ( 1 4 , 1 + 1 4 k 2,1 ) = -1 -1 4 k 2,1 = -15 16 4 k 4,1 = f 1 ( 1 2 , 0 + 1 2 k 3,3 ) = 1 2 k 3,3 = -15 32 k 4,3 = f 3 ( 1 2 , 1 + 1 2 k 3,1 ) = -1 -1 2 k 3,1 = -7 8 - y 1,1 = y 0,1 + 1 12 (k 1,1 + 2k 2,1 + 2k 3,1 + k 4,1 ) y 1,3 = y 0,3 + 1 12 (k 1,3 + 2k 2,3 + 2k 3,3 + k 4,3 ) = 1 + 1 12 (0 + -1 2 + -1 2 + -15 32 ) = 0 + 1 12 (-1 -2 -15 8 - 7

Remark. It is possible that a threshold block changes mode between two intermediate steps of the Runge-Kutta methods. As we tackle the threshold crossing immediately after, we apply Runge-Kutta method with fixed modes (the mode of time t i ).

Threshold crossing. Let σ 1 = next int (t i , σ 0 ). We determine if there is a mode change between t i and σ 1 . In this case, we compute the minimal instant of mode change between all threshold blocks, by doing a linear interpolation over all critical inputs. This linear interpolation is evaluated with a maximal precision given by ε V .

Computing new values. The following step

t i+1 = next z (t i , next int (t i , next s (t i ))
given by the procedures described above, we compute output of all blocks at this given time in the block order. For Integration blocks, as the ODE45 method has already given a guarantee over precision, computing their new value with RK4 is good enough. For other blocks, we apply operators with respect to the current mode.

Then, similarly to the exact semantics, we compute the new values at t i+1 + ε T to deduce the mode at t i+1 . If a mode change occurs in ]t i+1 , t i+1 + ε T [ then the result is fail. Approximation conjecture. This operational semantics is proposed with the aim to provide an accurate approximation of the exact semantics. This is formally expressed as follows:

Definition 56 (Approximation)

Let w be the trajectory of an SK-model M. Then, the set of values W produced by the operational semantics accurately approximates w if, for any ε, there exist values of ε T and ε V such that: This condition precludes classical problematic situations like those illustrated in Figure 3.26, where signals s 1 , s 2 and s 3 have an unsafe crossing of value v for t c = 1, 3 and 5 respectively.

∀t i , ∀k, ∀o, |W k,o [i] -w k,o (t i )| < ε.
We conclude this chapter with the conjecture: Overview. Cosmos is built around three main parts: • a Code Generator, which transforms the model files and the description of the HASL formula into C++ code that will be included in the Simulator; • the Simulator, which is built using a library (that we describe shortly after) and the C++ code generated by the Code Generator; • the Server, which runs in parallel several instances of the Simulator then aggregates their results and computes the statistical evaluation, based on the chosen statistical method.

Technical Description of Cosmos

General notations. In this chapter, code listings are used to present the data structures and the class templates. Functions will be presented using readable pseudo-code.

Code optimization.

In the actual code, several optimization techniques are implemented, some of them using C++ templates: in particular the curiously recurring template pattern1 which implements static polymorphism. We will not describe in detail the use of templates in the code, but will present the functional aspects of the classes.

For instance, we present the two classes SPN and SKModel as fulfilling the DEDS interface which communicates with the Simulator, as shown in Figure 4

Code Generator

In this section, we first present the Marking class used to store the current state of a model. Then, we describe the DEDS class that will be instantiated on code generation. We then show the version of this class used in the library for Petri nets, reducing the number of elements generated by the Code Generator. Usually, the generalized stochastic Petri nets are given to Cosmos in the .grml format, that can be produced using the GreatSPN Editor2 .

DEDS class.

As shown in algorithms 2 and 3, the Event Queue management is handled by a class corresponding to the DEDS interface (that we shorten to DEDS class). To this aim, two functions will be generated by the Code Generator: initialEventQueue which initialises the event queue at simulation start, and update which updates the event queue after the firing of a transition. This generic structure is currently instantiated for Petri nets and Simulink models.

The DEDS class handles the current state, stored into DEDState which contains, in the case of Petri nets, an instance of the Marking class. This state is updated by the fire function according to the specification of the transition. The initial value of the state is filled through the constructor of the class. Petri Net models. In the case of a Petri Net, the update and initialEventQueue functions are implemented as a part of the library, through SPNBase.cpp:

• the initialEventQueue iterates over all transitions, checking whether each transition is enabled and, in that case, adding it to the event queue with a random delay corresponding to the specified distribution (using generateEvent). • the update function, called after the firing of transition T , shown in Algorithm 1, updates the Event Queue with an incremental method:

first, it checks whether transition T is still enabled, and removes or updates it accordingly; -then, it checks each possibly enabled transition, and add new possible events.

It uses the PossiblyEnabled vector, generated by the Code Generator, that maps each transition t to the list of transitions t for which there exists a place p with t → p → t ; -then, it checks each possibly disabled transition and remove events that have been disabled. It uses the PossiblyDisabled vector that maps each transition t to the list of transitions t for which there exists a place p with p → t and p → t and t = t . -finally, each transition for which the time distribution depends on the value of the current marking (stored into the MarkingDependant vector) is checked and updated accordingly.

Example 37. In the net shown in Figure 4.4, we have:

t
• PossiblyEnabled(t 1 ) = {t 3 , t 4 }; • PossiblyEnabled(t 2 ) = {t 3 }; • PossiblyEnabled(t 3 ) = PossiblyEnabled(t 4 ) = ∅; • PossiblyDisabled(t 1 ) = PossiblyDisabled(t 2 ) = ∅; • PossiblyDisabled(t 3 ) = {t 4 };
• and PossiblyDisabled(t 4 ) = {t 3 }.

In the case if t 2 is fired, PossiblyEnabled(t 2 ) is checked and tells that we need to check if t 3 is enabled. Then, we check the number of tokens of places P 1 and P 2 . If there is at least one token in each of these places, we know that t 3 is enabled. We then check if this transition is already in the Event Queue. If not, this transition is added to the Event Queue (the time and weight are drawn through generateEvent).

There is a more clever way to perform this update function with a better incremental update where a value is associated with each transition giving the number of conditions that are not yet satisfied. However, simulations using this approach are only faster for ordinary Petri nets. The advantage of the chosen algorithm here is its flexibility and, in particular, its adaptability to high-level Petri nets.

Generating a model class. The generalized stochastic Petri nets are given to Cosmos in the .grml format, and then transformed into a C++ code that will have the behavior of the net. For example, the Petri net shown in Figure 4.4 would have the fire function shown in Figure 4.5. For each transition t, this function applies the token changes for the impacted places. Other model functions are generated in a similar way: IsEnabled would return the result of the test associated to the places of the transition (in this case, for t 3 , it would be Marking.P->_PL_P1 > 0 && Marking.P->_PL_P2 > 0). • CurrentLocation and CurrentTime are respectively an integer representing the current node of the automaton, and a double representing the current time; • FormulaVal is a vector of doubles containing the current value of each path expression; • FormulaValQual is a vector of booleans indicating if a path expression, used for binomial distribution evaluation, is currently satisfied.

It also contains a set of protected functions, which will be used by the class LHA_orig (described in the next paragraph) that will serve as an interface for the Simulator. Among them:

• DoElapsedTimeUpdate updates the variables of the LHA according to the given elapsed time and the current marking of the model; • GetEdgeEnablingTime gets the earliest time at which an edge of the LHA will be enabled (when its guard becomes true), by solving a linear equation system; • DoEdgeUpdates updates the LHA variables after the firing of an edge of the LHA, according to the current marking of the model.

Complete LHA class.

The final class for the Simulator is another class, called LHA_orig (see Figure 4.7), that inherits from the LHA class described in the previous paragraph. It contains mainly public functions:

• fireAutonomous firing an autonomous transition of the LHA;

• updateLHA which updates all the values of the LHA and HASL path expressions according to an elapsed time whose value is indicated by DeltaT and to the current state of the model; • synchroniseWith checks if there is an enabled synchronisation edge (with GetEnabled_S_Edge) and fires it, and returns whether there was such an edge. Some of them rely on protected functions. Among them:

• fireLHA fires a transition of the LHA, by using the DoEdgeUpdates function to update the LHA variables according to the current marking of the model, then changes the current state. • GetEnabled_S_Edge checks all LHA transitions starting from the current LHA state and which are labelled by the model transition; for each of these transitions, it then checks whether the guard of the destination and the edge are fulfilled, in which cases it returns the corresponding LHA transition.

Simulator

We describe here the Simulator class template which handles calls for model simulation and automaton synchronisation. The algorithms will then be presented, followed by the Event structure used in the EventQueue class.

Simulator class. The essential functions of the header file of the simulator are shown in Figure 4.8. As can be seen, the Simulator is a template that extends the SimulatorBase template. This SimulatorBase template has three generic types:

• the Simulator S itself;

• the structure EQT for an Event Queue;

• and the simulated model DEDS.

It implements several functions, the essential ones are: • SimulateOneStep which fires autonomous edges of the LHA until the next activation time of the model, and then performs the corresponding step. It returns whether the model has a next simulation step. • SimulateSinglePath which simulates the model until:

1. a final state of the LHA is reached; 2. a model transition cannot be synchronised with the LHA; 3. the Event Queue is empty. • interactiveSimulation which implements a command-line interface where the user can control the simulation by choosing to execute exactly one step, or execute the model for a chosen delay, or firing one of the enabled transitions. Note that unlike the model which needs a template as its structure might be quite generic, the LHA_orig class is fully instantiated.

Simulator algorithms. The SimulateOneStep function is described by Algorithm 2. This function first defines τ next the time of the next event (+∞ if there is no such event) and stores the next event of the model into a variable E. Then, while there is an autonomous edge of the LHA for which condition is fulfilled before τ next , it fires it and updates the LHA variables accordingly. In the case where a final state is reached during this loop, the simulation is accepted and the function returns false as there are no more possible simulation steps. Once this loop is finished, it checks whether there is a possible event in the model. If there is none, the simulation is rejected and the function returns false. Finally, it updates the LHA variables for their new value at τ next and fires the event E. If the transition of event E cannot be synchronised with the LHA, it rejects the simulation and returns false. Otherwise, the function performs the synchronisation and checks whether the LHA has reached a final state. In that case, it accepts the simulation and returns false. In the other case, it finally updates the Event Queue (using the model's update function) and returns true.

The SimulateSinglePath function is built directly using SimulateOneStep. It is described by Algorithm 3. It first calls the initialEventQueue of the model, initialising the Event Queue, then runs the SimulateOneStep function in a loop as long as the simulation can progress.

The procedure for interactive simulation is also based on SimulateOneStep, and is described by Algorithm 4. It is similar to the single path function, except that it maintains a variable minInteractive which contains the minimum time at which control is given back to the user. When the user has control, he may:

1. either fire a specific transition, which gives full priority to the event associated with this transition (executed at current time, with highest priority and lowest weight); 2. or execute a single step of the model; 3. or execute steps of the models until a given time (stored into minInteractive). Other features are implemented into the interactive simulation, such as waiting for a specific transition to be fired, or drawing the current Petri net state (through Graphviz's DOT language 3 ).

A verbose level can also be given when running Cosmos, which is useful for debugging. Depending on the verbose level, additional information will be provided about the current state of the model, the Event Queue, etc. This information will be printed by SimulateOneStep.

Event class. It is used to represent events of a model and consists of five variables: • an integer transition: the identifier of the transition; • a double time: the time at which the event will happen; • a double priority: the priority of the event (for example, the priority of the transition); • a double weight: the weight of the event; A fifth variable binding corresponds to the binding (of the variables of the transition) of the event. It is used in high-level Petri nets and its type depends on the structure of the net (this will be explained in section 4.2).

This class also implements several functions: • update functions setTime, setPriority, and setWeight that update the corresponding variables with the parameter of each function; • a function isPriorer which takes another event e 2 as parameter and returns true if e 2 should be fired before this event (see section 1.1.3)

EventQueue class. The event queue is defined as a class (see Figure 4.9) that maintains:

• a vector of events evtTbl which maps each possible (transition, binding) pair to an event; • a vector evtHeap, which is a heap of enabled (transition, bindings) pairs in the order of priority of their events; • and a vector evtHeapIndex which maps each possible (transition, binding) pair to the position of the corresponding event in the heap (or -1 if the event is not in the queue). Note that the heap (evtHeap) refers to each event through its (transition, binding) pair. The event is then found in the evtTbl vector. In order for the model to update its Event Queue, several standard functions are available: • insert which inserts a new Event in the Event Queue; • replace which replaces the Event of the same (transition,binding) pair by a new element, to update time or priority; • pause and restart which allows to pause the time of an event, and reschedule the event with the remaining time, in case of an age memory policy. • isScheduled which checks whether there is a scheduled Event of the (transition, binding) couple; • remove which removes of the queue the Event corresponding to the given (transition,binding) pair. These functions use the classical heap functions siftUp, siftDown and swapEvt described below. Event Queue Algorithms. The heap is built as a binary tree, implemented by the functions described in this section. The vectors are initialised through the class initialisation function: for that, it checks the number of possible bindings for each transition. The following functions are defined to ease further definitions: getLeftChildIndex(k) = 2k + 1, getRightChildIndex(k) = 2k + 2 and getParentIndex(k) = (k -1)/2 . We first define the basic functions: • swapEvt exchanges the positions of the events of nodes i and j of the heap, and modify the vector evtHeapIndex accordingly. • siftUp first checks if the node is the root node, in which case it does nothing.

Otherwise, it gets the parent Event and checks whether the current event (in the i th node) has a higher priority (with respect to the isPriorer function). In that case, it swaps both events (with swapEvt) and performs a recursive call on the new index (basically, (i -1)/2 ). • siftDown first finds the highest priority child (after checking whether there are zero, one, or two childs), then checks if the current event has a higher priority; In that case, it swaps both events and performs a recursive call on the new index.

The following functions are then defined: • isScheduled gets the transition and binding identifiers, checks the value in evtHeapIndex, and returns false if the corresponding value is -1 and true otherwise; • InPosition simply reads the current value in evtHeapIndex;

• getEvent takes a (transition, binding) pair and returns the corresponding event using evtTbl; • insert inserts an event in the queue by putting it at the highest index, and then using siftUp until stabilisation. Using the (transition, binding) pair as an identifier prohibits the multiple (or infinite) server policy, since in this case an additional indexing should be necessary to distinguish between several instances of the same pair.

TimeGen class. This class handles the generation of the delays for newly-enabled events, according to the distribution associated with their transition. It implements a function called GenerateTime which uses the library Boost and the Mersenne Twister random number generator.

Server

We now describe the functionalities of the server which include: Server. The server algorithm is shown in algorithm 5. It requires the number of parallel simulators to launch, the list of HASL formulas, and the statistical parameters. When KillClient is called, it also erases all the current batches made by running simulators.

Statistical methods.

The evaluateHASL function of the server will compute the new interval confidence for each HASL formula, using the Boost.Math library Possible Improvements. Cosmos is currently designed to work on a single computer, and can benefit of multiple cores of a processor by launching several clients in parallel. The current implementation can however be transformed for multi-machine simulation by: • transforming UNIX files (used for communication between processes) to sockets;

• broadcasting the Simulator binary to every machine. This raise the problem of having a similar infrastructure so that the generated binary can actually run on each of the machines. The current implementation is not flexible enough to support binaries generated on different infrastructures: there is no way to certify that the format of batch results will be the same on each machine.

Improving the Binding Mechanism for highlevel nets

The original implementation for high-level Petri nets has two limitations: it cannot handle infinite color classes (like the integers), and it was inefficient for large color classes. Here, we describe our first implementation contribution: an alternative data structure for markings on which we developed a more efficient search for enabled bindings, in particular in the case where color classes are very large or even infinite. This alternative data structure is toggled via a command-line option by the user.

Standard Implementation

Processes = {a, b, c, d, e} Ressources = {1, 2, 3} a c
Processes Waiting In this file, each basic (finite) color domain dom is represented with: • an enumeration of its colors (dom_color_Classe) with an element

Color_dom_dom_IC_i for each possible i a possible value of the color), Color_dom_Total used to end iterations and Color_dom_All a special value used to represent all colors of a domain in an edge; • a structure dom_Token with two variables, one for its class c0 and the other for its multiplicity mult. This structure will handle the iteration of the enumeration of color classes (with functions iter, next(int i), and islast()); • a structure dom_Domain with one variable mult which is a vector of the multiplicity of each possible color class. This structure will be used to handle the marking of each place. Finally, the class abstractMarkingImpl is defined, creating a C++ variable for each place (which contains the marking of the place) and each variable of the net (used to enumerate all the possible bindings). These C++ variables are typed with the Domain structure described above.

Example 38.

For the model shown in Figure 4.11, the initial marking is represented in the following way: • Processes: 1 0 1 0 0 ; • Ressources: 2 1 0 ; • Acquired: 0 . . . 0 (with 5 × 5 × 3 = 75 null entries).

The model of the net is similar to previously explained: • the InitialEventQueue iterates over all transitions, then over bindings of the variables of the transitions, then checks whether the transition is enabled using this marking. If the transition is enabled, it is added to the event queue; • the update function, shown in Algorithm 6 updates the event queue with the same incremental method as before, but this time iterating over bindings when needed.

Binding Iteration. As shown in Algorithm 6 describing the update function for high-level Petri nets, the function nextPossiblyEnabledBinding is called to iterate over bindings of a transition. This is done by iterating over the color domains of each variable. The Hospital model. The Petri net from [5] shown in Figure 4.13 describes the flow of patients in a hospital emergency department. This net has been modified to be able to scale the number of patients, and the diseases (and the patient priority) are now chosen through stochastic transitions (instead of a static relation between patient and disease priority). There is a single color class, representing the set of patients. All patients, represented by their tokens, start in the Healthy place. They can FallIll, and are assessed at their arrival. Depending on the priority of their illness, they go through different examinations and treatments.

Benchmarks for the Hospital model. We have run 5 000 trajectories of 200 time units with both implementations of the high-level models, for 20 or 100 patients. Results are shown in Table 4.1. The standard implementation runs this model 25% faster than the alternative implementation, due to having the same magnitude of tokens than of colors. In this case, the alternative implementation does not have significant impact on the memory used. Benchmarks for the SIA model. We have run 60 trajectories of the model until either the controlled vehicle exits the section of the motorway or a collision happens. We do these simulations with different values of cells (number of possible PosX values) and vehicles. Results are shown in Table 4.14. Note that the standard implementation is far slower than the alternative implementation, and not even running for high number of cells (because the size of the binding vector is higher than the available memory). Finally, the alternative implementation runs faster with more vehicles in the case of 500 cells, due to a higher collision rate. Moreover, this alternative implementation uses negligible memory space compared to the first color implementation in this example.

The SIA model. The Petri net shown in

Conclusion.

In the standard implementation, the memory space used to store a marking is exponential with respect to the number of class occurrences in the domain definition. For a transition, in the worst case, the binding enumeration time is exponential with respect to the size of the expression labelling the edge. Due to subtle optimisations of the rules, in some cases (as in the Hospital model), the binding can be performed in constant time.

In the alternative implementation, the memory space used to store a marking is linear with respect to the number of tokens. For a transition, in the worst case, the binding is polynomial in the number of tokens in each place.

Simulink Integration

The Model class has been introduced in section 4.1.1 as a generalization of the structure used for Petri nets. In this section, we describe how Cosmos transforms Simulink models into this structure, with the approximate semantics described in section 3.3.

The SKTime class. In order to avoid undesirable approximations when dealing with time operations, we chose to represent time by the SKTime class. This class consists of two integers: time and precision. A (t, p) SKTime would represent a precise time of t.10 -p . Usual arithmetic operations and relations have been implemented in this class. Furthermore, cast must be done explicitly using the getDouble function. This is particularly useful when dealing with latencies. Code generation. The simulink models, given in the default .slx format, are a set of XML files describing the models and the architecture. The transformation of these models into the C++ code implementing their behaviour through the functions described in this section is done via OCaml. The information from the XML tree is parsed using the standard Xml-light library and converted into a graph which is represented by a pair:

• a list of blocks (the edges), which are 4-uplets (a string blocktype, an integer blockid, a string name and a list of (string,string) pairs values; • a list of simulinkLink values which are 4-uplets of integers (fromblock, fromport, toblock, toport. The default block values, provided separately in the Simulink model, are joined to this graph representation as a list of (string blocktype,(string,string) values) pairs. They are then used to complete the values of each block, not overwriting any of the provided values of the block. The latency time of each block is then computed from the given block values, as the Simulink value can be in a relative format.

At this step, the blocks are usually only ordered by their identifier. It is now necessary to reorder them, using the Block Ordering described in section 3.2.2. It is done the following way, using a hash table c mapping each block ID to a number, and a set s of blocks that can be added to the final sorted list (decreasing corresponding c values and adding it to the list):

• first, fill this hash table by getting the number of incoming edges of each block;

• initialise s by adding all blocks with no incoming edges; • add blocks not using any of their inputs to the final list of blocks;

• order blocks by increasing latency, then add blocks with non-infinitesimal latency to the final list; • then process the infinitesimal latency blocks (which is the Integrator blocks, plus the DiscreteIntegrator blocks using the Forward-Euler method); • finally, perform the topologic sort by dealing with each element of s. The list of non-processed links between blocks is kept during this whole procedure, and an error is returned if this list is not empty at the end.

Finally, the code is generated from this final list of blocks using each block's definition.

Multi-model Simulation

The requirements for automotive transports has led to deal with two different formalisms: one for the probabilistic aspects used to model vehicular environments (a high-level Petri net), and the other for the control aspects (a Simulink model). This multi-model phenomenon tends to increases, which led us to to propose a generic approach for multi-formalism modelisation.

In this section, we describe how the simulation with two models is handled, including how communication is done between the two models. It is still a work in progress, with the multi-model simulation currently only possible between Petri nets and Simulink models.

Generic multi-model loop

Two models into one. In order to reuse most of the code handling the simulation in the case of two interacting models, we embed them into one single model. This, however, requires to dispatch the event handling between models. The MultiModel class has two variables m1 and m2, each corresponding to one of the two models. An additional shift parameter is added to the DEDS specification which should be given when creating an instance of the corresponding class. The DEDS functions are adapted as follows:

• initialEventQueue calls the function of both models;

• fire calls the appropriate fire function, by checking whether the transition ID tr is greater or equal to m2.shift, in which case it calls the m2 fire function with the appropriate ID (tr -m2.shift). In the other case, it calls the m1 fire function with the appropriate ID (tr); • update calls the appropriate update function with the same parameters, except for the transition ID. This way, the Event Queue can be shared for both models. The only modifications to model functions consist in dealing correctly with shift in calls of functions of the Event Queue. This MultiModel is defined as a library, for any kinds of models.

Interface for communication. In the implementation of the MultiModel class, we further add two functions to the MultiModel specification. These functions are specific to the pair of considered models. They are generated when parsing both models. These are: • synchronize_fire which, after a transition is fired into a model, updates the state of the other model considering as input the new state of the first model; • and synchronize_update which fires the appropriate synchronisation transition before calling the update function of the corresponding model. Both functions are called in the respective fire and update functions, after the operations described in the previous paragraph are performed.

Future works.

As seen in the previous paragraphs, while the simultaneous simulation of multiple models can be done in a generic way, the communication between these models is done in a ad hoc manner. One possible way to generalize the construction of communication between models would be introducing a formal language to describe the relation between each model. This language would be parsed and the appropriate code generated. Another improvement in the case of multiple models would be to use communication paradigms of distributed models, such as broadcast. One may also need to have informations shared between different models.

Simulink and Petri Net Simulation

As seen in Chapter 1, it would be interesting to combine stochastic high-level Petri nets and hybrid models since they present complementary features (randomness versus continuous time and space). Hence, they are good candidates for instanciating the generic multi-model loop presented in the previous section.

Interface transitions. In this paragraph, we detail the communication interface between Simulink models and stochastic high-level nets.

There are two kinds of interface transitions: SK -in transitions, directed from the net to the SK -model and SK -out transitions in the other direction.

The input arcs of a SK -in transition (see Figure 4.17(a)) are test arcs (explained with the firing) connected to places of the net. Any output arc is connected to the input of a Simulink block. An SK -in is enabled when the content of at least one of its input place is modified. The firing of such a transition proceeds as follows: a function is associated with each output arc, taking as parameters the contents of the input places. When the firing takes place, the function is evaluated. This function can be specified by a multiset of tokens as illustrated, or by a C code associated with the arcs.

The input arcs of an SK -out transition (see Figure 4.17 Parsing these transitions leads to the generation of the the synchronize_fire and synchronize_update functions. Target Simulation loop. We now instanciate the multi-model simulation loop for Simulink models and stochastic high-levels nets. All enabled transitions are stored into an event queue implemented as a binary heap, with their time of occurrence, their priority and weight. The next Simulink step is added as a possible event. At each simulation step, the earliest event is chosen. Among simultaneous events, the (decreasing) priority order is the following: 1. SK -out firings, 2. ordinary transition firings, 3. SK -in firings, 4. SK -event. In case of equal priorities, the choice is randomized according to the weights. Once an event is selected: • If it is an ordinary transition firing, the marking is updated, the associated C code is executed; transitions that are newly enabled trigger new events while events corresponding to disabled transitions are removed. • If it is an SK -in firing, the Simulink signals are updated and the time of the Simulink event is set to the current time.

• if it is the SK -event, all output signals are updated and the time of the Simulink event is updated as presented in Section 3.3. Finally, the SK -out transitions are added to the event queue with the current time. • if it is a SK -out firing, the contents of output places are updated.

To simulate a discrete event system, at each step, one only has to compute what the next event will be and increase the simulation time to the time of this event. This is how ordinary Petri net transitions are fired in Cosmos. This leads to an efficient simulation of such system as the time to compute a simulation depends on the number of events and not on the simulated time. Unfortunately, this property is lost when simulating hybrid systems: the SK -event is triggered at least at a fixed frequency (δ max ).

Implementation of the simulation loop. In practice, this simulation loop is done using the synchronize_fire and synchronize_update functions we have described in the previous paragraph. These functions are generated when processing the Petri net. The generated functions work as follows:

• whenever the Simulink transition is called (which is the only case a SK -out transition exists in the queue), the corresponding synchronize_fire updates the corresponding markings in the net. Then, the corresponding synchronize_update fires the SK -out transition described in the net before updating the Event Queue using the net update function; • the SK -in transitions appear in the net PossiblyEnabled vectors, which means that whenever a transition changes the entry place of a SK -in transition, then the SK -in transition is added to the event queue. The synchronize_fire corresponding to this SK -in transition then use the setInput function of the Simulink model, updating the corresponding Inport block. The synchronize_update then calls updateInput which advances the SimulinkTransition to the current time.

Use Case: Double Thermostat

Among the multiple systems that can be modeled using this enhanced version of Cosmos, we choose a well known toy example (which is presented in our Petri Nets 18 article [START_REF] Barbot | Integrating simulink models into the model checker cosmos[END_REF]): a device with two heaters prone to faults, and using bang-bang controllers to keep the temperature in a room between 20 Figure 4.20 shows a simulation of the system. In the first period there is only a small failure of heater 2, and we can observe the bang-bang behaviours of the system. In the second period both heaters fail at the same time while the outside temperature is low, thus the temperature quickly drops to 13 • C before the first heater is repaired.

T ext - + c ext T h 1 F 1 - + c 1 × × × Heater 1 T h 2 F 2 + - c 2 × × × Heater 2 + + + 1 s Ṫ T
We are interested in several performance indices. The first type concerns the reliability of the model measured by two indices: the minimal temperature observed along a trajectory (I 1 ) and the time spent in a state where the temperature is below 20 • C (I 2 ). The second type concerns the average behaviour: the average temperature (I 3 ), the average number of switched-on heaters (I 4 ), which is correlated with the energy consumption of the system, and the average time during which the repairman is idle (I 5 ). These indices are specified in HASL with an LHA (Fig. 4.21) which accepts the trajectories after S time time units. The LHA contains a hybrid variable tc with derivative 1 when the temperature is below 20 • C and 0 otherwise. The HASL expressions start with a probability operator: here AVG is used for all indices to specify the average value over all trajectories; then a path operator (Min, Last, Mean) which is defined along each path. Path operators take as parameters algebraic expressions over the Petri net places and the LHA variables. For example, I 1 specifies the minimal temperature along a trajectory and then the average value over all trajectories. The model as it is described above is referred to as M 0 . In order to study the overhead of integral computations over the stochastic simulation, we build two additional alternative models. The first one M 1 is a model where the integration block in the Simulink diagram has been replaced by a discrete time integrator. In the model M 2 , the simulink part is omitted, keeping only events that are transition firings. Each model was run for 500 000 simulations of 2 000 seconds, with ε V = 0.01 and δ max = 1. The sine wave frequency was 0.01 and oscillating between 5 • C and 25 • C, and the time step of the discrete-time integrator was δ max (1 second). We used 4.3. The left table reports the computed confidence interval for the different indices, the right one reports simulation and building times. Tool analysis. The build time is always less than the simulation time and becomes negligible when models include a Simulink part. The critical factors for simulation time are: (i) the speed of step firing, about 10 -7 sec. for net firing compared to 10 -6 sec. for Simulink steps, and (ii) the number of steps per trajectory, about 40 for M 2 vs. 2000 for M 1 . As expected, the use of a discrete-time Integrator yields a faster simulation, albeit still far longer than the net alone, while it affects the accuracy of the index values and more precisely triggers a larger variation of temperature over time.

Indices

T h 1 = 55 • C, T h 2 = 65 • C, c 1 = 0.
Property analysis. We focus on the most pertinent model M 0 , with two antagonist goals: minimizing the installation cost (depending on the parameters of heaters and repairman), and maximizing the comfort of the user (depending on the temperature evolution). With the current parameters, each heater is active about 1/4 of the time and the repairman is idle 92% of the time. The average temperature is about 22 • C, reaching the objective, while the minimal temperature is slightly above 18 • C.

Motorway Segment

Here, we consider a one-way motorway segment, with one controlled vehicle and other vehicles with random behavior. We want to evaluate the given controller with various performance indices, such as: • the collision rate: probability of a collision before the vehicle exits the section; • the distance to collision: expected travelled distance before collision, if any.

Model of the Motorway

The situation is depicted in Figure 5.1. The controlled vehicle, shown in red, and the environment vehicles, shown in white rectangles, are handled separately. Their position is characterized by the cell coordinates (i, j) representing position i on lane j. The lateral and longitudinal speeds are also tracked for all vehicles, while the accelerations on both axes are only kept for the controlled vehicle. Indeed, the behaviour of the other vehicles is modeled by a random choice of speed every time unit, while the acceleration of the controlled vehicle is produced by the controller and the new speed is then computed in a standard way. Vehicle Interaction. This part of the net contains the transitions removing the vehicles from the simulation in two situations:

• Transition exit (respectively end) is fired when the vehicle longitudinal position is higher than N pos -v max , and removes the environment (respectively controlled) vehicle from the simulation. • Transition collision is fired when some environment vehicle has the same (x, y) position as the controlled vehicle. Transition enter corresponds to the asynchronous part of the simulation, adding a vehicle in the first v max positions of a lane.

Specification of the Performance Indices

Observe that a trajectory ends when either end or collision is fired. The firing of end cannot be synchronised with any of the HASL automata (see below), and leads to an unsuccessful trajectory.

The HASL formula depicted in Figure 5.7 (a) is used to evaluate the probability of collision, where T denotes the set of all transitions of the Petri net in Figure 5.2. The synchronization implies that state s col is reached exactly when transition collision is fired in the simulation and the expression PROB means that the value returned is the probability of reaching s col .

The HASL formula depicted in Figure 5.7 (b) is used to evaluate the expected distance covered by the controlled vehicle before collision (if it occurs). With respect to the previous automaton, we add a new transition reaching state s end that will be synchronized with the transition end of the Petri net, and a variable x that will take the current value of the horizontal position when the transition is fired. The expression E(LAST(x)) (expected value of the last value of x) yields the desired index. Table 5.1: Evaluation of controllers. For each measure, controller and other vehicle behaviours, we provide a 99% confidence interval of the measure, the time in seconds required to perform the simulation and the number of simulated trajectories.

Evaluation

The first objective of this benchmark is to compare the performances of different controllers in various situations. As explained above, we measure the probability of a collision (I 1 ) and the expected covered distance before collision if any (I 2 ). The second objective is to quantify the influence of these parameters on the simulation time.

Scenarios.

To evaluate our approach we have modelled simple controllers. We first study a basic controller (see Figure 5.8) which only aims at monitoring the vehicle ahead (l.12-20) and aligns its velocity to avoid collision with it. When there is no vehicle ahead or it is too far away, this basic controller speeds up until it reaches a maximal velocity v max (l.22). Then we study a more advanced controller (see Figure 5.9) which is able to pass the vehicle ahead. In order to change its lane the controller monitors vehicles in the two adjacent lanes (l. [START_REF] Berthomieu | Modeling and Verification of Time Dependent Systems using Time Petri Nets[END_REF][START_REF] Bouissou | An operational semantics for simulink's simulation engine[END_REF][START_REF] Bérard | Collision avoidance in Intelligent Transport Systems: towards an application of control theory[END_REF]. The controller initiates a take over only when the speed of the vehicle in front of him is slower than the speed v max and when there is no vehicle in the two adjacent lanes (l.28-39).

We have evaluated the controllers in the environment described in section 5.1.1 with two kinds of behaviours for other vehicles. In one set of experiments, all speed changes are smooth: at each step of simulation vehicles either stay at the same speed, moderately accelerate (the speed is increased by 1) or moderately decelerate (the speed is decreased by 1). This is done through move_other function as shown in Figure 5.10. Those behaviours produce a highly predictive environment making the controller job easier. In contrast the second set of experiment involves vehicles with unpredicables speed changes: at each step of simulation, every vehicle randomly selects a new speed in a predefined interval. The move_other function of this scenario is shown in Figure 5.11. For these computations we simulate a motorway with 2 lanes and 1000 cells. The controlled vehicle is surrounded by 20 other vehicles for which speed varies between 3 and 5 cells per step. More precisely, all vehicles (including the controlled one) are generated randomly, with a uniform distribution on both lanes and cells (up to the 200 th ). If a vehicle is generated on a cell where there is another vehicle, we generate new values for both lane and cell. The I 1 performance index is estimated with a fixed confidence interval width of 0.01, and the I 2 performance with a relative confidence interval width of 5%. All simulations were performed on a single machine with 12 Intel Xeon cores using 8 parallel threads.

Controller evaluation.

Comparing the probability of collision between the two controllers, we observe that the lane changing controller is far better than the basic one: 2.3% probability of collision vs 83.8% in the smooth environment and 2.2% vs 84.7% in the more aggressive one. We notice that the behaviours of other vehicles do not have a major impact on the probability of collisions. Similarly, the average travelled distance before collision is highly different between the two controllers. In Table 5.1 we also observe that most collisions occurred at the beginning of the simulation. They occur around cell 104 for the lane changing controller and around cell 275 for the basic one. This shows that the simulation is biased by its initial state. More reliable results could be obtained by simulating a sliding window of the motorway around the controlled vehicle. Instead of simulating the controlled vehicle on the entire motorway we could make the controlled vehicle static in the middle of a short motorway and simulate all vehicles relatively to the controlled one.

Furthermore, in real life, the other vehicles do not behave randomly but follow an implicit control law thus reducing again this probability. In summary, this benchmark could be improved by adding these elements in order to obtain more accurate results (see section 5.2). Scalability evaluation. We finally observe in Table 5.1 that when the probability of collisions becomes small, the number of simulations required to obtain relevant estimates becomes large. This can be seen for instance when comparing the required time to estimate the average distance before collision for the two controllers (71s vs 4147s and 74s vs 3766s). Thus the simulation could be intractable when the probability of collisions becomes close to zero. This is a well-known difficulty in presence of rare events for which several theoretical techniques and practical tools have been developed in the literature [START_REF] Rubino | Rare Event Simulation Using Monte Carlo Methods[END_REF][START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF] and can be used to improve this simulation.

Possible Improvements. We have identified other limits of our simulation which should be overcome to produce more precise results. Compared to real vehicles which feature turn signals and braking lights, there is no communication between vehicles in our simulation. This prevents the controller from synchronising its action with those of other vehicles. The behaviours of vehicles on the motorway should be enhanced to include some communication.

Currently, a synchronous simulation step is done in linear time with respect to the number of vehicles. A modeling integrating such synchronisation between vehicles would induce a quadratic blowup of this synchronous step and would require the introduction of the sliding window to overcome this complexity.

Entrance Ramp

We now consider an entrance ramp into a motorway segment, again with one controlled vehicle. A given controller will be evaluated with various performance indices, such as: • the vehicle collision rate, i.e. the probability of a collision with another vehicle before the controlled vehicle exits the section;

• the barrier collision rate, i.e. the probability of the vehicle not entering the motorway segment at the end of the entrance ramp (and colliding with the barrier); • the severity of collision, related to the velocities during a collision.

Model

Our model is based on an official French document [START_REF]'Équipement, des Transports, du Logement, du Tourisme et de la Mer. Conception des accès sur Voies Rapides Urbaines de Type A (VRU A)[END_REF]; the situation is depicted in Figure 5.12, using a longitudinal cell size of 7 meters. The model is similar to the case described in the previous section, except that the lanes are cut into 5 lateral sections. The lateral granularity is finer since we need to take into account the progressive entrance of the vehicle. The choices of 40 pre-entrance cells and 40 entrance cells are based on a 100km/h motorway segment according to [START_REF]'Équipement, des Transports, du Logement, du Tourisme et de la Mer. Conception des accès sur Voies Rapides Urbaines de Type A (VRU A)[END_REF]. This speed is called the expected speed of the segment, which is below the maximal speed. We add 20 cells to the section to check whether the vehicle is correctly inserted in the section. We also suppose that the controlled vehicle has all traffic information and tracks vehicles with unique identifiers. This choice is justified with respect to the capacities of sensors and the length of the section. With this size of cell and time step, the minimum positive acceleration is already above the standard acceleration. Synchronisation of simulation steps. This part is identical to the one in the motorway segment model. Transitions sim_step1, sim_step2 and sim_step3 act as schedulers for the simulation parts: sim_step1 has a Dirac distribution with a parameter corresponding to the time between two steps, while sim_step2 and sim_step3 are immediate, thus enforcing the synchronous behaviour of the simulation. The priorities of the transitions enforce the order of the simulation. There must be a token in step1 for random_other to be fired, and this token will not leave step1 until all tokens of otherVehicles have been processed. Then, this token is consumed by sim_step2, producing a token in step2 which allows update_other to be fired. Similarly, the step2 token cannot leave the place until all tokens of inProcess have been processed. Then, this token is consumed by sim_step3 which produces a token in step3, allowing update_self then Controller to be fired in this order. Controlled vehicle decision. This single transition represents the call to the controller. We now consider that the vehicle can either see all the vehicles of the section via, for instance, a radar (as the section is fairly small) or get this information via communication with the motorway. The controller can obtain through the orderedVehicle function a vector of vehicles that are ordered by increasing lateral position, increasing longitudinal position, then incresing lateral and longitudinal speed. The controller function being the target of the evaluation, it is described in Section 5.2.3.

Vehicle Interaction. This part of the net contains the transitions removing the vehicles from the simulation in two situations:

• Transition exit (respectively end) is fired when the vehicle longitudinal position is higher than N pos -v max and removes the environment (respectively controlled) vehicle from the simulation; • Transition out is fired when the position of the controlled vehicle is outside the segment (in the case it does not leave the entrance ramp on time); • Transition collision is fired when any the environment vehicle has the same longitudinal position and either the same lateral position, or a lateral position that is one higher or lower.

Specification of the Performance Indices

A trajectory ends when either end, out or collision is fired. The three HASL formulas shown in Figure 5.15 are based on the single automaton, which sets the appropriate variable to true, according to the outcome of the simulation. CollForce is set to 0; • if the controlled vehicle goes outside the motorway segment, CollForce is set to 1 + v 2 x whwere v x is the longitudinal speed of the vehicle; • if the controlled vehicle collides with another vehicle w, CollForce is set to 1 + |w x -v x | 2 where v x and w x are their respective longitudinal speeds. This choice of values should give a decent approximation: the collision is not too dangerous at a low speed, but the higher the relative speed, the (quadratically) greater the damages are.

Evaluation

The first objective of this benchmark is to compare the performances of two controllers in various situations. As explained above, we measure the probability of a collision (I 1 ), the probability of leaving the entrance ramp unproperly (I 2 ), the probability the entering motorway properly (I 3 ) and the severity of a collision (I 4 ). We also measure the probability of a timeout (T O) where the controlled vehicle fails to enter the motorway into the 500 simulation steps. The second objective is to quantify the influence of some parameters on the simulation time, in order to check the scalability of our approach. These parameters are: the number of vehicles on the motorway, the presence of vehicles on the entrance ramp, and whether a new environment vehicle is inserted after an exit.

Controller.

Let us now detail the two controllers. The first one tracks two vehicles, V f ront and V behind , the vehicle directly behind V f ront on the target lane. Then, it compares their positions and speeds, before either choosing a new front vehicle or starting the entrance process. This entrance process also requires that the controlled vehicles accelerates or decelerates up to the front vehicle speed. It is initiated if the front vehicle is ahead of the controlled vehicle and the speed is reached. These tasks are described more thoroughly later on.

The second controller is built over the first one, trying to avoid collisions with other vehicles in the Entrance Ramp by accelerating if the vehicle behind is too fast and close, and decelerating if the vehicle ahead is too slow and close. Moreover, when the barrier is too close, it changes its strategy and avoids colliding with a barrier by slowing down to zero.

Environment vehicles. The behaviour of the environment vehicles is more refined here than in the first case study. We assign a mode to each of these vehicles. Let us describe these modes: 0. In this (default) mode, the vehicle adjusts its speed to the target vehicle (the vehicle ahead). If this target vehicle is far away, then it adjusts its speed to the expected speed of the segment. It only does so by unit increments or decrements. This models the expected behaviour on the motorway. 1. In this mode, the vehicle does not check anything and continues at its current speed. This models, for example, a micro-sleep of its driver. 2. In this mode, the vehicle decelerates if its speed is above two cells per time unit. This could model the driver being interrupted by a phone call. 3. In this mode, the vehicle changes lane to the right lane of the motorway segment. This is used in the two following cases: either the vehicle is on the left lane and is preparing for the next exit, or the vehicle is also in the entrance ramp and is preparing to enter the motorway.

The choice of the current mode is randomized (see Figure 5.16) in the following way:

• The starting mode is mode 0 (speed adjustment); • If the vehicle is currently in the entrance ramp at position x, then it switches to mode 3 (lane change) with probability (x-c 1 +1) c 2 -c 1 +1 (meaning that the further it is to the end of the entrance ramp, the higher is the chance of going into the motorway segment). We set c 1 to 40 (the first position where the vehicle may Figure 5.17: The two target vehicles 1. synchronising with the vehicules inside the segment; 2. entering into the motorway segment; 3. driving in the motorway segment. In order to perform these tasks, the controller maintains two variables: targetspeed the current target longitudinal speed, and targetway the current target lateral position. These tasks rely on finding two vehicles between which the controlled vehicle will enter the motorway segment (see Figure 5.17); the two target vehicles identifiers will be stored in targetfront and targetbehind. The procedure is described in Algorithm 7. The controller enumerates the vehicles with respect to the decreasing longitudinal position, until the vehicle A is behind the controlled vehicle and there is enough space between A and B the vehicle in front of A. The variables targetbehind and targetfront are respectively set to A and B. The different tasks are described in Algorithm 8. During the acceleration phase, the controller synchronises its longitudinal speed to the front target longitudinal speed. It only moves to the next task after reaching the longitudinal position p e . Simulation parameters. The two previously described controllers are evaluated through four different scenarios, presented in order of increasing difficulty for the controller:

• the first one has 10 vehicles in the motorway segment; • the second one has additional vehicles in the entrance ramp;

• in the third one, the enivronment vehicles are added back to the beginning of the motorway section when they have reached the end; • in the fourth one, there are 20 vehicles in the motorway segment that are still added back, and the vehicles on entrance ramp are kept.

In Table 5.2, these scenarios are represented using symbols: the number indicates the number of vehicles in the motorway segment, the symbol indicates that exiting vehicles are replaced and the symbol indicates that some environment vehicles are present in the entrance ramp.

Let us recall the performance indices: • I 1 the probability of collision with another vehicle; • I 2 the probability of collision with a barrier of the entrance ramp or the motorway segment (when the vehicle did not enter properly); • I 3 the probability of a successful entrance in the motorway; • I 4 the average severity of the collision as described in section 5.2.2; • and T O the probability of a timeout (in the case where the vehicle stays at null velocity). We have chosen to perform these evaluations with a 10% relative width, a 99% confidence interval, and a limit of 2 000 000 simulations for a single scenario, on a single machine with 12 Intel Xeon cores, using 8 parallel threads.

Controller evaluation. The results are shown in Table 5.2. On the easiest scenario, the first controller already fails to enter almost 8% of the time and often collides with the barrier. The examination of Algorithm 8 reveals that no check is performed in the synchronisation state on either the vehicle behind, or the closeness of the end of the entrance ramp. It only checks whether the front target is indeed in front of the vehicle.

The chance of collision increases as expected with respect to the scenario difficulty. More interestingly, the main factor is the number of the vehicles either in the main lane or the entrance ramp. For instance, the number of collisions more than doubles between the easiest and most difficult scenarios.

However, adding the possibility for vehicles to loop in the section does not increase the chances of a collision, except slightly increasing a barrier collision chance, because less place is available for the vehicle to enter. It is only when a lot more vehicles are added to the simulation that the chance of success drops to around 83% (17% of failure).

For the second controller, designed to avoid any barrier collision and reduce the collisions in the entrance ramp, the results are surprisingly worse than those of the first controller. Indeed, the controlled vehicle completely avoids barriers but this cost of a greatly increased collision probability. In addition, with a significant probability, the vehicle fails to enter the motorway segment within the 500 time units. Whatever the scenario, the probability of entrance success is now between 77% and 82%.

These results show that to avoid collisions with vehicles in the entrance ramp, or the barrier, requires more attention for the design of the controllers. Besides, the overall collision rates are far higher than the realistic ones. It remains to understand whether the scenarios are too pessimistic or the controllers are too basic. Scalability evaluation. We have chosen a relative width since we expected that some performance indices could be very small. In fact, this occurs in two cases, corresponding to the two hardest scenarios for the second controller. In the first case ( 10) we stop without reaching the 10% relative width because of the 2 000 000 simulation limit. This case illustrates the rare event phenomenon.

In the other case, the simulation stops when reaching the 10% relative width on the three other performance indices, because Cosmos entirely ignores, for its termination condition, a performance index whose simulated mean is zero. It is very likely that rare events are associated with this performance index: the lower bound is zero, since no simulation has timed out; however, this simulation cannot guarantee that it never happens. This would require a further study.

With only 10 environment vehicles, 90 simulations are performed per second, which is more than in the Motorway Segment case study. Progressively adding difficulties, it drops to 79 per second with up to two vehicles in the entrance ramp, then to 54 per second with the vehicles being allowed to loop. The worst case scenario has only 19 simulation per second, performed with 6 cores, which is only a third of the simulation speed of the first case study.

Possible Improvements. The chances of collision are really high for the case study: over such a segment of 120 cells (which is around 840 meters) the collision rate is approximately 10 -9 . None of the simulations have produced a better entrance rate than 93%. We have already dicsussed that the first controller was rough, not even trying to notice whether other vehicles are interacting in the entrance ramp. The discussion could be extended to the second controller, for which the additional checks to avoid barrier collisions lead to a greater vehicle collision rate than the sum of collision rates of the first controller. However, the design of controllers is outside the range of the thesis, and we must limit to some improvement ideas: a new controller could improve its estimation of the next positions of the target vehicles, and allow for acceleration during a manoeuver.

We have only roughly modeled the behaviour of the other vehicle: in real situation, a vehicle on the motorway might decide to slightly go on the left to ease the entrance of a vehicle in the entrance ramp. This, with the controller slightly enforcing its entrance before the end of the entrance ramp, could help reducing the number of collisions. approach, in the context of the autonomous vehicle. The first one is a motorway with heavy traffic, and the second one is an entrance ramp to a motorway.

Perspectives

Our perspectives are organised from short to long term. Short term perspectives. On the theoretical side, we plan to establish the approximation conjecture possibly with modifications of the current hypotheses. On the implementation side, while we have formally described Stateflow Charts in our Simulink semantics, they are not yet integrated into Cosmos. Concerning case studies, we plan to extend our models to other interesting simple situations: for instance, the presence of pedestrians or lane reductions.

Middle term perspectives.

We have defined the communication between models in an ad hoc way, which so far only applies to stochastic high-level nets and Simulink. We plan to define, or use, a pivot formalism to combine several models that have been designed using various formalisms. This could be based on the FMI 2 [START_REF]Functional Mock-up Interface[END_REF] standard of the Modelica Association, in a manner similar to its application to Uppaal [START_REF] Jensen | Integrating Tools: Co-simulation in UPPAAL using FMI-FMU[END_REF].

Cosmos has been previously used for the verification of DNA Walkers Circuits [START_REF] Barbot | On quantitative modelling and verification of dna walker circuits using stochastic petri nets[END_REF]. So the work performed during this thesis could be also used in bioinformatic for more involved case studies. It would require to introduce a dedicated formalism for designing, for example, regulation circuits, and instanciate them in the multi-model framework.

Long term perspective. In the IRT SystemX project SVA [START_REF]Simulation of Autonomous Vehicle Safety -IRT SystemX[END_REF], Wei Chen started in september 2016 a thesis called Formal Models for the Conceptualization and Caracterisation of Use Cases for the Autonomous Vehicle to provide a formalism to generate a large number of situations to be tested [START_REF] Wei | An Ontology-based Approach to Generate the Advanced Driver Assistance Use Cases of Highway Traffic[END_REF]. The creation of a tool for generating scenarios from this new formalism, and the analysis of such scenarios, could be the subject of another thesis.

Figure 1 . 2 :

 12 Figure 1.2: An example of Kripke structure : a coffee machine

Figure 1 . 3 :

 13 Figure 1.3: Another example of a Kripke structure The word (∅∅∅{p, q}) ω is generated by the Kripke structure shown in Figure 1.3.

Figure 1 . 4 :

 14 Figure 1.4: A realization of a discrete event stochastic process

Figure 1 . 5 :

 15 Figure 1.5: Tandem Queue

  Figure 1.6, target((n, m), e, d) is defined as follows for enabled(s): target((n, m), e, d) = 1, m) if e = in, (n -1, m + 1) if e = move, (n, m -1)

Figure 1 . 7 :

 17 Figure 1.7: The graph of a finite DTMC

Figure 1 . 8 :

 18 Figure 1.8: The graph of a finite CTMC

Figure 1 . 9 :

 19 Figure 1.9: A Petri net

  This firing is denoted by m t → m . Let σ = t 1 . . . t n ∈ T * be a sequence of transitions. Then σ is firable and and leads to marking m (written m σ → m ) if there exists a sequence of markings m 0 . . . m n such that m = m 0 and m = m n and ∀0 ≤ k < n, m k t k+1 → m k+1 . We denote by Reach(N , m 0 ) = {m|∃σ ∈ T * : m 0 σ → m} the set of reachable markings.

Example 8 .

 8 Using this notation, 2 • b + c and 2 • a + b + 3 • c are multisets over the set {a, b, c}.Definition 14 (High-level Petri net)A high-level Petri net is a 8-tuple N = (P, T, W -, W + , m 0 , C, Domain, G) with: • P and T the respective non-empty sets of places and transitions, with P ∩T = ∅;• C = {C 1 , . . . , C k } a setof finite color classes; • Domain : (P ∪ T ) → C * a function that associates with each place and transition a sequence of color classes, representing a cartesian product C i 1 × • • • × C i l ; • W -and W + are respectively forward and backward incidence matrices. For each p ∈ P and t ∈ T : W -(p, t) and W + (p, t) are mapping from Domain(t) to Bag(Domain(p)); • G the guard is a mapping from T such that G(t) is a subset of Domain(t); • m 0 the initial marking with m 0 (p) ∈ Bag(Domain(p)).

Figure 1 . 10 :

 110 Figure 1.10: A simple high-level Petri net A simple high-level Petri net is shown in Figure 1.10, with two color classes: Ressources and Processes. The domain of place P 3 is the cartesian product Processes × Resources. Firing transition Get corresponds to a process acquiring a resource (like a file). In the low part of the figure, a new marking is shown that is obtained after firing Get with x bound to process b and y bound to resource 2.

Figure 1 . 12 :

 112 Figure 1.12: Two examples of Büchi automata

Figure 1 . 13 :Figure 1 .

 1131 Figure 1.13: The beginning of an execution tree Figure 1.13 shows the beginning part of the execution tree for the Kripke structure of Figure 1.2. The nodes are labelled by the corresponding states of the Kripke structure.

Figure 1 . 14 :

 114 Figure 1.14: An example of HASL formula A HASL formula is described on Figure 1.14 over the tandem queue system, represented by the DESP shown in Figure 1.6, equipped with the indicator Size :(n, m) → n + m. This formula contains an automaton with three states ( 0 , 1 , 2 ) and an expression E(LAST (r)) which describes the proportion of time during which too many people are in the queues (over 10 simultaneous clients), using a overall duration of 100 time units. This automaton handles three variables: • x 1 representing the elapsed time; • x 2 representing the time during which the queues were too full; • and r representing the time proportion over which the queues were too full. The initial state 0 corresponds to the system running without overfilling. In this state, only x 1 evolves over time. In state 1 , the variable x 2 also evolves, measuring the time during which the queues are overfilled. The measure ends in state 2 , with the maximum simulation time. There are two types of transitions: those who are fired by a process event, shown in blue, and the autonomous transitions fired by a condition over variables (in this case, the execution time) represented in red. In the expression E(LAST (r)), the operator LAST is used to evaluate r at the end of an execution, and E the expectation (over successful runs).

Proposition 3 .

 3 A word is recognized by the Büchi automaton K ⊗ B if and only if it is recognized by both the Kripke structure K and the Büchi automaton B, id est: L(K ⊗ B) = L(K) ∩ L(B).

Figure 1 . 15 :

 115 Figure 1.15: Product of the Kripke structure and the Büchi automaton

Figure 1 . 16 :

 116 Figure 1.16: Does the realisation (x, y) belong in the quarter-circle?

  7 [46] is a tool for qualitative and quantitative analysis of generalised stochastic Petri nets, that relies on Interval Decision Diagrams to represent symbolically the state space of the Petri net. It also has a simulation engine for the model checking of PLTL formulas. Marcie has been mainly developed for the study of chemical reaction networks.ToolEvaluation. The experiments were run with Cosmos, Prism (version 4.0.2), Uppaal-smc (version 4.1.13), Plasma (version 1.1.4), Ymer (version 3.1) and Marcie (version 1178M). Two models were considered:

Figure 1 . 18 :

 118 Figure 1.18: Comparison of simulation time for probability measures

Figure 1 . 19 :

 119 Figure 1.19: Runtime for a probability measure of the TQS model for Sequential Testing

Contents 2 . 1

 21 Geometry-based approaches . . . . . . . . . . . . . . . . 40 2.1.

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Two controller models

Figure 2 . 3 :

 23 Figure 2.3: A motorway segmentA segment of the motorway consists of n positions (as shown in Figure2.3).A position is defined by its coordinates (x, y) with y ∈ 0, n -1 the lane and x ∈ 0, -1 the horizontal position.Physical parameters (such as the maximal speed vmax, maximal acceleration amax, minimum entry speed vmin and minimum delay between two entering vehicles) are assumed to be related to the motorway section. Moreover, it is supposed that state and position changes are synchronous. We give below a specification for relevant information of a vehicle and the configuration of a section:A vehicle present in the section is a tuple a = a.x, a.y, a.v, a.c where: (a.x, a.y) is the position of a; a.v ∈ 0, vmax is its current horizontal speed; and a.c ∈ B a boolean denoting whether the vehicle is controllable. A configuration of the section is a tuple s = d, A where A is a finite set of vehicles, and d is an integer array indexed by the lanes such that d[i] denotes the time elapsed since the last entrance on lane i when d[i] < dmin, and d[i] = dmin otherwise.A possible configuration is shown in Figure2.3. We consider a motorway segment of one kilometer, with two lanes and = 100 positions. There are at most 25 vehicles on each lane, due to the safety distance. The speed values are in 0, 4 , which corresponds to a maximum speed of 40m.s -1 . The possible acceleration values are -1, 0 and 1.

  game. A configuration of the game is a tuple s = d, A where A is a finite set of vehicles, and d is an integer array indexed by the lanes such that d[i] denotes the time elapsed since the last entrance on lane i when d[i] < dmin, and d[i] = dmin otherwise.

Figure 2 . 4 :

 24 Figure 2.4: Risky behaviours to be avoided

Figure 2 . 5 :

 25 Figure 2.5: Illustrating the saturation algorithmIn Figure2.5, we present a simplified version where the states are partitioned into two sets: the states where the environment can play (represented with squares), and those where the controller can play (represented by circles). The states that lead to a collision are represented in red. First, we have Sbad e = S f ail = {q crash } and Sbad c = ∅. During the first step, the controller state q 4 that can only lead to a crash is added to Sbad c . Then, the environment states q 2 and q 3 leading to q 4 ∈ Sbad c are added to Sbad e . Finally, the controller state q 1 which can only lead to Sbad e is added to Sbad c . The green

Figure 2 . 6 :

 26 Figure 2.6: CICAS-SSA hierarchical heterogeneous verification tree ([59], p. 99)

3 )

 3 Models M 21 , M 22 and M 23 are constructed from M 13 by enforcing the choice of the vehicle: in M 21 it has to go left, in M 22 to go straight, and in M 23 it has to go right. These are similar to the model shown in Figure 2.7 except we use only one of the three right branches, and enforce its variable name on the waiting state.

Figure 2 . 7 :

 27 Figure 2.7: Verification model M 13 ([59], p. 101)

Figure 2 .

 2 Figure 2.10: IntersectionUsing this decomposition, an automaton describing the behaviour of the vehicle is built (see Figure2.11), considering the need of a minimal speed for the crossing. The intersection is considered as a resource. The goal is to check that every reachable state satisfies a safety property, such as safety distance preservation.

Figure 2 .

 2 Figure 2.11: Vehicle Control Automaton

  block type is a tuple BT = (n, m, (op i ) 1≤i≤n , b c , b l , b i , b s , Param) where: • n and m are the numbers of output and input signals respectively; • (op i ) 1≤i≤n is a tuple of operators, one for each output signal; • b c is a boolean indicating if the block is discrete or continuous; • b l ∈ {imm, inf, pos} indicates if the block is immediate, with an infinitesimal or a positive latency; • b s is a boolean indicating if the block is a threshold block; • Param is a set of additional parameters depending on the block type, including a sampling delay δ for discrete synchronous blocks and a value r > 0 for blocks with a positive latency. A generic block type over a set G is a family (BT g ) g∈G of block types.

Figure 3 . 1 :

 31 Figure 3.1: Three classical block types

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: Blocks without inputs

  be extended in a similar way to Add.

  Figure 3.4: Sampling blocks

d) Continuous blocks i 1 (Figure 3 . 5 :

 135 Figure 3.5: Continuous blocks

Figure 3 . 6 :

 36 Figure 3.6: Threshold blocks Block Switch is an immediate continuous block with two parameters ∈ {>, ≥} and v ∈ R (with => and v = 0 on Figure 3.6(a)). It compares the value of the critical input i c with v according to the chosen comparison , then returns the value of i 1 or i 2 accordingly. Its operator is op(i 1 , i c , i 2 )(t) = if i c (t) v then i 1 (t) else i 2 (t).More generally, we define a generic block type Switch Tp 1 ,Tp 2 with Tp 1 the type of the input signals i 1 , i 2 and the output signal, and Tp 2 the type of the conditional input. If the output type Tp 1 is a boolean type, then the block is an asynchronous discrete block.Block Relay (Figure3.6(b)) is an immediate asynchronous discrete block with four parameters: v on , v off , s on , and s off . It represents a bang-bang controller which outputs v on when it is activated and v off otherwise. It has one input and one output and its current state can be viewed as an internal signal. Its unique operator is defined as follow:

Figure 3 . 7 :Definition 36 .

 3736 Figure 3.7: Block Chart A (Stateflow) Chart is a discrete asynchronous block describing a deterministic state machine, where transitions are equipped with guards and behaviours are triggered by events. An event is defined by an input signal and an activation mode: Definition 36. An event e is a pair (x, tr) composed of a boolean input signal x = src(e) and a trigger type tr = trig(e) ∈ {↑, ↓}. We denote by E the set of events.

Definition 37 .Figure 3 . 8 :

 3738 Figure 3.8: A signal and its ↑ event

Figure 3

 3 Figure 3.9: A ChartA Stateflow Chart is depicted in Figure3.9. It has two states Q = {q 1 , q 2 }, two variables X = {x, y}, two events E = {e, e } with trig(e) =↑ and trig(e ) =↓, and an initial action y ← 0. For this Stateflow Chart, no priority function is required as there is only one possible transition at each state. Finally, while event e is not linked to any transition, it can be used to activate the Stateflow Chart.

Figure 3 . 10 :

 310 Figure 3.10: Sequence of activation times

3. 1 . 4 Example 25 .Figure 3 . 11 :

 1425311 Figure 3.11: An SK -model In this figure (and later on) we use intuitive abbreviations: init 2 is the initial value

Figure 3 . 12 :

 312 Figure 3.12: Graph of the SK -model of Figure 3.11

Figure 3 . 13 :

 313 Figure 3.13: Illustrating the correctness of SK -models

Figure 3 .Example 29 .graph of block B 3 Figure 3 . 14 :

 3293314 Figure 3.14: A polynomial differential equation system encoded into a SK -model

Figure 3 .Figure 3 . 16 :

 3316 Figure 3.15: A model for a sinusoidal function For the model depicted in Figure 3.15, the equations are ẋ1 = x 3 , ẋ3 = -x 1 , which lead to the characteristic equation of a sinusoidal function ẍ1 = -x 1 .

Figure 3 Figure 3 . 7 xFigure 3 . 18 :Figure 3 . 19 :

 337318319 Figure 3.17: A model for which constraint at t = 1 cannot be solved Let T = [0, 2]; we have t 0 = 0. The mode at t 0 is mode 0 = {] -∞, 1]}. The differential equation corresponding to that mode is ẏ = 1. We choose t 1 = 1 based on this differential equation and the condition. The solution y(t) = t satisfies this mode for ]0, 1[; we can easily deduce the output values of each blocks on [0, 1[. On [1, 1 + ε T [, y(t) = t is still a solution of the differential equation ẏ(t) = 1. However, for any chosen value of t 2 , this solution does not satisfy mode {] -∞, 1]} over ]1, t 2 [. The model semantics is fail over [0, 2].

Figure 3 . 20 :

 320 Figure 3.20: The start of a unfinished trajectory of model Figure 3.18

Figure 3 . 22 :

 322 Figure 3.22: Searching for the next step t i+1

Figure 3 .

 3 Figure 3.23: A Simulink model computing an approximation of cos(t)

Figure 3 . 24 :

 324 Figure 3.24: Backward graphs of the Figure 3.23 Simulink model

Figure 3 . 25 :

 325 Figure 3.25: Applying RK4 method to Simulink model of Figure 3.23

3 Figure 3 . 26 :

 3326 Figure 3.26: Unsafe crossings of value v

Figure 4 . 1 :

 41 Figure 4.1: Overview of the architecture of Cosmos

Figure 4 . 2 :

 42 Figure 4.2: Relation of DEDS classes and the Simulator

1 class

 1 a b st r a ct M a rk i n gI m p l {

Figure 4 . 3 :

 43 Figure 4.3: A Marking structure example

Figure 4 . 4 :

 44 Figure 4.4: A toy Petri net

Figure 4 . 7 :

 47 Figure 4.7: A simplified version of LHA_orig.hpp

Figure 4 . 8 :

 48 Figure 4.8: A simplified version of Simulator.hpp

Figure 4 . 9 :

 49 Figure 4.9: A simplified version of EventQueue.hpp

  • launching simulators working in parallel; • aggregating the simulation results; • computing the values of HASL expressions; class BatchR { public : BatchR ( size_t , size_t ) ; unsigned long int I ; unsigned long int Isucc ; double simTime ; std :: vector < bool > IsBernoulli ; std :: vector < double > Mean ; std :: vector < double > M2 ; std :: vector < double > M3 ; std :: vector < double > M4 ; std :: vector < double > Min ; std :: vector < double > Max ; std :: vector < unsigned long int > bernVar ; void addSim ( const SimOutput &) ; void unionR ( const BatchR &) ; void outputR ( std :: ostream & f ) ; bool inputR ( FILE * f ) ; };

Figure 4 .

 4 Figure 4.10: A simplified version of BatchR.hpp

Figure 4 . 11 :

 411 Figure 4.11: An example of high-level Petri net

Example 39 .Figure 4 .

 394 Figure 4.12: A simplified version of EventQueueSet.hpp

Figure 4 . 2 :

 42 2 is used to model vehicles in a jammed motorway. It has been first presented in[START_REF] Barbot | Statistical Model-Checking for Autonomous Vehicle Safety Validation[END_REF] and will be described in more details in Chapter 5. It has six color classes: PosX, PosY (positions in both axes), VitX, VitY (speed in both axes), AccX and AccY (acceleration in both axes). The controlled vehicle is represented using the color domain PosX × PosY × VitX × VitY × AccX × AccY and the other vehicles using the color domain PosX × PosY × VitX × VitY. Note that the transition collision has a lot of different y, ẋ, ẏ, ẍ, ÿ > < x, y, ẋ, ẏ, ẍ, ÿ > y, ẋ, ẏ > < x, y, ẋ , ẏ , ẍ, ÿ > < x, y, ẋ, ẏ, ẍ, ÿ > exit x > x max < x, y, ẋ, ẏ > enter U N IF (a, b) The SIA model possible bindings. The behaviour of the different vehicles are modelled through C code added to the transitions.

Example 41 .

 41 The code shown in Figure4.15 represent time calculations, using a simple model with a static step (0.2) and a block of delay 0.4. The expected Model functions. As seen in section 4.1.1 that a Model class has to implement three functions: initialEventQueue, update and fire. A Simulink model comes with a unique transition, called SimulinkTransition, which remains in the Event Queue at a time corresponding to the next simulation step. These functions are defined as follows:• initialEventQueue inserts the SimulinkTransition in the Event Queue at time 0 with a minimal priority; • fire computes the current values of each block at given time t, according to the block ordering, then computes the value at t + ε T where ε T is the minimum increment of the SKTime value. • update which computes the next simulation step, starting with the default step size, and refining it due to the Runge-Kutta-Felhberg integration method and the threshold crossing search. These functions rely on basic ones: • executeBlock evaluates the value of the outputs signals of the given block at the given simulation step. It calls another function, findLatencyIndex which returns the index corresponding to the delay of blocks with latency; • computeBkwds computes the output values of the backward graph of the given block at the given simulation step; • estimateIntegrators performs the computation of the Runge-Kutta-Fehlberg integrator over the whole model, using the dedicated cells. It also generates a new approximation of step size if necessary.

  (b)) are output signals of an SK -model and the output arcs are overwriting arcs connected to places that can only be connected to ordinary transitions of the net by read arcs. Similarly to SK -in transitions, the output arcs are labelled by functions of the incoming signals. Such a transition is activated at every sampling time of the SK -model. Upon firing, it rewrites the contents of the output places according to the evaluation of the function.

  An SK -out transition

Figure 4 . 17 :

 417 Figure 4.17: Petri net/Simulink Interface transitions

Figure 4 . 18 :

 418 Figure 4.18: Petri net handling faults and repairs of a double heater

Figure 4 .

 4 Figure 4.19: A Simulink model computing differential equations for the double heater. It contains four parts: the two heater temperatures, the external temperature and the block completing the differential equation

2 Figure 4 .

 24 Figure 4.20: A trace of simulation: T and T ext represent the inside and outside temperatures, Op i corresponds to heater i being operational and On i corresponds to heater i being switched on. Gray areas highlight failure of at least one heater when T ext < T min .

I 1 : 20 l 2 :

 1202 AVG(Min(T)) I 2 : AVG(Last(tc)) I 3 : AVG(Mean(T)) I 4 : AVG(Mean(Active 1 + Active 2 )) I 5 : AVG(Mean(Idle)) (a) HASL formulas l 1 : ṫc = 0 T ≥

Figure 4 . 21 :

 421 Figure 4.21: HASL specification for performance indices

  02, c 2 = 0.013 and c ext = 0.04. Results are reported in Table

Figure 5 . 1 :Figure 5 . 2 :Figure 5 . 3 :

 515253 Figure 5.1: A motorway segment

Figure 5 . 4 :Figure 5 . 5 :

 5455 Figure 5.4: Function move_self

Figure 5 . 7 :

 57 Figure 5.7: Two HASL formulas (a) for the probability of collision and (b) for expected travelled distance

Figure 5 . 10 :Figure 5 . 11 :

 510511 Figure 5.10: Function move_other with smooth speed changes

Figure 5 . 12 :

 512 Figure 5.12: An entrance ramp to a motorway segment

Figure 5 . 13 :

 513 Figure 5.13: Entrance ramp model

Figure 5 . 14 :

 514 Figure 5.14: Function move_self for the Entrance Ramp

Figure 5 . 15 :

 515 Figure 5.15: The HASL formulas for the probabilities of each outcome These three performance indices only measure the kind of collision (if any) that has happened. A fourth performance index, E(LAST(CollForce)) is added, measuring the severity of the collision. The value of CollForce is defined as follows:• if the controlled vehicle reaches the end of the motorway section without collision,CollForce is set to 0; • if the controlled vehicle goes outside the motorway segment, CollForce is set to 1 + v 2 x whwere v x is the longitudinal speed of the vehicle; • if the controlled vehicle collides with another vehicle w, CollForce is set to 1 + |w x -v x | 2 where v x and w x are their respective longitudinal speeds. This choice of values should give a decent approximation: the collision is not too dangerous at a low speed, but the higher the relative speed, the (quadratically) greater the damages are.

Algorithm 7 : 6 if firstvehicle then 7 targetbehind ← vehicle 8 else 9 targetfront ← targetbehind ; 10 targetbehind

 7678910 The targetChoice procedure Data: the ordered vehicle vector vehicles Result: The two target vehicles targetfront and targetbehind 1 firstvehicle ← 1 ; 2 targetfront ← 0 ; 3 targetbehind ← 0 // Id of controlled vehicle 4 while there is another vehicle vehicle in vehicles do 5 if vehicle.py ∈ [5; 10] then // The vehicle is at least partly in target lane ← vehicle ; // We are checking the vehicles in decreasing px order: the new vehicle is behind the previous.

  • A ϕ ≡ ¬ E ♦ ¬ϕ. A state s statisfies this formula if ϕ always holds on all states 22CHAPTER 1. MODELING AND ANALYSIS OF CYBERPHYSICAL SYSTEMS of all paths starting from s ; • A ♦ ϕ ≡ A U ϕ. A state s statisfies this formula if ϕ holds eventually on all paths starting from s. The Kripke structure K satisfies a CTL formula ϕ, written K |= ϕ if its initial state satisfies ϕ, i.e. s 0 |= ϕ.

	Example 13.
	s 0	insert coin	s 1	choose	s 2
	∅		{p}		{p, c}	make coffee
		cancel
		give back change

Figure 1.2: An example of Kripke structure : a coffee machine

  .15.

		¬p, ¬q	¬p, ¬q	¬p, ¬q	p, q
		¬p, ¬q	¬p, ¬q	¬p, ¬q	p, q
	p				p, q
	¬q	¬p, ¬q	¬p, ¬q	¬p, ¬q

Table 1 .

 1 The results on hypothesis testing are reported in Table1.1. Each value is the mean over 100 experiments. The threshold value for the hypothesis is always very close to the

	10	0.17505	0.17	4.78	12.02	3.29	3.36
	20	0.33574	0.33	11.54	23.78	7.48	5.08
	40	0.56931	0.564	21.23	46.47 14.00	7.78
	80	0.81894	0.814	20.10	43.46 13.01	7.60
	200	0.98655	0.981	2.81	8.10	1.92	2.74

Time NumValue p ≥? Uppaal Prism Ymer Cosmos 1: Runtime comparison for the TQS for Sequential Testing the time bound increases. It is a consequence of the evolution of the satisfaction probability of φ T QS , which goes to 1 when T goes to infinity. Cosmos is again the fastest of the tools using a Chernoff-Hoerding method. When the time bound is T = 200, Plasma is approximately 2.6 times slower, Marcie 3.7 times slower, Uppaal 4.2 times slower, and Prism using APMC 6.5 times slower. Among the tools using sequential procedures, when the time bound is 40, Cosmos is still the fastest and the two versions of Prism have similar runtime, but 2.8 times slower than Cosmos.

Sequential hypothesis testing for Tandem Queue System.

  If there is a non-controlled vehicle a, it can be moved according to vehicle movement constraints. Let s 1 = d 1 , A 1 be a configuration (obtained after an environment step), then s = d , A can be reached in one controller transition by moving each controlled vehicle of A 1 according to vehicle movement constraints.Let s = d, A be a configuration. Then s = d , A can be reached by a complete transition if there exists an intermediate configuration s 1 = d 1 , A 1 such that an environment transition leads from s to s 1 and a controller transition leads from s 1 to s . Moreover, for every a 1 = a 2 ∈ A ∩ A :1. a 1 .x = a 2 .x ∨ a 1 .y = a 2 .y (two distinct vehicles do not share position); 2. If a 1 .x < a 2 .x, then either a 1 .x < a 2 .x or (a 2 .x ≤ a 1 .x and a 1 .y = a 1 .y and

dmin) is computed and a new vehicle (a.x = 0, a.y = i, a.v ∈ {vmin, . . . , vmax}, a.c = may be inserted if d 1 [i] = dmin; the delay must then be reinitialised: d 1 [i] = 0; 3.

3.1 Syntax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3

  Differential equations of an SK -model . . . . . . . . . . 68 3.2.2 Interval Partitioning . . . . . . . . . . . . . . . . . . . . 70 3.2.3 Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.

.1.1 Types, signals and operators . . . . . . . . . . . . . . . 56 3.1.2 Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.1.3 Block list . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.1.4 SK -models . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2 Exact semantics . . . . . . . . . . . . . . . . . . . . . . . 67 3.2.1

3 Operational semantics . . . . . . . . . . . . . . . . . . . . 78 55

  The differential equation associated with an integration block B, for which the variable is denoted by x B , is obtained by associating with each output o of each block B of G B -an expression y B ,o defined inductively over the blocks of G B -following the reverse topological order:• if B is a block without input, then y B ,o = op B ,o (); • if B is an integration block, then y B ,o = x B ; • if B is a positive latency block, then y B ,o = h (B ,o) ; • if B isan immediate block with inputs (B 1 , o 1 ), . . . , (B m , o m ) then y B ,o = op B ,o (y B 1 ,o 1 , . . . , y Bm,om ). The differential equation is then defined by ẋB = y B -,o-.

	Example 30.
	1
	s

  state( ) and state( )

	m>n,{x←x-1;n←m}

  An SK-model M satisfies the Safe Crossing condition if there exists d min and α > 0 such that, for any output signal s of M, for any threshold value v, if s(t c ) = v for some t c ∈ T, then:1. s is continuous over the interval [t c -α, t c + α]; ≥ d min (the right and left derivatives of the signal on the crossing are above the minimal value).

	2. For any β ≤ α, (s(t c -β) -v)(s(t c + β) -v) < 0 (the threshold value is
	effectively crossed);
	3. | ds dt (t + c )| ≥ d min and | ds dt (t -c )|
	Such a result would require reasonable hypotheses on the model. We now propose
	a semantical condition on signals, called Safe Crossing, under which we conjecture
	that the result holds:
	Definition 57 (Safe Crossing)

  .2.

	<<interface>>
	DEDS		Simulator
	initialEventQueue() fire(t: transition, b: binding,	1	SimulateOneStep(): boolean SimulateSinglePath()
	time: real)			interactiveSimulation()
	update(time: real, t: transition)
	SPN	SKModel
	marking : stateImpl	state : stateImpl
	initialEventQueue()	initialEventQueue()
	fire(t: transition, b: binding,	fire(t: transition, b: binding,
	time: real)	time: real)
	update(time: real, t: transition)	update(time: real, t: transition)

Table 4 .

 4 1: Benchmarks for the Hospital model

			Standard	Alternative
	Cells Vehicles Time Memory Time Memory
	50	10 418.9 s	2.60GB	4.09 s	4.06MB
	50	100 1 610 s	2.60GB	9.27 s	4.07MB
	500	10	Out of Memory	140.14 s 5.42MB
	500	100	Out of Memory	77.93 s	5.57MB
		Figure 4.14: Benchmarks for the SIA model

  • C and 25 • C. The system is modeled by a stochastic Petri net (Figure4.18) with randomized faults and repairs. The evolution of room temperature and heater behaviours are hybrid and thus are modeled in Simulink (Figure4.19). The fault transitions of the net have an exponential time distribution (with different rates). The repairman, initially at the Idle state, randomly chooses which (faulty) heater he will repair, then proceeds

	Op 1	Fault 1	TR 1	Start 1	Idle Start 2	TR 2	Fault 2 Op 2
	In 1		Active 1	Active 2	In 2
				Done 1	Done 2

11 if targetfront

  .pxtargetbehind.px ≥ 3 and targetbehind.px < px then return; // The distance is good enough for an insertion. [0.112, 0.121] [0.069, 0.078] [8.1 • 10 -4 , 9.3 • 10 -4 ] [0, 3.4 • 10 -4 ]

			Controller Version 1
		10	10	10	20
	I 1 [0.016, 0.019] [0.043, 0.047]	[0.041, 0.046]	[0.126, 0.133]
	I 2 [0.053, 0.057] [0.067, 0.073]	[0.070, 0.076]	[0.035, 0.039]
	I 3 [0.925, 0.929] [0.881, 0.887]	[0.880, 0.887]	[0.829, 0.838]
	I 4 [0.016, 0.019] [0.100, 0.112]	[0.100, 0.111]	[0.184, 0.197]
	T O
	Simulations	151 000	81 000	83 000	71 000
	Time	1 672s	1 014s	1 519s	3 601s
			Controller Version 2
		10	10	10	20
	I 1 [0.080, 0.089] [0.144, 0.155]	[0.182, 0.185]	[0.202, 0.220]
	I 2
	I 3 [0.793, 0.805] [0.770, 0.783]	[0.815, 0.817]	[0.780, 0.798]
	I 4 [0.259, 0.286] [0.466, 0.504]	[0.601, 0.607]	[0.486, 0.537]
	T O Simulations	37 000	34 000	2 000 000*	16 000
	Time	515s	486s	39 041s	828s
		Non meaningful	Zero by construction
	N

*

Reached maximum simulation number

Table 5 .

 5 2: Evaluation of controllers (with 99% confidence interval and 10% relative width)

Cosmos is available at http://cosmos.lacl.fr/.

Prism is available at https://www.prismmodelchecker.org/.

Uppaal is available at http://www.uppaal.org/

This extension is available at http://people.cs.aau.dk/~adavid/smc/

Plasma is available at https://project.inria.fr/plasma-lab/

Ymer is available through a git repository at https://github.com/hlsyounes/ymer

Marcie is available at http://www-dssz.informatik.tu-cottbus.de/DSSZ/Software/ Marcie

nearest by time-to-intersection, determined from positions and approach velocities

https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern

http://www.di.unito.it/~amparore/mc4cslta/editor.html

see https://www.graphviz.org/doc/info/lang.html

ACKNOWLEDGMENTS

Conclusion.

This approach has been validated by several case studies. It can be generalised to the case of multiple trajectories provided by the planner. In this case, the controller selects the one minimising intersection. [START_REF] Rizaldi | On time-memory trade-off for collision detection[END_REF] Objective. The objective is to provide an approach to reduce the computation time of collision detection in order to match human reaction delay.

Time-Memory Tradeoff for Collision Detection

Assumptions. There are two vehicles: the controlled one ego and another vehicle other of fixed sizes, both represented by rectangles located on a given area.

We suppose that the generated trajectory is represented by a sequence of connected segments of fixed length δ > 0. Moreover, for two adjacent segments, the difference between their orientations is assumed lower than φ max ∈ [0, π 2 [. Furthermore, the time evaluation interval [t 0 , t n ] is splitted in intervals [t k-1 , t k ] for k ∈ 1, n .

Database. The database consists of a single relation indicating whether the two rectangles intersect or not. The value domain is finite, due to two factors: the boundedness of the area, and the discretisation of the positions. As the size of the vehicle is supposed fixed, the fields of the relation are the center and orientation of each vehicle, and a boolean indicating if there is an intersection. Verification of trajectories. Two trajectories are considered: the trajectory of the controlled vehicle, and the center of the lane on which the other vehicle is located. The safety verification of the trajectory repeats the following steps, until t n :

1. compute reachability and occupation sets of the vehicles, taking into account their dynamics and the uncertainties over [t k-1 , t k ]; 2. cover the occupation set of the vehicles by rectangles; 3. checks if there is a collision using the database.

Conclusion.

A benchmark of the procedure has been done using a computer equipped with a 2.80GHz i5-4330M processor, 12GB of RAM and the following librairies and softwares: • for the construction of the database, MATLAB 2014a (and it required 23.37MB of memory); • for time measures, the Windows API (QueryPerformanceCounter),

• for collision detection algorithms, a C++ program compiled using Microsoft Visual C++ 2013 (with the O2 optimisations).

Algorithm 1: Function updating the Event Queue for a Petri Net Data: the Event Queue EQ, a Time Generator T G 1 Function update(ctime,lastT r,EQ,T G) Input : the current time ctime, the last activated transition lastT r, the Event Queue EQ and Time Generator T G • in the case of a E(f ) HASL formula, it computes the mean value by dividing Mean by the number of successful simulations. Then, it uses the M2 value to obtain the variance (V = m2 -mean 2 + 1 Succ ), using also the number of succesful simulations (needed in the Chows-Robbin algorithm). The width of the confidence interval is then obtained using the variance and the wanted confidence level; • in the case of a sequential hypothesis testing, it constructs a confidence interval (to be handled by the other functions of Cosmos). It uses a logarithmic version of the likelihood ratio of the Wald method, for probabilities v1 and v2. If this logarithmic value is negative, then returns the confidence interval [0, v2]. In the other case, [v1, 1] is returned.

Using these confidence intervals, an approximation on the progress of the simulation is made, using the ratio between the current width of the interval and the target width. Note that it is possible to use a relative width. Finally, continueSim checks whether the simulation should be continued:

• if the number of simulation exceeds the maximum number of simulations, then it returns false; • otherwise, if there is no target width and no sequential testing, then it returns true; • otherwise, if the the minimum value of all progress ratios is smaller than 1, then it returns true; • otherwise, it returns false.

Note that there is a slight difference in sequential hypothesis testing, as we collect results by batch of simulations instead of one simulation at a time. However, one may consider that this generated bias is small enough in most case, the size of a batch being small enough compared to the final number of simulations.

Algorithm 6: Function updating the Event Queue for a high-level Petri Net

Data: the Event Queue EQ, a Time Generator T G 1 Function update(ctime,lastT r,EQ,T G)

Input : the current time ctime, the last activated transition lastT r, the last associated binding lb, the Event Queue EQ and Time Generator T G Note that the number of possible bindings can increase exponentially with respect to the size of the expression labelling the edge. This was already observed and lead to some optimization when the color functions are simple (like single variables).

Alternative Implementation

Changes to the Marking class. Instead of using a map in dom_Domain, we now use an ordered set of (dom_Token,mult) pairs. In C++, the set data structure is implemented with red-black trees.

Example 40. For the model shown in Figure 4.11, the initial marking would now be represented this way:

In this new setting, the memory used to store the marking is linear with respect to the number of tokens.

The EventQueueSet class. This class (see Figure 4.12) is similar to the EventQueue class as it still maintains three tables: • evtTbl maps a transition and a binding hash to an Event; • evtHeap is now a vector of Events, organized as a heap; • evtHeapIndex maps a transition and a binding hash to the position of the corresponding event in the heap (and -1 if there is no corresponding event). Note that these definitions are different from those presented in section 4.1.2.

Binding Iteration. This time, the function nextPossiblyEnabledBinding tries to find values for each variable by iterating over the content of each input place. For remaining variables, it iterates over their color domain.

Benchmarks

In this section we describe two models, one of the patient flow in a hospital (with one token per patient, each having a different color) and the other of a small section of a motorway (with one token per vehicle). On these models, we run both implementations on one core of an AMD Turion(tm) II Ultra Dual-Core Mobile M640. This choice of a single core compared to a parallel implementation has no impact on the relative performances. The simulation parameters have been chosen so that all the simulation runs in less than 30 minutes. Moreover, the memory consumption has been tested on the simulator alone (using /usr/bin/time -v). behavior is that both time and target values are equal to respectively 1.2 and 0.8 after the 6 steps of the loop. However, running the executable generated from this code prints Greater or equal to 1.2 and Smaller than 0.8. This is explained by the fact 0.2 has no finite representation as a binary floating-point value.

Model state. The state of a Simulink model in Cosmos is represented by a class stateImpl using vectors indexed by the simulation step. There is one such vector for the time values of the simulation (in the SKTime class), and one vector by output signal for every block. The integer lastPrintEntry stores the current step of simulation. Note that this class can be used, like markingImpl, as an input for any function requiring a DEDState.

In order to deal with two particular features, these vectors have more cells than the number of simulation step. Cosmos has nine additional cells: • the buffer for application of the Runge-Kutta-Fehlberg method and its intermediate values requires six cells; • the computation of the time of threshold crossings requires three cells: current lower and upper bounds, and new candidate value.

CHAPTER 5

CASE STUDIES

Abstract. The increasing development of autonomous vehicles leads to safety challenges for the involved companies like Waymo and Tesla. The goal of the SVA 1 industrial project [START_REF]Simulation of Autonomous Vehicle Safety -IRT SystemX[END_REF] hosted by IRT SystemX is to provide methods, based on simulation, for assessing the safety of autonomous vehicles. The members of the project work in cooperation with major french car manufacturers such as Renault or PSA Peugeot Citroën, but also suppliers of vehicle components such as Continental or Valeo. In this section, we apply the earlier chapter developments to two relevant case studies. In the first one, we describe the study of a simple motorway segment, published in an industrial conference [START_REF] Barbot | Statistical Model-Checking for Autonomous Vehicle Safety Validation[END_REF]. The second case concerns an entrance ramp. It is more sophisticated since it requires a more complex decision-making algorithm for the autonomous vehicle.

In order to provide an accurate modeling, we have identified the quantitative and qualitative caracteristics of the scenarios: 1. the granularity of the representation of the lane(s), 2. the random behaviour of environment vehicles, 3. the information given to the controller, and 4. the controller algorithms to be analysed.

Contents

tokens of color domain OtherVehicle by x, y, ẋ, ẏ .

General structure. An ideal model would use the formalism of Cosmos as a specification language, all the features of the case study being specified only using high-level Petri nets. However, some features are currently not easily expressible with such nets, for instance the random generation over numbers inside a token. Furthermore, the controller algorithm is more naturally specified by a C code.

The model is shown in Figure 5. Synchronisation of simulation steps. Transitions sim_step1, sim_step2, sim_step3 act as schedulers for the simulation parts: sim_step1 has a Dirac distribution with a parameter corresponding to the time between two steps, while sim_step2 and sim_step3 are immediate, thus enforcing the synchronous behaviour of the simulation. The priorities of the transitions enforce the order of the simulation. There must be a token in step1 for random_other to be fired, and this token will not leave step1 until all tokens of otherVehicles have been processed. Then, this token is consumed by sim_step2, producing a token in step2 which allows update_other to be fired. Similarly, the step2 token cannot leave the place until all tokens of inProcess have been processed. Then, this token is consumed by sim_step3 which produces a token in step3, allowing update_self then Controller to be fired in this order.

Note that, as enter, exit, collision and end have a zero priority, those transitions cannot be fired during the synchronous step of the simulation, but only once all steps have been completed.

Movement of other vehicles. The other vehicles are described using tokens of color domain OtherVehicle. Due to the C function move_other, each firing of the transition update_other produces an updated token with new speed and position in place inProcess. Then, in the next simulation step, all tokens from place inProcess are moved back to otherVehicles. As the definition of the behaviour of other vehicles is part of the definition of the scenario, it is described in section 5.1.3. to the mode 0, only if it has reached the destination lane; • If the vehicle is in mode 1 (micro-sleep), there is at each time step a probability 25% to return in default mode. Similarly, in mode 2 (phone call), we assign a probability 20%. • If the vehicle is in mode 0, it can change its mode either to mode 1 (with probability 2.5%), mode 2 (with probability 0.5%), or mode 3 (with probability 1%, only if it is not yet in the middle of the right lane).

Moreover, we consider two possible situations for the other vehicles: • In the first situation, a vehicle leaving the motorway from the main lanes is replaced by another one entering it (with the same ID, same speed, and same lane). This keeps the number of vehicle constants. This situation is depicted with the pictogram in the tables; • In the second case, a vehicle leaving the motorway from the main lanes is simply removed.

Detailed controller specification. The controllers perform four consecutive tasks, the current one being stored in the currstate variable as described below: 0. accelerating in the entrance ramp before the entrance location;

Summary

Based on several case studies of autonomous vehicle situations, we have identified the need for an expressive formalism of probabilistic hybrid systems. Rather than building it from scratch, we decided to combine stochastic high-level Petri nets with Simulink. The interest of high-level stochastic Petri nets comes from the fact that it is supported by the performant tool Cosmos, while Simulink is a de facto standard for the design of vehicle controllers. Furthermore, the features of the two models are complementary. Generalising our objectives, we proposed to follow a multi-model approach to combine several formalisms. In order to empirically validate this approach by simulations over examples, this multi-model formalism had to be integrated into Cosmos. This required several steps:

• providing a formal syntax and exact semantics for (a significant fragment of) Simulink. The block-diagram models are transformed into a set of differential equations over subintervals of the simulation time, taking into account the modes related to threshold crossings and the delays of blocks with latencies; • providing an operational semantics intended to accurately approximate the exact semantics, when the trajectory exists. Our approach relies on several approximation methods, such as Runge-Kutta integration method, zero-crossing detection, and linear interpolation for latencies; • implementing three important upgrades in Cosmos:

adding a new binding mechanism for high-level Petri nets, that largely improves the performance for nets with large color classes but small number of tokens; -implementing the operational semantics of Simulink; -adding hooks for multi-model simulation giving a communication framework between different models. In the case of stochastic net and Simulink cosimulation, these hooks were instanciated as special transitions in the net. All these upgrades were illustrated and tested on several toy examples. We finally proposed two significant case studies for the validation of this Abstract : This thesis takes place in the context of autonomous vehicle design, and concerns more specifically the verification of controllers of such vehicles. Our contributions are the following : (1) give a syntax and a semantics for a hybrid system model, (2) extend the capacities of the model-checker Cosmos to that kind of models, and (3) empirically confirm the relevance of our approach on typical case studies handling autonomous vehicles. We chose to combine high-level stochastic Petri nets (which is the input formalism of Cosmos) with the input formalism of Simulink, to obtain an adequate expressive power. Indeed, Simulink is largely used in the automotive industry and numerous controllers have been specified using this tool. However, there is no formal semantics for Simulink, which lead us to define such a semantics in two steps : first, we propose an exact (but not operational) semantics, then we complete it by an approximate semantics that includes the targeted approximation level. In order to combine the discrete event model of Petri nets and the continous model specified in Simulink, we define a syntactic interface that relies on new transition types ; its semantics consists of an extension of the simulation loop. The evaluation of this new formalism has been entirely implemented into Cosmos. Using this new formalism, we have designed and studied the two following case studies : on one hand, a heavy traffic on a motorway segment, and on the other hand the insertion of a vehicle into a motorway. Our approach has been validated by the analysis of the corresponding models.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France