map : (a ↔ b) → ([a] ↔ [b]) λg.μf. ⎛ ⎝ [] ↔ [] h : t ↔ let x = g h in let y = f t in x : y ⎞ ⎠ , mapAccu : (a ⊗ b ↔ a ⊗ c) → (a ⊗ [b] ↔ a ⊗ [c]) λg.μf. ⎛ ⎜⎜⎝ x, [] ↔x, [] x, (h : t)↔ let y, h- = g x, h in let z, t- = f y, t in z, (h- : t - ) ⎞ ⎟⎟⎠ .

The three examples cnot\*, map and mapAccu uses fixpoints which are clearly terminating in any finite context. Indeed, the functions are structurally recursive. A formal definition of this notion for the reversible language is as follows.

**Definition 1.** Define a *structurally recursive type* as a type of the form [a] ⊗ b1⊗...⊗bn. Let ω = {vi ↔ ei | i ∈ I} be an iso such that f : a ↔ b ω ω : a ↔ c where a is a structurally recursive type. We say that μf.ω is *structurally recursive* provided that for each i ∈ I, the value vi is either of the form [], p1,...pn or of the form h : t, p1,...pn. In the former case, ei does not contain f as a free variable. In the latter case, ei is of the form C[ft, p 1,...,p n] where C is a context of the form C[−] ::= [−] | let p = C[−] in t | let p = t in C[−].

This definition will be critical for quantum loops in the next section.

### **5 From Reversible Isos to Quantum Control**

In the language presented so far, an iso ω : a ↔ b describes a bijection between the set Ba of closed values of type a and the set Bb of closed values of type b. If one regards Ba and Bb as the basis elements of some vector space a and b, the iso ω becomes a 0/1 matrix.

As an example, consider an iso ω defined using three clauses of the form { | v1 ↔ v 1 | v2 ↔ v 2 | v3 ↔ v 3 }. From the exhaustivity and non-overlapping conditions derives the fact that the space a can be split into the direct sum of the three subspaces av*i* (i = 1, 2, 3) generated by vi. Similarly, b is split into the direct sum of the subspaces bv- *i* generated by v i. One can therefore represent ω as the matrix ω in Fig. 4: The "1" in each column vi indicates to which subspace b an element of av*i* is sent to.

v*j* In Sect. 2.2 we discussed the fact that vi⊥vj when i = j. This notation hints at the fact that a and b could be seen as Hilbert spaces and the mapping ω as a unitary map from a to b. The purpose of this section is to extend and formalize precisely the correspondence between isos and unitary maps.

The definition of clauses is extended following this idea of seeing isos as unitaries, and not only bijections on basis elements of the

