We obtain the following theorem which is the basis of our algorithm.

**Theorem 4.** *For MDP* M *let* G*,* S?*, and* u *be as above. Assume* σ ∈ SM *and* k ≥ 0 *such that* σ ∈ arg maxσ-∈SMPrσ- (♦≤kG) + Prσ- (≤kS?) · u *and* Prσ s (≤kS?) < 1 *for all* s ∈ S?*. We have*

$$\begin{aligned} &\Pr^{\sigma}(\Diamond^{\leq k}G) + \Pr^{\sigma}(\Box^{\leq k}S\_{?}) \cdot \min\_{s \in S; \; 1 = \Pr^{\sigma}(\Box^{\leq k}S\_{?})} \frac{\Pr^{\sigma}(\Diamond^{\leq k}G)}{1 - \Pr^{\sigma}(\Box^{\leq k}S\_{?})} \\ &\leq \Pr^{\max}(\Diamond G) \leq \Pr^{\sigma}(\Diamond^{\leq k}G) + \Pr^{\sigma}(\Box^{\leq k}S\_{?}) \cdot u. \end{aligned}$$

Similar to the results for MCs it also holds that Prmax(♦G) ≤ maxσ∈SM uˆσ k with

$$
\hat{u}\_k^{\sigma} := \operatorname{Pr}^{\sigma}(\Diamond^{\leq k} G) + \operatorname{Pr}^{\sigma}(\Box^{\leq k} S\_?) \cdot \max\_{s \in S; \, } \frac{\operatorname{Pr}\_s^{\sigma}(\Diamond^{\leq k} G)}{1 - \operatorname{Pr}\_s^{\sigma}(\Box^{\leq k} S\_?)}.
$$

**Fig. 2.** Example MDP with corresponding step bounded probabilities.

However, this upper bound can not trivially be embedded in a value iteration based procedure. Intuitively, in order to compute the upper bound for iteration k, one can not necessarily build on the results for iteration k − 1.

*Example 6.* Consider the MDP M given in Fig. 2(a). Let G = {s3, s4} be the set of goal states. We therefore have S? = {s0, s1, s2}. In Fig. 2(b) we list step bounded probabilities with respect to the possible schedulers, where σα, σβα, and σββ refer to schedulers with σα(s0) = α and for γ ∈ {α, β}, σβγ(s0) = β and σβγ(s0βs0) = γ. Notice that the probability measures Prσ s1 and Prσ s2 are independent of the considered scheduler σ. For step bounds k ∈ {1, 2} we get

– maxσ∈SM uˆσ 1 = ˆuσα 1 =0+0.8 · max(0, 1, 0) = 0.8 and – maxσ∈SM uˆσ 2 = ˆuσββ 2 = 0.42 + 0.16 · max(0.5, 0.19, 1) = 0.5.

#### **4.2 Extending the Value Iteration Approach**

The idea of our algorithm is to compute the bounds for Prmax(♦G) as in Theorem 4 for increasing k ≥ 0. Algorithm 2 outlines the procedure. Similar to Algorithm 1 for MCs, vectors xk, yk ∈ R|S| store the step bounded probabilities Prσk s (♦≤kG) and Prσk s (≤kS?) for any s ∈ S. In addition, schedulers σk and upper bounds uk ≥ maxs∈S? Prmax s (♦G) are computed in a way that Theorem 4 is applicable.

**Lemma 3.** *After executing* k *iterations of Algorithm 2 we have for all* s ∈ S? *that* xk[s] = Prσk s (♦≤kG)*,* yk[s] = Prσk s (≤kS?)*, and* k ≤ Prmax s (♦G) ≤ uk*, where* σk ∈ arg maxσ∈SMPrσ s (♦≤kG) + Prσ s (≤kS?) · uk*.*

The lemma holds for k = 0 as x0, y0, and u0 are initialized accordingly. For k > 0 we assume that the claim holds after k − 1 iterations, i.e., for xk−1, yk−1, uk−1 and scheduler σk−1. The results of the kth iteration are obtained as follows.

The function *findAction* illustrated in Algorithm 3 determines the choices of a scheduler σk ∈ arg maxσ∈SMPrσ s (♦≤kG) + Prσ s (≤kS?) · uk−1 for s ∈ S?. The idea is to consider at state s an action σk(s) = α ∈ *Act*(s) that maximizes

$$\Pr\_s^{\sigma\_k}(\Diamond^{\leq k} G) + \Pr\_s^{\sigma\_k}(\Box^{\leq k} S\_?) \cdot u\_{k-1} = \sum\_{s' \in S} \mathbf{P}(s, \alpha, s') \cdot (x\_{k-1}[s'] + y\_{k-1}[s'] \cdot u\_{k-1}).$$

**Input :** MDP M = (S, *Act*, **P**, sI , ρ), absorbing states G ⊆ S, precision ε > 0 **Output :** r ∈ R with |r − Prmax(♦G)| < ε **1** S0 ← {s ∈ S | Prmax s (♦G)=0} **2** assert that M is contracting with respect to G ∪ S0 **3** S? ← S \ (S0 ∪ G) **4** initialize x0, y0 ∈ R|S| with x0[G] = 1, x0[S \ G] = 0, y0[S?] = 1, y0[S \ S?]=0 **5** 0 ← −∞; u0 ← +∞; d0 ← −∞ **6** k ← 0 **7 repeat 8** k ← k + 1 **9** initialize xk, yk ∈ R|S| with xk[G] = 1, xk[S0] = 0, yk[S \ S?]=0 **10** dk ← dk−1 **11 foreach** s ∈ S? **do 12** α ← *findAction*(xk−1, yk−1, s, uk−1) **13** dk ← max(dk, *decisionValue*(xk−1, yk−1, s, α)) **14** xk[s] ← s-∈S **P**(s, α, s- ) · xk−1[s- ] **15** yk[s] ← s-∈S **P**(s, α, s- ) · yk−1[s- ] **16 if** yk[s] < 1 for all s ∈ S? **then 17** k ← max(k−1, mins∈S? xk[s] 1−yk[s] ) **18** uk ← min(uk−1, max(dk, max,∈S? xk[s] 1−yk[s] )) **19 until** yk[sI ] · (uk − k) < 2 · ε **20 return** xk[sI ] + yk[sI ] · k+uk 2

#### **Algorithm 2:** Sound value iteration for MDPs

For the case where no real upper bound is known (i.e., uk−1 = ∞) we implicitly assume a sufficiently large value for uk−1 such that Prσ s (♦≤kG) becomes negligible. Upon leaving state s, σk mimics σk−1, i.e., we set σk(sαs1α1 ...sn) = σk−1(s1α1 ...sn). After executing Line 15 of Algorithm 2 we have xk[s] = Prσk s (♦≤kG) and yk[s] = Prσk s (≤kS?).

It remains to derive an upper bound uk. To ensure that Lemma 3 holds we require (i) uk ≥ maxs∈S? Prmax s (♦G) and (ii) uk ∈ Uk, where

$$U\_k = \{ u \in \mathbb{R} \mid \sigma\_k \in \operatorname\*{arg\,max}\_{\sigma \in \mathfrak{S}^{\mathcal{M}}} \text{Pr}\_s^\sigma(\diamond^{\leq k} G) + \text{Pr}\_s^\sigma(\square^{\leq k} S\_{\heartsuit}) \cdot u \text{ for all } s \in S\_{\heartsuit} \}.$$

Intuitively, the set Uk ⊆ R consists of all possible upper bounds u for which σk is still optimal. Uk ⊆ is convex as it can be represented as a conjunction of inequalities with U0 = R and u ∈ Uk if and only if u ∈ Uk−1 and for all s ∈ S? with σk(s) = α and for all β ∈ *Act*(s) \ {α}

$$\sum\_{s' \in S} \mathbf{P}(s, \alpha, s') \cdot (x\_{k-1}[s'] + y\_{k-1}[s'] \cdot u) \ge \sum\_{s' \in S} \mathbf{P}(s, \beta, s') \cdot (x\_{k-1}[s'] + y\_{k-1}[s'] \cdot u).$$

The algorithm maintains the so-called *decision value* dk which corresponds to the minimum of Uk (or −∞ if the minimum does not exist). Algorithm 4 outlines the

```
1 function findAction(x, y, s, u)
2 if u = ∞ then
3 return α ∈ arg maxα∈Act(s)

                                       s-
                                         ∈S P(s, α, s-

                                                    ) · (x[s-

                                                           ] + y[s-

                                                                  ] · u)
4 else
5 return α ∈ arg maxα∈Act(s)

                                       s-
                                         ∈S P(s, α, s-

                                                    ) · (y[s-

                                                           ])
```
**Algorithm 3:** Computation of optimal action.

```
1 function decisionValue(x, y, s, α)
2 d ← −∞
3 foreach β ∈ Act(s) \ {α} do
4 yΔ ←
                s-
                 ∈S(P(s, α, s-

                            ) − P(s, β, s-

                                       )) · y[s-

                                             ]
5 if yΔ > 0 then
6 xΔ ←
                   s-
                    ∈S(P(s, β, s-

                               ) − P(s, α, s-

                                          )) · x[s-

                                                ]
7 d ← max(d, xΔ/yΔ)
8 return d
```
**Algorithm 4:** Computation of decision value.

procedure to obtain the decision value at a given state. Our algorithm ensures that uk is only set to a value in [dk, uk−1] ⊆ Uk.

**Lemma 4.** *After executing Line 18 of Algorithm 2:* uk ≥ maxs∈S? Prmax s (♦G)*.*

To show that uk is a valid upper bound, let smax ∈ arg maxs∈S?Prmax s (♦G) and u∗ = Prmax smax (♦G). From Theorem 4, uk−1 ≥ u∗, and uk−1 ∈ Uk we get

$$\begin{split} &u^\* \leq \max\_{\sigma \in \mathfrak{S}^{\mathcal{M}}} \Pr\_{s\_{\max}}^{\sigma} (\Diamond^{\leq k} G) + \Pr\_{s\_{\max}}^{\sigma} (\Box^{\leq k} S\_{?}) \cdot u\_{k-1} \\ &= \Pr\_{s\_{\max}}^{\sigma\_k} (\Diamond^{\leq k} G) + \Pr\_{s\_{\max}}^{\sigma\_k} (\Box^{\leq k} S\_{?}) \cdot u\_{k-1} = x\_k [s\_{\max}] + y\_k [s\_{\max}] \cdot u\_{k-1} \end{split}$$

which yields a new upper bound xk[smax] + yk[smax] · uk−1 ≥ u∗. We repeat this scheme as follows. Let v0 := uk−1 and for i > 0 let vi := xk[smax]+yk[smax]·vi−1. We can show that vi−1 ∈ Uk implies vi ≥ u∗. Assuming yk[smax] < 1, the sequence v0, v1, v2,... converges to v∞ := limi→∞ vi = xk[smax] 1−yk[smax] . We distinguish three cases to show that uk = min(uk−1, max(dk, maxs∈S? xk[s] 1−yk[s])) ≥ u∗.

*Example 7.* Reconsider the MDP M from Fig. 2(a) and goal states G = {s3, s4}. The maximal reachability probability is attained for a scheduler that always chooses β at state s0, which results in Prmax(♦G)=0.5. We now illustrate how Algorithm 2 approximates this value by sketching the first two iterations. For the first iteration *findAction* yields action α at s0. We obtain:

$$\begin{aligned} x\_1[s\_0] &= 0, \; x\_1[s\_1] = 0.1, \; x\_1[s\_2] = 0.1, \; y\_1[s\_0] = 0.8, \; y\_1[s\_1] = 0.9, \; y\_1[s\_2] = 0, \\\ d\_1 &= 0.3/(0.8 - 0.4) = 0.75, \; \ell\_1 = \min(0, 1, 0) = 0, \; u\_1 = \max(0.75, 0, 1, 0) = 1. \end{aligned}$$

