{'fraction_non_alphanumeric': 0.04582871429582282, 'fraction_numerical': 0.010708304832889119, 'mean_word_length': 3.806484525563995, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Approximate probabilistic verification of hybrid systems

Hybrid systems whose mode dynamics are governed by non-linear ordinary differential equations (ODEs) are difficult to analyze. To address this, we develop an approximate probabilistic analysis method. We do so by approximating the mode transitions as stochastic events. More precisely, we assume that the probability of making a mode transition is proportional to the measure of the set of pairs of time points and value states at which the mode transition is enabled. To ensure a sound mathematical basis, we impose a natural continuity property on the non-linear ODEs. We also assume that the states of the system are observed at discrete time points but that the mode transitions may take place at any time between two successive discrete time points. This leads to a discrete time Markov chain as a probabilistic approximation of the hybrid system. We then show that for BLTL (bounded linear time temporal logic) specifications the hybrid system meets a specification iff its Markov chain approximation meets the same specification with probability 1. Based on this, we formulate a sequential hypothesis testing procedure for verifying -approximately- that the Markov chain meets a BLTL specification almost certainly. Our two case studies indicate that our scheme can be applied in a number of realistic settings.


INTRODUCTION
The analysis of hybrid systems is difficult. The high expressive power of the mixed dynamics renders even basic reachability questions undecidable [18]. Various lines of work have explored ways to mitigate this problem [15,13,9,3,2,19]. A common technique is to restrict the continuous dynamics. However, for many of the models arising in the domains of cyber physical systems and systems biology the mode dynamics will be governed by a system of (non-linear) ordinary differential equations (ODEs). To analyze such systems, we develop here a probabilistic scheme under which a rich class of hybrid systems with non-linear dynamics can be approximated as a discrete time Markov chain.
A key difficulty in analyzing a hybrid system's behavior arises from the interaction between the continuous mode dynamics and the discrete mode transitions governed by the guards. Specifically, the time points and value states at which a trajectory meets a guard will depend on the solutions to the ODE systems associated with the modes. For high-dimensional systems these solutions will not be available in closed form. To get around this, we assume that the mode transitions are stochastic events and that the probability of a mode transition is proportional to the measure of the value state and time point pairs at which this transition is enabled. More sophisticated hypotheses could be considered. For instance one could tie the mode transition probability to how long the guard has been continuously enabled or how deeply within a guard region the current state is. It could also be the case that some of the mode transitions are controllable, say, by a digital controller. This will, however, entail MDPs (Markov decision processes) based approximations, and we postpone the study of this more complex setting to the future.
To secure a sound mathematical basis for our approximation, we further assume: (i) The vector fields associated with the ODEs are C 1 functions.(ii) The states of the hybrid system are observable only at discrete time points. (iii) The set of initial states and the guard sets are bounded open sets.(iv) The hybrid dynamics is strictly non-Zeno in the sense there is uniform upper bound on the number of transitions that can take place in a unit time interval. For technical convenience we in fact assume that time discretization is so chosen that at most one mode transition takes place between two successive discrete time points. A rich class of engineered systems and biological processes can be modeled under these restrictions.
In this setting, we construct a Markov chain M as an approximation of the hybrid system H. To relate the behaviors of M and H, we use BLTL (bounded linear-time temporal logic [10]) to specify time bounded dynamic properties of H. We then show that H meets the specification ψ-i.e. every trajectory of H is a model of ψ-iff M meets the specification ψ with probability 1. This allows us to approximately verify interesting properties of the hybrid system using its Markov chain approximation. However, the transition probabilities of the Markov chain will depend on the solutions to the ODEs associated with the modes, which will not be available in a closed form. To cope with this, we design a statistical model checking procedure to approximately verify that the chain almost certainly meets the specification. We do so by randomly generating trajectories of H through numerical simulations in a way that corresponds to randomly sampling the paths of the Markov chain according to its underlying transition probabilities.
In establishing these results, we first assume that the atomic propositions in the specification are interpreted over the modes of the hybrid system. Consequently one can specify patterns of mode visitations while quantitative properties can be inferred only indirectly and in a limited fashion. As an extension we then introduce quantitative atomic propositions ("the current velocity of the vehicle is greater than 40 miles/hour") and show that our results go through -with a more refined construction-provided the dynamics is restricted to robust trajectories. Informally, under a natural topology, a robust trajectory is one that is contained in an open neighborhood of trajectories which satisfy the same specification. The procedure for sampling trajectories can be easily modified to sample for robust trajectories. Finally, we demonstrate the applicability of our method using a room heating system and a model of the activation dynamics of cardiac muscle cells.
Due to space limitations we present only proof sketches for a number of the results. more details can be found in the appendix.

Related work
Mode transitions have been approximated as random events in the literature. In [1] the dynamics of a hybrid system is approximated by substituting the guards with probabilistic barrier functions. Our transition probabilities are constructed using similar but simpler considerations. We have done so in order to be able to carry out temporal logic based verification based on simulations. Another related work is [11], which studies a network of hybrid automata that communicate with each other through input/output actions. Translated into our setting, the main idea roughly boils down to this: the time point in (0, 1), at which the decision about what the mode should be up to the next discrete time point, is determined by the uniform distribution over (0, 1). In our setting this probability is determined by the continuous mode dynamics and the guard sets.
An alternative approach to approximately verifying nonlinear hybrid systems is one based on δ-reals [14]. Here one verifies bounded reachability properties that are robust under small perturbations of the numerical values mentioned in the specification. Since the approximation involved is of a very different kind, it is difficult to compare this line of work with ours. However, it may be fruitful to combine the two approaches to verify a richer set of reachability properties.
The present work may be viewed as an extension of [23] where a single system of ODEs is considered. This method however, breaks down in the multi-mode hybrid setting and one needs to construct-as we do here-an entirely new machinery. Finally, a wealth of literature is available on the analysis of stochastic automata [8,6,21,4]. It will be interesting to explore if these methods can be transported to our setting.

HYBRID AUTOMATA
As usual we will use hybrid automata as a model of hybrid systems. Accordingly, we fix n real-valued variables {xi} n i=1 viewed as functions of time xi(t) with t ∈ R+, the set of non-negative reals. A valuation of {xi} n i=1 is v ∈ R n with v(i) ∈ R representing the value of xi. The language of guards is given by: (i) a < xi and xi < b are guards where a, b are rationals and i ∈ {1, 2, . . . , n}. (ii) If g and g are guards then so are g ∧ g and g ∨ g .
G denotes the set of guards. We define v |= g (i.e. v satisfies the guard g) via: v |= a < xi iff a < v(i) and similarly for xi < b. The clauses for conjunction and disjunction are standard. We let g = {v | v |= g}. We note that g is an open subset of R n for every guard g. We will abbreviate g as g.
• Q is a finite set of modes and qin ∈ Q is the initial mode.
We have not associated invariant conditions with the modes or reset conditions with the mode transitions. They can be introduced with some additional work.
Fixing a suitable unit time interval ∆, we discretize the time domain as t = 0, ∆, 2∆, . . .. We assume the states of the system are observed only at these discrete time points. Furthermore, we shall assume that only a bounded number of mode changes can take place between successive discrete time points. Both in engineered and biological processes this is a reasonable assumption. Given this, we shall in fact assume that ∆ is such that at most one mode change takes place within a ∆ time interval. We note that there can be multiple choices for ∆ that meet this requirement and in practice one must choose this parameter carefully. (Our method can be extended to handle a bounded number of mode transitions in a unit time interval but this will entail notational complications that will obscure the main ideas.) In what follows, for technical convenience we also assume the time scale has been normalized so that ∆ = 1. As a result, the discretized set of time points will be {0, 1, 2, . . .}.

Trajectories
In what follows, we fix a hybrid automaton H as defined above. The behavior of H will consist of its trajectories. To define this notion and for later use, we start with some preliminaries. We recall that a function f : R → R is C 1 if f , the derivative of f , exists everywhere and is continuous. This notion extends to R n in the obvious way. Next, the function F : where | · | is the standard Euclidean norm on R n .
We have assumed that for every mode q, the right hand side of the ODEs, Fq(x), is C 1 for each component. As a result, for each value v ∈ R n and in each mode q, the system of ODEs dx/dt = Fq(x) will have a unique solution Zq,v(t) [20] satisfying Zq,v(0) = v and dZq,v(t)/dt = F (Zq,v(t)). We are also guaranteed that Zq,v(t) is Lipschitz and hence measurable [20]. This fact will play a crucial role in what follows. It will be convenient to work with two sets of functions derived from solutions to the ODE systems.
The (unit interval) flow Φq : (0, 1) × R n → R n is given by Φq(t, v) = Zq,v(t). Φq will also be Lipschitz. Next we define the parametrized family of functions Φq,t : R n → R n given by Φq,t(v) = Φq(t, v). In addition to being Lipschitz, these functions will be bijective as well. Further, their inverses will also be bijective and Lipschitz.
We say that the trajectory τ as defined above starts from q0 and ends in q k . Further, its initial value state is v0 and its final value state is v k . We let T RJ denote the set of all finite trajectories that start from the initial mode qin and with an initial value state in INIT.

THE MARKOV CHAIN APPROXIMATION
A (finite) path in H is a sequence ρ = q0q1 . . . q k such that for 0 ≤ j < k, there exists a guard gj such that qj g j → qj+1. We say that this path starts from q0, ends at q k and is of length k + 1. We let paths H denote the set of all finite paths that start from qin.
In what follows µ will denote the standard Lebesgue measure over finite dimensional Euclidean spaces. We will construct MH = (Υ, ⇒), the Markov chain approximation of H inductively. Each state in Υ will be of the form (ρ, X, PX ) with ρ ∈ paths H , X an open subset of R n of non-zero, finite measure and PX a probability distribution over SA(X), the σ-algebra generated by X.
We start with (qin, INIT, PINIT) ∈ Υ. Clearly, INIT is an open set and it is of non-zero, finite measure since µ(INIT) = i (Ui − Li). For technical convenience we shall assume PINIT to be the uniform probability distribution. This reflects the assumption that each member of INIT is an equally likely initial value state. However we can handle other distributions over INIT as well.
Assume inductively that (ρ, X, PX ) is in Υ with X an open subset of R n of non-zero, finite measure and PX a probability distribution over SA(X). Suppose ρ ends in q and there are m outgoing transitions q g 1 → q1, . . . , q gm → qm from q in H ( Fig. 1 illustrates this inductive step). Then for 1 ≤ j ≤ m we define the triples (ρqj, Xj, PX j ) as follows. In doing so we will assume the required properties of the objects involved in this construction. We will then establish these properties and thus the soundness of the construction. For convenience, through the remaining parts of this section j will range over {1, 2, . . . , m}.
For each v ∈ X and each j we first define the set of time points Tj(v) ⊆ (0, 1) via Thus Tj(v) is the set of time points in (0, 1) at which the guard gj is satisfied if the system starts from v in mode q at time k and evolves according to dynamics of mode q up to time k + t. We next define Xj for each j as Thus Xj is the set of all value states obtained by starting from some v ∈ X at time k, evolving up to k +t according to the dynamics q, making an instantaneous mode switch to qj at this time point, and evolving up to time k + 1 according to dynamics of mode qj.
To complete the definition of the triples (ρqj, Xj, PX j ), we first denote by P T j (v) the uniform probability distribution over Tj(v). Our construction can be easily extended to handle other kinds of distributions as well. We now define the probability distributions PX j over SA(Xj) as follows. Suppose Y is a measurable subset of Xj. Then As usual 1Z is the indicator function of the set Z while dP T j (v) indicates that the inner integration over Tj(v) is w.r.t. the (uniform) probability measure P T j (v) and dPX indicates that the outer integration over X is w.r.t. the probability measure PX . Thus PX j (Y ) captures the probability that the value state Φq j (1−t, Φq(t, v)) lands in Y ⊆ Xj by taking the transition q g j → qj at some time point in Tj(v) given that one started with some value state in X.
Next we define the triples ((ρ, X, PX ), pj, (ρqj, Xj, PX j )), where pj is given by Thus pj captures the probability of taking the mode transition q g j → qj when starting from the value states in X and mode q. For every j we add the state (ρqj, Xj, PX j ) to Υ and the triple ((ρ, X, PX ), pj, (ρqj, Xj, PX j )) to ⇒ iff µ(Xj) > 0.
Finally, (qin, INIT, PINIT) is the initial state of MH . We can summarize the key properties of our construction as follows (while assuming the associated terminology and notations).
is an open set of finite measure for each v ∈ X and each j.

2.
Xj is open and is of finite measure for each j.
4. PX j is a probability distribution for each j.

MH = (Υ, ⇒)
is an infinite state Markov chain whose underlying graph is a finitely branching tree.
Proof. To prove the first part, suppose t ∈ Tj(v). Then Φq(t, v) = v ∈ gj and gj is open. Hence v will be contained in an open neighborhood U contained in gj. Since Φq is Lipschitz we can pick U such that Clearly each Ii has a finite measure. By the Lipschitz continuity of Φq we know that there exists a constant c such that µ(Φq(Ii)) < c · µ(Ii) for all i. We thus have Therefore Φq((0, 1), X) has a finite measure. By a similar argument we can show that Φq j ((0, 1), Φq((0, 1), X)) has a finite measure as well. Since Xj = t Φq j ,1−t(Φq,t(X )∩g) ⊆ Φq j ((0, 1), Φq((0, 1), X)), it must have a finite measure.
The remaining parts of the theorem follow easily from the definitions and basic measure theory. Figure 1: Markov chain construction.

RELATING THE BEHAVIORS OF H AND M H
We shall use bounded linear-time temporal logic (BLTL) [10] to specify time bounded properties and use it to relate the behaviors of H and MH . For convenience we shall write M instead of MH from now on.
We assume a finite set of atomic propositions AP and a valuation function Kr : Q → 2 AP . Formulas of BLTL are defined as: (i) Every atomic proposition as well as the constants true, false are formulas. (ii) If ψ, ψ are formulas then ¬ψ and ψ ∨ ψ are formulas. (iii) If ψ, ψ are formulas and is a positive integer then ψU ≤ ψ is a formula. The derived operators F ≤ and G ≤ are defined as usual: We shall assume through the rest of the paper that the behavior of the system is of interest only up to a maximum time point K > 0. This is guided by the fact that given a BLTL formula ψ there is a constant K ψ that depends only on ψ so that it is enough to evaluate an execution trace of length at most K ψ to determine whether ψ is satisfied [5].
Hence we assume that a sufficiently high K has been chosen to handle the specifications of interest. Having fixed K, we denote by T RJ K+1 the trajectories of length K +1, and view this set as representing the time bounded non-deterministic behavior of H of interest.
To develop the corresponding notion for M , we first define a finite path in M to be a sequence η0η1 . . . η k such that ηj ∈ Υ for 0 ≤ j ≤ k. Furthermore for 0 ≤ j < k there exists pj ∈ (0, 1] such that ηj p j ⇒ ηj+1. Such a path is said to start from η0 and its length is k + 1. We define paths M to be the set of finite paths that start from the initial state of M while paths K+1 M is the set of paths in paths M of length K + 1. The trajectory semantics. Let τ = (q0, v0) (q1, v1) . . .
(q k , v k ) be a finite trajectory, ψ a BLTL formula and 0 ≤ j ≤ K. Then τ, j |=H ψ is defined via: where A is an atomic proposition.
• ¬ and ∨ are interpreted in the usual way.
• The remaining clauses are defined just as in the case of |=H .

The main result
We wish to show that H meets the specification ψ iff Pr ≥1 (ψ). To this end let π = η0η1 . . . η k be a path in M with ηj = (q0q1 . . . qj, Xj, PX j ) for 0 ≤ j ≤ k and let τ = (q 0 , v0) (q 1 , v1) . . . (q k , v k ) be a trajectory. Then we say that π and τ are compatible iff k = k and qj = q j and vj ∈ Xj for 0 ≤ j ≤ k. The following three observations based on this notion will easily lead to the main result. Lemma 1.
2. Suppose π is a path in M. Then there exists a trajectory τ such that π and τ are compatible. Furthermore if π ∈ paths M then τ ∈ T RJ.
3. Suppose τ is a trajectory. Then there exists a path π in M such that τ and π are compatible. Furthermore if τ ∈ T RJ then π ∈ paths M .
The proof follows via a systematic application of the definitions, and details can be found in the Appendix. We can now establish the main result of this section.
Proof. Suppose H does not meet the specification ψ. Then there exists τ ∈ T RJ K+1 such that τ, 0 |=H ψ. By the third part of Lemma 1 there exists π ∈ paths K+1 M which is compatible with τ . By the first part of Lemma 1 we then have π / ∈ modelsM (ψ) which leads to P r<1(ψ).
Next suppose that P r<1(ψ). Then there exists π ∈ paths K+1 such that π, 0 |=M ψ. By the second part of Lemma 1 there exists τ ∈ T RJ K+1 which is compatible with π. By the first part of Lemma 1 this implies τ, 0 |=H ψ and this in turn implies that H does not meet the specification ψ.

THE SMC PROCEDURE
To verify whether H meets the specification ψ, we solve the equivalent problem whether Pr ≥1 (ψ) on M . However, M cannot be constructed explicitly since its transition probabilities are defined in terms of the solutions to the ODEs which will not be available. Therefore we shall use randomly generated trajectories to sample the paths of M and formulate a sequential hypothesis test to decide with bounded error rate whether Pr ≥1 (ψ) holds. Algorithm 1 describes our trajectory sampling procedure. Generate time points T := {t1, . . . , t J } uniformly in (0, 1).

4:
Simulate Let Tj := {t ∈ T : v j ∈ gj } be the time points where gj is enabled.

6:
Pick g randomly according to probabilities {pj

7:
Pick t uniformly at random from T .
The initial value v0 is sampled uniformly on INIT, and we start in mode qin, consistent with the initial state (qin,INIT,PINIT) of M . Inductively, suppose η = (ρ, X, PX ) is a state of M with ρ ending in q. Suppose η p j ⇒ ηj is a transition in M such that ηj = (ρqj, Xj, PX j ).
Proposition 1. Suppose, we obtain a sample v ∼ PX . The probability of choosing guard gj whose target mode is qj in Algorithm 1 tends to pj as J → ∞.
Proof. According to Algorithm 1, the probability of picking guard gj for a trajectory starting at v ∈ X is defined as | Tj|/ m i=1 | Ti|, which, by the law of large numbers tends to as J tends to ∞.
Now if v is randomly sampled according to PX , then the Figure 2: Propagating a single value v ∈ X to v ∈ Xj when taking the transition q → qj at time t ∈ Tj(v).
probability of picking guard j can be expressed as the ex- which by (4) is equal to pj, the corresponding transition probability in the Markov chain.
Similarly, picking the transition time t from Tj will approximate sampling t ∼ P T j (v) , for sufficiently high J. Next, assume that we have picked q g j → qj as the transition to take. We sample t ∼ P T j (v) , and obtain v by numerical simulation via: Proof. Clearly it suffices to show that for a measurable subset Y ⊆ Xj, Pr(v ∈ Y ) = PX j (Y ). We start with Integrating now over all possible choices of v with respect to PX we have Consequently, the trajectory being generated will now be in mode qj with v ∈ Xj and v distributed according to PX j , compatible with the state ηj = (ρqj, Xj, PX j ) of M . Inductively it is hence guaranteed that each subsequent iteration of Algorithm 1 will produce values compatible with a path of M .
Whether the generated trajectory of length K +1 (and hence the corresponding path of M ) is a model of ψ can be determined using a standard BLTL model checker [10]. In fact this can be done on the fly which will often avoid generating the whole trajectory. Based on this, we can test whether Pr ≥1 (ψ) on M by testing the following alternative pair of hypotheses: H0 : Pr ≥1 (ψ) and H1 : Pr <1−δ (ψ), where 0 < δ < 1 is a parameter chosen by the user marking the interval [1 − δ, 1) as an indifference region in which accepting either hypothesis is fine. In our setting, whenever we encounter a sample (i.e. a randomly generated trajectory) that does not satisfy ψ, we can reject H0 and accept H1. Thus we only have to deal with false positives (when H0 is accepted while H1 happens to be true).
This leads to Algorithm 2 that repeatedly generates a random trajectory (using Algorithm 1), and decides after a finite number of tries between H0 and H1. For doing so we also fix a user-defined false positive rate α. The accuracy of Algorithm 2 is captured by the next result.

Theorem 3
1. The probability of choosing H1 when H0 is true (false negative) is 0.
Proof. As observed earlier the first part is obvious. To prove the second part, if H1 is true, then we know that Pr <1−δ (ψ). The probability of N sampled trajectories all satisfying ψ (and thus returning H0, a false positive) is at most (1 − δ) N . Therefore we have α ≤ (1 − δ) N , and this leads to N ≥ log α/ log(1 − δ).

QUANTITATIVE SPECIFICATIONS
To specify quantitative properties we fix a finite set of atomic propositions APqt of the form xi < c or xi > c where c is a rational constant. In what follows we shall assume for convenience that all the atomic propositions that we encounter are members of APqt. It will be straightforward to extend our arguments to include qualitative atomic propositions as well.
We partition R n into hypercubes according to the constants mentioned in the quantitative atomic propositions in APqt.
(Actually one could just focus on the members of APqt that appear in a given specification but we wish to deal with specifications later). Accordingly, define Ci to be the set of rational constants so that c ∈ Ci iff an atomic proposition of the form xi < c or xi > c appears in APqt. We next define for each dimension i the set of intervals

This leads to the set of hypercubes H given by
Clearly H is a partition of R n . The states of the Markov chain Mqt we wish to define as the approximation of H will be the states of M defined previously but now refined using H. More precisely we define Mqt = (Υqt, ⇒qt) inductively as follows: ∈ Υqt and it is the initial state of Mqt. Every other state in Υqt will be of the form (ρ, X, h, PX ) where ρ is a path in H, X is an open subset of R n of finite non-zero measure, h ∈ H and PX is a probability distribution over X. Furthermore X ⊆ h.
As the first step, we form the tuples (qin, where p h is defined as: We add all such triples ( , p h , (qin, X h , h, P X h )) to ⇒qt. For each v ∈ X h and each j we first define the set of time points Tj(v) ⊆ (0, 1) as Tj(v) = {t | Φq(t, v) ∈ gj}.
We next define Xj as Finally for each h ∈ H define X h j to be Next we denote by P T j (v) the uniform probability over Tj(v). Then the distribution P X h j is given by: Suppose Y is a measurable subset of X h j . Then We add the state (ρqj, X h j , h , P X h j ) to Υqt provided µ(X h j ) > 0. For each such state we add the triple to ⇒qt where p j,h is given by Using arguments very similar to the ones used in Section 3, we can easily establish that Mqt is a Markov chain whose state space has the structure of a finitely branching tree.

The two semantics
For interpreting BLT L formulas over Mqt it will be convenient to assume the following syntax in which negation is immediately followed by a quantitative atomic proposition: Clearly, every BLTL formula can be transformed into an equivalent formula that has the above syntax. This can be achieved by pushing negation inwards using equivalences such as ¬(ϕ1∨ϕ2) ≡ ¬ϕ1∧¬ϕ2, ¬G ≤k ϕ ≡ F ≤k ¬ϕ, ¬(ϕ1U ≤k ϕ2) ≡ G ≤k ¬ϕ2 ∨ (¬ϕ2U ≤k (¬ϕ1 ∧ ¬ϕ2)) etc.
The trajectory semantics is defined along previous lines but the atomic propositions are handled as follows. Let τ = (q0, v0) (q1, v1) . . . (q k , v k ) be a finite trajectory and 0 ≤ ≤ k. Then τ, |=H,qt xi < c iff v (i) < c. On the other hand τ, |=H,qt ¬ xi < c iff τ, |=H xi < c . The clauses for the other cases are defined in the obvious way. As before τ is a (trajectory) model of ψ iff τ ∈ T RJ K+1 and τ, 0 |=H,qt ψ.
For v ∈ R n let v |= A denote the fact that v(i) < c in case A = xi < c and v(i) > c in case A = xi > c . Next suppose (ρ, X, h, PX ) is a state of Mqt and A ∈ APqt. Then X ⊆ h by construction. Furthermore it is easy to check that v |= A for every v ∈ h or v |= A for every v ∈ h. Thus the semantics defined above will be consistent in the sense it will be the case that either π, |=qt A or π, |=qt ¬A but not both.

The correspondence result
We shall relate the behavior of H to that Mqt using the notion of robust trajectories. To start with, for v ∈ R n we let hc(v) be the hypercube h in H such that v ∈ h. Since H is a partition of R n we have that hc(v) exists and is unique. In what follows we let range over {0, 1, . . . , K}. We now define the equivalence relation ≈⊆ T RJ K+1 as follows: Let τ, τ ∈ T RJ K+1 with τ ( ) = (q , v ) and τ ( ) = (q , v ). Then τ ≈ τ iff q = q and hc(v ) = hc(v ) for each . We let [τ ] denote the ≈-equivalence class containing τ .
It is easy to verify that [τ ]( ) is a measurable set (but perhaps with measure 0) for each .
The trajectory τ ∈ T RJ K+1 is said to be robust iff µ([τ ]( )) > 0 for every . We will say that H robustly satisfies the specification ψ-and this is denoted by H |=R ψ iff τ, 0 |=H ψ for every robust trajectory τ in T RJ K+1 . It is now straightforward to show (along the lines of the proof of 4.1) show: Thus in terms of the sub-dynamics consisting of robust trajectories there is again a strong relationship between the behaviors of H and Mqt. It also turns out that in measuretheoretic terms the non-robust trajectories can be ignored. More precisely if one starts with the discrete topology over Q K+1 and the usual topology over R n K+1 one can easily define a natural measure space over the product topology Q K+1 × R n K+1 . In this space for every non-robust trajectory τ the representation of [τ ] will be measurable but with measure 0. In this sense the contributions made by the nonrobust trajectories to the dynamics of H are negligible.
The SMC procedure. The SMC procedure remains the same. The only difference is we must sample robust trajectories. One can do so by discarding the value states that hit the constants mentioned in the quantitative atomic propositions and restarting the simulation from the previous time step. That this leads with probability 1 to the generation of a robust trajectory is easy to establish.

CASE STUDIES
We have implemented our method using MATLAB. The source code is available at [25]. The experiments were carried out on a PC with a 3.4GHz Intel Core i7 processor with 8GB RAM. We first evaluated our method on a room heating model adapted from [12]. We also applied our method a model of the electrical dynamics of the cardiac cell [7]. The ∆ time step parameter for the room heating system and the cardiac cell model were set to 0.025, and 0.5, respectively. The parameters used for the statistical model checking were δ = 0.01 and α = 0.01, and hence the maximal sample size was set to N = log 0.01/ log 0.99 = 459. We used J = 10 as the number of intermediate time steps for choosing mode transitions (an empirical justification of this value is given in Section 7.3).

Room heating system
Our system consists of 3 rooms and 2 heaters [12]. The temperature in each room depends on the temperature of the adjacent rooms, on the outside temperature and on whether a heater is in the room. Each room may have at most one heater and a heater is switched ON if the temperature is below a certain threshold, and OFF if it is above another threshold. When the temperature in a room falls below a certain level, it may get a heater from one of the adjacent rooms, provided that the temperature in that room is significantly higher.
The system dynamics is given bẏ where xi is the temperature in room i, hi is 1 when a heater is switched ON in room i (0 otherwise) and u denotes the outside temperature. The system has a total of 12 modes corresponding to the locations of the 2 heaters and whether each heater is ON or OFF. Additional details about model parameters can be found in the Appendix.
Property R1 Initially, heaters are in room 1 and room 2, and no heater will be moved within 5 days.
The property was verified to be true when the outside temperature is ≥ 16 • C. However, for an outside temperature of 4 • C, the property evaluates to false. The next property involves quantitative atomic propositions.

Property R2
The temperature in all three rooms stabilizes within 1 day between 18 and 22 degrees, and stays in that range for 4 days.
The property was verified to be true for an outside temperature of ≥ 16 • C. However, the property was verified to be false with an outside temperature of 4 • C.

Cardiac cell model
Cardiac muscle cells enable the heart to beat and pump blood throughout the body. Heart rhythm depends on the organized opening and closing of gates-called ion channelson the cell membrane, which govern the electrical activity of cardiac cells. Disordered electric wave propagation in heart muscle can cause cardiac abnormalities such as tachycardia and fibrillation. The dynamics of the electrical activity of a single cardiac cell has been modeled as a hybrid automaton [7,16] shown in Figure 3. The model contains 4 state variables and 26 parameters. The parameters of the model are summarized in Table 2 in the Appendix. An action potential (AP) is a change in the cell's transmembrane potential u, as a response to an external stimulus (current) . The flow of total currents is controlled by a fast channel gate v and two slow gates w and s.
In mode q0, the "Resting mode", the cell is waiting for stimulation. We assume an external stimulus equal to 1 mV lasting for 1 millisecond. The stimulation causes u to increase which may trigger a mode transition to mode q1. In mode q1, gate v starts closing and the decay rate of u changes. The system will jump to mode q2 if u ≥ θw. In mode q2, gate w is also closing. When u ≥ θv, mode q3 can be reached, which means a successful "AP initiation". In mode q3, u reaches its peak due to the fast opening of a sodium channel. The cardiac muscle then contracts and u starts decreasing.
Property C1 It is known that the cardiac cell can lose its excitability, which will lead to disorders such as ventricular tachycardia and fibrillation. We formulated the property for responding to stimulus by leaving the resting mode: The property was verified to be true under the healthy condition. However, under a disease condition (for example τo1 = 0.004 or τo2 = 0.1 [22]) the property was verified to be false no matter what stimulation value of was used. Consequently, a region of such unexcitable cells blocks the impulse conduction and can lead to cardiac disorders such as fibrillation. This is consistent with experimental results reported in [26].
Property C2 After successfully generating an AP with transmembrane potential above 1.2, the cardiac cell should return to transmembrane potential below 0.006 and wait for the next stimulation. The corresponding formula (which uses quantitative propositions) is The above query was verified to be true under the healthy condition and transient stimulation. However, if we change the stimulation profile from transient to sustained, i.e. assuming lasts for 500 milliseconds, the property was verified to be false-the cell doesn't return to and settle at a low transmembrane potential resting state. In ventricular tissue the stimulus can be delivered from neighboring cells [7]. Thus, our results suggest that the transient activation of a single cardiac cell depends on the stimulation profile of its neighboring cells.

Performance
The overall run times for verifying the properties are summarized in Table 1. Our current implementation is a simple sequential one. The performance can be considerably improved through an optimized parallel implementation. In our experiments, we used J = 10 as the number of intermediate time steps for choosing mode transitions. We investigated whether this choice is sufficient for accurate simulation. We simulated 1000 independent realizations of the cardiac cell system with J = 10 and J = 100, and compared the distributions of the modes that the system is in at a series of discrete time points. The Kolmogorov-Smirnov statistical test did not reject the hypothesis that the two distributions are the same (at confidence level 95%). This indicates that using J = 10 is adequate.

CONCLUSION
We have presented an approximate probabilistic verification method for analyzing the dynamics of a hybrid system H in terms of a Markov chain M . For bounded time properties, we have shown a strong correspondence between H and M , namely H meets the BLTL specification ψ iff M meets ψ with probability 1. We have also extended this result to handle quantitative atomic propositions and shown a similar correspondence result for the sub-dynamics consisting of robust trajectories. Thus the intractable verification problem for H can be solved approximately using its Markov chain approximation. Accordingly, we have devised a statistical model checking procedure to verify that M almost certainly meets a BLTL specification and then applied this procedure to two examples to demonstrate the applicability of our approximation scheme. A parallel implementation of the trajectory sampling procedure -which should be easy to work out-will considerably improve the performance and scalability of our method.
As an extension, one could consider more sophisticated stochastic assumptions regarding the time points and value states at which the mode transitions take place. Further, we have assumed here a closed system setting in which there is no distinction made between controllable and uncontrollable transitions. However, especially for cyber physical systems a subset of the mode transitions will be controllable. Hence it will be fruitful to extend our technique to this setting and study controller synthesis problems. It will however involve MDPs instead of Markov chains as approximations and hence will require new machinery.
It will be easy to extend our results to bounded PCTL [17] specifications. It is also easy to define the Markov chain approximation of H for infinite trajectories. However, relating the behavior of H to its Markov chain approximation when the specification logic is LTL or PCTL appears to be a hard problem.   Figure 3: The hybrid automaton model for the cardiac cell system [16].
We define the notion of compatibility as before. Let π = η0η1 . . . η k be a path in Mqt with ηj = (q0q1 . . . qj−1, X h j j , hj, Pr X h j j ) for 0 < j ≤ k, and η0 = . Let τ = (q 1 , v1)(q 2 , v1) . . . (q k , v k ) be a trajectory. Then we say that π and τ are compatible iff k = k and for 1 ≤ j ≤ k, qj = q j and vj ∈ X h j j . As it will turn out, if τ and π are compatible then τ will be robust.
In what follows we shall assume that our BLTL specifications involve only quantitative atomic propositions in APqt and the formulas obey the syntax in which negation is immediately followed by an atomic proposition. Further the semantic notions |=H and |=M qt (abbreviated as |=qt) are defined in the expected way.
2. Suppose π is a path in Mqt starting from . Then there exists a robust trajectory τ in T RJ such that π and τ are compatible.
3. Suppose τ is a robust trajectory in T RJ. Then there exists a path π in Mqt starting from such that τ and π are compatible.
Proof. 1. From the definitions it follows that if A ∈ APqt and h ∈ H then v |= A for every v ∈ h or v |= ¬A for every v ∈ h but not both. Since vj ∈ hj we then have τ, j |=H A iff π, j |=qt A and τ, j |=H ¬A iff π, j |=qt ¬A for every atomic proposition. The remaining cases now follow easily by structural induction on ψ.
2. Let π = η0η1 . . . η k in Mqt with η0 = and ηj = (q0q1 . . . qj−1, X h j j , hj, Pr X h j j ) for 0 < j ≤ k. For notational convenience we will write Xj instead of X h j j . Since µ(X k ) > 0 we can fix v k ∈ X k . Further h k being a product of open intervals in R with X k ⊆ h k , we can find an open set O k of non-zero measure such that From the construction of Mqt it follows there exists . Continuing this way we can find vj, Oj, Yj for 1 ≤ j ≤ k (with Y k = X k ) such that τ = (q1, v1)(q2, v2) . . . (q k , v k ) is a trajectory and vj ∈ Oj ⊆ Yj ⊆ hj for 1 ≤ j ≤ k. From the construction of Mqt it follows that Yj = [τ ](j) for 1 ≤ j ≤ k. From Lemma 3 it follows that π and τ are compatible. It is also clear due to Lemma 3 that τ is robust.
3. Suppose τ = (q1, v1)(q2, v1) . . . (q k , v k ) ∈ T RJ is robust. Then by Lemma 3 there exist open sets Oj of non zero measure and hj ∈ H such that vj ∈ Oj ⊆ [τ ](j) ⊆ hj for 1 ≤ j ≤ k. Let τ (j) denote the j-length prefix of τ for 1 ≤ j ≤ k. We now define Xj = [τ (j) ](j) for 1 ≤ j ≤ k. Then using the construction of Mqt it is easy to show that there exists distributions P rj over Xj such that π = η1η2 . . . η k is a path in Mqt with ηj = (qj, Xj, hj, P rj) for 1 ≤ j ≤ k and that π is compatible with τ .
We can now prove Theorem 4.
Proof. Suppose H |=R ψ. Then there exists τ ∈ T RJ such that τ is robust and τ, 0 |=H ψ. By Lemma 4, there exists a path π in Mqt which is compatible with τ . Hence again by Lemma 4 we then have π / ∈ modelsM qt (ψ) which leads to P r<1(ψ). Next suppose that P r<1(ψ). Then there exists a path π in Mqt such that π, 1 |=M qt ψ. By Lemma 4, there exists a robust trajectory τ which is compatible with π and τ, 0 |=H ψ. This implies H |=R ψ.
Finally, we wish to show that the number of non-robust trajectories are negligible compared with the robust ones. Hence they do not contribute much towards the dynamics of H. For that we need the following lemma.
Thus by the above lemma, if a trajectory τ ∈ T RJ K+1 is not robust then there exists a j ∈ {0, 1, . . . , K} such that µ(Yj) = 0. This implies that in the product topology of Q K+1 × R K+1 , [τ ] has measure 0. Thus, the contribution made by the non-robust trajectories to the dynamics of H is negligible.
Trajectory simulation for quantitative specifications Algorithm 3 gives the procedure for simulating robust trajectories for the verification of quatitative BLTL specifications. By Lemma 3 , a trajectory is robust iff it does not hit any of the constants mentioned in the atomic propositions. The procedure is the same as Algorithm 1 before, except that whenever a value state v k at any time step k hits a constant mentioned in any of the atomic propositions, we discard v k and start the simulation again from the value state of the previous time step. Generate time points T := {t1, . . . , t J } uniformly in (0, 1).

7:
Let Tj := {t ∈ T : v ∈ gj } be the time points where gj is enabled.

8:
Pick g randomly according to probabilities pj :=

9:
Pick t uniformly at random from T .

11:
until v (i) / ∈ Ci for any i 12: Set q k := q , v k := v , and extend τ := (q0, v0) . . . (q k , v k ). 13: end for 14: return τ To see that the algorithm terminates with probability 1, note that if v0 ∈ h and h(i) = {c} for some c ∈ Ci then µ(h) = 0. Thus Step 1 repeats with probability 0. As a result with probability 1 it will be repeated only a finite number of times. Similarly the repeat loop of Step 4-11 will terminate with probability 1.

Room heating system
The heaters are switched ON or OFF, and moved from one room to another based on temperature thresholds. For each room we define two thresholds oni and of fi. The heater in room i is ON if xi ≤ oni and OFF if xi ≥ of fi. A heater is moved from room j to an adjacent room i if the following holds: room i has no heater, room j has a heater, temperature xi ≤ geti, and the difference xj − xi ≥ difi.
We set the initial temperature to be uniformly distributed in the interval [20, 20.5] in all 3 rooms. We use the model parameters as they appear in [12], which are as follows. Further, of fi = 21, oni = 20, geti = 18, and dif fi = 1.