$$\left\{ \begin{array}{l} \mid \ v\_{1} \leftrightarrow a\_{11}v'\_{\clubsuit} + a\_{21}v'\_{\clubsuit} + a\_{31}v'\_{\clubsuit} \\ \mid \ v\_{2} \leftrightarrow a\_{12}v'\_{\clubsuit} + a\_{22}v'\_{\clubsuit} + a\_{23}v'\_{\clubsuit} \\ \mid \ v\_{3} \leftrightarrow a\_{31}v'\_{1} + a\_{32}v'\_{2} + a\_{33}v'\_{3} \end{array} \right\}$$

input space. We therefore essentially propose to generalize the clauses to complex, linear combinations of values on the right-hand-side, such as shown on the left, with the side conditions on

that the matrix of Fig. 5 is unitary. We define in Sect. 5.1 how this extends to second-order.

#### **5.1 Extending the Language to Linear Combinations of Terms**

The quantum unitary language extends the reversible language from the previous section by closing extended values and terms under complex, finite linear combinations. For example, if v1 and v2 are values and α and β are complex numbers, α · v1 + β · v2 is now an extended value.

Several approaches exist for performing such an extension. One can update the reduction strategy to be able to reduce these sums and scalar multiplications to normal forms [12,18], or one can instead consider terms modulo the usual algebraic equalities [13,18]: this is the strategy we follow for this paper.

When extending a language to linear combination of terms in a naive way, this added structure might generate inconsistencies in the presence of unconstrained fixpoints [12,13,18]. The weak condition on termination we imposed on fixpoints in the classical language was enough to guarantee reversibility. With the presence of linear combinations, we want the much stronger guarantee of unitarity. For this reason, we instead impose fixpoints to be *structurally recursive*.

The quantum unitary language is defined by allowing sums of terms and values and multiplications by complex numbers: if t and t are terms, so is α · t + t . Terms and values are taken modulo the equational theory of modules. We furthermore consider the value and term constructs −, −, let p = − in −, injl (−), injr (−) distributive over sum and scalar multiplication. We do *not* however take iso-constructions as distributive over sum and scalar multiplication: { | v1 ↔ αv2 + βv3 } is *not* the same thing as α { | v1 ↔ v2 } + β { | v1 ↔ v3 }. This is in the spirit of Lineal [11,12].

The typing rules for terms and extended values are updated as follows. We only allow linear combinations of terms and values of the same type and of the same free variables. Fixpoints are now required to be *structurally recursive*, as introduced in Definition 1. Finally, an iso is now not only performing an "identity" as in Fig. 4 but a true unitary operation:

$$\frac{\begin{array}{l} \Delta\_{1}; \Psi \vdash\_{v} v\_{1}: a & \dots & \Delta\_{n}; \Psi \vdash\_{v} v\_{n}: a \\ \Delta\_{1}; \Psi \vdash\_{v} e\_{1}: b & \dots & \Delta\_{n}; \Psi \vdash\_{v} e\_{n}: b \\ \text{OD}\_{a} \{v\_{1}, \dots, v\_{n}\} & \text{OD}\_{b}^{ext} \{e\_{1}, \dots, e\_{n}\} & \begin{pmatrix} a\_{11} & \dots & a\_{1n} \\ \vdots & & \vdots \\ a\_{n1} & \dots & a\_{nn} \end{pmatrix} \text{ is unitary} \end{array}}{\begin{array}{l} \text{OD}\_{a} \{v\_{1}, \dots, v\_{n}\} \qquad \begin{pmatrix} v\_{1} \dashv \cdots \dashv a\_{1n} \\ \vdots \\ a\_{n1} \dashv \cdots \dashv a\_{nn} \end{pmatrix} \text{ is unitary} \end{array}}$$

The reduction relation is updated in a way that it remains deterministic in this extended setting. It is split into two parts: the reduction of pure terms, i.e. non-extended terms or values, and linear combinations thereof. Pure terms and values reduce using the reduction rules found in Table 3. We do not extend applicative contexts to linear combinations. For linear combinations of pure terms, we simply ask that *all* pure terms that are not normal forms in the combination are reduced. This makes the extended reduction relation deterministic.

*Example 1.* This allows one to define an iso behaving as the Hadamard gate, or a slightly more complex iso conditionally applying another iso, whose behavior as a matrix is shown in Fig. 6.

$$\begin{array}{llll} \mathsf{Had}: \mathbb{B} \leftrightarrow \mathbb{B} & \mathsf{Gate}: \mathbb{B} \otimes \mathbb{B} \leftrightarrow \mathbb{B} \otimes \mathbb{B} \\ \begin{pmatrix} \mathtt{tt} \leftrightarrow \frac{1}{\sqrt{2}} \mathtt{tt} + \frac{1}{\sqrt{2}} \mathtt{ff} \\ \mathtt{ff} \leftrightarrow \frac{1}{\sqrt{2}} \mathtt{tt} - \frac{1}{\sqrt{2}} \mathtt{ff} \end{pmatrix}, \quad \begin{pmatrix} \langle \mathtt{tt}, x \rangle \leftrightarrow \ \mathtt{let} \ x \ \mathtt{in} \ \frac{1}{\sqrt{2}} \langle \mathtt{tt}, y \rangle + \frac{1}{\sqrt{2}} \langle \mathtt{ff}, y \rangle \\ \langle \mathtt{ff}, x \rangle \leftrightarrow \ \mathtt{let} \ y = \mathtt{Id} \ x \ \mathtt{in} \ \frac{1}{\sqrt{2}} \langle \mathtt{tt}, y \rangle - \frac{1}{\sqrt{2}} \langle \mathtt{ff}, y \rangle \end{pmatrix}. \end{array}$$

With this extension to linear combinations of terms, one can characterize normal forms as follows.

**Lemma 1 (Structure of the Normal Forms).** *Let* ω *be such that* ω ω : a ↔ b*. For all closed values* v *of type* a*, the term* ω v *rewrites to a normal form* N i=1 αi · wi *where* N < ∞*, each* wi *is a closed value of type* b *and* i |αi| = 1*.*

*Proof.* The fact that ω v converges to a normal form is a corollary of the fact that we impose structural recursion on fixpoints. The property of the structure of the normal form is then proven by induction on the maximal number of steps it takes to reach it. It uses the restriction on the introduction of sums in the typing rule for clauses in isos and the determinism of the reduction.

In the classical setting, isos describe bijections between sets of closed values: it was proven by considering the behavior of an iso against its inverse. In the presence of linear combinations of terms, we claim that isos describe more than bijections: they describe unitary maps. In the next section, we discuss how types can be understood as Hilbert spaces (Sect. 5.2) and isos as unitary maps (Sects. 5.3 and 5.4).

#### **5.2 Modeling Types as Hilbert Spaces**

By allowing complex linear combinations of terms, closed normal forms of finite types such as **B** or **B** ⊗ **B** can be regarded as complex vector spaces with basis consisting of closed values. For example, **B** is associated with -**B** = {α · tt + β · ff | α, β ∈ **C**} ≡ **C**2. We can consider this space as a complex Hilbert space where the scalar product is defined on basis elements in the obvious way: v|v = 1 and v|w = 0 if v = w. The map Had of Example 1 is then effectively a unitary map on the space -**B**.

The problem comes from lists: the type [**1**] is inhabited by an infinite number of closed values: [], [()], [(),()], [(),(),()],. . . To account for this case, we need to consider infinitely dimensional complex Hilbert spaces. In general, a complex Hilbert space [19] is a complex vector space endowed with a scalar product that is complete with respect the distance induced by the scalar product. The completeness requirement implies for example that the infinite linear combination []+ 1 2 ·[()]+ 1 4 [(),()]+ 1 8 [(),(),()]+··· needs to be an element of -[**B**]. To account for these limit elements, we propose to use the standard [19] Hilbert space 2 of infinite sequences.

**Definition 2.** Let a be a value type. As before, we write Ba for the set of closed values of type a, that is, Ba = {v | v v : a}. The *span of a* is defined as the Hilbert space a = 2(Ba) consisting of sequences (φv)v∈B*a* of complex numbers indexed by Ba such that v∈B*a* |φv| 2 < ∞. The scalar product on this space is defined as (φv)v∈B*a* |(ψv)v∈B*a* = v∈B*a* φvψv.

We shall use the following conventions. A closed value v of a is identified with the sequence (δv,v- )v-∈B*a* where δv,v = 1 and δv,v- = 0 if v = v . An element (φv)v∈B*a* of a is also written as the infinite, formal sum v∈B*a* φv · v.

#### **5.3 Modeling Isos as Bounded Linear Maps**

We can now define what is the linear map associated to an iso.

**Definition 3.** For each closed iso ω ω : a ↔ b we define ω as the linear map from a to b sending the closed value v : a to the normal form of ω v : b under the rewrite system.

In general, the fact that ω is well-defined is not trivial. If it is formally stated in Theorem 3, we can first try to understand what could go wrong. The problem comes from the fact that the space a is not finite in general. Consider the iso map Had : [**B**] ↔ [**B**]. Any closed value v : [**B**] is a list and the term (map Had) v rewrites to a normal form consisting of a linear combination of lists. Denote the linear combination associated to v with Lv. An element of -[**B**] is a sequence φ = (φv)v∈B[**B**] . From Definition 3, the map ω sends the element φ ∈ -[**B**] to v∈B[**B**] φv ·Lv. This is an infinite sum of sums of complex numbers: we need to make sure that it is well-defined: this is the purpose of the next result. Because of the constraints on the language, we can even show that it is a *bounded* linear map.

In the case of the map map Had, we can understand why it works as follows. The space -[**B**] can be decomposed as the direct sum ∞ i=0 Ei, where Ei is generated with all the lists in **B** of size i. The map map Had is acting locally on each finitely-dimensional subspace Ei. It is therefore well-defined. Because of the unitarity constraint on the linear combinations appearing in Had, the operation performed by map Had sends elements of norm 1 to elements of norm 1. This idea can be formalized and yield the following theorem.

**Theorem 3.** *For each closed iso* ω ω : a ↔ b *the linear map* ω : a → b *is well-defined and bounded.*

#### **5.4 Modeling Isos as Unitary Maps**

In this section, we show that not only closed isos can be modeled as bounded linear maps, but that these linear maps are in fact unitary maps. The problem comes from fixpoints. We first consider the case of isos written without fixpoints, and then the case with fixpoints.

*Without recursion.* The case without recursion is relatively easy to treat, as the linear map modeling the iso can be compositionally constructed out of elementary unitary maps.

**Theorem 4.** *Given a closed iso* ω ω : a ↔ b *defined without the use of recursion, the linear map* π : a → b *is unitary.*

The proof of the theorem relies on the fact that to each closed iso ω ω : a ↔ b one can associate an operationally equivalent iso ω ω : a ↔ b that does not use iso-variables nor lambda-abstractions. We can define a notion of *depth* of an iso as the number of nested isos. The proof is done by induction on this depth of the iso ω: it is possible to construct a unitary map for ω using the unitary maps for each ωij as elementary building blocks.

As an illustration, the semantics of Gate of Example 1 is given in Fig. 6.

*Isos with structural recursion.* When considering fixpoints, we cannot rely anymore on this finite compositional construction: the space a cannot anymore be regarded as a *finite* sum of subspaces described by each clause.

We therefore need to rely on the formal definition of unitary maps in general, infinite Hilbert spaces. On top of being bounded linear, a map ω : a → b is unitary if (1) it preserves the scalar product: ω(e)|ω(f) = e|f for all e and f in a and (2) it is surjective.

**Theorem 5.** *Given a closed iso* ω ω : a ↔ b *that can use structural recursion, the linear map* π : a → b *is unitary.*

The proof uses the idea highlighted in Sect. 5.4: for a structurally recursive iso of type [a] ⊗ b ↔ c, the Hilbert space -[a] ⊗ b can be split into a canonical decomposition E0 ⊕E1 ⊕E2 ⊕··· , where Ei contains only the values of the form [x1 ...xi], y, containing the lists of size i. On each Ei, the iso is equivalent to an iso without structural recursion.

### **6 Conclusion**

In this paper, we proposed a reversible language amenable to quantum superpositions of values. The language features a weak form of higher-order that is nonetheless expressible enough to get interesting maps such as generalized Toffoli operators. We sketched how this language effectively encodes bijections in the classical case and unitary operations in the quantum case. It would be interesting to see how this relates to join inverse categories [14,15].

In the vectorial extension of the language we have the same control as in the classical, reversible language. Tests are captured by clauses, and naturally yield quantum tests: this is similar to what can be found in QML [5,6], yet more general since the QML approach is restricted to if-then-else constructs. The novel aspect of quantum control that we are able to capture here is a notion of *quantum loops*. These loops were believed to be hard, if not impossible. What makes it work in our approach is the fact that we are firmly within a closed quantum system, without measurements. This makes it possible to only consider unitary maps and frees us from the L¨ower order on positive matrices [6]. As we restrict fixpoints to structural recursion, valid isos are regular enough to capture unitarity. Ying [7] also proposes a framework for quantum while-loops that is similar in spirit to our approach at the level of denotations: in his approach the control part of the loops is modeled using an external systems of "coins" which, in our case, correspond to conventional lists. Reducing the manipulation of this external coin system to iteration on lists allowed us to give a simple operational semantics for the language.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# Quantitative Models

# **The Complexity of Graph-Based Reductions for Reachability in Markov Decision Processes**

St´ephane Le Roux1(B) and Guillermo A. P´erez2

1 Department of Mathematics, Technische Universit¨at Darmstadt, Darmstadt, Germany<EMAIL_ADDRESS>2 Departement d'Informatique, Universit´e libre de Bruxelles, Brussels, Belgium

<EMAIL_ADDRESS>
**Abstract.** We study the never-worse relation (NWR) for Markov decision processes with an infinite-horizon reachability objective. A state *q* is never worse than a state *p* if the maximal probability of reaching the target set of states from *p* is at most the same value from *q*, regardless of the probabilities labelling the transitions. Extremal-probability states, end components, and essential states are all special cases of the equivalence relation induced by the NWR. Using the NWR, states in the same equivalence class can be collapsed. Then, actions leading to sub-optimal states can be removed. We show that the natural decision problem associated to computing the NWR is coNP-complete. Finally, we extend a previously known incomplete polynomial-time iterative algorithm to under-approximate the NWR.

### **1 Introduction**

Markov decision processes (MDPs) are a useful model for decision-making in the presence of a stochastic environment. They are used in several fields, including robotics, automated control, economics, manufacturing and in particular planning [20], model-based reinforcement learning [22], and formal verification [1]. We elaborate on the use of MDPs and the need for graph-based reductions thereof in verification and reinforcement learning applications below.

Several verification problems for MDPs reduce to reachability [1,5]. For instance, MDPs can be model checked against linear-time objectives (expressed in, say, LTL) by constructing an omega-automaton recognizing the set of runs that satisfy the objective and considering the product of the automaton with the original MDP [6]. In this product MDP, accepting end components—a generalization of strongly connected components—are identified and selected as target components. The question of maximizing the probability that the MDP behaviours satisfy the linear-time objective is thus reduced to maximizing the probability of reaching the target components.

The maximal reachability probability is computable in polynomial time by reduction to linear programming [1,6]. In practice, however, most model checkers c The Author(s) 2018

use value iteration to compute this value [9,17]. The worst-case time complexity of value iteration is pseudo-polynomial. Hence, when implementing model checkers it is usual for a graph-based pre-processing step to remove as many unnecessary states and transitions as possible while preserving the maximal reachability probability. Well-known reductions include the identification of extremal-probability states and maximal end components [1,5]. The intended outcome of this pre-processing step is a reduced amount of transition probability values that need to be considered when computing the number of iterations required by value iteration.

The main idea behind MDP reduction heuristics is to identify subsets of states from which the maximal probability of reaching the target set of states is the same. Such states are in fact redundant and can be "collapsed". Figure 1 depicts an MDP with actions and probabilities omitted for clarity. From p and q there are strategies to ensure that s is reached with probability 1. The same holds for t. For instance, from p, to get to t almost surely, one plays to go to the distribution directly below q; from q, to the distribution above q. Since from the state p, there is no strategy to ensure that q is reached with probability 1, p and q do not form an *end component*. In fact, to the best of our knowledge, no known MDP reduction heuristic captures this example (i.e., recognizes that p and q have the same maximal reachability probability for all possible values of the transition probabilities).

**Fig. 1.** An MDP with states depicted as circles and distributions as squares. The maximal reachability probability values from *p* and *q* are the same since, from both, one can enforce to reach *s* with probability 1, or *t* with probability 1, using different strategies.

In reinforcement learning the actual probabilities labelling the transitions of an MDP are not assumed to be known in advance. Thus, they have to be estimated by experimenting with different actions in different states and collecting statistics about the observed outcomes [14]. In order for the statistics to be good approximations, the number of experiments has to be high enough. In particular, when the approximations are required to be *probably approximately correct* [23], the necessary and sufficient number of experiments is pseudo-polynomial [13]. Furthermore, the expected number of steps before reaching a particular state even once may already be exponential (even if all the probabilities are fixed). The fact that an excessive amount of experiments is required is a known drawback of reinforcement learning [15,19].

A natural and key question to ask in this context is whether the maximal reachability probability does indeed depend on the actual value of the probability labelling a particular transition of the MDP. If this is not the case, then it need not be learnt. One natural way to remove transition probabilities which do not affect the maximal reachability value is to apply model checking MDP reduction techniques.

*Contributions and Structure of the Paper.* We view the directed graph underlying an MDP as a directed bipartite graph. Vertices in this graph are controlled by players *Protagonist* and *Nature*. Nature is only allowed to choose full-support probability distributions for each one of her vertices, thus instantiating an MDP from the graph; Protagonist has strategies just as he would in an MDP. Hence, we consider infinite families of MDPs with the same support. In the game played between Protagonist and Nature, and for vertices u and v, we are interested in knowing whether the maximal reachability probability from u is never (in any of the MDPs with the game as its underlying directed graph) worse than the same value from v.

In Sect. 2 we give the required definitions. We formalize the *never-worse relation* in Sect. 3. We also show that we can "collapse" sets of equivalent vertices with respect to the NWR (Theorem 1) and remove sub-optimal edges according to the NWR (Theorem 2). Finally, we also argue that the NWR generalizes most known heuristics to reduce MDP size before applying linear programming or value iteration. Then, in Sect. 4 we give a graph-based characterization of the relation (Theorem 3), which in turn gives us a coNP upper bound on its complexity. A matching lower bound is presented in Sect. 5 (Theorem 4). To conclude, we recall and extend an iterative algorithm to efficiently (in polynomial time) under-approximate the never-worse relation from [2].

*Previous and Related Work.* Reductions for MDP model checking were considered in [5,7]. From the reductions studied in both papers, extremal-probability states, essential states, and end components are computable using only graphbased algorithms. In [3], learning-based techniques are proposed to obtain approximations of the maximal reachability probability in MDPs. Their algorithms, however, do rely on the actual probability values of the MDP.

This work is also related to the widely studied model of interval MDPs, where the transition probabilities are given as intervals meant to model the uncertainty of the numerical values. Numberless MDPs [11] are a particular case of the latter in which values are only known to be zero or non-zero. In the context of numberless MDPs, a special case of the question we study can be simply rephrased as the comparison of the maximal reachability values of two given states.

In [2] a preliminary version of the iterative algorithm we give in Sect. 6 was described, implemented, and shown to be efficient in practice. Proposition 1 was first stated therein. In contrast with [2], we focus chiefly on characterizing the never-worse relation and determining its computational complexity.

### **2 Preliminaries**

We use set-theoretic notation to indicate whether a letter b ∈ Σ *occurs* in a word α = a0 ...ak ∈ Σ∗, i.e. b ∈ α if and only if b = ai for some 0 ≤ i ≤ k.

Consider a directed graph G = (V,E) and a vertex u ∈ V . We write uE for the set of *successors* of u. That is to say, uE := {v ∈ V | (u, v) ∈ E}. We say that a path π = u0 ...uk ∈ V ∗ in G *visits* a vertex v if v ∈ π. We also say that π is a v–T path, for T ⊆ V , if u0 = v and uk ∈ T.

#### **2.1 Stochastic Models**

Let S be a finite set. We denote by D(S) the set of all *(rational) probabilistic distributions* on S, i.e. the set of all functions f : S → Q≥0 such that - s∈S f(s) = 1. A probabilistic distribution f ∈ D(S) has *full support* if f(s) > 0 for all s ∈ S.

**Definition 1 (Markov chains).** *A* Markov chain C *is a tuple* (Q, δ) *where* Q *is a finite set of states and* δ *is a probabilistic transition function* δ : Q → D(Q)*.*

A *run* of a Markov chain is a finite non-empty word = p0 ...pn over Q. We say *reaches* q if q = pi for some 0 ≤ i ≤ n. The *probability of the run* is 0≤i<n δ(pi, pi+1).

Let T ⊆ Q be a set of states. The *probability of (eventually) reaching* T in C from q0, which will be denoted by Pq0 C [♦T], is the measure of the runs of C that start at q0 and reach T. For convenience, let us first define the *probability of staying in states from* S ⊆ Q *until* T *is reached* 1, written Pq0 C [S U T], as 1 if q0 ∈ T and otherwise

$$\sum \left\{ \prod\_{0 \le i < n} \delta(q\_i, q\_{i+1}) \left| q\_0 \dots q\_n \in (S \nmid T)^\* T \text{ for } n \ge 1 \right. \right\}.$$

We then define Pq0 C [♦T] := Pq0 C [Q U T].

When all runs from q0 to T reach some set U ⊆ Q before, the probability of reaching T can be decomposed into a finite sum as in the lemma below.

**Lemma 1.** *Consider a Markov chain* C = (Q, δ)*, sets of states* U, T ⊆ Q*, and a state* q0 ∈ Q \ U*. If* Pq0 C [(Q \ U) U T]=0*, then*

$$\mathbb{P}\_{\mathcal{C}}^{q\_0}[\lozenge T] = \sum\_{u \in U} \mathbb{P}\_{\mathcal{C}}^{q\_0}[(Q \nmid U) \bullet u] \mathbb{P}\_{\mathcal{C}}^u[\lozenge T].$$

**Definition 2 (Markov decision processes).** *A* (finite, discrete-time) Markov decision process M*, MDP for short, is a tuple* (Q, A, δ, T) *where* Q *is a finite set of states,* A *a finite set of actions,* δ : Q × A → D(Q) *a* probabilistic transition function*, and* T ⊆ Q *a set of* target *states.*

For convenience, we write δ(q|p, a) instead of δ(p, a)(q).

1 *S* U *T* should be read as "*S until T*" and not understood as a set union.

**Definition 3 (Strategies).** *A* (memoryless deterministic) strategy σ *in an MDP* M = (Q, A, δ, T) *is a function* σ : Q → A*.*

Note that we have deliberately defined only memoryless deterministic strategies. This is at no loss of generality since, in this work, we focus on maximizing the probability of reaching a set of states. It is known that for this type of objective, memoryless deterministic strategies suffice [18].

*From MDPs to Chains.* An MDP M = (Q, A, δ, T) and a strategy σ induce the Markov chain Mσ = (Q, μ) where μ(q) = δ(q, σ(q)) for all q ∈ Q.

**Fig. 2.** On the left we have an MDP with actions {*a, b*}. On the right we have the Markov chain induced by the left MDP and the strategy {*p* -→ *a, q* -→ *b*}.

*Example 1.* Figure 2 depicts an MDP on the left. Circles represent states; doublecircles, target states; and squares, distributions. The labels on arrows from states to distributions are actions; those on arrows from distributions to states, probabilities.

Consider the strategy σ that plays from p the action a and from q the action b, i.e. σ(p) = a and σ(q) = b. The Markov chain on the right is the chain induced by σ and the MDP on the left. Note that we no longer have action labels.

The probability of reaching a target state from q under σ is easily seen to be 3/4. In other words, if we write M for the MDP and T for the set of target states then Pq Mσ [♦T] = 3 4 .

#### **2.2 Reachability Games Against Nature**

We will speak about families of MDPs whose probabilistic transition functions have the same support. To do so, we abstract away the probabilities and focus on a game played on a graph. That is, given an MDP M = (Q, A, δ, T) we consider its *underlying directed graph* GM = (V,E) where V := Q∪(Q×A) and E := {(q,q, a ) ∈ Q × (Q × A)}∪{(p, a , q) | δ(q|p, a) > 0}. In GM, *Nature* controls the vertices Q × A. We formalize the game and the *arena* it is played on below.

**Definition 4 (Target arena).** *A* target arena A *is a tuple* (V,VP ,E,T) *such that* (VP , VN := V \VP , E) *is a bipartite directed graph,* T ⊆ VP *is a set of* target *vertices, and* uE = ∅ *for all* u ∈ VN *.*

Informally, there are two agents in a target arena: *Nature*, who controls the vertices in VN , and *Protagonist*, who controls the vertices in VP .

*From Arenas to MDPs.* A target arena A = (V,VP ,E,T) together with a family of probability distributions μ = (μu ∈ D(uE))u∈VN induce an MDP. Formally, let Aμ be the MDP (Q, A, δ, T) where Q = VP {⊥}, A = VN , δ(q|p, a) is μa(q) if (p, a),(a, q) ∈ E and 0 otherwise, for all p ∈ VP ∪ {⊥} and a ∈ A we have δ(⊥|p, a) = 1 if (p, a) ∈ E.

*The Value of a Vertex.* Consider a target arena A = (V,VP ,E,T) and a vertex v ∈ VP . We define its *(maximal reachability probability) value* with respect to a family of full-support probability distributions μ as Valμ(v):= maxσ Pv Aσ μ [♦T]. For u ∈ VN we set Valμ(u) := -{μu(v)Valμ(v) | v ∈ uE}.

### **3 The Never-Worse Relation**

We are now in a position to define the relation that we study in this work. Let us fix a target arena A = (V,VP ,E,T).

**Definition 5 (The never-worse relation (NWR)).** *A subset* W ⊆ V *of vertices is* never worse *than a vertex* v ∈ V *, written* v W*, if and only if*

$$\forall \mu = (\mu\_u \in \mathbb{D}(uE))\_{u \in V\_N}, \exists w \in W : \text{Val}^\mu(v) \le \text{Val}^\mu(w).$$

*where all the* μu *have full support. We write* v ∼ w *if* v {w} *and* w {v}*.*

It should be clear from the definition that ∼ is an equivalence relation. For u ∈ V let us denote by ˜u the set of vertices that are ∼-equivalent and belong to the same owner, i.e. ˜u is {v ∈ VP | v ∼ u} if u ∈ VP and {v ∈ VN | v ∼ u} otherwise.

**Fig. 3.** Two target arenas with *T* = {*fin*} are shown. Round vertices are elements from *VP* ; square vertices, from *VN* . In the left target arena we have that *p* {*q*} and *q* {*p*} since any path from either vertex visits *t* before *T*—see Lemma 1. In the right target arena we have that *t* {*p*}—see Proposition 1.

*Example 2.* Consider the left target arena depicted in Fig. 3. Using Lemma 1, it is easy to show that neither p nor q is ever worse than the other since t is visited before *fin* by all paths starting from p or q.

The literature contains various heuristics which consist in computing sets of states and "collapsing" them to reduce the size of the MDP without affecting the maximal reachability probability of the remaining states. We now show that we can collapse equivalence classes and, further, remove sub-optimal distributions using the NWR.

#### **3.1 The Usefulness of the NWR**

We will now formalize the idea of "collapsing" equivalent vertices with respect to the NWR. For convenience, we will also remove self-loops while doing so.

Consider a target arena A = (V,VP ,E,T). We denote by A/∼ its ∼*-quotient*. That is, A/∼ is the target arena (S, SP , R, U) where SP = {v˜ | ∃v ∈ VP }, S = {v˜ | ∃v ∈ VN } ∪ SP , U = {t ˜ | ∃t ∈ T}, and

$$\begin{aligned} R &= \{ (\tilde{u}, \tilde{v}) \mid \exists (u, v) \in (V\_P \times V\_N) \cap E : vE \mid \tilde{u} \neq \mathcal{Q} \}, \\ &\cup \{ (\tilde{u}, \tilde{v}) \mid \exists (u, v) \in (V\_N \times V\_P) \cap E \}. \end{aligned}$$

For a family μ = (μu ∈ D(uE))u∈VN of full-support distributions we denote by μ/∼ the family ν = (νu˜ ∈ D(˜uR))u˜∈SN defined as follows. For all ˜u ∈ SN and all v˜ ∈ uR˜ we have νu˜(˜v) = - w∈v˜ μu(w), where u is any element of ˜u.

The following property of the ∼-quotient follows from the fact that all the vertices in ˜v have the same maximal probability of reaching the target vertices.

**Theorem 1.** *Consider a target arena* A = (V,VP ,E,T)*. For all families* μ = (μu ∈ D(uE))u∈VN *of full-support probability distributions and all* v ∈ VP *we have*

$$\max\_{\sigma} \mathbb{P}^{\nu}\_{\mathcal{A}^{\sigma}\_{\mu}}[\Diamond T] = \max\_{\sigma'} \mathbb{P}^{\tilde{v}}\_{\mathcal{B}^{\sigma'}\_{\nu}}[\Diamond U],$$