In the second iteration *findAction* yields again α for s0 and we get:

x2[s0]=0.08, x2[s1]=0.19, x2[s2]=0.1, y2[s0]=0.72, y2[s1]=0, y2[s2]=0, d2 = 0.75, 2 = min(0.29, 0.19, 0.1) = 0.1, u2 = max(0.75, 0.29, 0.19, 0.1) = 0.75.

Due to the decision value we do not set the upper bound u2 to 0.29 < Prmax(♦G).

**Theorem 5.** *Algorithm 2 terminates for any MDP* M*, goal states* G *and precision* ε > 0*. The returned value* r *satisfies* |r − Prmax(♦G)| ≤ ε*.*

The correctness of the algorithm follows from Theorem 4 and Lemma 3. Termination follows since M is contracting with respect to S0 ∪ G, implying limk→∞ Prσ(≤kS?) = 0. The optimizations for Algorithm 1 mentioned in Sect. 3.4 can be applied to Algorithm 2 as well.

**Fig. 3.** Comparison of sound value iteration (x-axis) and interval iteration (y-axis).

#### **5 Experimental Evaluation**

**Implementation.** We implemented sound value iteration for MCs and MDPs into the model checker Storm [8]. The implementation computes reachability probabilities and expected rewards using explicit data structures such as sparse matrices and vectors. Moreover, Multi-objective model checking is supported, where we straightforwardly extend the value iteration-based approach of [22] to sound value iteration. We also implemented the optimizations given in Sect. 3.4.

The implementation is available at www.stormchecker.org.

#### **Experimental Results.** We considered a wide range of case studies including

In total, 130 model and property instances were considered. For CTMCs and Markov automata we computed (untimed) reachability probabilities or expected rewards on the underlying MC and the underlying MDP, respectively. In all experiments the precision parameter was given by ε = 10−6.

We compare sound value iteration (SVI) with interval iteration (II) as presented in [12,13]. We consider the Gauss-Seidel variant of the approaches and compute initial bounds 0 and u0 as in [12]. For a better comparison we consider the implementation of II in Storm. [21] gives a comparison with the implementation of II in PRISM. The experiments were run on a single core (2GHz) of an HP BL685C G7 with 192GB of available memory. However, almost all experiments required less than 4GB. We measured model checking times and required iterations. All logfiles and considered benchmarks are available at [25].

Figure 3(a) depicts the model checking times for SVI (x-axis) and II (y-axis). For better readability, the benchmarks are divided into four plots with different scales. Triangles () and circles (•) indicate MC and MDP benchmarks, respectively. Similarly, Fig. 3(b) shows the required iterations of the approaches. We observe that SVI converged faster and required fewer iterations for almost all MCs and MDPs. SVI performed particularly well on the challenging instances where many iterations are required. Similar observations were made when comparing the topological variants of SVI and II. Both approaches were still competitive if no a priori bounds are given to SVI. More details are given in [21].

Figure 4 indicates the model checking times of SVI and II as well as their topological variants. For reference, we also consider standard (unsound) value iteration (VI). The x-axis depicts the number of instances that have been solved by the corresponding approach within the time limit indicated on the y-axis. Hence, a point (x, y) means that for x instances the model checking time was less or equal than y. We observe that the topological variant of SVI yielded the best run times among all sound approaches and even competes with (unsound) VI.

**Fig. 4.** Runtime comparison between different approaches.

#### **6 Conclusion**

In this paper we presented a sound variant of the value iteration algorithm which safely approximates reachability probabilities and expected rewards in MCs and MDPs. Experiments on a large set of benchmarks indicate that our approach is a reasonable alternative to the recently proposed interval iteration algorithm.

#### **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.

### **Safety-Aware Apprenticeship Learning**

Weichao Zhou and Wenchao Li(B)

Department of Electrical and Computer Engineering, Boston University, Boston, USA<EMAIL_ADDRESS>
**Abstract.** Apprenticeship learning (AL) is a kind of Learning from Demonstration techniques where the reward function of a Markov Decision Process (MDP) is unknown to the learning agent and the agent has to derive a good policy by observing an expert's demonstrations. In this paper, we study the problem of how to make AL algorithms inherently safe while still meeting its learning objective. We consider a setting where the unknown reward function is assumed to be a linear combination of a set of state features, and the safety property is specified in Probabilistic Computation Tree Logic (PCTL). By embedding probabilistic model checking inside AL, we propose a novel *counterexample-guided* approach that can ensure safety while retaining performance of the learnt policy. We demonstrate the effectiveness of our approach on several challenging AL scenarios where safety is essential.

#### **1 Introduction**

The rapid progress of artificial intelligence (AI) comes with a growing concern over its safety when deployed in real-life systems and situations. As highlighted in [3], if the objective function of an AI agent is wrongly specified, then maximizing that objective function may lead to harmful results. In addition, the objective function or the training data may focus only on accomplishing a specific task and ignore other aspects, such as safety constraints, of the environment. In this paper, we propose a novel framework that combines explicit safety specification with learning from data. We consider safety specification expressed in Probabilistic Computation Tree Logic (PCTL) and show how probabilistic model checking can be used to ensure safety and retain performance of a learning algorithm known as *apprenticeship learning* (AL).

We consider the formulation of apprenticeship learning by Abbeel and Ng [1]. The concept of AL is closely related to *reinforcement learning* (RL) where an agent learns what actions to take in an environment (known as a policy) by maximizing some notion of long-term reward. In AL, however, the agent is not given the reward function, but instead has to first estimate it from a set of expert demonstrations via a technique called *inverse reinforcement learning* [18]. The formulation assumes that the reward function is expressible as a linear combination of *known state features*. An expert demonstrates the task by maximizing this reward function and the agent tries to derive a policy that can match the feature expectations of the expert's demonstrations. Apprenticeship learning can also be

c The Author(s) 2018

H. Chockler and G. Weissenbacher (Eds.): CAV 2018, LNCS 10981, pp. 662–680, 2018. https://doi.org/10.1007/978-3-319-96145-3\_38

viewed as an instance of the class of techniques known as Learning from Demonstration (LfD). One issue with LfD is that *the expert often can only demonstrate how the task works but not how the task may fail*. This is because failure may cause irrecoverable damages to the system such as crashing a vehicle. In general, the lack of "negative examples" can cause a heavy bias in how the learning agent constructs the reward estimate. In fact, *even if all the demonstrations are safe, the agent may still end up learning an unsafe policy*.

The key idea of this paper is to incorporate formal verification in apprenticeship learning. We are inspired by the line of work on formal inductive synthesis [10] and counterexample-guided inductive synthesis [22]. Our approach is also similar in spirit to the recent work on safety-constrained reinforcement learning [11]. However, our approach uses the results of model checking in a novel way. We consider safety specification expressed in probabilistic computation tree logic (PCTL). We employ a verification-in-the-loop approach by embedding PCTL model checking as a safety checking mechanism inside the learning phase of AL. In particular, when a learnt policy does not satisfy the PCTL formula, we leverage counterexamples generated by the model checker to steer the policy search in AL. In essence, counterexample generation can be viewed as supplementing negative examples for the learner. Thus, the learner will try to find a policy that not only imitates the expert's demonstrations but also stays away from the failure scenarios as captured by the counterexamples.

In summary, we make the following contributions in this paper.

The rest of the paper is organized as follows. Section 2 reviews background information on apprenticeship learning and PCTL model checking. Section 3 defines the safety-aware apprenticeship learning problem and gives an overview of our approach. Section 4 illustrates the counterexample-guided learning framework. Section 5 describes the proposed algorithm in detail. Section 6 presents a set of experimental results demonstrating the effectiveness of our approach. Section 7 discusses related work. Section 8 concludes and offers future directions.

#### **2 Preliminaries**

#### **2.1 Markov Decision Process and Discrete-Time Markov Chain**

Markov Decision Process (MDP) is a tuple M = (S, A, T, γ, s0, R), where S is a finite set of states; A is a set of actions; T : S × A × S → [0, 1] is a transition function describing the probability of transitioning from one state s ∈ S to another state by taking action a ∈ A in state s; R : S → R is a reward function which maps each state s ∈ S to a real number indicating the reward of being in state s; s0 ∈ S is the initial state; γ ∈ [0, 1) is a discount factor which describes how future rewards attenuate when a sequence of transitions is made. A deterministic and stationary (or memoryless) policy π : S → A for an MDP M is a mapping from states to actions, i.e. the policy deterministically selects what action to take solely based on the current state. In this paper, we restrict ourselves to deterministic and stationary policy. A policy π for an MDP M induces a Discrete-Time Markov Chain (DTMC) Mπ = (S, Tπ, s0), where Tπ : S × S → [0, 1] is the probability of transitioning from a state s to another state in one step. A trajectory τ = s0 Tπ(s0,s1)>0 −−−−−−−−→ s1 Tπ(s1,s2)>0 −−−−−−−−→ s2,..., is a sequence of states where si ∈ S. The accumulated reward of τ is - ∞ i=0 γi R(si). The value function Vπ : S → R measures the expectation of accumulated reward E[ - ∞ i=0 γi R(si)] starting from a state s and following policy π. An *optimal policy* π for MDP M is a policy that maximizes the value function [4].

#### **2.2 Apprenticeship Learning via Inverse Reinforcement Learning**

*Inverse reinforcement learning (IRL)* aims at recovering the reward function R of M\R = (S, A, T, γ, s0) from a set of m trajectories ΓE = {τ0, τ1,...,τm−1} demonstrated by an expert. *Apprenticeship learning (AL)* [1] assumes that the reward function is a linear combination of state features, i.e. R(s) = ωT f(s) where f : S → [0, 1]k is a vector of known features over states S and ω ∈ Rk is an unknown weight vector that satisfies ||ω||2 ≤ 1. The expected features of a policy π are the expected values of the cumulative discounted state features f(s) by following π on M, i.e. μπ = E[ -∞ t=0 γt f(st)|π]. Let μE denote the expected features of the unknown expert's policy πE. μE can be approximated by the expected features of expert's m demonstrated trajectories ˆμE = 1 m - τ∈ΓE - ∞ t=0 γt f(st) if m

is large enough. With a slight abuse of notations, we use μΓ to also denote the expected features of a set of paths Γ. Given an error bound , a policy π∗ is defined to be *-close* to πE if its expected features μπ∗ satisfies ||μE −μπ∗ ||2 ≤ . The expected features of a policy can be calculated by using Monte Carlo method, value iteration or linear programming [1,4].

The algorithm proposed by Abbeel and Ng [1] starts with a random policy π0 and its expected features μπ0 . Assuming that in iteration i, a set of i candidate policies Π = {π0, π1,...,πi−1} and their corresponding expected features {μπ|π ∈ Π} have been found, the algorithm solves the following optimization problem.

$$\delta = \max\_{\omega} \min\_{\pi \in \Pi} \omega^T (\hat{\mu}\_E - \mu\_\pi) \qquad s.t. \ ||\omega||\_2 \le 1 \tag{1}$$

The optimal ω is used to find the corresponding optimal policy πi and the expected features μπi . If δ ≤ , then the algorithm terminates and πi is produced as the output. Otherwise, μπi is added to the set of features for the candidate policy set Π and the algorithm continues to the next iteration.

#### **2.3 PCTL Model Checking**

Probabilistic model checking can be used to verify properties of a stochastic system such as "is the probability that the agent reaches the unsafe area within 10 steps smaller than 5%?". *Probabilistic Computation Tree Logic* (PCTL) [7] allows for probabilistic quantification of properties. The syntax of PCTL includes state formulas and path formulas [13]. A state formula φ asserts property of a single state s ∈ S whereas a path formula ψ asserts property of a trajectory.

$$\phi ::= \text{true} \mid l\_i \mid \neg \phi\_i \mid \phi\_i \land \phi\_j \mid P\_{\mathsf{M}p^\*} [\psi] \tag{2}$$

$$
\psi ::= \mathbf{X}\phi \mid \phi\_1 \mathbf{U}^{\leq k} \phi\_2 \mid \phi\_1 \mathbf{U} \phi\_2 \tag{3}
$$

where li is atomic proposition and φi, φj are state formulas; - ∈ {≤, ≥, <, >}; Pp∗ [ψ] means that the probability of generating a trajectory that satisfies formula ψ is p∗. **X**φ asserts that the next state after initial state in the trajectory satisfies φ; φ1 **U**≤k φ2 asserts that φ2 is satisfied in at most k transitions and all preceding states satisfy φ1; φ1 **U** φ2 asserts that φ2 will be eventually satisfied and all preceding states satisfy φ1. The semantics of PCTL is defined by a satisfaction relation |= as follows.

$$s \mid = \text{true} \quad \text{iff state } s \in S \tag{4}$$

$$s \mid \quad \phi \quad \text{iff state s satisfies the state formula } \phi \tag{5}$$

$$
\tau \vdash \quad \psi \quad \text{iff trajectory } \tau \text{ satisfies the path formula } \psi. \tag{6}
$$

Additionally, |=min denotes the minimal satisfaction relation [6] between τ and ψ. Defining pref(τ ) as the set of all prefixes of trajectory τ including τ itself, then τ |=min ψ iff (τ |= ψ) ∧ (∀τ ∈ pref(τ )\τ,τ ψ). For instance, if ψ = φ1 **U**≤k φ2, then for any finite trajectory τ |=min φ1**U**≤kφ2, only the final state in τ satisfies φ2. Let P(τ ) be the probability of transitioning along a trajectory τ and let Γψ be the set of all finite trajectories that satisfy τ |=min ψ, the value of PCTL property ψ is defined as P=?|s0 [ψ] = - τ∈Γψ P(τ ). For a DTMC

Mπ and a state formula φ = P≤p∗ [ψ], Mπ |= φ iff P=?|s0 [ψ] ≤ p∗.

A *counterexample* of φ is a set cex ⊆ Γψ that satisfies τ∈cex P(τ ) > p∗. Let P(Γ) = - τ∈Γ P(τ ) be the sum of probabilities of all trajectories in a set Γ. Let CEXφ ⊆ 2Γψ be the set of all counterexamples for a formula φ such that (∀cex ∈ CEXφ, P(cex) > p∗) and (∀Γ ∈ 2Γψ \CEXφ, P(Γ) ≤ p∗). A *minimal counterexample* is a set cex ∈ CEXφ such that ∀cex ∈ CEXφ, |cex|≤|cex |. By converting DTMC Mπ into a weighted directed graph, counterexample can be found by solving a k-shortest paths (KSP) problem or a hop-constrained

KSP (HKSP) problem [6]. Alternatively, counterexamples can be found by using Satisfiability Modulo Theory solving or mixed integer linear programming to determine the minimal critical subsystems that capture the counterexamples in Mπ [23].

A policy can also be synthesized by solving the objective minπ P=?[ψ] for an MDP M. This problem can be solved by linear programming or policy iteration (and value iteration for step-bounded reachability) [14].

#### **3 Problem Formulation and Overview**

Suppose there are some unsafe states in an MDP\RM = (S, A, T, γ, s0). A safety issue in apprenticeship learning means that an agent following the learnt policy would have a higher probability of entering those unsafe states than it should. There are multiple reasons that can give rise to this issue. First, it is possible that the expert policy πE itself has a high probability of reaching the unsafe states. Second, human experts often tend to perform only successful demonstrations that do not highlight the unwanted situations [21]. This *lack of negative examples* in the training set can cause the learning agent to be unaware of the existence of those unsafe states.

**Fig. 1.** The 8 × 8 grid-world. (a) Lighter grid cells have higher rewards than the darker ones.  The grid cells enclosed by red lines are considered *unsafe*. (b) The blue line is an example trajectory demonstrated by the expert. (c) Only the goal states are assigned high rewards and there is little difference between the unsafe states and their nearby states. As a result, the learnt policy has a high probability of passing through the unsafe states as indicated by the cyan line. (d) p∗ = 20%. The learnt policy is optimal to a reward function that correctly assigns low rewards to the unsafe states. (Color figure online)

We use a 8 × 8 grid-world navigation example as shown in Fig. 1 to illustrate this problem. An agent starts from the upper-left corner and moves from cell to cell until it reaches the lower-right corner. The 'unsafe' cells are enclosed by the red lines. These represent regions that the agent should avoid. In each step, the agent can choose to stay in current cell or move to an adjacent cell but with 20% chance of moving randomly instead of following its decision. The goal area, the unsafe area and the reward mapping for all states are shown in Fig. 1(a). For each state s ∈ S, its feature vector consists of 4 radial basis feature functions with respect to the squared Euclidean distances between s and the 4 states with the highest or lowest rewards as shown in Fig. 1(a). In addition, a specification Φ formalized in PCTL is used to capture the safety requirement. In (7), p∗ is the required upper bound of the probability of reaching an unsafe state within t = 64 steps.

$$\Phi ::= P\_{\leq p^\*} [\texttt{true } \mathbf{U}^{\leq t} \text{ шызабе}] \tag{7}$$

Let πE be the optimal policy under the reward map shown in Fig. 1(a). The probability of entering an unsafe region within 64 steps by following πE is 24.6%. Now consider the scenario where the expert performs a number of demonstrations by following πE. *All demonstrated trajectories in this case successfully reach the goal areas without ever passing through any of the unsafe regions.* Figure 1(b) shows a representative trajectory (in blue) among 10, 000 such demonstrated trajectories. The resulting reward map by running the AL algorithm on these 10,000 demonstrations is shown in Fig. 1(c). Observe that only the goal area has been learnt whereas the agent is oblivious to the unsafe regions (treating them in the same way as other dark cells). In fact, the probability of reaching an unsafe state within 64 steps with this policy turns out to be 82.6% (thus violating the safety requirement by a large margin). To make matters worse, the value of p∗ may be decided or revised after a policy has been learnt. In those cases, even the original expert policy πE may be unsafe, e.g., when p∗ = 20%. Thus, we need to adapt the original AL algorithm so that it will take into account of such safety requirement. Figure 1(d) shows the resulting reward map learned using our proposed algorithm (to be described in detail later) for p∗ = 20%. It clearly matches well with the color differentiation in the original reward map and captures both the goal states and the unsafe regions. This policy has an unsafe probability of 19.0%. We are now ready to state our problem.

**Definition 1.** *The safety-aware apprenticeship learning (SafeAL) problem is, given an* MDP\R*, a set of* m *trajectories* {τ0, τ1,...,τm−1} *demonstrated by an expert, and a specification* Φ*, to learn a policy* π *that satisfies* Φ *and is -close to the expert policy* πE*.*

*Remark 1.* We note that a solution may not always exist for the SafeAL problem. While the decision problem of checking whether a solution exists is of theoretical interest, in this paper, we focus on tackling the problem of finding a policy π that satisfies a PCTL formula Φ (if Φ is satisfiable) and whose performance is as close to that of the expert's as possible, i.e. we relax the condition on μπ being -close to μE.

#### **4 A Framework for Safety-Aware Learning**

In this section, we describe a general framework for safety-aware learning. This novel framework utilizes information from both the expert demonstrations and a verifier. The proposed framework is illustrated in Fig. 2. Similar to the *counterexample-guided inductive synthesis* (CEGIS) paradigm [22], our framework consists of a *verifier* and a *learner*. The verifier checks if a candidate policy satisfies the safety specification Φ. In case Φ is not satisfied, the verifier generates a counterexample for Φ. The main difference from CEGIS is that our framework considers not only functional correctness, e.g., safety, but also performance (as captured by the learning objective). Starting from an initial policy π0, each time the learner learns a new policy, the verifier checks if the specification is satisfied. If true, then this policy is added to the candidate set, otherwise the verifier will generate a (minimal) counterexample and add it to the counterexample set. During the learning phase, the learner uses both the counterexample set and candidate set to find a policy that is close to the (unknown) expert policy and far away from the counterexamples. The goal is to find a policy that is -close to the expert policy and satisfies the specification. For the grid-world example introduced in Sect. 3, when p∗ = 5% (thus presenting a stricter safety requirement compared to the expert policy πE), our approach produces a policy with only 4.2% of reaching an unsafe state within 64 steps (with the correspondingly inferred reward mapping shown in Fig. 1(d)).

**Fig. 2.** Our safety-aware learning framework. Given an initial policy π0, a specification Φ and a learning objective (as captured by ), the framework iterates between a *verifier* and a *learner* to search for a policy π∗ that satisfies both Φ and . One invariant that this framework maintains is that all the πi's in the candidate policy set satisfy Φ.

Learning from a (minimal) counterexample cexπ of a policy π is similar to learning from expert demonstrations. The basic principle of the AL algorithm proposed in [1] is to find a weight vector ω under which the expected reward of πE maximally outperforms any mixture of the policies in the candidate policy set Π = {π0, π1, π2,...}. Thus, ω can be viewed as the normal vector of the hyperplane ωT (μ−μE) = 0 that has the maximal distance to the convex hull of the set {μπ | π ∈ Π} as illustrated in the 2D feature space in Fig. 3(a). It can be shown

**Fig. 3.** (a) Learn from expert. (b) Learn from both expert demonstrations and counterexamples.

that ωT μπ ≥ ωT μπ for all previously found π s. Intuitively, this helps to move the candidate μπ closer to μE. Similarly, we can apply the same max-margin separation principle to maximize the distance between the candidate policies and the counterexamples (in the μ space). Let CEX = {cex0, cex1, cex2, ...} denote the set of counterexamples of the policies that do not satisfy the specification Φ. Maximizing the distance between the convex hulls of the sets {μcex|cex ∈ CEX} and {μπ | π ∈ Π} is equivalent to maximizing the distance between the parallel supporting hyperplanes of the two convex hulls as shown in Fig. 3(b). The corresponding optimization function is given in Eq. (8).

$$\delta = \max\_{\omega} \min\_{\pi \in \Pi, \epsilon ex \in CEX} \omega^T (\mu\_{\pi} - \mu\_{\epsilon ex}) \qquad s.t. \ ||\omega||\_2 \le 1 \tag{8}$$

To attain good performance similar to that of the expert, we still want to learn from μE. Thus, the overall problem can be formulated as a multi-objective optimization problem that combines (1) and (8) into (9).

$$\max\_{\omega} \min\_{\pi \in \Pi, \tilde{\pi} \in \Pi, cex \in CEX} \left( \omega^T (\mu\_E - \mu\_\pi), \ \omega^T (\mu\_{\tilde{\pi}} - \mu\_{cex}) \right) \qquad s.t. \ ||\omega||\_2 \le 1 \tag{9}$$

#### **5 Counterexample-Guided Apprenticeship Learning**

In this section, we introduce the CounterExample Guided Apprenticeship Learning (CEGAL) algorithm to solve the SafeAL problem. It can be viewed as a special case of the safety-aware learning framework described in the previous section. In addition to combining policy verification, counterexample generation and AL, our approach uses an adaptive weighting scheme to weight the separation from μE with the separation from μcex.