*where* B = A/∼*,* ν = μ/∼*, and* U = {t ˜ | ∃t ∈ T}*.*

We can further remove edges that lead to sub-optimal Nature vertices. When this is done after ∼-quotienting the maximal reachability probabilities are preserved.

**Theorem 2.** *Consider a target arena* A = (V,VP ,E,T) *such that* A/∼ = A*. For all families* μ = (μu ∈ D(uE))u∈VN *of full-support probability distributions, for all* (w, x) ∈ E ∩(VP × VN ) *such that* x(wE \ {x})*, and all* v ∈ VP *we have*

$$\max\_{\sigma} \mathbb{P}^{\upsilon}\_{\mathcal{A}^{\sigma}\_{\mu}}[\Diamond T] = \max\_{\sigma'} \mathbb{P}^{\upsilon}\_{\mathcal{B}^{\sigma'}\_{\mu}}[\Diamond T],$$

*where* B = (V,VP , E \ {(w, x)}, T)*.*

#### **3.2 Known Efficiently-Computable Special Cases**

We now recall the definitions of the set of extremal-probability states, end components, and essential states. Then, we observe that for all these sets of states their maximal probability reachability coincide and their definitions are independent of the probabilities labelling the transitions of the MDP. Hence, they are subsets of the set of the equivalence classes induced by ∼.

**Extremal-Probability States.** The set of *extremal-probability states* of an MDP M = (Q, A, δ, T) consists of the set of states with maximal probability reachability 0 and 1. Both sets can be computed in polynomial time [1,4]. We give below a game-based definition of both sets inspired by the classical polynomialtime algorithm to compute them (see, e.g., [1]). Let us fix a target arena A = (V,VP ,E,T) for the sequel.

For a set T ⊆ V , let us write **Z**T := {v ∈ V | T is not reachable from v}.

*(Almost-Surely Winning) Strategies.* A strategy for Protagonist in a target arena is a function σ : VP → VN . We then say that a path v0 ...vn ∈ V ∗ is *consistent with* σ if vi ∈ VP =⇒ σ(vi) = vi+1 for all 0 ≤ i<n. Let **Reach**(v0, σ) denote the set of vertices reachable from v0 under σ, i.e. **Reach**(v0, σ) := {vk | v0 ...vk is a path consistent with σ}.

We say that a strategy σ for Protagonist is *almost-surely winning from* u0 ∈ V *to* T ⊆ VP if, after modifying the arena to make all t ∈ T into sinks, for all v0 ∈ **Reach**(u0, σ) we have **Reach**(v0, σ)∩T = ∅. We denote the set of all such strategies by **Win**v0 T .

The following properties regarding almost-surely winning strategies in a target arena follow from the correctness of the graph-based algorithm used to compute extremal-probability states in an MDP [1, Lemma 10.108].

**Lemma 2 (From** [1]**).** *Consider a target arena* A = (V,VP ,E,T)*. For all families* μ = (μu ∈ D(uE))u∈VN *of full-support probability distributions, for all* v ∈ VP *the following hold.*

*(i)* maxσ Pv Aσ μ [♦T]=0 ⇐⇒ v ∈ **Z**T *(ii)* ∀σ : σ ∈ **Win**v T ⇐⇒ Pv Aσ μ [♦T]=1

**End Components.** Let us consider an MDP M = (Q, A, δ, T). A set S ⊆ Q of states is an *end component* in M if for all pairs of states p, q ∈ S there exists a strategy σ such that Pp Mσ [S U q] = 1.

*Example 3.* Let us consider the MDP shown on the left in Fig. 2. The set {p, q} is an end component since, by playing a from both states, one can ensure to reach either state from the other with probability 1.

It follows immediately from the definition of end component that the maximal probability of reaching T from states in the same end component is the same.

**Lemma 3.** *Let* S ⊆ Q *be an end component in* M*. For all* p, q ∈ S *we have that* maxσ Pp Mσ [♦T] = maxσ Pq Mσ [♦T]*.*

We say an end component is *maximal* if it is maximal with respect to set inclusion. Furthermore, from the definition of end components in MDPs and Lemma 2 it follows that we can lift the notion of end component to target arenas. More precisely, a set S ⊆ VP is an end component in A if and only if for some family of full-support probability distributions μ we have that S is an end component in Aμ (if and only if for all μ the set S is an end component in Aμ-).

The set of all maximal end components of a target arena can be computed in polynomial time using an algorithm based on the strongly connected components of the graph [1,8].

**Essential States.** Consider a target arena A = (V,VP ,E,T) and let be the smallest relation satisfying the following. For all u ∈ VP we have u u. For all u0, v ∈ VP \ **Z**T such that u0 = v we have u0 v if for all paths u0u1u2 we have that u2 v and there is at least one such path. Intuitively, u v holds whenever all paths starting from u reach v. In [7], the maximal vertices according to are called *essential states*2.

**Lemma 4 (From** [7]**).** *Consider a target arena* A = (V,VP ,E,T)*. For all families* μ = (μu ∈ D(uE))u∈VN *of full-support probability distributions, for all* v ∈ VP *and all essential states* w*, if* v w *then* maxσ Pv Aσ μ [♦T] = maxσ- Pw Aσ- μ [♦T]*.*

Note that, in the left arena in Fig. 3, p t does not hold since there is a cycle between p and q which does not visit t.

It was also shown in [7] that the relation is computable in polynomial time.

#### **4 Graph-Based Characterization of the NWR**

In this section we give a characterization of the NWR that is reminiscent of the topological-based value iteration proposed in [5]. The main intuition behind our characterization is as follows. If v W does not hold, then for all 0 <ε< 1 there is some family μ of full-support distributions such that Valμ(v) is at least 1 − ε, while Valμ(w) is at most ε for all w ∈ W. In turn, this must mean that there is a path from v to T which can be assigned a high probability by μ while, from W, all paths go with high probability to **Z**T .

We capture the idea of separating a "good" v–T path from all paths starting from W by using partitioning of V into layers Si ⊆ V . Intuitively, we would like it to be easy to construct a family μ of probability distributions such that from all vertices in Si+1 all paths going to vertices outside of Si+1 end up, with high probability, in lower layers, i.e. some Sk with k<i. A formal definition follows.

**Definition 6 (Drift partition and vertices).** *Consider a target arena* A = (V,VP ,E,T) *and a partition* (Si)0≤i≤k *of* V *. For all* 0 ≤ i ≤ k*, let* S+ i := ∪i<jSj *and* S− i := ∪j<iSj *, and let* Di := {v ∈ Si ∩ VN | vE ∩ S− i = ∅}*. We define the set* D := ∪0<i<kDi *of* drift vertices*. The partition is called a* drift partition *if the following hold.*

2 This is not the usual notion of essential states from classical Markov chain theory.

Using drift partitions, we can now formalize our characterization of the negation of the NWR.

**Theorem 3.** *Consider a target arena* A = (V,VP ,E,T)*, a non-empty set of vertices* W ⊆ V *, and a vertex* v ∈ V *. The following are equivalent*

Before proving Theorem 3 we need an additional definition and two intermediate results.

**Definition 7 (Value-monotone paths).** *Let* A = (V,VP ,E,T) *be a target arena and consider a family of full-support probability distributions* μ = (μu ∈ D(uE))u∈VN *. A path* v0 ...vk *is* μ-non-increasing *if and only if* Valμ(vi+1) ≤ Valμ(vi) *for all* 0 ≤ i<k*; it is* μ-non-decreasing *if and only if* Valμ(vi) ≤ Valμ(vi+1) *for all* 0 ≤ i<k*.*

It can be shown that from any path in a target arena ending in T one can obtain a simple non-decreasing one.

**Lemma 5.** *Consider a target arena* A = (V,VP ,E,T) *and a family of fullsupport probability distributions* μ = (μu ∈ D(uE))u∈VN *. If there is a path from some* v ∈ V *to* T*, there is also a simple* μ*-non-decreasing one.*

Additionally, we will make use of the following properties regarding vertexvalues. They formalize the relation between the value of a vertex, its owner, and the values of its successors.

**Lemma 6.** *Consider a target arena* A = (V,VP ,E,T) *and a family of fullsupport probability distributions* μ = (μu ∈ D(uE))u∈VN *.*

*(i) For all* u ∈ VP *, for all successors* v ∈ uE *it holds that* Valμ(v) ≤ Valμ(u)*. (ii) For all* u ∈ VN *it holds that*

$$(\exists v \in uE : \text{Val}^{\mu}(u) < \text{Val}^{\mu}(v)) \implies (\exists w \in uE : \text{Val}^{\mu}(w) < \text{Val}^{\mu}(u)).$$