$$\begin{cases} \max\_{\omega} \min\_{\pi \in \Pi\_S, \tilde{\pi} \in \Pi\_S, \operatorname{exex} \in CEX} \omega^T (k(\mu\_E - \mu\_\pi) + (1 - k)(\mu\_\tilde{\pi} - \mu\_{\text{ex}})) & (10) \\ s.t. \ ||\omega||\_2 \le 1, \ k \in [0, 1] \\ \omega^T(\mu\_E - \mu\_\pi) \le \omega^T(\mu\_E - \mu\_{\pi'}), \ \forall \pi' \in \Pi\_S \\ \omega^T(\mu\_\tilde{\pi} - \mu\_{\text{ex}}) \le \omega^T(\mu\_{\tilde{\pi}'} - \mu\_{\text{ex}}, \epsilon), \ \forall \tilde{\pi}' \in \Pi\_S, \forall \epsilon \text{ex}' \in CEX \end{cases} (10)$$

In essence, we take a weighted-sum approach for solving the multi-objective optimization problem (9). Assuming that ΠS = {π1, π2, π3,...} is a set of candidate policies that all satisfy Φ, CEX = {cex1, cex2, cex3,...} is a set of counterexamples. We introduce a parameter k and change (9) into a weighted sum optimization problem (10). Note that π and ˜π can be different. The optimal ω solved from (10) can be used to generate a new policy πω by using algorithms such as policy iteration. We use a probabilistic model checker, such as PRISM [13], to check if πω satisfies Φ. If it does, then it will be added to ΠS. Otherwise, a counterexample generator, such as COMICS [9], is used to generate a (minimal) counterexample cexπω , which will be added to CEX.

**Algorithm 1.** Counterexample-Guided Apprenticeship Learning (CEGAL)

1: **Input**: 2: M ← A partially known MDP\R; f ← A vector of feature functions 3: μE ← The expected features of expert trajectories {τ0, τ1,...,τm} 4: Φ ← Specification; ← Error bound for the expected features; 5: σ, α ∈ (0, 1) ← Error bound σ and step length α for the parameter k; 6: **Initialization**: 7: **If** ||μE − μπ0 ||2 ≤ , **then return** π0 π0 is the **initial safe policy** 8: ΠS ← {π0}, CEX ← ∅ Initialize candidate and counterexample set 9: inf ← 0, sup ← 1, k ← sup Initialize multi-optimization parameter k 10: π1 ← Policy learnt from μE via apprenticeship learning 11: **Iteration** i (i ≥ 1): 12: **Verifier:** 13: status ← Model Checker(M, πi, Φ) 14: **If** status = SAT, **then go to Learner** 15: **If** status = UNSAT 16: cexπi ← Counterexample Generator(M, πi, Φ) 17: Add cexπi to CEX and solve μcexπi , **go to Learner** 18: **Learner:** 19: **If** status = SAT 20: **If** ||μE − μπi ||2 ≤ , **then return** π∗ ← πi 21: Terminate. πi is -close to πE 22: Add πi to ΠS, inf ← k, k ← sup Update ΠS, inf and reset k 23: **If** status = UNSAT 24: **If** |k − inf| ≤ σ, **then return** π∗ ← argmin π∈ΠS ||μE − μπ||2 25: Terminate. k is too close to its lower bound. 26: k ← α · inf + (1 − α)k Decrease k to learn for safety 27: ωi+1 ← argmax ω min π∈ΠS ,π˜∈ΠS ,cex∈CEX ωT (k(μE − μπ) + (1 − k)(μπ˜ − μcex)) 28: Note that the multi-objective optimization function recovers AL when k = 1 29: πi+1, μπi+1 ← Compute the optimal policy πi+1 and its expected features μπi+1 for the MDP M with reward R(s) = ωT i+1f(s) 30: **Go to next iteration**

Algorithm 1 describes CEGAL in detail. With a constant sup = 1 and a variable inf ∈ [0, sup] for the upper and lower bounds respectively, the learner determines the value of k within [inf, sup] in each iteration depending on the outcome of the verifier and uses k in solving (10) in line 27. Like most nonlinear optimization algorithms, this algorithm requires an initial guess, which is an initial safe policy π0 to make ΠS nonempty. A good initial candidate would be the maximally safe policy for example obtained using PRISM-games [15]. Without loss of generality, we assume this policy satisfies Φ. Suppose in iteration i, an intermediate policy πi learnt by the learner in iteration i − 1 is verified to satisfy Φ, then we increase inf to inf = k and reset k to k = sup as shown in line 22. If πi does not satisfy Φ, then we reduce k to k = α · inf + (1 − α)k as shown in line 26 where α ∈ (0, 1) is a step length parameter. If |k − inf| ≤ σ and πi still does not satisfy Φ, the algorithm chooses from ΠS a best safe policy π∗ which has the smallest margin to πE as shown in line 24. If πi satisfies Φ and is -close to πE, the algorithm outputs πi as show in line 19. For the occasions when πi satisfies Φ and inf = sup = k = 1, solving (10) is equivalent to solving (1) as in the original AL algorithm.

*Remark 2.* The initial policy π0 does not have to be maximally safe, although such a policy can be used to verify if Φ is satisfiable at all. Naively safe policies often suffice for obtaining a safe and performant output at the end. Such a policy can be obtained easily in many settings, e.g., in the grid-world example one safe policy is simply staying in the initial cell. In both cases, π0 typically has very low performance since satisfying Φ is the only requirement for it.

**Theorem 1.** *Given an initial policy* π0 *that satisfies* Φ*, Algorithm 1 is guaranteed to output a policy* π∗*, such that (1)* π∗ *satisfies* Φ*, and (2) the performance of* π∗ *is at least as good as that of* π0 *when compared to* πE*, i.e.* μE − μπ∗ 2 ≤ μE − μπ0 2*.*

*Proof Sketch.* The first part of the guarantee can be proven by case splitting. Algorithm 1 outputs π∗ either when π∗ satisfies Φ and is -close to πE, or when |k − inf| ≤ σ in some iteration. In the first case, π∗ clearly satisfies Φ. In the second case, π∗ is selected from the set ΠS which contains all the policies that have been found to satisfy Φ so far, so π∗ satisfies Φ. For the second part of the guarantee, the initial policy π0 is the final output π∗ if π0 satisfies Φ and is close to πE. Otherwise, π0 is added to ΠS if it satisfies Φ. During the iteration, if |k−inf| ≤ σ in some iteration, then the final output is π∗ = argmin π∈ΠS ||μE −μπ||2, so it must satisfy μE −μπ∗ 2 ≤ μE −μπ0 2. If a learnt policy π∗ satisfies Φ and is -close to πE, then Algorithm 1 outputs π∗ without adding it to ΠS. Obviously μE − μπ 2 > , ∀π ∈ ΠS, so μE − μπ∗ 2 ≤ μE − μπ0 2.

*Discussion.* In the worst case, CEGAL will return the initial safe policy. However, this can be because a policy that simultaneously satisfies Φ and is -close to the expert's demonstrations does not exist. Comparing to AL which offers no safety guarantee and finding the maximally safe policy which has very poor performance, CEGAL provides a principled way of guaranteeing safety while retaining performance.

*Convergence.* Algorithm 1 is guaranteed to terminate. Let inft be the t th assigned value of inf. After inft is given, k is decreased from k0 = sup iteratively by ki = α · inft + (1 − α)ki−1 until either |ki − inft| ≤ σ in line 24 or a new safe policy is found in line 18. The update of k satisfies the following equality.

$$\frac{|k\_{i+1} - \imath inf\_t|}{|k\_i - \imath inf\_t|} = \frac{\alpha \cdot \imath inf\_t + (1 - \alpha)k\_i - \imath inf\_t}{k\_i - \imath inf\_t} = 1 - \alpha \tag{11}$$

Thus, it takes no more than 1 + log1−α σ sup−inft iterations for either the algorithm to terminate in line 24 or a new safe policy to be found in line 18. If a new safe policy is found in line 18, inf will be assigned in line 22 by the current value of k as inft+1 = k which obviously satisfies inft+1 −inft ≥ (1−α)σ. After the assignment of inft+1, the iterative update of k resumes. Since sup−inft ≤ 1, the following inequality holds.

$$\frac{|\inf\_{t+1} - \sup|}{|\inf\_t - \sup|} \le \frac{\sup - \inf\_t - (1 - \alpha)\sigma}{\sup - \inf\_t} \le 1 - (1 - \alpha)\sigma \tag{12}$$

Obviously, starting from an initial inf = inf0 < sup, with the alternating update of inf and k, inf will keep getting close to sup unless the algorithm terminates as in line 24 or a safe policy -close to πE is found as in line 19. The extreme case is that finally inf = sup after no more than sup−inf0 (1−α)σ updates on inf. Then, the problem becomes AL. Therefore, the worst case of this algorithm can have two phases. In the first phase, inf increases from inf = 0 to inf = sup. Between each two consecutive updates (t, t + 1) on inf, there are no more than log1−α (1−α)σ sup−inft updates on k before inf is increased from inft to inft+1. Overall, this phase takes no more than

$$\sum\_{0 \le i < \frac{\sup - \inf\_{\Omega}}{(1 - \alpha)\sigma}} \log\_{1 - \alpha} \frac{(1 - \alpha)\sigma}{\sup - \inf\_{\Omega} - i \cdot (1 - \alpha)\sigma} = \sum\_{0 \le i < \frac{1}{(1 - \alpha)\sigma}} \log\_{1 - \alpha} \frac{(1 - \alpha)\sigma}{1 - i \cdot (1 - \alpha)\sigma} \tag{13}$$

iterations to reduce the multi-objective optimization problem to original apprenticeship learning and then the second phase begins. Since k = sup, the iteration will stop immediately when an unsafe policy is learnt as in line 24. This phase will not take more iterations than original AL algorithm does to converge and the convergence result of AL is given in [1].

In each iteration, the algorithm first solves a second-order cone programming (SOCP) problem (10) to learn a policy. SOCP problems can be solved in polynomial time by interior-point (IP) methods [12]. PCTL model checking for DTMCs can be solved in time linear in the size of the formula and polynomial in the size of the state space [7]. Counterexample generation can be done either by enumerating paths using the k-shortest path algorithm or determining a critical subsystem using either a SMT formulation or mixed integer linear programming (MILP) [23]. For the k-shortest path-based algorithm, it can be computationally expensive sometimes to enumerate a large amount of paths (i.e. a large k) when p∗ is large. This can be alleviated by using a smaller p∗ during calculation, which is equivalent to considering only paths that have high probabilities.

#### **6 Experiments**

We evaluate our algorithm on three case studies: (1) grid-world, (2) cart-pole, and (3) mountain-car. The cart-pole environment1 and the mountain-car environment2 are obtained from OpenAI Gym. All experiments are carried out on a quad-core i7-7700K processor running at 3.6 GHz with 16 GB of memory. Our prototype tool was implemented in Python3. The parameters are γ = 0.99, = 10, σ = 10−5, α = 0.5 and the maximum number of iterations is 50. For the OpenAI-gym experiments, in each step, the agent sends an action to the OpenAI environment and the environment returns an observation and a reward (0 or 1). We show that our algorithm can guarantee safety while retaining the performance of the learnt policy compared with using AL alone.

#### **6.1 Grid World**

We first evaluate the scalability of our tool using the grid-world example. Table 1 shows the average runtime (per iteration) for the individual components of our tool as the size of the grid-world increases. The first and second columns indicate the size of the grid world and the resulting state space. The third column shows the average runtime that policy iteration takes to compute an optimal policy π for a known reward function. The forth column indicates the average runtime that policy iteration takes to compute the expected features μ for a known policy. The fifth column indicates the average runtime of verifying the PCTL formula using PRISM. The last column indicates the average runtime that generating a counterexample using COMICS.

**Table 1.** Average runtime per iteration in seconds.

#### **6.2 Cart-Pole from OpenAI Gym**

In the cart-pole environment as shown in Fig. 4(a), the goal is to keep the pole on a cart from falling over as long as possible by moving the cart either to the left or to the right in each time step. The maximum step length is t = 200. The

1 https://github.com/openai/gym/wiki/CartPole-v0.

2 https://github.com/openai/gym/wiki/MountainCar-v0.

3 https://github.com/zwc662/CAV2018.

position, velocity and angle of the cart and the pole are continuous values and observable, but the actual dynamics of the system are unknown4.

**Fig. 4.** (a) The cart-pole environment. (b) The cart is at −0.3 and pole angle is −20◦. (c) The cart is at 0.3 and pole angle is 20◦.

A maneuver is deemed *unsafe* if the pole angle is larger than ±20◦ while the cart's horizontal position is more than ±0.3 as shown in Fig. 4(b) and (c). We formalize the safety requirement in PCTL as (14).

$$\Phi ::= P\_{\leq p^\*} \left[ true \text{ } \mathbf{U}^{\leq t} \text{ } (angle \leq -20^\circ \land position \leq -0.3) \right]$$

$$\forall (angle \geq 20^\circ \land position \geq 0.3) \text{ } \tag{14}$$

**Table 2.** In the cart-pole environment, *higher* average steps mean better performance. The safest policy is synthesized using PRISM-games.

We used 2000 demonstrations for which the pole is held upright without violating any of the safety conditions for all 200 steps in each demonstration. The safest policy synthesized by PRISM-games is used as the initial safe policy. We also compare the different policies learned by CEGAL for different safety threshold p∗s. In Table 2, the policies are compared in terms of model checking results

4 The MDP is built from sampled data. The feature vector in each state contains 30 radial basis functions which depend on the squared Euclidean distances between current state and other 30 states which are uniformly distributed in the state space.

('MC Result') on the PCTL property in (14) using the constructed MDP, the average steps ('Avg. Steps') that a policy (executed in the OpenAI environment) can hold across 5000 rounds (the higher the better), and the number of iterations ('Num. of Iters') it takes for the algorithm to terminate (either converge to an -close policy, or terminate due to σ, or terminate after 50 iterations). The policy in the first row is the result of using AL alone, which has the best performance but also a 49.1% probability of violating the safety requirement. The safest policy as shown in the second row is always safe has almost no performance at all. This policy simply letts the pole fall and thus does not risk moving the cart out of the range [−0.3, 0.3]. On the other hand, it is clear that the policies learnt using CEGAL always satisfy the safety requirement. From p∗ = 30% to 10%, the performance of the learnt policy is comparable to that of the AL policy. However, when the safety threshold becomes very low, e.g., p∗ = 5%, the performance of the learnt policy drops significantly. This reflects the phenomenon that the tighter the safety condition is the less room for the agent to maneuver to achieve a good performance.

#### **6.3 Mountain-Car from OpenAI Gym**

Our third experiment uses the mountain-car environment from OpenAI Gym. As shown in Fig. 5(a), a car starts from the bottom of the valley and tries to reach the mountaintop on the right as quickly as possible. In each time step the car can perform one of the three actions, accelerating to the left, coasting, and accelerating to the right. The agent fails if the step length reaches the maximum (t = 66). The velocity and position of the car are continuous values and observable while the exact dynamics are unknown5. In this game setting, the car cannot reach the right mountaintop by simply accelerating to the right. It has to accumulate momentum first by moving back and forth in the valley. The safety rules we enforce are shown in Fig. 5(b). They correspond to speed limits when the car is close to the left mountaintop or to the right mountaintop (in case it is a cliff on the other side of the mountaintop). Similar to the previous experiments, we considered 2000 expert demonstrations for which all of them successfully reach the right mountaintop without violating any of the safety conditions. The average number of steps for the expert to drive the car to the right mountaintop is 40. We formalize the safety requirement in PCTL as (15).

$$\Phi ::= P\_{\leq p\*} \left[ true \text{ } \mathbf{U}^{\leq t} \text{ } (speed \leq -0.04 \land position \leq -1.1) \text{ } \right.\tag{15}$$

$$\lor (speed \geq 0.04 \land position \geq 0.5) \text{ } \right.\tag{15}$$

We compare the different policies using the same set of categories as in the cart-pole example. The numbers are averaged over 5000 runs. As shown in the

5 The MDP is built from sampled data. The feature vector for each state contains 2 exponential functions and 18 radial basis functions which respectively depend on the squared Euclidean distances between the current state and other 18 states which are uniformly distributed in the state space.

**Fig. 5.** (a) The original mountain-car environment. (b) The mountain-car environment with traffic rules: when the distance from the car to the left edge or the right edge is shorter than 0.1, the speed of the car should be lower than 0.04.

first row, the policy learnt via AL6 has the highest probability of going over the speed limits. We observed that this policy made the car speed up all the way to the left mountaintop to maximize its potential energy. The safest policy corresponds to simply staying in the bottom of the valley. The policies learnt via CEGAL for safety threshold p∗ ranging from 60% to 50% not only have lower probability of violating the speed limits but also achieve comparable performance. As the safety threshold p∗ decreases further, the agent becomes more conservative and it takes more time for the car to finish the task. For p∗ = 20%, the agent never succeeds in reaching the top within 66 steps (Table 3).

**Table 3.** In the mountain-car environment, *lower* average steps mean better performance. The safest policy is synthesized via PRISM-games.

#### **7 Related Work**

A taxonomy of AI safety problems is given in [3] where the issues of misspecified objective or reward and insufficient or poorly curated training data are highlighted. There have been several attempts to address these issues from different angles. The problem of *safe exploration* is studied in [8,17]. In particular, the latter work proposes to add a safety constraint, which is evaluated by amount

6 AL did not converge to an -close policy in 50 iterations in this case.

of damage, to the optimization problem so that the optimal policy can maximize the return without violating the limit on the expected damage. An obvious shortcoming of this approach is that actual failures will have to occur to properly assess damage.

Formal methods have been applied to the problem of AI safety. In [5], the authors propose to combine machine learning and reachability analysis for dynamical models to achieve high performance and guarantee safety. In this work, we focus on probabilistic models which are natural in many modern machine learning methods. In [20], the authors propose to use formal specification to synthesize a control policy for reinforcement learning. They consider formal specifications captured in Linear Temporal Logic, whereas we consider PCTL which matches better with the underlying probabilistic model. Recently, the problem of *safe reinforcement learning* was explored in [2] where a monitor (called shield) is used to enforce temporal logic properties either during the learning phase or execution phase of the reinforcement learning algorithm. The shield provides a list of safe actions each time the agent makes a decision so that the temporal property is preserved. In [11], the authors also propose an approach for controller synthesis in reinforcement learning. In this case, an SMT-solver is used to find a scheduler (policy) for the synchronous product of an MDP and a DTMC so that it satisfies both a probabilistic reachability property and an expected cost property. Another approach that leverages PCTL model checking is proposed in [16]. A so-called abstract Markov decision process (AMDP) model of the environment is first built and PCTL model checking is then used to check the satisfiability of safety specification. Our work is similar to these in spirit in the application of formal methods. However, while the concept of AL is closely related to reinforcement learning, an agent in the AL paradigm needs to learn a policy from demonstrations without knowing the reward function a priori.

A distinguishing characteristic of our method is the tight integration of formal verification with learning from data (apprenticeship learning in particular). Among imitation or apprenticeship learning methods, margin based algorithms [1,18,19] try to maximize the margin between the expert's policy and all learnt policies until the one with the smallest margin is produced. The apprenticeship learning algorithm proposed by Abbeel and Ng [1] was largely motivated by the support vector machine (SVM) in that features of expert demonstration is maximally separately from all features of all other candidate policies. Our algorithm makes use of this observation when using counterexamples to steer the policy search process. Recently, the idea of learning from failed demonstrations started to emerge. In [21], the authors propose an IRL algorithm that can learn from both successful and failed demonstrations. It is done by reformulating maximum entropy algorithm in [24] to find a policy that maximally deviates from the failed demonstrations while approaching the successful ones as much as possible. However, this entropy-based method requires obtaining many failed demonstrations and can be very costly in practice.

Finally, our approach is inspired by the work on formal inductive synthesis [10] and counterexample-guided inductive synthesis (CEGIS) [22]. These frameworks typically combine a constraint-based synthesizer with a verification oracle. In each iteration, the agent refines her hypothesis (i.e. generates a new candidate solution) based on counterexamples provided by the oracle. Our approach can be viewed as an extension of CEGIS where the objective is not just functional correctness but also meeting certain learning criteria.

### **8 Conclusion and Future Work**

We propose a counterexample-guided approach for combining probabilistic model checking with apprenticeship learning to ensure safety of the apprenticehsip learning outcome. Our approach makes novel use of counterexamples to steer the policy search process by reformulating the feature matching problem into a multi-objective optimization problem that additionally takes safety into account. Our experiments indicate that the proposed approach can guarantee safety and retain performance for a set of benchmarks including examples drawn from OpenAI Gym. In the future, we would like to explore other imitation or apprenticeship learning algorithms and extend our techniques to those settings.

**Acknowledgement.** This work is funded in part by the DARPA BRASS program under agreement number FA8750-16-C-0043 and NSF grant CCF-1646497.

#### **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.

### **Deciding Probabilistic Bisimilarity Distance One for Labelled Markov Chains**

Qiyi Tang(B) and Franck van Breugel

DisCoVeri Group, York University, Toronto, Canada<EMAIL_ADDRESS>
**Abstract.** Probabilistic bisimilarity is an equivalence relation that captures which states of a labelled Markov chain behave the same. Since this behavioural equivalence only identifies states that transition to states that behave exactly the same with exactly the same probability, this notion of equivalence is not robust. Probabilistic bisimilarity distances provide a quantitative generalization of probabilistic bisimilarity. The distance of states captures the similarity of their behaviour. The smaller the distance, the more alike the states behave. In particular, states are probabilistic bisimilar if and only if their distance is zero. This quantitative notion is robust in that small changes in the transition probabilities result in small changes in the distances.

During the last decade, several algorithms have been proposed to approximate and compute the probabilistic bisimilarity distances. The main result of this paper is an algorithm that decides distance one in O(n2 + m2), where n is the number of states and m is the number of transitions of the labelled Markov chain. The algorithm is the key new ingredient of our algorithm to compute the distances. The state of the art algorithm can compute distances for labelled Markov chains up to 150 states. For one such labelled Markov chain, that algorithm takes more than 49 h. In contrast, our new algorithm only takes 13 ms. Furthermore, our algorithm can compute distances for labelled Markov chains with more than 10,000 states in less than 50 min.

**Keywords:** Labelled Markov chain · Probabilistic bisimilarity Probabilistic bisimilarity distance

#### **1 Introduction**

A *behavioural equivalence* captures which states of a model give rise to the same behaviour. Bisimilarity, due to Milner [22] and Park [25], is one of the best known behavioural equivalences. Verifying that an implementation satisfies a specification boils down to checking that the model of the implementation gives rise to the same behaviour as the model of the specification, that is, the models are behavioural equivalent (see [1, Chap. 3]).

In this paper, we focus on models of probabilistic systems. These models can capture randomized algorithms, probabilistic protocols, biological systems and c The Author(s) 2018

many other systems in which probabilities play a central role. In particular, we consider *labelled Markov chains*, that is, Markov chains the states of which are labelled.

The above example shows how the behaviour of rolling a die can be mimicked by flipping a coin, an example due to Knuth and Yao [19]. Six of the states are labelled with the values of a die and the other states are labelled zero. In this example, we are interested in the labels representing the value of a die. As the reader can easily verify, the states with these labels are each reached with probability 1 6 from the initial, top most, state. In general, labels are used to identify particular states that have properties of interest. As a consequence, states with different labels are not behaviourally equivalent.

*Probabilistic bisimilarity*, due to Larsen and Skou [21], is a key behavioural equivalence for labelled Markov chains. As shown by Katoen et al. [16], minimizing a labelled Markov chain by identifying those states that are probabilistic bisimilar speeds up model checking. Probabilistic bisimilarity only identifies those states that behave exactly the same with exactly the same probability. If, for example, we replace the fair coin in the above example with a biased one, then none of the states labelled with zero in the original model with the fair coin are behaviourally equivalent to any of the states labelled with zero in the model with the biased coin. Behavioural equivalences like probabilistic bisimilarity rely on the transition probabilities and, as a result, are sensitive to minor changes of those probabilities. That is, such behavioural equivalences are not robust, as first observed by Giacalone et al. [12].

The *probabilistic bisimilarity distances* that we study in this paper were first defined by Desharnais et al. in [11]. Each pair of states of a labelled Markov chain is assigned a distance, a real number in the unit interval [0, 1]. This distance captures the similarity of the behaviour of the states. The smaller the distance, the more alike the states behave. In particular, states have distance zero if and only if they are probabilistic bisimilar. This provides a quantitative generalization of probabilistic bisimilarity that is robust in that small changes in the transition probabilities give rise to small changes in the distances. For example, we can model a biased die by using a biased coin instead of a fair coin in the above example. Let us assume that the odds of heads of the biased coin, that is, going to the left, is 51 100 . A state labelled zero in the model of the fair die has a *non-trivial* distance, that is, a distance greater than zero and smaller than one, to the corresponding state in the model of the biased die. For example, the initial states have distance about 0.036. We refer the reader to [7] for a more detailed discussion of a similar example.

As we already mentioned earlier, behavioural equivalences can be used to verify that an implementation satisfies a specification. Similarly, the distances can be used to check how similar an implementation is to a specification. We also mentioned that probabilistic bisimilarity can be used to speed up model checking. The distances can be used in a similar way, by identifying those states that behave almost the same, that is, have a small distance (see [3,23,26]).

We focus in this paper on computing the probabilistic bisimilarity distances. In particular, we present a *decision procedure* for *distance one*. That is, we compute the set of pairs of states that have distance one. Recall that distance one is the maximal distance and, therefore, captures that states behave very differently. States with different labels have distance one. However, also states with the same label can have distance one, as the next example illustrates.

Instead of computing the set of state pairs that have distance one, we compute the complement, that is, the set of state pairs with distance smaller than one. Obviously, the set of state pairs with distance zero is included in this set. First, we decide distance zero. As we mentioned earlier, distance zero coincides with probabilistic bisimilarity. The first decision procedure for probabilistic bisimilarity was provided by Baier [4]. More efficient decision procedures were subsequently proposed by Derisavi et al. [10] and also by Valmari and Franceschinis [30]. The latter two both run in O(m log n), where n and m are the number of states and transitions of the labelled Markov chain. Subsequently, we use a traversal of a directed graph derived from the labelled Markov chain. This traversal takes O(n2 + m2).

The decision procedures for distance zero and one can be used to compute or approximate probabilistic bisimilarity distances as indicated below.

Once we have computed the sets D0 and D1 of state pairs that have distance zero or one, we can easily compute the number of state pairs with non-trivial distances. If the number of non-trivial distances is small, then we can use the *simple policy iteration* (SPI) algorithm due to Bacci et al. [2] to compute those distances. Otherwise, we can either compute all distances smaller than a chosen ε > 0 or we can approximate the distances up to some chosen accuracy α > 0. In the former case, we first compute a query set Q of state pairs that contains all state pairs the distances of which are at most ε. Subsequently, we apply the *simple partial policy iteration* (SPPI) algorithm due to Bacci et al. [2] to compute the distances for all state pairs in Q. In the latter case, we start with a pair of distance functions, one being a lower-bound and the other being an upper-bound of the probabilistic bisimilarity distances, and iteratively improve the accuracy of those until they are α close. We call this new approximation algorithm *distance iteration* (DI) as it is similar in spirit to Bellman's value iteration [5].

Chen et al. [8] presented an algorithm to compute the distances by means of Khachiyan's ellipsoid method [17]. Though the algorithm is polynomial time, in practice it is not as efficient as the policy iteration algorithms (see the examples in [28, Sect. 8]). The state of the art algorithm to compute the probabilistic bisimilarity distances consists of two components: D0 and SPI. To compare this algorithm with our new algorithm consisting of the components D0, D1 and SPI, we implemented all the components in Java and ran both implementations on several labelled Markov chains. These labelled Markov chains model randomized algorithms and probabilistic protocols that are part of the distribution of probabilistic model checkers such as PRISM [20]. Whereas the original state of the art algorithm can handle labelled Markov chains with up to 150 states, our new algorithm can handle more than 10,000 states. Furthermore, for one such labelled Markov chain with 150 states, the original algorithm takes more than 49 h, whereas our new algorithm takes only 13 ms. Also, the new algorithm consisting of the components D0, D1, Q and SPPI to compute only small distances along with the new algorithm consisting of the components D0, D1 and DI to approximate the distances give rise to even less execution times for a number of the labelled Markov chains.

The main contributions of this paper are

Furthermore, by means of experiments we have shown that these three new algorithms are very effective, improving significantly on the state of the art.

#### **2 Labelled Markov Chains and Probabilistic Bisimilarity Distances**

We start by reviewing the model of interest, labelled Markov chains, its most well known behavioural equivalence, probabilistic bisimilarity due to Larsen and Skou [21], and the probabilistic bisimilarity pseudometric due to Desharnais et al. [11]. We denote the set of rational probability distributions on a set S by Distr(S). For μ ∈ Distr(S), its support is defined by support(μ) = { s ∈ S | μ(s) > 0 }. Instead of S × S, we often write S2.

**Definition 1.** *A labelled Markov chain is a tuple* S, L, τ, *consisting of*

For the remainder of this section, we fix such a labelled Markov chain S, L, τ, .

**Definition 2.** *Let* μ*,* ν ∈ Distr(S)*. The set* Ω(μ, ν) *of couplings of* μ *and* ν *is defined by*

$$\Omega(\mu, \nu) = \left\{ \omega \in \text{Distr}(S^2) \; \middle| \; \begin{array}{l} \forall s \in S: \sum\_{t \in S} \omega(s, t) = \mu(s) \land \\ \forall t \in S: \sum\_{s \in S} \omega(s, t) = \nu(t) \end{array} \right\}.$$

Note that ω ∈ Ω(μ, ν) is a joint probability distribution with marginals μ and ν. The following proposition will be used to prove Proposition 5.

**Proposition 1.** *For all* μ*,* ν ∈ Distr(S) *and* X ⊆ S2*,*

∀ω ∈ Ω(μ, ν) : support(ω) ⊆ X *if and only if* support(μ) × support(ν) ⊆ X.

**Definition 3.** *An equivalence relation* R ⊆ S2 *is a probabilistic bisimulation if for all* (s, t) ∈ R*,* (s) = (t) *and there exists* ω ∈ Ω(τ (s), τ (t)) *such that* support(ω) ⊆ R*. Probabilistic bisimilarity, denoted* ∼*, is the largest probabilistic bisimulation.*

The probabilistic bisimilarity pseudometric of Desharnais et al. [11] maps each pair of states of a labelled Markov chain to a distance, an element of the unit interval [0, 1]. Hence, the pseudometric is a function from S2 to [0, 1], that is, an element of [0, 1]S2 . As we will discuss below, it can be defined as a fixed point of the following function.

**Definition 4.** *The function* Δ : [0, 1]S2 → [0, 1]S2 *is defined by*

$$\Delta(d)(s,t) = \begin{cases} 1 & \text{if } \ell(s) \neq \ell(t) \\ \min\_{\omega \in \Omega(\tau(s), \tau(t))} \sum\_{u,v \in S} \omega(u,v) \, d(u,v) \, \text{otherwise} \end{cases}$$

Since a concave function on a convex polytope attains its minimum (see [18, p. 260]), the above minimum exists. We will use this fact in Proposition 4, one of the key technical results in this paper. We endow the set [0, 1]S2 of functions from S2 to [0, 1] with the following partial order: d e if d(s, t) ≤ e(s, t) for all s, t ∈ S. The set [0, 1]S2 together with the order form a complete lattice (see [9, Chap. 2]). The function Δ is monotone (see [6, Sect. 3]). According to the Knaster-Tarski fixed point theorem [29, Theorem 1], a monotone function on a complete lattice has a least fixed point. Hence, Δ has a least fixed point, which we denote by *μ*(Δ). This fixed point assigns to each pair of states their probabilistic bisimilarity distance.

Given that *μ*(Δ) captures the probabilistic bisimilarity distances, we define the following sets.

$$\begin{aligned} D\_0 &= \{ (s, t) \in S^2 \mid \mu(\Delta)(s, t) = 0 \} \\ D\_1 &= \{ (s, t) \in S^2 \mid \mu(\Delta)(s, t) = 1 \} \end{aligned}$$

The probabilistic bisimilarity pseudometric *μ*(Δ) provides a quantitative generalization of probabilistic bisimilarity as captured by the following result by Desharnais et al. [11, Theorem 1].

**Theorem 1.** D0 = { (s, t) ∈ S2 | s ∼ t }*.*

#### **3 Distance One**

We concluded the previous section with the characterization of D0 as the set of state pairs that are probabilistic bisimilar. In this section we present a characterization of D1 as a fixed point of the function introduced in Definition 5.

Let us consider the case that the probabilistic bisimilarity distance of states s and t is one, that is, *μ*(Δ)(s, t) = 1. Then Δ(*μ*(Δ))(s, t) = 1. From the definition of Δ, we can conclude that either (s) = (t), or for all couplings ω ∈ Ω(τ (s), τ (t)) we have support(ω) ⊆ D1.

We partition the set S2 of state pairs into

$$\begin{aligned} S\_0^2 &= \{ (s, t) \in S^2 \mid s \sim t \} \\ S\_1^2 &= \{ (s, t) \in S^2 \mid \ell(s) \neq \ell(t) \} \\ S\_?^2 &= S^2 \backslash (S\_0^2 \cup S\_1^2) \end{aligned}$$

Hence, if *μ*(Δ)(s, t) = 1, then either (s, t) ∈ S2 1 , or (s, t) ∈ S2 ? and for all couplings ω ∈ Ω(τ (s), τ (t)) we have support(ω) ⊆ D1. This leads us to the following function.

**Definition 5.** *The function* Γ : 2S2 → 2S2 *is defined by*

$$\Gamma(X) = S\_1^2 \cup \{ (s, t) \in S\_?^2 \mid \forall \omega \in \Omega(\tau(s), \tau(t)): \text{support}(\omega) \subseteq X \}. $$

**Proposition 2.** *The function* Γ *is monotone.*

Since the set 2S2 of subsets of S2 endowed with the order ⊆ is a complete lattice (see [9, Example 2.6(2)]) and the function Γ is monotone, we can conclude from the Knaster-Tarski fixed point theorem that Γ has a greatest fixed point, which we denote by *ν*(Γ). Next, we show that D1 is a fixed point of Γ.

#### **Proposition 3.** D1 = Γ(D1)*.*

Since we have already seen that D1 is a fixed point of Γ, we have that D1 ⊆ *ν*(Γ). To conclude that D1 is the greatest fixed point of Γ, it remains to show that *ν*(Γ) ⊆ D1, which is equivalent to the following.

**Proposition 4.** *ν*(Γ) \ D1 = ∅*.*

*Proof.* Towards a contradiction, assume that *ν*(Γ) \ D1 = ∅. Let

Since *ν*(Γ) \ D1 = ∅, we have that M = ∅. Furthermore,

$$M \subseteq \nu(\Gamma) \backslash D\_1. \tag{1}$$

Since *ν*(Γ) \ D1 ⊆ *ν*(Γ), we have

$$M \subseteq \nu(\varGamma) = \varGamma(\nu(\varGamma)) \subseteq S\_1^2 \cup S\_?^2. \tag{2}$$

For all (s, t) ∈ M,

$$\begin{aligned} (s,t) &\in \nu(\varGamma) \land (s,t) \notin D\_1 \quad [(1)] \\ &\Rightarrow (s,t) \in \Gamma(\nu(\varGamma)) \land (s,t) \notin S\_1^2 \\ &\Rightarrow \forall \omega \in \Omega(\tau(s), \tau(t)) : \operatorname{support}(\omega) \subseteq \nu(\varGamma). \end{aligned}$$

For each (s, t) ∈ M, let

$$\omega\_{s,t} = \underset{\omega \in \Omega(\tau(s), \tau(t))}{\text{argmin}} \sum\_{u,v \in S} \omega(u,v) \,\mu(\Delta)(u,v). \tag{4}$$

We distinguish the following two cases.

– Assume that there exists (s, t) ∈ M such that support(ωs,t) ∩ D1 = ∅. Let

$$p = \sum\_{(u,v)\in\nu(\varGamma)\cap D\_1} \omega\_{s,t}(u,v).$$

By (3), we have that support(ωs,t) ⊆ *ν*(Γ). Since support(ωs,t) ∩ D1 = ∅ by assumption, we can conclude that p > 0. Again using the fact that support(ωs,t) ⊆ *ν*(Γ), we have that

$$\sum\_{(u,v)\in\nu(I)\backslash D\_1} \omega\_{s,t}(u,v) = 1 - p.s.\tag{5}$$

Furthermore,

$$\begin{split} m &= \mu(\Delta)(s,t) \\ &= \Delta(\mu(\Delta))(s,t) \\ &= \min\_{\omega \in \Omega(\Gamma(s),\tau(t))} \sum\_{u,v \in S} \omega(u,v)\,\mu(\Delta)(u,v) \\ &= \sum\_{u,v \in S} \omega\_{s,t}(u,v)\,\mu(\Delta)(u,v) \quad [(4)] \\ &= \sum\_{(u,v) \in \nu(\varGamma)} \omega\_{s,t}(u,v)\,\mu(\Delta)(u,v) \quad [(3)] \\ &= \sum\_{(u,v) \in \nu(\varGamma) \cap D\_{1}} \omega\_{s,t}(u,v)\,\mu(\Delta)(u,v) + \sum\_{(u,v) \in \nu(\varGamma) \backslash D\_{1}} \omega\_{s,t}(u,v)\,\mu(\Delta)(u,v) \\ &= p + \sum\_{(u,v) \in \nu(\varGamma) \backslash D\_{1}} \omega\_{s,t}(u,v)\,\mu(\Delta)(u,v) \\ &\ge p + (1-p)m. \end{split}$$

The last step follows from (5) and the fact that *μ*(Δ)(u, v) ≥ m for all (u, v) ∈ *ν*(Γ) \ D1. From the facts that p > 0 and m ≥ p + (1 − p)m we can conclude that m ≥ 1. This contradicts (1).

– Otherwise, support(ωs,t) ∩ D1 = ∅ for all (s, t) ∈ M. Next, we will show that M is a probabilistic bisimulation under this assumption. From the fact that M is a probabilistic bisimulation, we can conclude from Theorem 1 that *μ*(Δ)(s, t) = 0 for all (s, t) ∈ M. Hence, since M = ∅ we have that M ∩S2 0 = ∅ which contradicts (2).

Next, we prove that M is a probabilistic bisimulation. Let (s, t) ∈ M. Since M ⊆ *ν*(Γ) \ D1 by (1), we have that (s, t) ∈ D1 and, hence, Δ(*μ*(Δ))(s, t) = *μ*(Δ)(s, t) < 1. From the definition of Δ, we can conclude that (s) = (t). Since

$$\begin{aligned} m &= \mu(\Delta)(s, t) \\ &= \sum\_{(u, v) \in \nu(I) \backslash D\_1} \omega\_{s, t}(u, v) \,\mu(\Delta)(u, v) \quad \text{[as above]} \end{aligned}$$

and *μ*(Δ)(u, v) ≥ m for all (u, v) ∈ *ν*(Γ) \ D1, we can conclude that *μ*(Δ)(u, v) = m for all (u, v) ∈ support(ωs,t). Hence, support(ωs,t) ⊆ M. Therefore, M is a probabilistic bisimulation.

#### **Theorem 2.** D1 = *ν*(Γ)*.*

*Proof.* Immediate consequence of Proposition 3 and 4.

We have shown that D1 can be characterized as the greatest fixed point of Γ. Next, we will show that D1 can be decided in polynomial time.

**Theorem 3.** *Distance one can be decided in* O(n2 + m2)*.*

*Proof.* As we will show in Theorem 5, distance smaller than one can be decided in O(n2 + m2). Hence, distance one can be decided in O(n2 + m2) as well.

#### **4 Distance Smaller Than One**

To compute the set of state pairs which have distance one, we can first compute the set of state pairs which have distance less than one. The latter set we denote by D<1. We can then obtain D1 by taking the complement of D<1. As we will discuss below, D<1 can be characterized as the least fixed point of the following function.

**Definition 6.** *The function* Γ: 2S2 → 2S2 *is defined by*

$$
\Pi(X) = S^2 \backslash \Gamma(S^2 \backslash X).
$$

The next theorem follows from Theorem 2.

#### **Theorem 4.** D<1 = *μ*( Γ)*.*

Next, we show that the computation of D<1 can be formulated as a reachability problem on a directed graph which is induced by the labelled Markov chain. Thus, we can use standard search algorithms, for example, breadth-first search, on the induced graph.

Next, we present the graph induced by the labelled Markov chain.

**Definition 7.** *The directed graph* G = (V,E) *is defined by*

$$\begin{array}{l} V = S\_0^2 \cup S\_? \\ E = \left\{ \langle (u,v), (s,t) \rangle \mid \tau(s)(u) > 0 \land \tau(t)(v) > 0 \right\} \end{array}$$

We are left to show that in the graph G defined above, a vertex (s, t) is reachable from some vertex in S2 0 if and only if the state pair (s, t) in the labelled Markov chain has distance less than one.

As we have discussed earlier, if a state pair (s, t) has distance one, either s and t have different labels, or for all couplings ω ∈ Ω(τ (s), τ (t)) we have that support(ω) ⊆ D1. To avoid the universal quantification over couplings, we will use Proposition 1 in the proof of following proposition.

**Proposition 5.** *μ*( Γ) = { (s, t) | (s, t) *is reachable from some* (u, v) ∈ S2 0 }*.*

**Theorem 5.** *Distance smaller than one can be decided in* O(n2 + m2)*.*

*Proof.* Distance smaller than one can be decided as follows.

By Theorem 4 and Proposition 5, we have that s and t have distance smaller than one if and only if (s, t) is reachable in the directed graph G from some (u, v) such that u and v have distance zero. These reachable state pairs can be computed using breadth-first search, with the queue initially containing S2 0 .

Distance zero, that is, probabilistic bisimilarity, can be decided in O(m log n) as shown by Derisavi et al. in [10]. The directed graph G has n2 vertices and m2 edges. Hence, breadth-first search takes O(n2 + m2).

#### **5 Number of Non-trivial Distances**

As we have already discussed earlier, distance zero captures that states behave exactly the same, that is, they are probabilistic bisimilar, and distance one indicates that states behave very differently. The remaining distances, that is, those greater than zero and smaller than one, we call non-trivial. Being able to determine quickly the number of non-trivial distances of a labelled Markov chain allows us to decide whether computing all these non-trivial distances (using some policy iteration algorithm) is feasible.

To determine the number of non-trivial distances of a labelled Markov chain, we use the following algorithm.

As first proved by Baier [4], distance zero, that is, probabilistic bisimilarity, can be decided in polynomial time. As we proved in Theorem 3, distance one can be decided in polynomial time as well. Hence, we can compute the number of non-trivial distances in polynomial time.

To decide distance zero, we implemented the algorithm to decide probabilistic bisimilarity due to Derisavi et al. [10] in Java. We also implemented our algorithm to decide distance one, described in the proof of Theorems 3 and 5.

We applied our implementation to labelled Markov chains that model randomized algorithms and probabilistic protocols. These labelled Markov chains have been obtained from the verification tool PRISM [20]. We compute the number of non-trivial distances for two models: the randomized self-stabilising algorithm due to Herman [14] and the bounded retransmission protocol by Helmink et al. [13].

For the randomized self-stabilising algorithm, the size of the labelled Markov chain grows exponentially in the numbers of processes, N. The results for the randomized self-stabilising algorithm are shown in the table below. As we can see from the table, for systems up to 128 states, the algorithm runs for less than a second. For the system with 512 states, the algorithm terminates within seven minutes. For the case N = 3, there are only 12 non-trivial distances. The size is so small that we can easily compute all the non-trivial distances. Section 6 will use the simple policy iteration algorithm as the next step to compute them. The same applies to the case N = 5. For N = 7 or 9, the number of non-trivial distances is around 11,000 and 200,000, respectively. This makes computing all of them infeasible. Thus, instead of computing all of them, we need to find alternative ways to handle systems with a large number of non-trivial distances. We will discuss two alternative ways in Sects. 7 and 8. Moreover, in this example, as |D1| = |S2 1 |, we know that all the state pairs with distance one are those that have different labels.

In the bounded retransmission protocol, there are two parameters: N denotes the number of chunks and M the maximum allowed number of retransmissions of each chunk. The results are shown in the table below. The algorithm can handle systems up to 3,526 states within 11 min. In this example, there are no non-trivial distances. As a consequence, deciding distance zero and one suffices to compute all the distances in this case.

#### **6 All Distances**

To compute all distances of a labelled Markov chain, we augment the existing state of the art algorithm, which is based on algorithms due to Derisavi et al. [10] (step 1) and Bacci et al. [2] (step 3), by incorporating our decision procedure (step 2) as follows.

Given that we not only decide distance zero, but also distance one, before running simple policy iteration, the correctness of the simple policy iteration algorithm in the augmented setting needs an adjusted proof.

As we already discussed in the previous section, step 1 and 2 are polynomial time. However, step 3 may take at least exponential time in the worst case, as we have shown in [27]. Hence, the overall algorithm is exponential time.

The first example we consider here is the synchronous leader election protocol of Itai and Rodeh [15] which is taken from PRISM. The protocol takes the number of processors, N, and a constant K as parameters. We compare the running time of our new algorithm with the state of the art algorithm, that combines algorithms due to Derisavi et al. and due to Bacci et al. The results are shown in the table below. In this protocol, the number of non-trivial distances is zero. Thus, our new algorithm terminates without running step 3 which is the simple policy iteration algorithm. On the other hand, the original simple policy iteration algorithm computes the distances of all the elements in the set D1 \S2 1 , the size of which is huge as can be seen from the last two columns of the table.

The simple policy iteration algorithm can only handle a limited number of states. For the labelled Markov chain with 26 states (N = 3 and K = 2) the simple policy iteration algorithm takes four seconds, while our new algorithm takes one millisecond. The speed-up is more than 4,000 times. For the labelled Markov chain with 61 states (N = 4 and K = 2), the simple policy iteration algorithm runs in 812 s, while our new algorithm takes three milliseconds. The speed-up of the new algorithm is 30,000 times. The biggest system the simple policy iteration algorithm can handle is the one with 147 states (N = 3 and K = 4) and it takes more than 49 h. In contrast, our new algorithm terminates within 13 ms. That makes the new algorithm seven orders of magnitude faster than the state of the art algorithm. This example also shows that the new algorithm can handle systems with at least 12,400 states.

In the second example, we model two dies, one using a fair coin and the other one using a biased coin. The goal is to compute the probabilistic bisimilarity distance between these two dies. An implementation of the die algorithm is part of PRISM. The resulting labelled Markov chain has 20 states.

As there are only 30 non-trivial distances, we run the simple policy iteration algorithm as step 3. The new algorithm is about 46 times faster than the original algorithm.

#### **7 Small Distances**

As we have discussed in Sect. 5, for systems of which the number of non-trivial distances is so large that computing all of them is infeasible, we have to find alternative ways. In practice, as we only identify the state pairs with small distances, we can cut down the number of non-trivial distances by only computing those with small distances.

To compute the non-trivial distances smaller than a positive number, ε, we use the following algorithm.

$$Q = \{ (s, t) \in S^2 \mid (D\_0 \cup D\_1) \mid \Delta(d)(s, t) \le \varepsilon \} $$

where

$$d(s,t) = \begin{cases} 1 \text{ if } (s,t) \in D\_1 \\ 0 \text{ otherwise} \end{cases}$$

#### 4. Simple partial policy iteration for Q.

The first two steps remain the same. In step 3, we compute a query set Q that contains all state pairs with distances no greater than ε, as shown in Proposition 6. In step 4, we use this set as the query set to run the simple partial policy iteration algorithm by Bacci et al. [2].

**Proposition 6.** *Let* d *be the distance function defined in step 3. For all* (s, t) ∈ S2 \ (D0 ∪ D1)*, if μ*(Δ)(s, t) ≤ ε*, then* Δ(d)(s, t) ≤ ε*.*

Given that we not only decide distance zero, but also distance one, before running simple partial policy iteration, the correctness of the simple partial policy iteration algorithm in the augmented setting needs an adjusted proof.

As we have seen before, step 1 and 2 take polynomial time. In step 3, computing Δ(d) corresponds to solving a minimum cost network flow problem. Such a problem can be solved in polynomial time using, for example, Orlin's network simplex algorithm [24]. As we have shown in [28], step 4 takes at least exponential time in the worst case. Therefore, the overall algorithm is exponential time.

We consider the randomized quicksort algorithm, an implementation of which is part of jpf-probabilistic [31]. The input of the algorithm is the list to be sorted. The list of size 6 gives rise to a labelled Markov chain with 82 states. We compare the running time of the new algorithm for small distances (D0+D1+Q + SPPI) to the original algorithm (D0 + SPI) and the new algorithm presented in Sect. 6 (D0 + D1 + SPI). The original algorithm (D0 + SPI) takes about 14 h, the new algorithm which incorporates the decision procedure of distance one takes less than 7 h. For ε = 0.1, the new algorithm for small distances takes 57 min. This makes it about 7 times faster than the algorithm presented in Sect. 6 and about 15 times faster than the original simple policy iteration algorithm. For ε = 0.01, the new algorithm for small distances takes even less time, namely 41 min. As can be seen in the table below, the total number of non-trivial distances is 2,300. The simple partial policy iteration algorithm starts with the query set Q but may have to compute the distances of other state pairs as well. The total number of state pairs considered by the simple partial policy iteration algorithm can be found in the column labelled Total.

#### **8 Approximation Algorithm**

We propose another solution to deal with a large number of non-trivial distances by approximating the distances rather than computing the exact values. To approximate the distances such that the approximate values differ from the exact ones by at most α, a positive number, we use the following algorithm.

$$\begin{aligned} \text{3. } \begin{cases} \text{3. } l(s,t) = \begin{cases} 1 \text{ if } (s,t) \in D\_1 \\ 0 \text{ otherwise} \end{cases} \\ u(s,t) = \begin{cases} 0 \text{ if } (s,t) \in D\_0 \\ 1 \text{ otherwise} \end{cases} \\ \text{re-peak} \\ \begin{cases} \text{or} \\ \text{for} \\ \end{cases} \text{each } (s,t) \in S^2 \\ \text{if} \quad l(s,t) \neq u(s,t) \\ l(s,t) = \Delta(l)(s,t) \\ u(s,t) = \Delta(u)(s,t) \\ \text{until } ||l-u|| \le \alpha \end{cases} \end{aligned}$$

Again, the first two steps remain the same. Step 3 contains the new approximation algorithm called *distance iteration* (DI). In this step, we define two distance functions, a lower-bound l and an upper-bound u. We repeatedly apply Δ to these two functions until the difference of the non-trivial distances in these two functions is smaller than the threshold α. For each state pair we end up with an interval of at most size α in which their distance lies. To prove the algorithm correct, we modify the function Δ defining the probabilistic bisimilarity distances slightly as follows.

**Definition 8.** *The function* Δ0 : [0, 1]S2 → [0, 1]S2 *is defined by*

$$\Delta\_0(d)(s,t) = \begin{cases} 0 & if \,(s,t) \in D\_0\\ \Delta(d)(s,t) \,\,otherwise \end{cases}$$

Some properties of Δ0, which are key to the correctness proof of the above algorithm, are collected in the following theorem.

#### **Theorem 6.**

*(a) The function* Δ0 *is monotone. (b) The function* Δ0 *is nonexpansive. (c) μ*(Δ0) = *μ*(Δ)*. (d) μ*(Δ0) = *ν*(Δ0)*. (e) μ*(Δ0) = supm∈N Δm 0 (d0)*, where* d0(s, t)=0 *for all* s, t ∈ S*. (f ) ν*(Δ0) = infn∈N Δn 0 (d1)*, where* d1(s, t)=1 *for all* s, t ∈ S*.*

Let us use randomized quicksort introduced in Sect. 7 and the randomized self-stabilising algorithm due to Herman [14] introduced in Sect. 5 as examples. Recall that for the randomized self-stabilising algorithm, when N = 7, the number of non-trivial distances is 11,032, which we are not able to handle using the simple policy iteration algorithm. We apply the approximation algorithm to this model and the randomized quicksort example with 82 states and present the results below. The accuracy α is set to be 0.01.

The approximation algorithm for randomized quicksort runs for about 14 min, which is about 3 to 4 times faster than the algorithm for small distances in Sect. 7. For the randomized self-stabilising algorithm with 128 states, the approximation algorithm terminates in about 54 h. Although the number of non-trivial distances for the randomized self-stabilising algorithm is about 5 times of that of the randomized quicksort, the running time is more than 200 times slower. It is unknown whether this approximation algorithm has exponential running time.

#### **9 Conclusion**

In this paper, we have presented a decision procedure for probabilistic bisimilarity distance one. This decision procedure provides the basis for three new algorithms to compute and approximate the probabilistic bisimilarity distances of a labelled Markov chain. The first algorithm decides distance zero, then decides distance one, and finally uses simple policy iteration to compute the remaining distances. As shown experimentally, this new algorithm significantly improves the state of the art algorithm that only decides distance zero and then uses simple policy iteration. The second algorithm computes all probabilistic bisimilarity distances that are smaller than some given upper bound, by deciding distance zero, deciding distance one, computing a query set, and running simple partial policy iteration for that query set. This second algorithm can handle labelled Markov chains that have considerably more non-trivial distances than our first algorithm. The third algorithm approximates the probabilistic bisimilarity distances up to a given accuracy, deciding distance zero, deciding distance one and running distance iteration. Also this third algorithm can handle labelled Markov chains that have considerably more non-trivial distances than our first algorithm. Whereas we know that the first two algorithms take at least exponential time in the worst case, the analysis of the running time of the third algorithm has not yet been determined. Moreover, if we are only interested in the probabilistic bisimilarity distances for a few state pairs, with pre-computation of distance zero and one we can exclude the state pairs with trivial distances. We can add the remaining state pairs to a query set and run simple partial policy iteration to get the distances. Alternatively, we can modify the distance iteration algorithm to approximate the distances for the predefined state pairs. The details of these new algorithms will be studied in the future.

**Acknowledgements.** The authors would like to thank Daniela Petrisan, Eric Ruppert and Dana Scott for discussions related to this research. The authors are also grateful to the referees for their constructive feedback.

#### **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.

#### Author Index

Abate, Alessandro I-270 Akshay, S. I-251 Albarghouthi, Aws I-327 Albert, Elvira II-392 Anderson, Greg I-407 Argyros, George I-427 Arndt, Hannah II-3

Backes, John II-20 Bansal, Suguman I-367, II-99 Bardin, Sébastien II-294 Barrett, Clark II-236 Bartocci, Ezio I-449, I-547 Bauer, Matthew S. II-117 Becchi, Anna I-230 Berzish, Murphy II-45 Biere, Armin I-587 Bloem, Roderick I-547 Blondin, Michael I-604 Blotsky, Dmitry II-45 Bonichon, Richard II-294 Bønneland, Frederik M.  I-307

Gacek, Andrew II-20 Ganesh, Vijay II-45, II-275 Gao, Pengfei II-157 Gao, Sicun II-219 Ghassabani, Elaheh II-20 Giacobazzi, Roberto II-75 Giacobbe, Mirco I-468 Goel, Shubham I-251 Gómez-Zamalloa, Miguel II-392 Goubault, Eric II-523 Grishchenko, Ilya I-51 Gu, Ronghui II-317 Gupta, Aarti I-124, I-164, II-136

Hahn, Christopher I-144, I-289 Hassan, Mostafa II-12 He, Jinlong II-487 Henzinger, Monika II-178 Henzinger, Thomas A.  I-191, I-407 Meggendorfer, Tobias I-567 Mertens, Eric II-430 Meyer, Philipp J. I-578 Mitra, Sayan I-347 Mora, Federico II-45 Mrazek, Vojtech I-612 Mullen, Eric II-430 Müller, Peter I-596, II-12, II-55 Münger, Severin II-55 Muñiz, Marco I-527 Mutluergil, Suha Orhun II-336

Namjoshi, Kedar S. I-367 Nguyen, Huyen T. T. II-354 Nickovic, Dejan I-547

Noll, Thomas II-3, II-447 Oraee, Simin II-178 Petrucci, Laure II-354 Pick, Lauren I-164 Pike, Lee II-413 Polgreen, Elizabeth I-270 Potet, Marie-Laure II-294 Prasad Sistla, A. II-117 Preiner, Mathias I-587, II-236 Pu, Geguang II-37 Püschel, Markus I-211 Putot, Sylvie II-523

Niemetz, Aina I-587, II-236

Qadeer, Shaz I-79, II-372 Quatmann, Tim I-643

Rabe, Markus N. II-256 Rakotonirina, Itsaka II-28 Ranzato, Francesco II-75 Rasmussen, Cameron II-256 Reynolds, Andrew II-236 Robere, Robert II-275 Rodríguez, César II-354 Roeck, Franz I-547 Rozier, Kristin Yvonne II-37 Rubio, Albert II-392

Sa'ar, Yaniv I-367 Sahlmann, Lorenz II-523 Satake, Yuki I-105 Schemmel, Daniel II-447 Schneidewind, Clara I-51 Schrammel, Peter I-183 Sekanina, Lukas I-612 Seshia, Sanjit A.