*Proof (of Theorem* 3*).* Recall that, by definition, (i) holds if and only if there exists a family μ = (μu ∈ D(uE))u∈VN of full-support probability distributions such that ∀w ∈ W : Valμ(w) < Valμ(v).

Let us prove (i) =⇒ (ii). Let x0 < x1 <... be the finitely many (i.e. at most |V |) values that occur in the MDP Aμ, and let k be such that Valμ(v) = xk. For all 0 ≤ i<k let Si := {u ∈ V | Valμ(u) = xi}, and let Sk := V \ ∪i<kSi. Let us show below that the Si form a drift partition.

We have that Valμ(w) < Valμ(v) = xk for all w ∈ W, by assumption, so W ⊆ S− k by construction. By Lemma 5 there exists a simple μ-non-decreasing path π from v to T, so all the vertices occurring in π have values at least Valμ(v), so π ⊆ Sk.

We will prove (ii) =⇒ (i) by defining some full-support distribution family μ. The definition will be partial only, first on π ∩ VN , and then on the drift vertices in V \ Sk. Let 0 <ε< 1, which is meant to be small enough. Let us write π = v0 ...vn so that v0 = v and vn ∈ T. Let us define μ on π ∩ VN as follows: for all i<n, if vi ∈ VN let μvi (vi+1) := 1 − ε. Let σ be an arbitrary Protagonist strategy such that for all i<n, if vi ∈ VP then σ(vi) := vi+1. Therefore

$$\begin{aligned} (1 - \varepsilon)^{|V|} &\le (1 - \varepsilon)^n &\text{since } \pi \text{ is simple} \\ &\le \prod\_{i < n, v\_i \in S\_N} \mu\_{v\_i}(v\_{i+1}) &\text{by definition of } \mu \\ &\le \mathbb{P}\_{\mathcal{A}\_\mu^\sigma}^v[\Diamond T] \\ &\le \max\_{\sigma'} \mathbb{P}\_{\mathcal{A}\_\mu^\sigma}^v[\Diamond T] = \text{Val}^\mu(v). \end{aligned}$$

So, for 0 <ε< 1 − 1 |V√| 2 , we have 1 2 < (1 − ε)|V | ≤ Valμ(v). Below we will further define μ such that Valμ(w) ≤ 1 − (1 − ε)|V | < 1 2 for all w ∈ W and all 0 <ε< 1 − 1 |V√| 2 , which will prove (ii) =⇒ (i). However, the last part of the proof is more difficult.

For all 1 ≤ i ≤ k, for all drift vertices u ∈ Si, let (u) be a successor of u in S− i . Such a (u) exists by definition of the drift vertices. Then let μu((u)) := 1 − ε. We then claim that

$$\forall u \in D: (1 - \varepsilon)(1 - \mathbb{P}^{\varrho(u)}\_{\mathcal{A}\_{\mu}^{\sigma}}[\Diamond T]) \le 1 - \mathbb{P}^{u}\_{\mathcal{A}\_{\mu}^{\sigma}}[\Diamond T]. \tag{2}$$

Indeed, 1 − Pu Aσ μ [♦T] is the probability that, starting at u and following σ, T is never reached; and (1 − ε)(1 − P(u) Aσ μ [♦T]) is the probability that, starting at u and following σ, the second vertex is (u) and T is never reached.

Now let σ be an arbitrary strategy, and let us prove the following by induction on j.

$$\forall 0 \le j < k, \forall w \in S\_j \cup S\_j^- : \mathbb{P}\_{\mathcal{A}\_\mu^\sigma}^w[\Diamond T] \le 1 - (1 - \varepsilon)^j$$

Base case, j = 0: by assumption W is non-empty and included in S− k , so 0 < k. Also by assumption T ⊆ Sk, so T ∩ S0 = ∅. By definition of a drift partition, there are no edges going out of S0, regardless of whether the starting vertex is in VP or VN . So there is no path from w to T, which implies Valμ(w)=0 for all w ∈ S0, and the claim holds for the base case. Inductive case, let w ∈ Sj , let D := D ∩ (Sj ∪ S− j ) and let us argue that every path π from w to T must at some point leave Sj ∪ S− j to reach a vertex with higher index, i.e. there is some edge (πi, πi+1) from πi ∈ Sj ∪ S− j to some πi+1 ∈ S with j< . By definition of a drift partition, πi must also be a drift vertex, i.e. πi ∈ D . Thus, if we let F := VP \ D , Lemma 1 implies that Pw Aσ μ [♦T] = - u∈D- Pw Aσ μ [F U u] Pu Aσ μ [♦T]. Now, since

$$\begin{split} & \sum\_{u \in D'} \mathbb{P}\_{\mathcal{A}\_{\mu}^{u}}^{u}[\Diamond T] \\ &= \sum\_{u \in D \cap S\_{j}^{-}} \mathbb{P}\_{\mathcal{A}\_{\mu}^{u}}^{u}[\Diamond T] + \sum\_{u \in D\_{\mu}} \mathbb{P}\_{\mathcal{A}\_{\mu}^{u}}^{u}[\Diamond T] \qquad \text{by splitting the sum} \\ & \leq \sum\_{u \in D \cap S\_{j}^{-}} \mathbb{P}\_{\mathcal{A}\_{\mu}^{u}}^{u}[\Diamond T] + \sum\_{u \in D\_{j}} (1 - (1 - \varepsilon)(1 - \mathbb{P}\_{\mathcal{A}\_{\mu}^{\mathcal{O}}}^{\varrho(u)}[\Diamond T])) \qquad \text{by (2)} \\ & \leq \sum\_{u \in D \cap S\_{j}^{-}} (1 - (1 - \varepsilon)^{j-1}) + \\ & \sum\_{u \in D \cap S\_{j}^{-}} (1 - (1 - \varepsilon)(1 - \varepsilon)^{j-1}) \qquad \forall x \in D\_{j} : \varrho(x) \in S\_{j}^{-} \\ & \leq \sum\_{u \in D'} (1 - (1 - \varepsilon)^{j}) \qquad \quad (1 - \varepsilon)^{j} \leq (1 - \varepsilon)^{j-1} \end{split}$$

and - u∈D- Pw Aσ μ [F U u] ≤ 1, we have that Pw Aσ μ [♦T] ≤ 1−(1−ε)j . The induction is thus complete. Since σ is arbitrary in the calculations above, and since j < k ≤ |V |, we find that Valμ(w) ≤ 1 − (1 − ε)|V | for all w ∈ W ⊆ S− k .

For 0 <ε< 1 − 1 |V√| 2 we have 1 2 < (1 − ε)|V | , as mentioned after (1), so Valμ(w) ≤ 1 − (1 − ε)|V | < 1 2 .

### **5 Intractability of the NWR**

It follows from Theorem 3 that we can decide whether a vertex is sometimes worse than a set of vertices by guessing a partition of the vertices and verifying that it is a drift partition. The verification can clearly be done in polynomial time.

**Corollary 1.** *Given a target arena* A = (V,VP ,E,T)*, a non-empty set* W ⊆ V *, and a vertex* v ∈ V *, determining whether* v W *is decidable and in* coNP*.*

We will now show that the problem is in fact coNP-complete already for Markov chains.

**Theorem 4.** *Given a target arena* A = (V,VP ,E,T)*, a non-empty vertex set* W ⊆ V *, and a vertex* v ∈ V *, determining whether* v W *is* coNP*-complete even if* |uE| = 1 *for all* u ∈ VP *.*

The idea is to reduce the 2-Disjoint Paths problem (2DP) to the existence of a drift partition witnessing that v {w} does not hold, for some v ∈ V . Recall that 2DP asks, given a directed graph G = (V,E) and vertex pairs (s1, t1),(s2, t2) ∈ V × V , whether there exists an s1–t1 path π1 and an s2–t2 path π2 such that π1 and π2 are vertex disjoint, i.e. π1 ∩ π2 = ∅. The problem is known to be NP-complete [10,12]. In the sequel, we assume without loss of generality that (a) t1 and t2 are reachable from all s ∈ V \ {t1, t2}; and (b) t1 and t2 are the only sinks G.

*Proof (of Theorem* 4*).* From the 2DP input instance, we construct the target arena A = (S, SP , R, T) with S := V ∪ E, R := {(u,u, v ),(u, v , v) ∈ S × S | (u, v) ∈ E or u = v ∈ {t1, t2}}, SP := V × V , and T := {t1, t1 }. We will show there are vertex-disjoint s1–t1 and s2–t2 paths in G if and only if there is a drift partition (Si)0≤i≤k and a simple s1–t1 path π such that π ⊆ Sk and s2 ∈ S− k . The result will then follow from Theorem 3.

Suppose we have a drift partition (Si)0≤i≤k with s2 ∈ S− k and a simple path π = v0v0, v1 ...vn−1, vn vn with v0 = s1, vn = t1. Since the set {t2,t2, t2 } is *trapping* in A, i.e. all paths from vertices in the set visit only vertices from it, we can assume that S0 = {t2,t2, t2 }. (Indeed, for any drift partition, one can obtain a new drift partition by moving any trapping set to a new lowest layer.) Now, using the assumption that t2 is reachable from all s ∈ V \ {t1, t2} one can show by induction that for all 0 ≤ j<k and for all = u0 ∈ Sj there is a path u0 ...um in G with um = t2 and ⊆ S− j+1. This implies that there is a s2–t2 path π2 in G such that π2 ⊆ S− k . It follows that π2 is vertex disjoint with the s1–t1 path v0 ...vn in G.

Now, let us suppose that we have s1–t1 and s2–t2 vertex disjoint paths π1 = u0 ...un and π2 = v0 ...vm. Clearly, we can assume both π1, π2 are simple. We will construct a partition (Si)0≤i≤m+1 and show that it is indeed a drift partition, that u0u0, u1 ...un−1, un un ⊆ Sm+1, and s2 = v0 ∈ S− m+1. Let us set S0 := {vm−1, vm , vm,t2, t2 }, Si := {vm−i−1, vm−i , vm−i} for all 0 < i ≤ m, and Sm+1 := S \ ∪0≤i≤mSi. Since π2 is simple, (Si)0≤i≤m+1 is a partition of V . Furthermore, we have that s2 = v0 ∈ S− m+1, and u0u0, u1 ...un−1, un un ⊆ Sm+1 since π1 and π2 are vertex disjoint. Thus, it only remains for us to argue that for all 0 ≤ i ≤ m + 1: for all w ∈ Si ∩SN we have wR ∩S+ i = ∅, and for all w ∈ Si ∩ VN we have wR ∩ S+ i = ∅ =⇒ wR ∩ S− i = ∅. By construction of the Si, we have that eR ⊆ Si for all 0 ≤ i ≤ m and all e ∈ Si ∩ SP . Furthermore, for all 0 < i ≤ m, for all x ∈ Si ∩ SN = {vm−i}, there exists y ∈ Si−1 ∩ SP = {vm−i, vm−i+1 } such that (x, y) ∈ R—induced by (vm−i, vm−1+1) ∈ E from π2. To conclude, we observe that since S0 = {vm−1, vm , vm = t2,t2, t2 } and {t2,t2, t2 } is trapping in A, the set t2R is contained in S0.

#### **6 Efficiently Under-Approximating the NWR**

Although the full NWR cannot be efficiently computed for a given MDP, we can hope for "under-approximations" that are accurate and efficiently computable.

**Definition 8 (Under-approximation of the NWR).** *Let* A = (V,VP ,E,T) *be a target arena and consider a relation* : V × P(V )*. The relation is an* under-approximation *of the NWR if and only if* ⊆ *.*

We denote by ∗ the *pseudo transitive closure* of . That is, ∗ is the smallest relation such that ⊆∗ and for all u ∈ V,X ⊆ V if there exists W ⊆ V such that u ∗ W and w ∗ X for all w ∈ W, then u ∗ X.

*Remark 1.* The empty set is an under-approximation of the NWR. For all underapproximations of the NWR, the pseudo transitive closure ∗ of is also an under-approximation of the NWR.

In [2], efficiently-decidable sufficient conditions for the NWR were given. In particular, those conditions suffice to infer relations such as those in the right MDP from Fig. 3. We recall (Proposition 1) and extend (Proposition 2) these conditions below.

**Proposition 1 (From** [2]**).** *Consider a target arena* A = (V,VP ,E,T) *and an under-approximation of the NWR. For all vertices* v0 ∈ V *, and sets* W ⊆ V *the following hold.*

*Proof (Sketch).* The main idea of the proof of item (i) is to note that S is visited before T. The desired result then follows from Lemma 1. For item (ii), we intuitively have that there is a strategy to visit T with some probability or visit W, where the chances of visiting T are worse than before. We then show that it is never worse to start from v0 to have better odds of visiting T.

The above "rules" give an iterative algorithm to obtain increasingly better under-approximations of the NWR: from i apply the rules and obtain a new under-approximation i+1 by adding the new pairs and taking the pseudo transitive closure; then repeat until convergence. Using the special cases from Sect. 3.2 we can obtain a nontrivial initial under-approximation 0 of the NWR in polynomial time.

The main problem is how to avoid testing all subsets W ⊆ V in every iteration. One natural way to ensure we do not consider all subsets of vertices in every iteration is to apply the rules from Proposition 1 only on the successors of Protagonist vertices.

In the same spirit of the iterative algorithm described above, we now give two new rules to infer NWR pairs.

**Proposition 2.** *Consider a target arena* A = (V,VP ,E,T) *and an underapproximation of the NWR.*

*Proof (Sketch).* Item (i) follows immediately from the definition of Val. For item (ii) one can use the Bellman optimality equations for infinite-horizon reachability in MDPs to show that since the successors of v are never worse than the non-dominated successors of u, we must have u {v}.

**Fig. 4.** Two target arenas with *T* = {*fin*} are shown. Using Propositions 1 and 2 one can conclude that *p* ∼ *q* in both target arenas.

The rules stated in Proposition 2 can be used to infer relations like those depicted in Fig. 4 and are clearly seen to be computable in polynomial time as they speak only of successors of vertices.

#### **7 Conclusions**

We have shown that the never-worse relation is, unfortunately, not computable in polynomial time. On the bright side, we have extended the iterative polynomialtime algorithm from [2] to under-approximate the relation. In that paper, a prototype implementation of the algorithm was used to empirically show that interesting MDPs (from the set of benchmarks included in PRISM [17]) can be drastically reduced.

As future work, we believe it would be interesting to implement an exact algorithm to compute the NWR using SMT solvers. Symbolic implementations of the iterative algorithms should also be tested in practice. In a more theoretical direction, we observe that the planning community has also studied maximizing the probability of reaching a target set of states under the name of MAXPROB (see, e.g., [16,21]). There, online approximations of the NWR would make more sense than the under-approximation we have proposed here. Finally, one could define a notion of never-worse for finite-horizon or quantitative objectives.

**Acknowledgements.** The research leading to these results was supported by the ERC Starting grant 279499: inVEST. Guillermo A. P´erez is an F.R.S.-FNRS Aspirant and FWA postdoc fellow.

We thank Nathana¨el Fijalkow for pointing out the relation between this work and the study of interval MDPs and numberless MDPs. We also thank Shaull Almagor, Micha¨el Cadilhac, Filip Mazowiecki, and Jean-Fran¸cois Raskin for useful comments on earlier drafts of this paper.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **A Hierarchy of Scheduler Classes for Stochastic Automata**

Pedro R. D'Argenio1,2,3, Marcus Gerhold4 , Arnd Hartmanns4(B) , and Sean Sedwards5

> 1 Universidad Nacional de C´ordoba, C´ordoba, Argentina<EMAIL_ADDRESS>2 CONICET, C´ordoba, Argentina 3 Saarland University, Saarbr¨ucken, Germany 4 University of Twente, Enschede, The Netherlands<EMAIL_ADDRESS>5 University of Waterloo, Waterloo, Canada<EMAIL_ADDRESS>
**Abstract.** Stochastic automata are a formal compositional model for concurrent stochastic timed systems, with general distributions and nondeterministic choices. Measures of interest are defined over *schedulers* that resolve the nondeterminism. In this paper we investigate the power of various theoretically and practically motivated classes of schedulers, considering the classic complete-information view and a restriction to non-prophetic schedulers. We prove a hierarchy of scheduler classes w.r.t. unbounded probabilistic reachability. We find that, unlike Markovian formalisms, stochastic automata distinguish most classes even in this basic setting. Verification and strategy synthesis methods thus face a tradeoff between powerful and efficient classes. Using lightweight scheduler sampling, we explore this tradeoff and demonstrate the concept of a useful approximative verification technique for stochastic automata.

### **1 Introduction**

The need to analyse continuous-time stochastic models arises in many practical contexts, including critical infrastructures [4], railway engineering [36], space mission planning [7], and security [28]. This has led to a number of discrete event simulation tools, such as those for networking [34,35,42], whose probabilistic semantics is founded on generalised semi-Markov processes (GSMP [21,33]). Nondeterminism arises through inherent concurrency of independent processes [11], but may also be deliberate underspecification. Modelling such uncertainty with probability is convenient for simulation, but not always adequate [3,29]. Various models and formalisms have thus been proposed to extend continuous-time

This work is supported by the 3TU.BSR, NWO BEAT (602.001.303) and JST ERATO HASUO Metamathematics for Systems Design (JPMJER1603) projects, by ERC grant 695614 (POWVER), and by SeCyT-UNC projects 05/BP12, 05/B497.

c The Author(s) 2018

C. Baier and U. Dal Lago (Eds.): FOSSACS 2018, LNCS 10803, pp. 384–402, 2018. https://doi.org/10.1007/978-3-319-89366-2\_21

stochastic processes with nondeterminism [8,10,19,23,27,38]. It is then possible to *verify* such systems by considering the extremal probabilities of a property. These are the supremum and infimum of the probabilities of the property in the purely stochastic systems induced by classes of *schedulers* (also called *strategies*, *policies* or *adversaries*) that resolve all nondeterminism. If the nondeterminism is considered controllable, one may alternatively be interested in the *planning* problem of synthesising a scheduler that satisfies certain probability bounds.

We consider closed systems of stochastic automata (SA [16]), which extend GSMP and feature both generally distributed stochastic delays as well as discrete nondeterministic choices. The latter may arise from non-continuous distributions (e.g. deterministic delays), urgent edges, and edges waiting on multiple clocks. Numerical verification algorithms exist for very limited subclasses of SA only: Buchholz et al. [13] restrict to phase-type or matrix-exponential distributions, such that nondeterminism cannot arise (as each edge is guarded by a single clock). Bryans et al. [12] propose two algorithms that require an a priori fixed scheduler, continuous bounded distributions, and that all active clocks be reset when a location is entered. The latter forces regeneration on every edge, making it impossible to use clocks as memory between locations. Regeneration is central to the work of Ballarini et al. [6], however they again exclude nondeterminism. The only approach that handles nondeterminism is the region-based approximation scheme of Kwiatkowska et al. [30] for a model closely related to SA, but restricted to bounded continuous distributions. Without that restriction [22], error bounds and convergence guarantees are lost.

Evidently, the combination of nondeterminism and continuous probability distributions is a particularly challenging one. With this paper, we take on the underlying problem from a fundamental perspective: we investigate the power of, and relationships between, different classes of schedulers for SA. Our motivation is, on the one hand, that a clear understanding of scheduler classes is crucial to design verification algorithms. For example, Markov decision process (MDP) model checking works well because memoryless schedulers suffice for reachability, and the efficient time-bounded analysis of continuous-time MDP (CTMDP) exploits a relationship between two scheduler classes that are sufficiently simple, but on their own do not realise the desired extremal probabilities [14]. When it comes to planning problems, on the other hand, practitioners desire *simple* solutions, i.e. schedulers that need little information and limited memory, so as to be explainable and suitable for implementation on e.g. resource-constrained embedded systems. Understanding the capabilities of scheduler classes helps decide on the tradeoff between simplicity and the ability to attain optimal results.

We use two perspectives on schedulers from the literature: the classic complete-information *residual lifetimes* semantics [9], where optimality is defined via history-dependent schedulers that see the entire current state, and *nonprophetic* schedulers [25] that cannot observe the timing of *future* events. Within each perspective, we define classes of schedulers whose views of the state and history are variously restricted (Sect. 3). We prove their relative ordering w.r.t. achieving optimal reachability probabilities (Sect. 4). We find that SA distinguish most classes. In particular, memoryless schedulers suffice in the completeinformation setting (as is implicit in the method of Kwiatkowska et al. [30]), but turn out to be suboptimal in the more realistic non-prophetic case. Considering only the relative order of clock expiration times, as suggested by the first algorithm of Bryans et al. [12], surprisingly leads to partly suboptimal, partly incomparable classes. Our distinguishing SA are small and employ a common nondeterministic gadget. They precisely pinpoint the crucial differences and how schedulers interact with the various features of SA, providing deep insights into the formalism itself.

Our study furthermore forms the basis for the application of *lightweight scheduler sampling* (LSS) to SA. LSS is a technique to use Monte Carlo simulation/statistical model checking with nondeterministic models. On every LSS simulation step, a pseudo-random number generator (PRNG) is re-seeded with a hash of the identifier of the current scheduler and the (restricted) information about the current state (and previous states, for history-dependent schedulers) that the scheduler's class may observe. The PRNG's first iterate then determines the scheduler's action deterministically. LSS has been successfully applied to MDP [18,31,32] and probabilistic timed automata [15,26]. Using only constant memory, LSS samples schedulers uniformly from a selected scheduler class to find "near-optimal" schedulers that conservatively approximate the true extremal probabilities. Its principal advantage is that it is largely indifferent to the size of the state space and of the scheduler space; in general, sampling efficiency depends only on the likelihood of selecting near-optimal schedulers. However, the mass of *near*-optimal schedulers in a scheduler class that also includes the optimal scheduler may be *less* than the mass in a class that does *not* include it. Given that the mass of optimal schedulers may be vanishingly small, it may be advantageous to sample from a class of less powerful schedulers. We explore these tradeoffs and demonstrate the concept of LSS for SA in Sect. 5.

**Other Related Work.** Alur et al. first mention nondeterministic stochastic systems similar to SA in [2]. Markov automata (MA [19]), interactive Markov chains (IMC [27]) and CTMDP are special cases of SA restricted to exponential distributions. Song et al. [37] look into partial information distributed schedulers for MA, combining earlier works of de Alfaro [1] and Giro and D'Argenio [20] for MDP. Their focus is on information flow and hiding in parallel specifications. Wolf et al. [39] investigate the power of classic (time-abstract, deterministic and memoryless) scheduler classes for IMC. They establish (non-strict) subset relationships for almost all classes w.r.t. trace distribution equivalence, a very strong measure. Wolovick and Johr [41] show that the class of measurable schedulers for CTMDP is complete and sufficient for reachability problems.

### **2 Preliminaries**

For a given set S, its power set is P(S). We denote by R, R+, and R+ 0 the sets of real numbers, positive real numbers and non-negative real numbers, respectively. A (discrete) *probability distribution* over a set Ω is a function μ: Ω → [0, 1], such that support(μ) def = { ω ∈ Ω | μ(ω) > 0 } is countable and - ω∈support(μ) μ(ω) = 1. Dist(Ω) is the set of probability distributions over Ω. We write D(ω) for the *Dirac* distribution for ω, defined by D(ω)(ω) = 1. Ω is *measurable* if it is endowed with a σ-algebra σ(Ω): a collection of *measurable* subsets of Ω. A (continuous) *probability measure* over Ω is a function μ: σ(Ω) → [0, 1], such that μ(Ω)=1 and μ(∪i∈I Bi) = - i∈I μ(Bi) for any countable index set I and pairwise disjoint measurable sets Bi ⊆ Ω. Prob(Ω) is the set of probability measures over Ω. Each μ ∈ Dist(Ω) induces a probability measure. Given probability measures μ1 and μ2, we denote by μ1 ⊗ μ2 the *product measure*: the unique probability measure such that (μ1 ⊗ μ2)(B1 × B2) = μ1(B1) · μ2(B2), for all measurable B1 and B2. For a collection of measures (μi)i∈I , we analogously denote the product measure by i∈I μi. Let *Val* def = V → R+ 0 be the set of valuations for an (implicit) set V of (non-negative real-valued) variables. **0** ∈ *Val* assigns value zero to all variables. Given X ⊆ V and v ∈ *Val*, we write v[X] for the valuation defined by v[X](x) = 0 if x ∈ X and v[X](y) = v(y) otherwise. For t ∈ R+ 0 , v + t is the valuation defined by (v + t)(x) = v(x) + t for all x ∈ V .

**Stochastic Automata** [16] extend labelled transition systems with stochastic *clocks*: real-valued variables that increase synchronously with rate 1 over time and expire some random amount of time after having been *restarted*. Formally:

**Definition 1.** *A* stochastic automaton *(SA) is a tuple Loc*, C, A, E, F, *init, where Loc is a countable set of* locations*,* C *is a finite set of* clocks*,* A *is the finite* action alphabet*, and* E : *Loc* → P(P(C) × A × P(C) × Dist(*Loc*)) *is the* edge function*, which maps each location to a finite set of edges that in turn consist of a* guard set *of clocks, a label, a* restart set *of clocks and a distribution over target locations.* F : C → Prob(R+ 0 ) *is the* delay measure function *that maps each clock to a probability measure, and init* ∈ *Loc is the* initial location*.*

We also write G,a,R −−−−→E μ for G, a, R, μ ∈ E(). W.l.o.g. we restrict to SA where edges are fully characterised by source state and action label, i.e. whenever G1,a,R1 −−−−−→E μ1 and G2,a,R2 −−−−−→E μ2, then G1 = G2, R1 = R2 and μ1 = μ2.

Intuitively, an SA starts in *init* with all clocks expired. An edge G,a,R −−−−→E μ may be taken only if all clocks in G are expired. If any edge is enabled, some edge must be taken (i.e. all actions are *urgent* and thus the SA is *closed*). When an edge is taken, its action is a, all clocks in R are restarted, other expired clocks remain expired, and we move to successor location with probability μ( ). There, another edge may be taken immediately or we may need to wait until some further clocks expire, and so on. When a clock c is restarted, the time until it expires is chosen randomly according to the probability measure F(c).

*Example 1.* We show an example SA, M0, in Fig. 1. Its initial location is 0. It has two clocks, x and y, with F(x) and F(y) both being the continuous uniform distribution over the interval [0, 1]. No time can pass in locations 0 and 1, since they have outgoing edges with empty guard sets. We omit action labels and assume every edge to have a unique label. On entering 1, both clocks are restarted. The choice of going to either 2 or 3 from 1 is nondeterministic, since

**Fig. 1.** Example SA M0 **Fig. 2.** Excerpt of the TPTS semantics of M0

the two edges are always enabled at the same time. In 2, we have to wait until the first of the two clocks expires. If that is x, we have to move to location ✓; if it is y, we have to move to ✗. The probability that both expire at the same time is zero. Location 3 behaves analogously, but with the target states interchanged.

**Timed Probabilistic Transition Systems** form the semantics of SA. They are finitely-nondeterministic uncountable-state transition systems:

**Definition 2.** *A (finitely nondeterministic)* timed probabilistic transition system *(TPTS) is a tuple* S, A ,T,s*init.* S *is a measurable set of states.* A = R+ A *is the* alphabet*, partitioned into* delays *in* R+ *and* jumps *in* A*.* T : S → P(A × Prob(S)) *is the* transition function*, which maps each state to a finite set of transitions, each consisting of a label in* A *and a measure over target states. The initial state is* s*init* ∈ S*. For all* s ∈ S*, we require* |T(s)| = 1 *if* ∃ t, μ ∈ T(s): t ∈ R+*, i.e. states admitting delays are deterministic.*

We also write s a −→T μ for a, μ ∈ T(s). A *run* is an infinite alternating sequence s0a0s1a1... ∈ (S×A )ω, with s0 = s*init*. A *history* is a finite prefix of a run ending in a state, i.e. an element of (S × A )∗ × S. Runs resolve all nondeterministic and probabilistic choices. A *scheduler* resolves only the nondeterminism:

**Definition 3.** *A measurable function* s : (S ×A )∗ ×S → Dist(A × Prob(S)) *is a* scheduler *if, for all histories* h ∈ (S × A )∗ × S*,* a, μ ∈ support(s(h)) *implies lst*h a −→T μ*, where lst*h *is the last state of* h*.*

Once a scheduler has chosen si a −→T μ, the successor state si+1 is picked randomly according to μ. Every scheduler s defines a probability measure Ps on the space of all runs. For a formal definition, see [40]. As is usual, we restrict to *non-Zeno* schedulers that make time diverge with probability one: we require Ps(Π∞) = 1, where Π∞ is the set of runs where the sum of delays is ∞. In the remainder of this paper we consider extremal probabilities of reaching a set of goal locations G:

**Definition 4.** *For* G ⊆ *Loc, let* JG def = { , v, e ∈ S | ∈ G }*. Let* S *be a class of schedulers. Then* PS min(G) *and* PS max(G) *are the minimum and maximum* reachability probabilities *for* G *under* S*, defined as* PS min(G) = infs∈S Ps(ΠJ*G* ) *and* PS max(G) = sups∈S Ps(ΠJ*G* )*, respectively.*

**Semantics of Stochastic Automata.** We present here the residual lifetimes semantics of [9], simplified for closed SA: any delay step must be of the minimum delay that makes some edge become enabled.

**Definition 5.** *The semantics of an SA* M = *Loc*, C, A, E, F, *init is the TPTS*

$$\mathbb{E}\left[M\right] = \langle Loc \times Val \times Val, A \uplus \mathbb{R}^+, T\_M, \langle \ell\_{init}, \mathbf{0}, \mathbf{0} \rangle \rangle$$

*where the states are triples* , v, e *of the current location , a valuation* v *assigning to each clock its current value, and a valuation* e *keeping track of all clocks' expiration times.* TM *is the smallest transition function satisfying inference rules*

$$\frac{\ell \xrightarrow{G, a, R} \Sigma \text{ } \mu \quad \text{En}(G, v, e)}{\langle \ell, v, e \rangle \xrightarrow{a}\_{T\_M} \mu \otimes \mathcal{D}(v[R]) \otimes \text{Sample}\_e^R}$$

$$\frac{\ell \in \mathbb{R}^+ \quad \exists \ell \xrightleftharpoons^{G, a, R} \mu \colon \text{En}(G, v + t, e) \quad \forall \, t' \in [0, t), \ell \xrightarrow{G, a, R} \nu \colon \text{En}(G, v + t', e)}{\langle \ell, v, e \rangle \xrightarrow{t \to\_M} \mathcal{D}(\langle \ell, v + t, e \rangle)}$$

*with* En(G, v, e) def = ∀ x ∈ G: v(x) ≥ e(x) *characterising the enabled edges and*

$$\text{Sample}\_e^R \overset{\text{def}}{=} \bigotimes\_{c \in \mathcal{C}} \begin{cases} F(c) & if \, c \in R \\ \mathcal{D}(e(c)) & if \, c \notin R. \end{cases}$$

The second rule creates *delay* steps of t time units if no edge is enabled from now until just before t time units have elapsed (third premise) but then, after exactly t time units, some edge becomes enabled (second premise). The first rule applies if an edge G,a,R −−−−→E μ is enabled: a transition is taken with the edge's label, the successor state's location is chosen by μ, v is updated by resetting the clocks in R to zero, and the expiration times for the restarted clocks are resampled. All other expiration times remain unchanged. Notice that [[M]] is also a nondeterministic labelled Markov process [40] (a proof can be found in [17]).

*Example 2.* Figure 2 outlines the semantics of M0. The first step from 0 to all the states in 1 is a single transition. Its probability measure is the product of F(x) and F(y), sampling the expiration times of the two clocks. We exemplify the behaviour of all of these states by showing it for the case of expiration times e(x) and e(y), with e(x) < e(y). In this case, to maximise the probability of reaching ✓, we should take the transition to the state in 2. If a scheduler s can see the expiration times, noting that only their order matters here, it can always make the optimal choice and achieve P{s} max({ ✓ }) = 1.

### **3 Classes of Schedulers**

We now define classes of schedulers for SA with restricted information, hiding in various combinations the history and parts of states such as clock values and expiration times. All definitions consider TPTS as in Definition 5 with states , v, e and we require for all s that a, μ ∈ support(s(h)) ⇒ *lst*h a −→T μ, as in Definition 3.

#### **3.1 Classic Schedulers**

We first consider the "classic" complete-information setting where schedulers can in particular see expiration times. We start with restricted classes of historydependent schedulers. Our first restriction hides the values of all clocks, only revealing the total time since the start of the history. This is inspired by the stepcounting or time-tracking schedulers needed to obtain optimal step-bounded or time-bounded reachability probabilities on MDP or Markov automata:

**Definition 6.** *A classic history-dependent* global-time *scheduler is a measurable function* s : (S|,t,e × A )∗ × S|,t,e → Dist(A × Prob(S))*, where* S|,t,e def = *Loc* × R+ 0 ×*Val with the second component being the total time* t *elapsed since the start of the history. We write* S*hist* ,t,e *for the set of all such schedulers.*

We next hide the values of all clocks, revealing only their expiration times:

**Definition 7.** *A classic history-dependent* location-based *scheduler is a measurable function* s : (S|,e × A )∗ × S|,e → Dist(A × Prob(S))*, where* S|,e def = *Loc* × *Val , with the second component being the clock expiration times* e*. We write* S*hist* ,e *for the set of all such schedulers.*

Having defined three classes of classic history-dependent schedulers, S*hist* ,v,e, S*hist* ,t,e and S*hist* ,e , noting that S*hist* ,v,e denotes all schedulers of Definition 3, we also consider them with the restriction that they only see the relative order of clock expiration, instead of the exact expiration times: for each pair of clocks c1, c2, these schedulers see the relation ∼∈{<, =, >} in e(c1) − v(c1) ∼ e(c2) − v(c2). E.g. in 1 of Example 2, the scheduler would not see e(x) and e(y), but only whether e(x) < e(y) or vice-versa (since v(x) = v(y) = 0, and equality has probability 0 here). We consider this case because the expiration order is sufficient for the first algorithm of Bryans et al. [12], and would allow optimal decisions in M0 of Fig. 1. We denote the relative order information by o, and the corresponding scheduler classes by S*hist* ,v,o, S*hist* ,t,o and S*hist* ,o . We now define memoryless schedulers, which only see the current state and are at the core of e.g. MDP model checking. On most formalisms, they suffice to obtain optimal reachability probabilities.

**Definition 8.** *A classic* memoryless *scheduler is a measurable function* s : S → Dist(A × Prob(S))*. We write* S*ml* ,v,e *for the set of all such schedulers.*

We apply the same restrictions as for history-dependent schedulers:

**Definition 9.** *A classic memoryless global-time scheduler is a measurable function* s : S|,t,e → Dist(A × Prob(S))*, with* S|,t,e *as in Definition 6. We write* S*ml* ,t,e *for the set of all such schedulers.*

**Definition 10.** *A classic memoryless location-based scheduler is a measurable function* s : S|,e → Dist(A × Prob(S))*, with* S|,e *as in Definition 7. We write* S*ml* ,e *for the set of all such schedulers.*

Again, we also consider memoryless schedulers that only see the expiration order, so we have memoryless scheduler classes S*ml* ,v,e, S*ml* ,t,e, S*ml* ,e, S*ml* ,v,o, S*ml* ,t,o and S*ml* ,o. Class S*ml* ,o is particularly attractive because it has a compact finite domain.

#### **3.2 Non-prophetic Schedulers**

Consider the SA M0 in Fig. 1. No matter which of the previously defined scheduler classes we choose, we always find a scheduler that achieves probability 1 to reach ✓, and a scheduler that achieves probability 0. This is because they can all see the expiration times or expiration order of x and y when in 1. When in 1, x and y have not yet expired—this will only happen later, in 2 or 3—yet the schedulers already know which clock will "win". The classic schedulers can thus be seen to make decisions based on the timing of *future* events. This *prophetic* scheduling has already been observed in [9], where a "fix" in the form of the *spent lifetimes* semantics was proposed. Hartmanns et al. [25] have shown that this not only still permits prophetic scheduling, but even admits *divine* scheduling, where a scheduler can *change* the future. The authors propose a complex *non-prophetic* semantics that provably removes all prophetic and divine behaviour.

Much of the complication of the non-prophetic semantics of [25] is due to it being specified for open SA that include delayable actions. For the closed SA setting of this paper, prophetic scheduling can be more easily excluded by hiding from the schedulers all information about what will happen in the future of the system's evolution. This information is only contained in the expiration times e or the expiration order o. We can thus keep the semantics of Sect. 2 and modify the definition of schedulers to exclude prophetic behaviour by construction.

In what follows, we thus also consider all scheduler classes of Sect. 3.1 with the added constraint that the expiration times, resp. the expiration order, are not visible, resulting in the *non-prophetic* classes S*hist* ,v , S*hist* ,t , S*hist* , S*ml* ,v, S*ml* ,t and S*ml* . Any non-prophetic scheduler can only reach ✓ of M0 with probability 1 2 .

### **4 The Power of Schedulers**

Now that we have defined a number of classes of schedulers, we need to determine what the effect of the restrictions is on our ability to optimally control an SA. We thus evaluate the power of scheduler classes w.r.t. unbounded reachability probabilities (Definition 4) on the semantics of SA. We will see that this simple setting already suffices to reveal interesting differences between scheduler classes.

For two scheduler classes S1 and S2, we write S1 - S2 if, for all SA and all sets of goal locations G, PS1 min(G) ≤ PS2 min(G) and PS1 max(G) ≥ PS2 max(G). We write S1 S2 if additionally there exists at least one SA and set G where PS1 min(G ) < PS2 min(G ) or PS1 max(G ) > PS2 max(G ). Finally, we write S1 ≈ S2 for S1 - S2 ∧ S2 - S1, and S1 ≈ S2, i.e. the classes are incomparable, for S1 - S2 ∧ S2 - S1. Unless noted otherwise, we omit proofs for S1 - S2 when it is obvious that the information available to S1 includes the information available to S2. All our distinguishing examples are based on the resolution of a single nondeterministic choice between two actions to eventually reach one of two locations. We therefore prove only w.r.t. the maximum probability, pmax, for these examples since the minimum probability is given by 1 − pmax and an analogous proof for pmin can be made by relabelling locations. We may write Pmax(Sy x) for PS*y x* max({ ✓ }) to improve readability.

$$\begin{array}{ccccccccc}\mathfrak{S}^{ml}\_{\ell,o} & \prec & \mathfrak{S}^{ml}\_{\ell,t,o} & \prec & \mathfrak{S}^{hist}\_{\ell,o} & \prec & \mathfrak{S}^{hist}\_{\ell,o} \\ \mathfrak{R} & & & \nwarrow & \lambda & & \ll & \varnothing \\ \mathfrak{S}^{ml}\_{\ell,e} & \prec & \mathfrak{S}^{ml}\_{\ell,t,e} & \prec & \mathfrak{S}^{ml}\_{\ell,v,e} & \succ & \mathfrak{S}^{hist}\_{\ell,t,o} \\ \lambda & & \lambda & & \ll & \varnothing & & \ll \\ \mathfrak{S}^{hist}\_{\ell,e} & \simeq & \mathfrak{S}^{hist}\_{\ell,t,e} & \simeq & \mathfrak{S}^{hist}\_{\ell,v,e} & \succ & \mathfrak{S}^{hist}\_{\ell,v,o} \end{array}$$

$$\begin{array}{ccccc} \mathfrak{S}\_{\ell}^{ml} & \prec & \mathfrak{S}\_{\ell,t}^{ml} & \prec & \mathfrak{S}\_{\ell,v}^{ml} \\ \lambda & & \lambda & & \lambda \\ \mathfrak{S}\_{\ell}^{hist} & \approx & \mathfrak{S}\_{\ell,t}^{hist} & \approx & \mathfrak{S}\_{\ell,v}^{hist} \end{array}$$

**Fig. 3.** Hierarchy of classic scheduler classes **Fig. 4.** Non-prophetic classes

#### **4.1 The Classic Hierarchy**

We first establish that all classic history-dependent scheduler classes are equivalent:

**Proposition 1.** S*hist* ,v,e ≈ S*hist* ,t,e ≈ S*hist* ,e *.*

*Proof.* From the transition labels in A = A R+ in the history (S × A )∗, with S ∈ { S, S|,t,e, S|,e } depending on the scheduler class, we can reconstruct the total elapsed time as well as the values of all clocks: to obtain the total elapsed time, sum the labels in R+ up to each state; to obtain the values of all clocks, do the same per clock and perform the resets of the edges identified by the actions.

The same argument applies among the expiration-order history-dependent classes:

#### **Proposition 2.** S*hist* ,v,o ≈ S*hist* ,t,o ≈ S*hist* ,o *.*

However, the expiration-order history-dependent schedulers are strictly less powerful than the classic history-dependent ones:

#### **Proposition 3.** S*hist* ,v,e S*hist* ,v,o*.*

*Proof.* Consider the SA M1 in Fig. 5. Note that the history does not provide any information for making the choice in 1: we always arrive after having spent zero time in 0 and then having taken the single edge to 1. We can analytically determine that Pmax(S*hist* ,v,e) = 3 4 by going from 1 to 2 if e(x) ≤ 1 2 and to 3 otherwise. We would obtain a probability equal to 1 2 by always going to either 2 or 3 or by picking either edge with equal probability. This is the best we can do if e is not visible, and thus Pmax(S*hist* ,v,o) = 1 2 : in 1, v(x) = v(y) = 0 and the expiration order is always "y before x" because y has not yet been started.

Just like for MDP and unbounded reachability probabilities, the classic historydependent and memoryless schedulers with complete information are equivalent:

**Proposition 4.** S*hist* ,v,e ≈ S*ml* ,v,e*.*

*Proof sketch.* Our definition of TPTS only allows finite nondeterministic choices, i.e. we have a very restricted form of continuous-space MDP. We can thus adapt the argument of the corresponding proof for MDP [5, Lemma 10.102]: For each state (of possibly countably many), we construct a notional optimal memoryless (and deterministic) scheduler in the same way, replacing the summation by an integration for the continuous measures in the transition function. It remains to show that this scheduler is indeed measurable. For TPTS that are the semantics of SA, this follows from the way clock values are used in the guard sets so that optimal decisions are constant over intervals of clock values and expiration times (see e.g. the arguments in [12] or [30]).

On the other hand, when restricting schedulers to see the expiration order only, history-dependent and memoryless schedulers are no longer equivalent:

**Proposition 5.** S*hist* ,v,o S*ml* ,v,o*.*

*Proof.* Consider the SA M2 in Fig. 6. Let s *opt ml*(l,v,o) be the (unknown) optimal scheduler in S*ml* ,v,o w.r.t. the max. probability of reaching ✓. Define s*better hist*(l,v,o) ∈ S*hist* ,v,o as: when in 2 and the last edge in the history is the left one (i.e. x is expired), go to 3; otherwise, behave like s *opt ml*(l,v,o). This scheduler distinguishes S*hist* ,v,o and S*ml* ,v,o (by achieving a strictly higher max. probability than s *opt ml*(l,v,o)) if and only if there are some combinations of clock values (aspect v) and expiration orders (aspect o) in 2 that can be reached with positive probability via the left edge into 2, for which s *opt ml*(l,v,o) must nevertheless decide to go to 4.

All possible clock valuations in 2 can be achieved via either the left or the right edge, but taking the left edge implies that x expires before z in 2. It is thus sufficient to show that s *opt ml*(l,v,o) must go to 4 in *some* cases where x expires before z. The general form of schedulers in S*ml* ,v,o in 2 is "go to 3 iff (a) x expires before z and v(x) ∈ S1 or (b) z expires before x and v(x) ∈ S2" where the Si are measurable subsets of [0, 8]. S2 is in fact *irrelevant*: whatever s *opt ml*(l,v,o) does when (b) is satisfied will be mimicked by s*better hist*(l,v,o) because z can only expire before x when coming via the right edge into 2. Conditions (a) and (b) are independent.

With S1 = [0, 8], the max. probability is 77 96 = 0.80208¯3. Since this is the only scheduler in S*ml* ,v,o that is *relevant* for our proof and never goes to l4 when x expires before z, it remains to show that the max. probability under s *opt ml*(l,v,o) is > 77 96 . With S1 = [0, 35 12 ), we have a max. probability of 7561 9216 ≈ 0.820421. Thus s *opt ml*(l,v,o) must sometimes go to l4 even when the left edge was taken, so s*better hist*(l,v,o) achieves a higher probability and thus distinguishes the classes.

Knowing only the global elapsed time is less powerful than knowing the full history or the values of all clocks:

#### **Proposition 6.** S*hist* ,t,e S*ml* ,t,e *and* S*ml* ,v,e S*ml* ,t,e*.*

*Proof sketch.* Consider the SA M3 in Fig. 7. We have Pmax(S*hist* ,t,e) = 1: when in 3, the scheduler sees from the history which of the two incoming edges was used, and thus knows whether x or y is already expired. It can then make the optimal choice: go to 4 if x is already expired, or to 5 otherwise. We also have Pmax(S*ml* ,v,e) = 1: the scheduler sees that either v(x) = 0 or v(y) = 0, which implies that the other clock is already expired, and the argument above applies. However, Pmax(S*ml* ,t,e) < 1: the distribution of elapsed time t on entering 3 is itself independent of which edge is taken. With probability 1 4 , exactly one of e(x) and e(y) is below t in 3, which implies that that clock has just expired and thus the scheduler can decide optimally. Yet with probability 3 4 , the expiration times are not useful: they are both positive and drawn from the same distribution, but one unknown clock is expired. The wait for x in 1 ensures that comparing t with the expiration times in e does not reveal further information in this case.

In the case of MDP, knowing the total elapsed time (i.e. steps) does not make a difference for unbounded reachability. Only for step-bounded properties is that extra knowledge necessary to achieve optimal probabilities. With SA, however, it makes a difference even in the unbounded case:

#### **Proposition 7.** S*ml* ,t,e S*ml* ,e*.*

*Proof.* Consider SA M4 in Fig. 8. We have Pmax(S*ml* ,t,e) = 1: in 2, the remaining time until y expires is e(y) and the remaining time until x expires is e(x) − t for the global time value t as 2 is entered. The scheduler can observe all of these quantities and thus optimally go to 3 if x will expire first, or to 4 otherwise. However, Pmax(S*ml* ,e) < 1: e(x) only contains the absolute expiration time of x, but without knowing t or the expiration time of z in 1, and thus the current value v(x), this scheduler cannot know with certainty which of the clocks will expire first and is therefore unable to make an optimal choice in 2.

Finally, we need to compare the memoryless schedulers that see the clock expiration times with memoryless schedulers that see the expiration order. As noted in Sect. 3.1, these two views of the current state are incomparable unless we also see the clock values:

#### **Proposition 8.** S*ml* ,v,e S*ml* ,v,o*.*

*Proof.* S*ml* ,v,e S*ml* ,v,o follows from the same argument as in the proof of Proposition 3. S*ml* ,v,e - S*ml* ,v,o is because knowing the current clock values v and the expiration times e is equivalent to knowing the expiration order, since that is precisely the order of the differences e(c) − v(c) for all clocks c.

#### **Proposition 9.** S*ml* ,t,e ≈ S*ml* ,t,o*.*

*Proof.* S*ml* ,t,e S*ml* ,t,o follows from the same argument as in the proof of Proposition 3. For S*ml* ,t,e - S*ml* ,t,o, consider the SA M3 of Fig. 7. We know from the proof of Proposition 6 that Pmax(S*ml* ,t,e) < 1. However, if the scheduler knows the order in which the clocks will expire, it knows which one has already expired (the first one in the order), and can thus make the optimal choice in 3 to achieve Pmax(S*ml* ,t,o) = 1.

#### **Proposition 10.** S*ml* ,e ≈ S*ml* ,o*.*

*Proof.* The argument of Proposition 9 applies by observing that, in M3 of Fig. 7, we also have Pmax(S*ml* ,e) < 1 via the same argument as for S*ml* ,t,e in the proof of Proposition 6.

Among the expiration-order schedulers, the hierarchy is as expected:

**Proposition 11.** S*ml* ,v,o S*ml* ,t,o S*ml* ,o*.* *Proof sketch.* Consider M5 of Fig. 9. To maximise the probability, in 3 we should go to 4 whenever x is already expired or close to expiring, for which the amount of time spent in 2 is an indicator. S*ml* ,o only knows that x may have expired when the expiration order is "x before y", but definitely has not expired when it is "y before x". Schedulers in S*ml* ,t,o can do better: They also see the amount of time spent in 2. Thus S*ml* ,t,o S*ml* ,o. If we modify M5 by adding an initial delay on x from a new 0 to 1 as in M3, then the same argument can be used to prove S*ml* ,v,o S*ml* ,t,o: the extra delay makes knowing the elapsed time t useless with positive probability, but the exact time spent in l2 is visible to S*ml* ,v,o as v(x).

We have thus established the hierarchy of classic schedulers shown in Fig. 3, noting that some of the relationships follow from the propositions by transitivity.

#### **4.2 The Non-prophetic Hierarchy**

Each non-prophetic scheduler class is clearly dominated by the classic and expiration-order scheduler classes that otherwise have the same information, for example S*hist* ,v,e S*hist* ,v (with very simple distinguishing SA). We show that the non-prophetic hierarchy follows the shape of the classic case, including the difference between global-time and pure memoryless schedulers, with the notable exception of memoryless schedulers being weaker than history-dependent ones.

**Proposition 12.** S*hist* ,v ≈ S*hist* ,t ≈ S*hist .*

*Proof.* This follows from the argument of Proposition 1.

#### **Proposition 13.** S*hist* ,v S*ml* ,v*.*

*Proof.* Consider the SA M6 in Fig. 10. It is similar to M4 of Fig. 8, and our arguments are thus similar to the proof of Proposition 7. On M6, we have Pmax(S*hist* ,v ) = 1: in 2, the history reveals which of the two incoming edges was used, i.e. which clock is already expired, thus the scheduler can make the optimal choice. However, if neither the history nor e is available, we get Pmax(S*ml* ,v) = 1 2 : the only information that can be used in 2 are the values of the clocks, but v(x) = v(y), so there is no basis for an informed choice.

**Proposition 14.** S*hist* ,t S*ml* ,t *and* S*ml* ,v S*ml* ,t *.*

*Proof.* Consider the SA M3 in Fig. 7. We have Pmax(S*hist* ,t )=Pmax(S*ml* ,v) = 1, but Pmax(S*ml* ,t ) = 1 2 by the same arguments as in the proof of Proposition 6.

#### **Proposition 15.** S*ml* ,t S*ml .*

*Proof.* Consider the SA M4 in Fig. 8. The schedulers in S*ml* have no information but the current location, so they cannot make an informed choice in 2. This and the simple loop-free structure of M4 make it possible to analytically calculate the resulting probability: Pmax(S*ml* ) = 17 24 = 0.7083. If information about the global elapsed time t in 2 is available, however, the value of x is revealed. This allows making a better choice, e.g. going to 3 when t ≤ 1 2 and to 4 otherwise, resulting in Pmax(S*ml* ,t ) ≈ 0.771 (statistically estimated with high confidence).

We have thus established the hierarchy of non-prophetic schedulers shown in Fig. 4, where some relationships follow from the propositions by transitivity.

### **5 Experiments**

We have built a prototype implementation of lightweight scheduler sampling for SA by extending the Modest Toolset's [24] modes simulator, which already supports deterministic stochastic timed automata (STA [8]). With some care, SA can be encoded into STA. Using the original algorithm for MDP of [18], our prototype works by providing to the schedulers a discretised view of the continuous components of the SA's semantics, which, we recall, is a continuousspace MDP. The currently implemented discretisation is simple: for each realvalued quantity (the value v(c) of clock c, its expiration time e(c), and the global elapsed time t), it identifies all values that lie within the same interval [ i n , i+1 n ), for integers i, n. We note that better static discretisations are almost certainly possible, e.g. a region construction for the clock values as in [30].

We have modelled M1 through M6 as STA in Modest. For each scheduler class and model in the proof of a proposition, and discretisation factors n ∈ { 1, 2, 4 }, we sampled 10 000 schedulers and performed statistical model checking for each of them in the lightweight manner. In Fig. 11 we report the min. and max. estimates, (ˆpmin, pˆmax)..., over all sampled schedulers. Where different discretisations lead to different estimates, we report the most extremal values. The subscript denotes the discretisation factors that achieved the reported estimates. The analysis for each sampled scheduler was performed with a number of simulation runs sufficient for the overall max./min. estimates to be within ± 0.01 of the true maxima/minima of the *sampled* set of schedulers with probability ≥0.95 [18]. Note that ˆpmin is an upper bound on the true minimum probability and ˆpmax is a lower bound on the true maximum probability.

Increasing the discretisation factor or increasing the scheduler power generally increases the number of decisions the schedulers *can* make. This may also increase the number of *critical* decisions a scheduler *must* make to achieve the extremal probability. Hence, the sets of discretisation factors associated to specific experiments may be informally interpreted in the following way:

**Fig. 11.** Results from the prototype of lightweight scheduler sampling for SA

The results in Fig. 11 respect and differentiate our hierarchy. In most cases, we found schedulers whose estimates were within the statistical error of calculated optima or of high confidence estimates achieved by alternative statistical techniques. The exceptions involve M3 and M4. We note that M4 makes use of an additional clock, increasing the dimensionality of the problem and potentially making near-optimal schedulers rarer. The best result for M3 and class S*ml* l,t,e was obtained using discretisation factor n = 2: a compromise between nearness to optimality and rarity. A greater compromise was necessary for M4 and classes S*ml* l,t,e, S*ml* l,e , where we found near-optimal schedulers to be very rare and achieved best results using discretisation factor n = 1.

The experiments demonstrate that lightweight scheduler sampling can produce useful and informative results with SA. The present theoretical results will allow us to develop better abstractions for SA and thus to construct a refinement algorithm for efficient lightweight verification of SA that will be applicable to realistically sized case studies. As is, they already demonstrate the importance of selecting a proper scheduler class for efficient verification, and that restricted classes are useful in planning scenarios.

### **6 Conclusion**

We have shown that the various notions of information available to a scheduler class, such as history, clock order, expiration times or overall elapsed time, almost all make distinct contributions to the power of the class in SA. Our choice of notions was based on classic scheduler classes relevant for other stochastic models, previous literature on the character of nondeterminism in and verification of SA, and the need to synthesise simple schedulers in planning. Our distinguishing examples clearly expose how to exploit each notion to improve the probability of reaching a goal. For verification of SA, we have demonstrated the feasibility of lightweight scheduler sampling, where the different notions may be used to finely control the power of the lightweight schedulers. To solve stochastic timed planning problems defined via SA, our analysis helps in the case-by-case selection of an appropriate scheduler class that achieves the desired tradeoff between optimal probabilities and ease of implementation of the resulting plan.

We expect the arguments of this paper to extend to steady-state/frequency measures (by adding loops back from absorbing to initial states in our examples), and that our results for classic schedulers transfer to SA with delayable actions. We propose to use the results to develop better abstractions for SA, the next goal being a refinement algorithm for efficient lightweight verification of SA.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Symbolically Quantifying Response Time in Stochastic Models Using Moments and Semirings**

Hugo Bazille1, Eric Fabre1, and Blaise Genest2(B)

1 Univ Rennes, Inria, SUMO Team, Rennes, France 2 Univ Rennes, CNRS, IRISA, Rennes, France<EMAIL_ADDRESS>
**Abstract.** We study quantitative properties of the response time in stochastic models. For instance, we are interested in quantifying bounds such that a high percentage of the runs answers a query within these bounds. To study such problems, computing probabilities on a statespace blown-up by a factor depending on the bound could be used, but this solution is not satisfactory when the bound is large.

In this paper, we propose a new *symbolic* method to quantify bounds on the response time, using the moments of the distribution of simple stochastic systems. We prove that the distribution (and hence the bounds) is uniquely defined given its moments. We provide *optimal* bounds for the response time over all distributions having a pair of these moments. We explain how to *symbolically* compute in polynomial time any moment of the distribution of response times using adequatelydefined semirings. This allows us to compute optimal bounds in parametric models and to reduce complexity for computing optimal bounds in hierarchical models.

### **1 Introduction**

Response time has been considered lately as an important property of systems [8,15,21]. In this context, one does not simply want a query to be answered eventually, but to be answered in a reasonable amount of time. In the model-checking community, problems on response time have been studied mainly *qualitatively*, in the context of (pure, that is non stochastic) two-player games [8,21]. There, one looks for a strategy ensuring that the lim-sup of response time is finite. It ensures that under this strategy, there will be a bound on the response time to any query. This has been extended in [15] to a quantitative setting, where one wants to optimize the mean response time in a pure two-player game.

In this paper, we consider stochastic systems. In such systems, the response time is a random variable, unlikely to be bounded as even a single probabilistic loop on a reachable state will make the response time longer than T for a set of runs of small but positive probability, no matter T. Instead, we propose to quantify such response times. One way to do that is to obtain the distribution of response times. Another way is to compute, for a probability 0 <p< 1, the bound T that is satisfied (by a set of runs) with probability at least 1 − p. In this paper, we tackle both problems. For that, we use the concept of *moments* of the distribution of response times, as described next.

The *moment of order* r of a probability distribution δ over R or R+ is defined as the integral of xrδ(x) over the support of δ, when defined (that is if xrδ(x) is measurable and the integral is defined). For instance, the moment of order 1 is the expected value of δ, while the moment of order 2 allows one to compute the standard deviation of δ. Inspired by the computation of entropy for automata [10] (see also [1] for the computation of entropy for (non-Zeno) timed-automata), we design new semirings in which each moment corresponds to the sum of weights of runs reaching a state. This construction can be applied to probabilistic automata (that is, labeled discrete time Markov chains), as well as labeled *continuous time Markov chains*, where time is continuous and is drawn according to some rate. Adapting the Floyd-Warshall algorithm provides a *symbolic* way to perform the computation of the n first moments in time cubic in the number of states of the Markov Chain, and quadratic in n. For any n, we can thus compute the value of the first n moments. In some sense, we extend the approach of [12,16] from computing probabilities to computing any moments. This allows us to evaluate the distribution of response times in two ways:

Firstly, thanks to the symbolic expression of moments, we prove that there is a unique distribution having the moments of a distribution of response times of a probabilistic automaton. We can then build a sequence of distributions matching the first n moments, for instance the maximal entropy one [11]. Here, maximal entropy means assuming the least information besides these moments. This sequence of distributions is then ensured to converge in law towards the distribution of response times.

Secondly, we study optimal symbolic bounds on the time to answer a high percentage of queries, obtained from moments. The Tchebychev inequality provides optimal symbolic bounds when considering the space of distributions having one given moment, of any order i. We obtain bounds optimal in the space of distributions having two given moments, of any orders i, j. We show how this improves Tchebychev bounds on some example. Having symbolic methods allows for instance to deal with parametric systems where the parameters represent uncertain probabilities. In this case, we can compute optimal bounds satisfying all valuations of parameters. For hierarchical systems [3], which are compact representations of large systems, our symbolic method allows to design a much more efficient algorithm (e.g. it does not consider twice the same component) to compute the moments, and thus the bounds. Missing proofs can be found in [5].

*Related Work*: Response times in stochastic systems have been studied for a long time by the perf.eval. community under the name "first passage times", e.g. in [22]. Techniques used in this community to compute moments of Markov chains are mostly based on numerical methods, e.g. [13]. While [13] has the same complexity as our symbolic technique, it is very efficient on explicit models. However, these numerical methods are less adaptable than our symbolic algorithm, in particular concerning parametric or hierarchical systems.

Concerning the determinacy of the distribution given moments, it is known [20] that phase-type distributions of order n are determined by their first 2n − 1 moments. First passage distribution time in Markov chains with n states are phase type distribution of order n. However, [20] does not help characterizing bounds as it does not ensure that a non-phase type distribution cannot have the exact same moments as a phase type distribution, unlike our result.

Bounding the response time has also been studied in the perf.eval. community. Again, methods used there are mostly numerical [6,19]. In [19] (pp. 68–69), a symbolic bound is also provided in the particular case of moments of order 1, 2 and 3. In [2], it is shown how to use the two first moments of response time across various components to compute general bounds, using techniques close to ours, but restricted to moments of order 1 and 2. In our paper, we provide *optimal* bounds for any order (i, j) ∈ N2. Taking into account moments of order i, j > 3 is important when the proportion of runs to answer is close to 1.

Last, computing moments find other applications. For instance, in [4,7,14], complex functions describing the evolution of molecular species are approximated using the first k moments, for some k.

### **2 Probabilistic Automata**

We first introduce a simple class of models, namely *probabilistic automata* (also called *labeled discrete time Markov chains*), on which we can demonstrate our techniques. Later, we will extend our results to handle continuous time, considering Continuous-Time Markov Chains (CTMC), as well as parametric and hierarchical systems.

**Definition 1.** *A probabilistic automaton* A *over a finite alphabet* Σ *is a tuple* (S, P r, δ0) *where:*

*Example 1.* For instance, the model depicted on Fig. 1 is a probabilistic automaton with 3 states {1, 2, 3}. There is a transition between 1 and 2 labeled **query** with probability 1. From state 2, with probability .9 we stay in state 2 with a transition labeled **wait**, and with probability .1 we go to state 3 with a transition labeled **response**. We loop in state 3 with probability 1.

**Fig. 1.** A simple example of a query-response model

A finite sequence π = s0, a1, s1,...,an, sn ∈ (SΣ)nS is called a *finite path* starting from s0 and ending in sn, and a transition t ∈ π if t = siai+1si+1 for some i. We denote |π| = n the length of the path π. For a path π1 ending in sn and a path π2 starting from sn, we can define the concatenated path π1 · π2 where the last node of π1 and the first node of π2 are merged. A path π1 is a *prefix* of π if there exists a path π2 such that π1 · π2 = π.

For a path π starting in a state s0, we define P(π) = t∈π P r(t) the probability

that a path with prefix π is executed from s0. A path π is realizable if P(π) > 0. Let s be a state, and Π be a set of finite paths starting from s such that no path in Π is a prefix of another path in Π. Then the probability that a path starting from s has a prefix in Π is P(Π) = - P(ρ). We say that Π is

ρ∈Π *disjoint* if no path ρ of Π is a prefix of another path ρ = ρ of Π or similarly, Cyl(ρ) ∩ Cyl(ρ ) = ∅ with Cyl(ρ) = {π, ρ prefix of π}.

Some labels of an automaton will be of particular interest concerning response time.  For simplicity, we will assume that there is a unique query type ΣQ = {q} and a unique response type ΣR = {r}, with q = r. We will also assume that there is no path with two (similar) queries q. To handle cases with several query/response types, it suffices for each type to consider only queries and answers of that type and disregard other types.

**Problem Statement:** We are interested in quantifying the time between queries and responses, called the *response time*, which is a random variable. A way to quantify it is to produce the distribution of response times, either for each transition labeled by a query, or averaged on these transitions, weighted by the probability to see each of these transitions. Another way is to answer modelchecking questions such as: what is the smallest delay T such that the mass of paths unanswered after T units of time is smaller than some probability p?

To compute both the distribution and the delay T, we will use the so called *moments of the distribution of response times*. The moment of order 1 is the mean value, and the moment of order 2 allows to compute the standard deviation.

### **3** *Symbolically* **Computing Moments Using Semirings**

In this section, we define moments and explain how to compute them *symbolically* using appropriately-defined semirings.

Let X be the random variable of the response time. If all queries are answered, then X takes values in Nmax, else X takes values in Nmax ∪ {∞}. Let p(x) be the probability that the response is obtained x units of time after the query, that is, the probability that X = x. Variable p is a distribution over response time, with - x p(x) = 1.

#### **3.1 Semirings Associated with Moments**

We will compute moments of the distribution of response times by considering each query individually. We can then take e.g. the average over all queries (as we assumed that there are no two queries on the same path). Thus, we first fix a state q, target of a transition labeled by a query. State q symbolizes that a query has just been asked. We then let R be the set of target states of transitions labeled by a response. A state is in R if a response to this query has just been given. For instance, on Fig. 1, we have q = 2 and R = {3}.

We introduce a set of semirings that will allow us to compute symbolically the moment of order n of the distribution of response times to the query associated with state q, for all n ∈ N. We will compute the moment inductively on a disjoint subset - Π of paths of A from q to R. For an integer n, we denote μn(Π) = ρ∈Π P(ρ)|ρ| n. Let **Path**R q be the set of paths in the automaton A between q and

the first occurrence of R. Notice that **Path**R q is disjoint. Thus, we have that μn(**Path**R q ) is the moment of order n of the distribution of response times to the query associated with state q. To avoid some heavy notations, when R is reduced to one state t, let μn(**Path**t s) be the set of paths between s to the first occurrence of t and we denote μn(s, t) = μn(**Path**t s).

We now give some properties of μ. Let Π1 be a set of paths ending in some state s and let Π2 be a set of paths starting from s. We denote by Π1 · Π2 the set of paths ρ1ρ2 with ρ1 ∈ Π1 and ρ2 ∈ Π2.

$$\textbf{Proposition 1.}\text{ For all }n\text{, we have }\mu\_n(\boldsymbol{\Pi}\_1\cdot\boldsymbol{\Pi}\_2) = \sum\_{i=0}^n \binom{n}{i} \mu\_i(\boldsymbol{\Pi}\_1)\cdot\mu\_{n-i}(\boldsymbol{\Pi}\_2)$$

This property hints to a set of semirings (R, ⊕n, ⊗n, 0n, 1n) with good properties to compute moments. For (n + 1)-tuples (x0,...,xn) and (y0,...,yn), we define operations ⊕n and ⊗n:

$$\begin{aligned} \text{- } (x\_0, \dots, x\_n) \oplus\_n (y\_0, \dots, y\_n) &= (x\_0 + y\_0, \dots, x\_n + y\_n) \\ \text{- } (x\_0, \dots, x\_n) \otimes\_n (y\_0, \dots, y\_n) &= (z\_0, \dots, z\_n) \text{ with } z\_i = \sum\_{j=0}^i \binom{i}{j} x\_j y\_{i-j} \end{aligned}$$

The neutral element for ⊕n is 0n = (0,..., 0). 0n is an annihilator for ⊗n. The neutral element for ⊗n is 1n = (1, 0,..., 0). In the following, we will denote the different laws and elements by ⊕, ⊗, 0 and 1.

### **Proposition 2.** *For* n ≥ 0*,* (Rn+1, ⊕, ⊗, 0, 1) *defines a commutative semiring.*