In each state of the CSG, each player chooses between enabled commands of the corresponding modules; if no command is enabled, the player idles. The players move simultaneously so transitions are labelled with *lists* of action labels [a1,..., a*n*]. So the guarded command notation is extended accordingly: note how the channel's behaviour depends on which actions the two users take (the same principle applies when specifying reward structures; see send2). Furthermore, variable updates within a command can now be dependent on the updated values of other variables, provided there are no cyclic dependencies. See for example (s1'=c'?0:1), which updates s1 depending on whether there was a channel collision (reflected in c', the updated value of c). We use this mechanism to model interference on the channel: module channel specifies a joint probability distribution which is used to update variables s1 and s2 simultaneously.

#### **2.2 Property Specification**

PRISM-games 3.0 also extends the language used to specify properties for verification and strategy synthesis. The previous version already supported *zero-sum* queries for TSGs using the logic rPATL, which combines the game logic ATL with reward-based extensions of the probabilistic logic PCTL. Again, for the new real-time models, it is relatively easy to combine the existing rPATL notation with real-valued time bounds. So, we focus here on the case of CSGs, and in particular *equilibria-based* properties.

We compute values or synthesise strategies which are *social-welfare optimal Nash equilibria (SWNE)*, i.e., which maximise (or minimise) the sum of the values associated to the objectives for each player, but from which there is no incentive for any of them to unilaterally deviate in any state of the game. We express such properties by adding to rPATL the + operator, which is then used to denote the sum of the values associated to both *bounded* and *unbounded* objectives.

When using the rewards operator in equilibria-based properties, we can reason about *cumulative* (Ck), *instantaneous* (I=k) and *expected reachability* (F) objectives. For properties with the probability operator, we support bounded and unbounded reachability using the temporal operators *next* (X), *eventually* (F) and *until* (U). In order to express zero-sum properties for CSGs, we have implemented all the previous temporal operators for probabilistic queries and a subset of the rPATL operators reported in [8] for reward-based queries, adding to that the instantaneous reward operator.

Finally, following the style of rPATL we separate players into *coalitions* with the syntax -coalition, in order to specify the player or association of players for which we seek to maximise or minimise the values for a given zero-sum property. For equilibria-based properties, given that we maximise/minimise the sum, we use the same operator to separate players in different coalitions using a colon, while players in the same coalition are separated by a comma.

The following are examples of both zero-sum and equilibria-based properties for the medium access CSG model described in Fig. 1.

#### **3 Verification and Strategy Synthesis Algorithms**

#### **3.1 Zero-Sum Properties for CSGs**

When verifying zero-sum properties of CSGs, PRISM-games makes use of the model checking algorithms described in [19], which were based on the methods formulated in [2,3]. We rely on *value iteration* and classical convergence criteria to approximate/compute the values for all states of the game under study, and on solving a *linear program* to compute a *minimax* strategy at each state. This corresponds to solving a *matrix game*, which represents a *one-shot zero-sum* game for the actions of each player in a state. For unbounded properties, the solutions of the matrix games are used to synthesise an optimal (memoryless and randomised) strategy for each player. Prior to this numerical solution phase, we find and remove the states for which the optimal expected reward values are infinite by using the qualitative algorithms developed in [1].

Our current implementation uses the LPsolve [24] library to solve the matrix games at each state. CSGs are built and stored in a explicit-state fashion using an extension of PRISM's Java-implemented *explicit* (sparse-matrix based) engine.

#### **3.2 Equilibria-Based Properties for CSGs**

For equilibria-based properties of CSGs, PRISM-games implements the methods described in [20]. We rely on value iteration and *backwards induction* to approximate/compute values and synthesise strategies that are SWNE. For unbounded properties, we can only compute values that are ε-Nash equilibria, since Nash equilibria are not guaranteed to exist. At each state, we solve a *bimatrix game*, which is a representation of a *one-shot nonzero-sum* game and is a linear complementarity problem. We solve these games via *labelled polytopes*, finding all equilibria values through an SMT-based implementation, for which we use thirdparty SMT solvers Z3 [12] and Yices [13]. We make use of a precomputation step of finding and removing *dominated strategies* in order to minimise the number of calls to the solver.

Unlike zero-sum properties, the synthesised strategies for bounded and unbounded equilibria-based properties require (finite) memory. This is needed due to the fact that a player's choices may change once their objectives have been satisfied. We synthesise strategies by combining the strategy vectors computed for each bimatrix game and the strategy generated by computing optimal values for the MDP resulting from playing the game after either goal has been met. As we use value iteration to approximate values for infinite-horizon properties, we can only synthesise ε-Nash strategy profiles.

#### **3.3 Turn-Based Probabilistic Timed Games**

Verification and strategy synthesis of TPTGs relies on the algorithms from [18], which use the *digital clocks* approach that has been a developed for a variety of real-time models. A translation, at the level of the PRISM-games modelling languages, automatically converts the problem of analysing a TPTG into one of solving a (discrete-time) TSG, for which PRISM-games's existing engines can be used. Time-bounded properties are handled by automatically integrating a timing clock into the model prior to translation. As in the rest of PRISM-games, TSGs are also built and solved using the Java-based *explicit* engine.

## **4 Case Studies and Experimental Results**

The features added in PRISM-games 3.0 have been used for over 10 new case studies across a wide range of application domains, including computer security (intrusion detection, radio jamming, non-repudiation), communication protocols (medium access control, Aloha), incentive schemes for cooperative networking, multi-robot navigation problems and processor task scheduling. Details can be found in [18–20] and on the case studies section of the PRISM-games website [35]. Supporting material is at [33]. In this section, we showcase four selected case studies that demonstrate the benefits of the tool's new functionality. We also include a discussion of the scalability and performance of the tool.

**Future Markets Investor.** This example models two investors playing against the stock market. Investors choose when to invest or to cash in, and the stock market can decide to bar investments at certain points; fluctuations in share values are modelled stochastically. PRISM-games can, for example, synthesise optimal strategies for the two investors to maximise their expected joint profit over time, acting against the stock market which aims to minimise it.

(a) **Future markets investor:** avoiding unrealistic strategy choices using CSGs

(b) **Robot coordination:** using equilibria for mutually beneficial navigation plans

(c) **Non-repudiation:** Attack & defence strategies in a timed, randomised protocol (d) **Public good game:** Tuning incentive parameter *f* by synthesising equilibria

**Fig. 2.** Results illustrating the benefits of the new verification and strategy synthesis techniques implemented in PRISM-games 3.0; see Sect. 4 for details. (Color figure online)

Figure 2(a) shows the results obtained for this property using both a *turnbased* stochastic game (TSG) and a *concurrent* stochastic game (CSG). The former leads to unrealistic modelling as the market can see the choices made by the investors and gain an unfair advantage: the values in the blue plot in Fig. 2(a) are artificially low. In the CSG model, using PRISM 3.0, decisions are taken simultaneously, yielding the correct strategies and values (red plot).

**Robot Coordination.** Our next example models two robots navigating in opposite directions across a 10-by-10 grid as a CSG. Obstacles which hinder the robots as they move from location to location are modelled stochastically; and if the robots collide, both of them fail in their attempt to reach their goal. We use PRISM-games to find navigation strategies for the two robots, where each robot does not know the choice being made by the other at each step.

The objective for each robot is to navigate successfully, so we maximise the average probability (across the two robots) of success. Figure 2(b) shows the best value that can be achieved within a fixed period of k moves across the grid. One robot aiming single-handedly to achieve this goal performs reasonably well (blue plot), but we can achieve better collective performance by using PRISM-games to synthesise a (social welfare Nash) *equilibrium* strategy (red plot).

**Non-repudiation.** Next we consider a non-repudiation protocol [25], which permits an originator O to transfer information to a recipient R while guaranteeing non-repudiation, i.e., that neither O nor R can deny that they participated in the transfer. Here, both *probability* (the protocol is randomised) and *time* (the protocol relies on acknowledgement time-outs) are essential ingredients for checking correctness. Furthermore, we model the two participants of the protocol as opposing players, resulting in a TPTG model.

To verify the protocol, we check the worst-case probability that a malicious recipient R can obtain the information being transferred within time T. This can be done with a PTA model (as in [28]) but, with a timed game model, we can also analyse counter-strategies of the honest participant. The results (see Fig. 2(c)) show that, while it is not possible to prevent the information being received, it *is* possible to delay it (the red plot shows lower probabilities for higher times). Note that the bound T is an actual time bound, unlike the examples above, where step-bounded properties measure the number of steps or rounds.

**Public Good Game.** Lastly, we show a new case study modelling a *public good game*, a well studied model of social choice in economics where participants repeatedly decide how much of an endowment to keep for themselves or to share it with the other players. The total shared by the players is boosted by a factor f in order to incentivise sharing and then divided equally between the players.

Figure 2(d) shows results from a 2-player game, modelled as a CSG. Player choices are necessarily *concurrent*, to avoid cheating. We also need to use *equilibria* since the players have distinct individual goals (maximising personal expected profit). Figure 2(d) shows the values for each player in a synthesised optimal (social welfare Nash) equilibrium for varying f. Changes in f affect both the resulting profit *and* potential inequalities between players in equilibria, indicating the subtleties involved when tuning parameters in an incentive mechanism and the usefulness of analysing this with PRISM-games.

**Scalability and Performance.** Finally, we show some experimental results for a representative selection of larger examples, to give an indication of the scalability and performance of PRISM-games 3.0. Table 1 shows a range of models (the first 4 are CSGs; the last is a TPTG), the statistics for each one (number of

**Table 1.** Model statistics for some of the case studies.

players, states, transitions) and the time taken to build and verify the model for some example properties on a 2.10 GHz Intel Xeon with 8 GB of JVM memory.

Verification of CSGs is more computationally expensive than for TSGs supported in earlier versions of the tool, but PRISM-games 3.0 is able to build and analyse CSGs with more than 3 million states on relatively modest hardware. The majority of the time is spent solving (bi)matrix games, which is done repeatedly for all states of the model. Hence, the number of choices per state, which dictates the size of these games, has a greater impact on performance than for TSGs. Unsurprisingly, equilibria properties are slower than zero-sum ones. For both types of property, the number of players in the game does not have a major impact since they are grouped into coalitions yielding a 2-player game to solve. For TPTGs, the digital clocks translation is fast since it is done syntactically, and then a TSG is solved whose size depends on several factors, primarily the number of locations and the magnitude of any time bound in the property.

#### **5 Conclusions**

We have presented PRISM-games 3.0, which adds three major new features: (i) concurrent stochastic games; (ii) synthesis of equilibria; and (iii) timed probabilistic games. The usefulness of these has been illustrated on several newly created or extended applications.

CSGs are considerably more expensive to solve than their turn-based counterparts and a key challenge is efficiently solving the matrix game at each state, which is itself a non-trivial optimisation problem. For equilibria, the main difficulty is finding an optimal equilibrium, which currently relies on iteratively restricting the solution search space. Both problems are sensitive to the limitations and issues of floating-point arithmetic, particularly equilibria computation, and might benefit from arbitrary precision representations. Recent research has also pointed out the shortcomings of only using a lower bound approximation as a stopping criterion for value iteration, as it can lead to inaccuracies [4,14,17]. The impact of similar issues on model checking for games is still to be studied.

A range of further challenges exist for future work. These include providing support for *multi-coalitional* properties and implementing other techniques for equilibria computation. For timed games, we plan to investigate concurrent variants, and also zone-based solution techniques. More broadly speaking, partial information variants of games would be a useful addition.

**Acknowledgements.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 834115) and the EPSRC Programme Grant on Mobile Autonomy (EP/M019918/1).

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

## **Optimistic Value Iteration**

Arnd Hartmanns1(B) and Benjamin Lucien Kaminski2

1 University of Twente, Enschede, The Netherlands<EMAIL_ADDRESS>2 University College London, London, UK<EMAIL_ADDRESS>
**Abstract.** Markov decision processes are widely used for planning and verification in settings that combine controllable or adversarial choices with probabilistic behaviour. The standard analysis algorithm, value iteration, only provides lower bounds on infinite-horizon probabilities and rewards. Two "sound" variations, which also deliver an upper bound, have recently appeared. In this paper, we present a new sound approach that leverages value iteration's ability to usually deliver good lower bounds: we obtain a lower bound via standard value iteration, use the result to "guess" an upper bound, and prove the latter's correctness. We present this optimistic value iteration approach for computing reachability probabilities as well as expected rewards. It is easy to implement and performs well, as we show via an extensive experimental evaluation using our implementation within the mcsta model checker of the Modest Toolset.

#### **1 Introduction**

Markov decision processes (MDP, [30]) are a widely-used formalism to represent discrete-state and -time systems in which *probabilistic* effects meet controllable *nondeterministic* decisions. The former may arise from an environment or agent whose behaviour is only known statistically (e.g. message loss in wireless communication or statistical user profiles), or it may be intentional as part of a randomised algorithm (such as exponential backoff in Ethernet). The latter may be under the control of the system—then we are in a planning setting and typically look for a *scheduler* (or strategy, policy) that minimises the probability of unsafe behaviour or maximises a reward—or it may be considered adversarial, which is the standard assumption in verification: we want to establish that the maximum probability of unsafe behaviour is below, or that the minimum reward is above, a specified threshold. Extensions of MDP cover continuous time [11,26],

The authors are listed alphabetically. This work was partly performed while author B. L. Kaminski was at RWTH Aachen University, Aachen, Germany. This work was supported by ERC Advanced Grant 787914 (FRAPPANT), DFG Research Training Group 2236 (UnRAVeL), and NWO VENI grant no. 639.021.754.

S. K. Lahiri and C. Wang (Eds.): CAV 2020, LNCS 12225, pp. 488–511, 2020. https://doi.org/10.1007/978-3-030-53291-8\_26

and the analysis of complex formalisms such as stochastic hybrid automata [13] can be reduced to the analysis of MDP abstractions.

The standard algorithm to compute optimal (maximum or minimum) probabilities or reward values on MDP is *value iteration* (VI). It implicitly computes the corresponding optimal scheduler, too. It keeps track of a value for every state of the MDP, locally improves the values iteratively until a "convergence" criterion is met, and then reports the final value for the initial state as the overall result. The initial values are chosen to be an underapproximation of the true values (e.g. 0 for all states in case of probabilities or non-negative rewards). The final values are then an improved underapproximation of the true values. For unbounded (infinite-horizon) properties, there is unfortunately no (known and practical) convergence criterion that could guarantee a predefined error on the final result. Still, probabilistic model checkers such as Prism [24] report the final result obtained via simple relative or absolute global error criteria as the definitive probability. This is because, on *most* case studies considered so far, value iteration in fact converges fast enough that the (relative or absolute) difference between the reported and the true value approximately meets the error specified for the convergence criterion. Only relatively recently has this problem of soundness come to the attention of the probabilistic verification and planning communities [7,14,28]. First highlighted on hand-crafted counterexamples, it has by now been found to affect benchmarks and real-life case studies, too [3].

The first proposal to compute sound reachability probabilities was to use *interval iteration* (II [15], first presented in [14]). The idea is to perform two iterations concurrently, one starting from 0 as before, and one starting from 1. The latter improves an overapproximation of the true values, and the process can be stopped once the (relative or absolute) difference between the two values for the initial state is below the specified -, or at any earlier time with a correspondingly larger but known error. Baier et al. extended interval iteration to expected accumulated reward values [3]; here, the complication is to find initial values that are guaranteed to be an overapproximation. The proposed graph-based (i.e. not numerical) algorithm in practice tends to compute conservative initial values from which many iterations are needed until convergence. More recently, *sound value iteration* (SVI) [31] improved upon interval iteration by computing upper bounds on-the-fly and performing larger value improvements per iteration, for both probabilities and expected rewards. However, we found SVI tricky to implement correctly; some edge cases not considered by the algorithm as presented in [31] initially caused our implementation to deliver incorrect results or diverge on very few benchmarks. Both II and SVI fundamentally depend on the MDP being *contracting*; this must be ensured by appropriate structural transformations, e.g. by collapsing end components, a priori. These transformations additionally complicate implementations, and increase memory requirements.

*Our Contribution.* We present (in Sect. 4) a new algorithm to compute sound reachability probabilities and expected rewards that is both simple and practically efficient. We first (1) perform standard value iteration until "convergence", resulting in a lower bound on the value for every state. To this we (2) apply specific heuristics to "guess", for every state, a candidate upper bound value. Further value iterations (3) then confirm (if all values decrease) or disprove (if all values increase, or lower and upper bounds cross) the soundness of the upper bounds. In the latter case, we perform more lower bound iterations with reduced before retrying from step 2. We combine classic results from domain theory with specific properties of value iteration to show that our algorithm terminates. In problematic cases, many retries may be needed before termination, and performance may be worse than interval or sound value iteration. However, on many existing case studies, value iteration already worked well, and our approach attaches a soundness proof to its result with moderate overhead. We thus refer to it as *optimistic value iteration* (OVI). In contrast to II and SVI, it also works well for non-contracting MDP, albeit without a general termination guarantee. Our experimental evaluation in Sect. 5 uses all applicable models from the Quantitative Verification Benchmark Set [21] to confirm that OVI indeed performs as expected. It uses our publicly available implementations of II, SVI, and now OVI in the mcsta model checker of the Modest Toolset [20].

*Related Work.* In parallel to [15], the core idea behind II was also presented in [7] (later improved in [2]), embedded in a learning-based framework that manages to alleviate the state space explosion problem in models with a particular structure. In this approach, end components are statistically detected and collapsed on-thefly. II has recently been extended to stochastic games in [23], offering *deflating* as a new alternative to collapsing end components in MDP. Deflating does not require a structural transformation, but rather extra computation steps in each iteration applied to the states of all (a priori identified) end components.

The only known convergence criterion for pure VI was presented in [9, Sect. 3.5]: if we run VI until the absolute error between two iterations is less than a certain value α, then the computed values at that point are within α of the true values, and can in fact be rounded to the exact true values (as implemented in the *rational search* approach [5]). However, α cannot be freely chosen; it is a fixed number that depends on the size of the MDP and the largest denominator of the (rational) transition probabilities. The number of iterations needed is exponential in the size and the denominators. While not very useful in practice, this establishes an exponential upper bound on the number of iterations needed in unbounded-horizon VI. Additionally, Balaji et al. [4] recently showed the computations in finite-horizon value iteration to be EXPTIME-complete.

As an alternative to the iterative numeric road, guaranteed correct results (modulo implementation errors) can be obtained by using precise rational arithmetic. It does not combine too well with iterative methods like II or SVI due to the increasingly small differences between the values and the actual solution. The probabilistic model checker Storm [10] thus combines topological decomposition, policy iteration, and exact solvers for linear equation systems based on Gaussian elimination when asked to use rational arithmetic [22, Section 7.4.8]. The disadvantage is the significant runtime cost for performing the unlimitedprecision calculations, limiting such methods to relatively smaller MDP.

The only experimental evaluations using large sets of benchmarks that we are aware of compared VI with II to study the overhead needed to obtain sound

**Table 1.** VI and OVI example on Me

results via II [3], and II with SVI to show the performance improvements of SVI [31]. The learning-based method with deflation of [2] does not compete against II and SVI; its aim is rather in dealing with state space explosion (i.e. memory usage). Its performance was evaluated on 16 selected small (<400 k states) benchmark instances in [2], showing absolute errors on the order of 10−4 on many benchmarks with a 30-min timeout. SVI thus appears the most competitive technique in runtime and precision so far. Consequently, in our evaluation in Sect. 5, we compare OVI with SVI, and II for reference, using the default relative error of 10−6, including large and excluding clearly acyclic benchmarks (since they are trivial even for VI), with a 10-min timeout which is rarely hit.

#### **2 Preliminaries**

R+ 0 is the set of all non-negative real numbers. We write { x1 -→ y1,... } to denote the function that maps all xi to yi, and if necessary in the respective context, implicitly maps to 0 all x for which no explicit mapping is specified. Given a set S, its powerset is 2S. A (discrete) *probability distribution* over S is a function μ ∈ S → [0, 1] with countable *support spt*(μ) def = { s ∈ S | μ(s) > 0 } and - s∈*spt*(μ) μ(s) = 1. *Dist*(S) is the set of all probability distributions over S.

*Markov Decision Processes* (MDP) combine nondeterministic choices as in labelled transition systems with discrete probabilistic decisions as in discrete-time Markov chains (DTMC). We define them formally and describe their semantics.

**Definition 1.** *A* Markov decision process *(MDP) is a triple* M = S, sI , T *where* S *is a finite set of* states *with* initial state sI ∈ S *and* T : *S* → 2*Dist*(R+ 0 ×S) *is the* transition function*.* T(s) *must be finite and non-empty for all* s ∈ S*.*

For s ∈ S, an element of T(s) is a *transition*, and a pair r, s ∈ *spt*(T(s)) is a *branch* to successor state s with *reward* r and probability T(s)(r, s ). Let M(s- *I* ) be M but with initial state s I , and M0 be M with all rewards set to zero. *Example 1.* Figure 1 shows our example MDP Me. We draw transitions as lines to an intermediate node from which branches labelled with probability and reward (if not zero) lead to successor states. We omit the intermediate node and probability 1 for transitions with a single branch, and label some transitions to refer to them in the text. Me has 5 states, 7 transitions, and 10 branches.

In practice, higher-level modelling languages like Modest [17] are used to specify MDP. The semantics of an MDP is captured by its *paths*. A path represents a concrete resolution of all nondeterministic and probabilistic choices. Formally:

**Definition 2.** *A* finite path *is a sequence* πfin = s0 μ0 r0 s1 μ1 r1 ...μn−1rn−1sn *where* si ∈ S *for all* i ∈ { 0,...,n } *and* ∃ μi ∈ T(si): ri, si+1 ∈ *spt*(μi) *for all* i ∈ { 0,...,n − 1 }*. Let* |πfin| def = n*,* last(πfin) def = sn*, and* rew(πfin) def = n−1 i=0 ri*.* Π*fin is the set of all finite paths starting in* sI *. A* path *is an analogous infinite sequence* π*, and* Π *is the set of all paths starting in* sI *. We write* s ∈ π *if* ∃ i: s = si*, and* π→G *for the shortest prefix of* π *that contains a state in* G ⊆ S*, or* ⊥ *if* π *contains no such state. Let* rew(⊥) def = ∞*.*

A scheduler (or *adversary*, *policy* or *strategy*) only resolves the nondeterministic choices of M. For this paper, memoryless deterministic schedulers suffice [6].

**Definition 3.** *A function* s : S → *Dist*(R+ 0 × S) *is a* scheduler *if, for all* s ∈ S*, we have* s(s) ∈ T(s)*. The set of all schedulers of* M *is* S(M)*.*

Given an MDP M as above, let M|s = S, sI , T|s with T|s(s) = { s(s) } be the DTMC induced by s. Via the standard cylinder set construction [12, Sect. 2.2] on M|s, a scheduler induces a probability measure PM s on measurable sets of paths starting in sI . For goal state g ∈ S, the maximum and minimum **probability of reaching** *g* is defined as PM max( g) = sups∈S PM s ({ π ∈ Π | g ∈ π }) and PM min( g) = infs∈S PM s ({ π ∈ Π | g ∈ π }), respectively. The definition extends to sets G of goal states. Let RM G : Π → R+ 0 be the random variable defined by RM G (π) = rew(π→G) and let EM s (G) be the expected value of RM G under PM s . Then the maximum and minimum **expected reward to reach** *G* is defined as EM max(G) = sups EM s (G) and EM min(G) = infs EM s (G), respectively. We omit the superscripts for M when they are clear from the context. From now on, whenever we have an MDP with a set of goal states G, we assume that they have been made absorbing, i.e. for all g ∈ G we only have a self-loop: T(g) = {{0, g -→ 1 } }.

**Definition 4.** *An* end component *of* M *as above is a (sub-)MDP* S , T , s I *where* S ⊆ S*,* T (s) ⊆ T(s) *for all* s ∈ S *, if* μ ∈ T (s) *for some* s ∈ S *and* r, s ∈ *spt*(μ) *then* r = 0*, and the directed graph with vertex set* S *and edge set* { s, s |∃ μ ∈ T (s): 0, s ∈ *spt*(μ) } *is strongly connected.*

#### **3 Value Iteration**

The standard algorithm to compute reachability probabilities and expected rewards is *value iteration* (VI) [30]. In this section, we recall its theoretical foundations and its limitations regarding convergence.

```
1 function GSVI(M = -
                 S, sI , T, S?, v, α, diff )
2 repeat
3 error := 0
4 foreach s ∈ S? do
5 vnew := Φ(v)(s) // iterate lower bound
6 if vnew > 0 then error := max(error, diff (v(s), vnew ))
7 v(s) := vnew
8 until error ≤ α
```
**Algorithm 1.** Gauss-Seidel value iteration

#### **3.1 Theoretical Foundations**

Let V = { v | v : S → R+ 0 ∪ {∞} } be a space of vectors of values. It can easily be shown that V, with

$$v \preceq w \qquad \text{if and only if} \qquad \forall s \in S \colon v(s) \le w(s)$$

forms a complete lattice, i.e. every subset V ⊆ V has a supremum (and an infimum) in V with respect to . We write v ≺ w for v w ∧ v = w and v ∼ w for ¬(v w ∨ w v).

Minimum and maximum reachability probabilities and expected rewards can be expressed as the *least fixed point* of the *Bellman operator* Φ: V → V given by

$$\Phi(v) \stackrel{\text{def}}{=} \lambda \, s. \begin{cases} opt\_{\mu \in T(s)} \sum\_{\langle r, s' \rangle \in spt(\mu)} \mu(s') \cdot (r + v(s')) & \text{if } s \in S\_{\text{?}}\\ d & \text{if } s \notin S\_{\text{?}} \end{cases}$$

where *opt* ∈ { max, min } and the choice of both S? ⊆ S and d depends on whether we wish to compute reachability probabilities or expected rewards. In any case, the Bellman operator Φ can be shown to be Scott-continuous [1], i.e. in our case: for any subset V ⊆ V, we have Φ(sup V ) = supΦ(V ).

The Kleene fixed point theorem for Scott-continuous self-maps on complete lattices [1,27] guarantees that lfpΦ, the least fixed point of Φ, indeed exists. Note that Φ can still have more than one fixed point. In addition to mere existence of lfp Φ, the Kleene fixed point theorem states that lfp Φ can be expressed by

$$\text{lfp } \Phi = \lim\_{n \to \infty} \Phi^n(\bar{0}) \tag{1}$$

where ¯0 ∈ V is the zero vector and Φn(v) denotes n-fold application of Φ to v. Equation 1 is the basis of VI: the algorithm iteratively constructs a sequence of vectors

v0 = ¯0 and vi+1 = Φ(vi),

which converges to the sought-after least fixed point. This convergence is *monotonic*: for every n ∈ N, we have Φn(¯0) Φn+1(¯0) and hence Φn(¯0) lfp Φ. In particular, Φn(¯0)(sI ) is an *under*approximation of the sought-after quantity for every n. Note that iterating Φ on *any* underapproximation v lfp Φ (instead of ¯0) will still converge to lfp Φ and Φn(v) lfp Φ will hold for any n.

*Gauss-Seidel Value Iteration.* Algorithm 1 shows the pseudocode of a VI implementation that uses the so-called *Gauss-Seidel optimisation*: Whereas standard VI needs to store two vectors vi and vi+1, Gauss-Seidel VI stores only a single vector v and performs updates in place. This does not affect the correctness of VI, but may speed up convergence depending on the order in which the loop in line 4 considers the states in S?. The error metric *diff* is used to check for convergence.

*VI for Probabilities.* For determining reachability probabilities, we operate on M0 and set S? = S \ G and d = 1. Then the corresponding Bellman operator satisfies

$$(\mathfrak{fp}\,\Phi)(s) = \mathcal{P}\_{opt}^{M^{(s)}}(\diamond G),$$

and VI will iteratively approximate this quantity *from below*. The corresponding call to Algorithm 1 is GSVI(M0, S \ G, { s -→ 0 | s ∈ S \ G }∪{ s -→ 1 | s ∈ G }, α, *diff* ).

*VI for Expected Rewards.* For determining the expected reward EM(*s*) *opt* (G), we operate on M and first have to determine the set S∞ of states from which the minimum (if *opt* = max) or maximum (if *opt* = min) probability to reach G is less than 1.1 If sI ∈ S∞, then the result is ∞ due to the definition of rew(⊥). Otherwise, we choose S? = S \ S∞ and d = ∞. Then, for *opt* = max, the least fixed point of the corresponding Bellman operator satisfies

$$(\mathfrak{fp}\,\Phi)(s) = \mathrm{E}\_{opt}^{M^{(s)}}(G).$$

Again, VI underapproximates this quantity. The same holds for *opt* = min if M does not have end components containing states other than those in G and S∞. The corresponding call to Algorithm 1 is GSVI(M, S \ S∞, { s -→ 0 | s ∈ S \ S∞ }∪{ s -→∞| s ∈ S∞ }, α, *diff* ).

#### **3.2 Uniqueness of Fixed Points**

lfp Φ may not be unique for two reasons: states that cannot reach G under the optimal scheduler may take any value (causing fixed points greater than lfpΦ for Pmin and Pmax), and states in end components may take values higher than lfpΦ. The latter affects Pmax (higher fixed points) and Emin (lower fixed points).

*Example 2.* In Me of Fig. 1, s1 and s2 and the two transitions in-between form an end component. For PM*e* max( { s+ }), v = { s -→ 1 } is a non-least fixed point for the corresponding Bellman operator; with appropriate values for s1 and s2, we can obtain fixed points with any v(s0) > 0.5 of our choice. Similarly, we have EM min({ s+, s− })=0.6 (by scheduling b in s0), but due to the end component (with only zero-reward transitions by definition), the fixed point is s.t. v(s0) = 0.

1 This can be done via Algs. 2 (for S1 min) and 4 (for S1 max) of [12], respectively. These algorithms do not consider the probabilities, but only whether there is a transition and branch (with positive probability) from one state to another or not. We thus call them graph-based algorithms, as opposed to numeric algorithms like VI itself.

VI works for Pmin, Pmax, and Emax with multiple fixed points: we anyway seek lfpΦ and start from a (trivial) underapproximation. For Emin, (zero-reward) end components need to be collapsed: we determine the maximal end components using algorithms similar to [15, Alg. 1], then replace each of them by a single state, keeping all transitions leading out of the end component. We refer to this as the *ECC* transformation. However, such end components rarely occur in case studies for Emin since they indicate Zeno behaviour w.r.t. to the reward. As rewards are often associated to time progress, such behaviour would be unrealistic.

To make the fixed points unique, for Emax and Emin we fix the values of all states in G to 0. For Pmin, we precompute the set S0 min of states that reach G with minimum probability 0 using Alg. 1 of [12], then fix their values to 0. For Pmax, we analogously use S0 max via Alg. 3 of [12]. For Pmax and Emin, we additionally need to remove end components via ECC. In contrast to the precomputations, ECC changes the structure of the MDP and is thus more memory-intensive.

#### **3.3 Convergence**

VI and GSVI will not *reach* a fixed point in general, except for special cases such as acyclic MDP. It is thus standard to use a convergence criterion based on the difference between two consecutive iterations (lines 6 and 8) to make GSVI terminate: we either check the *absolute error*, i.e.

$$\|\vec{df} = \vec{df}\_{abs} \stackrel{\text{def}}{=} \lambda \left< v\_{old}, v\_{new} \right> , \ v\_{new} - v\_{old},$$

or the *relative error*, i.e.

$$diff = diff\_{rel} \stackrel{\text{def}}{=} \lambda \left< v\_{old}, v\_{new} \right> . (v\_{new} - v\_{old}) / v\_{new} \dots$$

By default, probabilistic model checkers like Prism and Storm use *diffrel* and α = 10−6. Upon termination of GSVI, v is then closer to the least fixed point, but remains an underapproximation. In particular, α has, in general, no relation to the final difference between v(sI ) and P*opt*( G) or E*opt*(G), respectively.

*Example 3.* Consider MDP Me of Fig. 1 again with G = { s+ }. The first four rows in the body of Table 1 show the values for v after the i-th iteration of the outer loop of a call to GSVI(M0 e , { s0, s1, s2 }, max, { s+ -→ 1 }∪{ s -→ 0 | s = s+ }, 0.05, *diffabs* ). After the fourth iteration, GSVI terminates since the error is less than α = 0.05; at this point, we have Pmax( s+) − v(s0)=0.08 > α.

To obtain a value within a prescribed error of the true value, we can compute an upper bound in addition to the lower bound provided by VI. Interval iteration (II) [3,15] does so by performing, in parallel, a second value iteration on a second vector u that starts from a known overapproximation. For probabilities, the vector ¯1 = { s -→ 1 } is a trivial overapproximation; for rewards, more involved graph-based algorithms need to be used to precompute (a very conservative) one [3]. II terminates when *diff* (v(sI ), u(sI )) ≤ 2-

**Table 2.** Preprocessing requirements of value iteration variants

and returns v*II* = 1 2 (u(sI ) + v(sI )). With v*true* = P*opt*( G), II thus guarantees that v*II* ∈ [v*true* − - · v*true* , v*true* + - · v*true* ] and analogously for expected rewards. However, to ensure termination, II requires a unique fixed point: u converges from above to the greatest fixed point gfp Φ, thus for every MDP where *diff* ((lfp Φ)(sI ),(gfp Φ)(sI )) > 2-, II diverges. For Pmax, we have gfp Φ(s*ec*)=1 for all s*ec* in end components, thus II tends to diverge when there is an end component. Sound value iteration (SVI) [31] is similar, but uses a different approach to derive upper bounds that makes it perform better overall, and that eliminates the need to precompute an initial overapproximation for expected rewards. However, SVI still requires unique fixed points.

We summarise the preprocessing requirements of VI, II, and SVI in Table 2. With unique fixed points, we can transform Pmin into Pmax by making S0 min states absorbing and setting G to S0 min, and Pmax into Emax by a similar transformation adding reward 1 to entering G. Most of the literature on VI variants works in such a setting and describes the Pmax or Emax case only. Since OVI also works with multiple fixed points, we have to consider all four cases individually.

#### **4 Optimistic Value Iteration**

We now present a new, practical solution to the convergence problem for unbounded reachability and expected rewards. It exploits the empirical observation that on many case studies VI delivers results which are roughly α-close to the true value—it only lacks the ability to prove it. Our approach, *optimistic value iteration* (OVI), extends standard VI with the ability to deliver such a proof.

The key idea is to exploit a property of the Bellman operator Φ and its Gauss-Seidel variant as in Algorithm 1 to determine whether a candidate vector is a lower bound, an upper bound, or neither. The foundation is basic domain theory: by Scott-continuity of Φ it follows that Φ is monotonic, meaning v w implies Φ(v) Φ(w). A principle called *Park induction* [29] for monotonic self-maps on complete lattices yields the following induction rules: For any u ∈ V,

 **function** OVI(M = -S, sI , T, S?, v, , α, diff ) GSVI(M,S?, v, α, diff ) // perform standard value iteration u := { s → diff +(s) | s ∈ S? }, viters := 0 // guess candidate upper bound **while** viters < 1 α **do** // start verification phase up∀ := true, down∀ := true, viters := viters + 1, error := 0 **foreach** s ∈ S? **do** v*new* := Φ(v)(s), u*new* := Φ(u)(s) // iterate both bounds **if** v*new* > 0 **then** error := max { error, diff (v(s), v*new* ) } **if** u*new* < u(s) **then** // upper value decreased: u(s) := u*new* , up∀ := false // update u with new lower u*new* **else if** u*new* > u(s) **then** // upper value increased: down∀ := false // discard new higher u*new* v(s) := v*new* // update v with new value v*new* **if** v(s) > u(s) **then goto** line 17 // lower bound crossed u **if** down∀ **then return** 1 (u(sI ) + v(sI )) // u is inductive upper bound **else if** up∀ **then goto** line 17 // u is inductive lower bound **return** OVI(M,S?, v, , *error* , diff ) // retry with reduced α

**Algorithm 2.** Optimistic value iteration

$$
\Phi(u) \preceq u \qquad \text{implies} \qquad \mathsf{lfp} \Phi \preceq u. \tag{2}
$$

$$\text{and} \qquad u \preceq \Phi(u) \qquad \text{implies} \qquad u \preceq \mathbf{gfp} \,\Phi. \tag{3}$$

Thus, if we can construct a candidate vector u s.t. Φ(u) u, then u is in fact an upper bound on the sought-after lfpΦ.  Optimistic value iteration uses this insight and can be summarised as follows:

The resulting procedure in more detail is shown as Algorithm 2. Starting from the same initial vectors v as for VI, we first perform standard Gauss-Seidel value iteration (in line 2). We refer to this as the *iteration phase* of OVI. After that, vector v is an improved underapproximation of the actual probabilities or reward values. We then "guess" a vector u of *upper values* from the *lower values* in v (line 3). The guessing heuristics depends on *diff* : if *diff* = *diffabs* , then we use

$$\operatorname{dif}^+(s) = \begin{cases} 0 & \text{if } v(s) = 0\\ v(s) + \epsilon & \text{otherwise} \end{cases}$$

if *diff* = *diffrel* , then

$$\operatorname{dif}^+(s) = v(s) \cdot (1 + \epsilon).$$

We cap the result at 1 for Pmin and Pmax. These heuristics have three important properties: (**H1**) v(s) = 0 implies *diff* +(s) = 0, (**H2**) *diff* (v(s), *diff* +(s)) ≤ 2-, and (**H3**) *diff* (v(s), *diff* +(s)) > 0 unless v(s) = 0 or v(s) = 1 for Pmin and Pmax.

Then the *verification phase* starts in line 4: we perform value iteration on the lower values v and upper values u at the same time, keeping track of the direction in which the upper values move. For u, line 7 and the conditions around line 10 mean that we actually use operator Φmin(u) = λ s. min(Φ(u)(s), u(s)). This may shorten the verification phases, and is crucial for our termination argument. A state s is *blocked* if Φ(u)(s) > Φmin(u)(s) and *unblocked* if Φ(u)(s) < u(s) here.

If, in some iteration, no state was blocked (line 15), then we had Φ(u) u before the start of the iteration. We thus know by Eq. 2 that the current u is an inductive upper bound for the values of all states, and the true value must be in the interval [v(sI ), u(sI )]. By property H2, our use of Φmin for u, and the monotonicity of Φ as used on v, we also know that *diff* (v(sI ), u(sI )) ≤ 2-, so we immediately terminate and return the interval's centre vI = 1 2 (u(sI ) + v(sI )). The true value v*true* = (lfp Φ)(sI ) must then be in [vI − - · v*true* , vI + -· v*true* ].

If, in some iteration, no state was unblocked (line 16), then again by Park induction we know that u gfpΦ. If we are in a situation of unique fixed points, this also means u lfp Φ, thus the current u is no upper bound: we cancel verification and go back to the iteration phase to further improve v before trying again. We do the same if v crosses u: then u(s) < v(s) ≤ (lfp Φ)(s) for some s, so this u was just another bad guess, too.

Otherwise, we do not yet know the relationship between u and lfp Φ, so we remain in the verification phase until we encounter one of the cases above, or until we exceed the verification budget of 1 α iterations (as checked by the loop condition in line 4). This budget is a technical measure to ensure termination.

*Optimisation.* In case the fixed point of Φ is *unique*, by Park induction (via Eq. 3) we know that u Φ(u) implies that u is a lower bound on lfp Φ. In such situations of single fixed points, we can—as an optimisation—additionally replace v by u before the *goto* in line 16.

*Heuristics.* OVI relies on heuristics to gain an advantage over alternative methods such as II or SVI; it cannot be better on *all* MDP. Concretely, we can choose

Algorithm 2 shows the choices made by our implementation. We employ the standard stopping criteria used by probabilistic model checkers for VI, and the "weakest" guessing heuristics that satisfies properties H1, H2, and H3 (i.e. guessing any higher values would violate one of these properties). The only arbitrary choice is how to reduce α, which we at least halve on every retry. We experimentally found this to be a good compromise on benchmarks that we consider in Sect. 5, where

*Example 4.* We now use the version of Φ to compute Pmax and call

OVI(M0 e , { s0, s1, s2 }, { s+ -→ 1 }∪{ s -→ 0 | s = s+ }, 0.05, 0.05, *diffabs* ).

Table 1 shows the values in v and u during this run, assuming that we use non-Gauss-Seidel iterations. The first iteration phase lasts from i = 0 to 4. At this point, u is initialised with the values shown in italics. The first verification phase needs only one iteration to realise that u is actually a lower bound (to a fixed point which is not the least fixed point, due to the uncollapsed end component). Blocked states are marked with a bar; unblocked states have a lower u-value than in the previous iteration. We resume GSVI from i = 6. The error in GSVI is again below α, which had been reduced to 0.008, during iteration i = 9. We thus start another verification phase, which immediately (in one iteration) finds the newly guessed vector u to be an upper bound, with *diff* (v(s0), u(s0)) < 2-.

#### **4.1 Termination of OVI**

We showed above that OVI returns an --correct result when it terminates. We now show that it terminates in all cases except for Pmax with multiple fixed points. Note that this is a stronger result than what II and SVI can achieve.

Let us first consider the situations where lfp Φ is the unique fixed point of Φ. First, GSVI terminates by Eq. 1. Let us now write vi and ui for the vectors u and v as they are at the beginning of verification phase iteration i. We know that v0 u0. We distinguish three cases relating the initial guess u0 to lfp Φ.

1. u0 ∼ lfpΦ or u0 ≺ lfpΦ, i.e. there is a state s with u0(s) < (lfpΦ)(s). Since we use Φmin on the upper values, it follows ui(s) ≤ u0(s) < (lfpΦ)(s) for all i. By Eq. 1, there must thus be a j such that vj (s) > uj (s), triggering a retry with reduced α in line 14. Such a retry could also be triggered earlier in line 16. Due to the reduction of α and Eq. 1, every call to GSVI will further increase some values in v or reach v = lfpΦ (in special cases), and for some subsequent guess u we must have u0(s) < u(s). Consequently, after some repetitions of this case 1, we must eventually guess a u with lfp Φ u.

ECC

**Fig. 2.** DTMC Md

e without

**Table 3.** Nontermination of OVI on M

2. lfp Φ ≺ u0. Observe that operators Φ and Φmin are *local* [9], i.e. a state's value can only change if a direct successor's value changes. In particular, a state's value can only decrease (increase) if a direct successor's value decreases (increases). If ui(s) < ui−1(s), then s cannot be blocked again in any later iteration j>i: for it to become blocked, a successor's upper value would have to increase, but Φmin ensures non-increasing upper values for all states. Analogously to Eq. 1, we know that [3, Lemma 3.3 (c)]

$$\text{lfp } \Phi \preceq u \quad \text{implies} \quad \lim\_{n \to \infty} \Phi\_{\text{min}}^n(u) = \text{lfp } \Phi$$

(for the unique fixpoint case, since [3] assumes contracting MDP as usual). Thus, for all states s, there must be an i such that ui(s) < ui−1(s); in consequence, there is also an iteration j where no state is blocked any more. Then the condition in line 15 will be true and OVI terminates.

3. lfp Φ u0 but not lfp Φ ≺ u0, i.e. there is a state s with u0(s)=(lfp Φ)(s). If there is an i where no state, including s, is blocked, then OVI terminates as above. For Pmin and Pmax, if u0(s) = 1, s cannot be blocked, so we can w.l.o.g. exclude such s. For other s not to be blocked in iteration i, we must have ui(s )=(lfp Φ)(s ) for all states s reachable from s under the optimal scheduler, i.e. all of those states must *reach* the fixed point. This cannot be guaranteed on general MDP. Since this case is a very particular situation unlikely to be encountered in practice with our heuristics, OVI adopts a pragmatic solution: it bounds the number of iterations in every verification phase (cf. line 4). Due to property H3 of our heuristics, u0(s)=(lfp Φ)(s) requires v0(s) < (lfp Φ)(s), thus some subsequent guess u will have u(s) > u0(s), and eventually we must get a u with lfp Φ ≺ u, which is case 2. Since we strictly increase the iteration bound on every retry, we will eventually encounter case 2 with a sufficiently high bound for termination.

Three of the four situations with multiple fixed points reduce to the corresponding unique fixed point situation due to property H1 of our guessing heuristics:

1. For Pmin, recall from Sect. 3.2 that the fixed point is unique if we fix the values of all S0 min states to 0. In OVI without preprocessing, such states are in S?, thus they initially have value 0. Φ will not increase their values, neither will guessing due to H1, and neither will Φmin. Thus OVI here operates on a sublattice of V, where the fixed point of Φ is unique.

The only case where OVI may not terminate is for Pmax without ECC. Here, end components may cause states to be permanently blocked. However, we did not encounter this on any benchmark used in Sect. 5, so in contrast to e.g. II, OVI is still *practically* useful in this case despite the lack of a termination guarantee.

*Example 5.* We turn Me of Fig. 1 into M e by replacing the c-labelled transition from s2 by transition { 0, s2 -→ 1 2 ,0, s+ -→ 1 4 ,1, s− -→ 1 4 }, i.e. we can now go from s2 back to s2 with probability 1 2 and to each of s+, s− with probability 1 4 . The probability-1 transition from s2 to s1 remains. Then Table 3 shows a run of OVI for Pmax with *diffabs* and α = 0.1. s0 is forever blocked from iteration 6 on.

#### **4.2 Variants of OVI**

While the core idea of OVI rests on classic results from domain theory, Algorithm 2 includes several particular choices that work together to achieve good performance and ensure termination. We sketch two variants to motivate these choices.

First, let us use Φ instead of Φmin for the upper values, i.e. move the assignment u(s) := u*new* down into line 13. Then we cannot prove termination because the arguments of case 2 for lfp Φ ≺ u0 no longer hold. Consider DTMC Md of Fig. 2 and Pmax( s+)=Pmin( s+). Let

$$u = \{ \begin{array}{c} s\_0 \mapsto 0.2, s\_1 \mapsto 1, s\_+ \mapsto 1, s\_- \mapsto 0 \end{array} \} \succ \{ \begin{array}{c} s\_0 \mapsto \frac{1}{9}, s\_1 \mapsto \frac{1}{9}, \dots \end{array} \} = \mathsf{lfp} \,\Phi.$$

Iterating Φ, we then get the following sequence of pairs u(s0), u(s1):

$$
\langle 0.2, 1 \rangle, \langle 1, 0.12 \rangle, \langle 0.12, 0.2 \rangle, \langle 0.2, 0.112 \rangle, \langle 0.112, 0.12 \rangle, \langle 0.12, 0.1112 \rangle, \dots
$$

Observe how the value of s0 increases iff s1 decreases and vice-versa. Thus we never encounter an inductive upper or lower bound. In Algorithm 2, we use Gauss-Seidel VI, which would not show the same effect on this model; however, if we insert another state between s0 and s1 that is updated last, Algorithm 2 would behave in the same alternating way. This particular u is contrived, but we could have guessed one with a similar relationship of the values leading to similar behaviour.

An alternative that allows us to use Φ instead of Φmin is to change the conditions that lead to retrying and termination: We separately store the initial guess of a verification phase as u0, and then compare each newly calculated u with u0. If u u0, then we know that there is an i such that u = Φi (u) u0. Φi retains all properties of Φ needed for Park induction, so this would also be a proof of lfp Φ u. The other conditions and the termination proofs can be adapted analogously. However, this variant needs ≈50 % more memory (to store an additional vector of values), and we found it to be significantly slower than Algorthm 2 and the first variant on almost all benchmark instances of Sect. 5.

#### **5 Experimental Evaluation**

We have implemented interval iteration (II) (using the "variant 2" approach of [3] to compute initial overapproximations for expected rewards), sound value iteration (SVI), and now optimistic value iteration (OVI) precisely as described in the previous section, in the mcsta model checker of the Modest Toolset [20], which is publicly available at modestchecker.net. It is cross-platform, implemented in C#, and built around the Modest [17] high-level modelling language. Via support for the Jani format [8], mcsta can exchange models with other tools like Epmc [18] and Storm [10]. Its performance is competitive with Storm and Prism [16]. We tried to spend equal effort performance-tuning our VI, II, SVI, and OVI implementations to avoid unfairly comparing highly-optimised OVI code with na¨ıve implementations of the competing algorithms.

In the following, we report on our experimental evaluation of OVI using mcsta on all applicable models of the Quantitative Verification Benchmark Set (QVBS) [21]. All models in the QVBS are available in Jani and can thus be used by mcsta. Most are parameterised, and come with multiple properties of different types. Aside from MDP models, the QVBS also includes DTMCs (which are a special case of MDP), continuous-time Markov chains (CTMC, for which the analysis of unbounded properties reduces to checking the embedded DTMC), Markov automata (MA [11], on which the embedded MDP suffices for unbounded properties), and probabilistic timed automata (PTA [26], some of which can be converted into MDP via the digital clocks semantics [25]). We use all of these model types. The QVBS thus gives rise to a large number of benchmark *instances*: combinations of a model, a parameter valuation, and a property to check. For every model, we chose one instance per probabilistic reachability and expected-reward property such that state space exploration did not run out of memory and VI took at least 10 s where possible. We only excluded

**Fig. 3.** OVI runtime and iteration count compared to VI (probabilistic reachability)

As a result, we considered 38 instances with probabilistic reachability and 41 instances with expected-reward properties, many comprising several million states.

We ran all experiments on an Intel Core i7-4790 workstation (3.6–4.0 GHz) with 8 GB of memory and 64-bit Ubuntu Linux 18.04. By default, we request a relative half-width of - = 10−6 for the result probability or reward value, and configure OVI to use the relative-error criterion with α = 10−6 in the iteration phase. We use a 600 s timeout ("TO"). Due to the number of instances, we show most results as scatter plots like in Fig. 3. Each such plot compares two methods in terms of runtime or number of iterations. Every point x, y corresponds to an instance and indicates that the method noted on the x-axis took x seconds or iterations to solve this instance while the method noted on the y-axis took y seconds or iterations. Thus points above the solid diagonal line correspond to instances where the x-axis method was faster (or needed fewer iterations); points above (below) the upper (lower) dotted diagonal line are where the x-axis method took less than half (more than twice) as long or as many iterations.

#### **5.1 Comparison with VI**

All methods except VI delivered correct results up to -. VI offers low runtime at the cost of occasional incorrect results, and in general the absence of any guarantee about the result. We thus compare with VI separately to judge the overhead caused by performing additional verification, and possibly iteration, phases. This is similar to the comparison done for II in [3]. Figures 3 and 4 show the results. The unfilled shapes indicate instances where VI produced an incorrect result. In terms of runtime, we see that OVI does not often take more than twice as long as VI, and frequently requires less than 50% extra time. On several instances where OVI incurs most overhead, VI produces an incorrect result, indicating

**Fig. 4.** OVI runtime and iteration count compared to VI (expected rewards)

that they are "hard" instances for value iteration. The unfilled CTMCs where OVI takes much longer to compute probabilities are all instances of the *embedded* model; the DTMC on the x-axis is *haddad-monmege*, an adversarial model built to highlight the convergence problem of VI in [14]. The problematic cases for expected rewards include most MA instances, the two expected-reward instances of the *embedded* CTMC, and again *haddad-monmege*. In terms of iterations, the overhead of OVI is even less than in runtime.

#### **5.2 Comparison with II and SVI**

We compare the runtime of OVI with the runtime of II and that of SVI separately for reachability probabilities (shown in Fig. 5) and expected rewards (shown in Fig. 6). As shown in Table 2, OVI has almost the same requirements on precomputations as VI, while II and SVI require extra precomputations and ECC for reachability probabilities. The precomputations and ECC need extra runtime (which turned out to be negligible in some cases but significant enough to cause a timeout in others) prior to the numeric iterations. However, doing the precomputations can reduce the size of the set S?, and ECC can reduce the size of the MDP itself. Both can thus reduce the runtime needed for the numeric iterations. For the overall runtime, we found that none of these effects dominates the other over all models. Thus sometimes it may be better to perform only the required precomputations and transformations, while on other models performing all applicable ones may lead to lower total runtime. For reachability probabilities, we thus compare OVI, II, and SVI in two scenarios: once in the default ("std") setting of mcsta that uses only required preprocessing steps

**Fig. 5.** OVI runtime compared to II and SVI (probabilities)

(without ECC for OVI; we report the total runtime for preprocessing and iterations), and once with all of them enabled ("pre", where we report only the runtime for numeric iterations, plus the computation of initial upper bounds in case of II).

For probabilistic reachability, we see in Fig. 5 that there is no clear winner among the three methods in the "std" setting (top plots). In some cases, the extra precomputations take long enough to give an advantage to OVI, while in others they speed up II and SVI significantly, compensating for their overhead. The "pre" setting (bottom), in which all three algorithms operate on exactly the same input w.r.t. to MDP M and set S?, however, shows a clearer picture: now OVI is faster, sometimes significantly so, than II and SVI on most instances.

**Fig. 6.** OVI runtime compared to II and SVI (expected rewards)

Expected-reward properties were more challenging for all three methods (as well as for VI, which produced more errors here than for probabilities). The plots in Fig. 6 paint a very clear picture of OVI being significantly faster for expected rewards than II (which suffers from the need to precompute initial upper bounds that then turn out to be rather conservative), and faster (though by a lesser margin and with few exceptions) than SVI.

In Fig. 7, we give a summary view combining the data from Figs. 3 to 6. For each algorithm, we plot the instances sorted by runtime, i.e. a point x, y on the line for algorithm z means that some instance took y seconds to solve via z, and there are x instances that z solves in less time. Note in particular that the times are *not* cumulative. The right-hand plot zooms into the left-hand one. We clearly see the speedup offered by OVI over SVI and especially II. Where the scatter plots merely show that OVI often does not obtain more than a 2× speedup compared to SVI, these plots provide an explanation: the VI line is a rough

**Fig. 7.** Summary comparison to VI, II, and SVI, instances ordered by runtime

**Fig. 8.** Influence of /α on runtime (expected rewards, relative error)

**Fig. 9.** Runtime comparison with absolute error (expected rewards)

bound on the performance that any *extension* of VI can deliver. Comparing the SVI and VI lines, over much of the plot's range, OVI thus cannot take less than half the runtime of SVI without outperforming VI itself.

#### **5.3 On the Effect of** **and** *α*

We also compared the four algorithms for different values of and, where applicable, α. We show a selection of the results in Fig. 8. The axis labels are of the form "algorithm, -/α". On the left, we see that the runtime of OVI changes if we set α to values different from -, however there is no clear trend: some instances are checked faster, some slower. We obtained similar plots for other combinations of α values, with only a slight tendency towards longer runtimes as α>-. mcsta thus uses α = as a default that can be changed by the user.

In the middle, we study the impact of reducing the desired precision by setting to 10−3. This allows OVI to speed up by factors mostly between 1 and 2; the same comparison for SVI and II resulted in similar plots, however VI was able to more consistently achieve higher speedups. When we compare the right plot with the right-hand plot of Fig. 6, we consequently see that the overall result of our comparison between OVI and SVI does not change significantly with the lower precision, although OVI does gain slightly more than SVI.

#### **5.4 Comparing Relative and Absolute Error**

In Fig. 9, we show comparison plots for the runtime when using *diffabs* instead of *diffrel* . Requiring absolute-error-correct results may make instances with low result values much easier and instances with high results much harder. We chose - = 10−2 as a compromise, and the leftmost plot confirms that we indeed chose an that keeps the expected-reward benchmarks on average roughly as hard as with 10−6 relative error. In the middle and right plots, we again see OVI compared with II and SVI. Compared to Fig. 6, both II and SVI gain a little, but there are no significant differences overall. Our experiments thus confirm that the relative performance of OVI is stable under varying precision requirements.

#### **5.5 Verification Phases**

On the right, we show histograms of the number of verification phases started (top, from 1 phase on the left to 20 on the right) and the percentage of iterations that are done in verification phases (bottom) over all benchmark instances (probabilities and rewards). We see that, in the vast majority of cases, we need few verification attempts, with many succeeding in the first attempt, and most iterations are performed in the iteration phases.

#### **6 Conclusion**

We have presented *optimistic value iteration* (OVI), a new approach to making non-exact probabilistic model checking via iterative numeric algorithms sound in the sense of delivering results within a prescribed interval around the true value (modulo floating-point and implementation errors). Compared to interval (II) and sound value iteration (SVI), OVI has slightly stronger termination guarantees in presence of multiple fixed points, and works in practice for max. probabilities without collapsing end components despite the lack of a guarantee. Like II, it can be combined with alternative methods for dealing with end components such as the new *deflating* technique of [23]. OVI is a *simple* algorithm that is *easy* to add to any tool that already implements value iteration, and it is *fast*, further closing the performance gap between VI and sound methods.

**Acknowledgments.** The authors thank Tim Quatmann (RWTH Aachen) for fruitful discussions when the idea of OVI initially came up in late 2018, and for his help in implementing and optimising the SVI implementation in mcsta.

**Data Availability.** A dataset to replicate our experimental evaluation is archived and available at DOI 10.4121/uuid:3df859e6-edc6-4e2d-92f3-93e478bbe8dc [19].

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

## PrIC3: Property Directed Reachability for MDPs

Kevin Batz1(B) , Sebastian Junges2 , Benjamin Lucien Kaminski3 , Joost-Pieter Katoen1 , Christoph Matheja4 , and Philipp Schröer1

	- 2 University of California, Berkeley, USA
	- 3 University College London, London, UK 4 ETH Zürich, Zürich, Switzerland

Abstract. IC3 has been a leap forward in symbolic model checking. This paper proposes PrIC3 (pronounced pricy-three), a conservative extension of IC3 to symbolic model checking of MDPs. Our main focus is to develop the theory underlying PrIC3. Alongside, we present a first implementation of PrIC3 including the key ingredients from IC3 such as generalization, repushing, and propagation.

#### 1 Introduction

IC3. Also known as property-directed reachability (PDR) [23], IC3 [13] is a symbolic approach for verifying finite transition systems (TSs) against safety properties like "*bad states are unreachable*". It combines bounded model checking (BMC) [12] and inductive invariant generation. Put shortly, IC3 either proves that a set B of bad states is *un*reachable by finding a set of non-B states closed under reachability—called an *inductive invariant*—or refutes reachability of B by a *counterexample* path reaching B. Rather than unrolling the transition relation (as in BMC), IC3 attempts to incrementally strengthen the invariant "no state in B is reachable" into an inductive one. In addition, it applies aggressive abstraction to the explored state space, so-called generalization [36]. These aspects together with the enormous advances in modern SAT solvers have led to IC3's success. IC3 has been extended [27,38] and adapted to software verification [19,44]. This paper develops a *quantitative* IC3 framework for probabilistic models.

*MDPs.* Markov decision processes (MDPs) extend TSs with discrete probabilistic choices. They are central in planning, AI as well as in modeling randomized distributed algorithms. A key question in verifying MDPs is *quantitative* reachability: "*is the (maximal) probability to reach* B *at most* λ*?*". Quantitative reachability [5,6]

This work has been supported by the ERC Advanced Grant 787914 (FRAPPANT), NSF grants 1545126 (VeHICaL) and 1646208, the DARPA Assured Autonomy program, Berkeley Deep Drive, and by Toyota under the iCyPhy center.

c The Author(s) 2020

S. K. Lahiri and C. Wang (Eds.): CAV 2020, LNCS 12225, pp. 512–538, 2020. https://doi.org/10.1007/978-3-030-53291-8\_27

reduces to solving linear programs (LPs). Various tools support MDP model checking, e.g., Prism [43], Storm [22], modest [34], and EPMC [31]. The LPs are mostly solved using (variants of) value iteration [8,28,35,51]. Symbolic BDD-based MDP model checking originated two decades ago [4] and is rather successful.

*Towards* IC3 *for MDPs.* Despite the success of BDD-based symbolic methods in tools like Prism, IC3 has not penetrated probabilistic model checking yet. The success of IC3 and the importance of quantitative reachability in probabilistic model checking raises the question *whether and how* IC3 *can be adapted—not just utilized—to reason about quantitative reachability in MDPs*. This paper addresses the challenges of answering this question. It extends IC3 in several dimensions to overcome these hurdles, making PrIC3—to our knowledge—*the first* IC3 *framework for quantitative reachability in MDPs*1. Notably, PrIC3 is conservative: For a threshold λ = 0, PrIC3 solves the same qualitative problem *and behaves (almost) the same as standard* IC3. Our main contribution is developing the theory underlying PrIC3, which is accompanied by a proof-of-concept implementation.

*Challenge 1 (Leaving the Boolean domain).* IC3 iteratively computes *frames*, which are over-approximations of sets of states that can reach B in a bounded number of steps. For MDPs, Boolean reachability becomes a *quantitative reachability probability*. This requires a shift: frames become real-valued functions rather than sets of states. Thus, there are infinitely many possible frames—even for finite-state MDPs—just as for infinite-state software [19,44] and hybrid systems [54]. Additionally, whereas in TSs a state reachable within k steps remains reachable on increasing k, the reachability probability in MDPs may increase. This complicates ensuring termination of an IC3 algorithm for MDPs. -

*Challenge 2 (Counterexamples* = *single paths).* For TSs, a single cycle-free path2 to B suffices to refute that "B *is not reachable*". This is not true in the probabilistic setting [32]. Instead, proving that the probability of reaching B exceeds the threshold λ requires *a set of possibly cyclic paths*—e.g., represented as a sub-MDP [15]—whose probability mass exceeds λ. Handling sets of paths as counterexamples in the context of IC3 is new. -

*Challenge 3 (Strengthening).* This key IC3 technique intuitively turns a proof obligation of type (i) "state s is unreachable from the initial state s*I* " into type (ii) "s's *predecessors* are unreachable from s*I* ". A first issue is that in the quantitative setting, the standard characterization of reachability probabilities in MDPs (the Bellman equations) inherently *reverses* the direction of reasoning (cf. "reverse" IC3 [53]): Hence, strengthening turns (i) "s cannot reach B" into (ii) "s's *successors* cannot reach B".

A much more challenging issue, however, is that in the quantitative setting obligations of type (i) read "s is reachable *with at most probability* δ". However,

1 Recently, (standard) IC3 for TSs was *utilized* in model checking Markov chains [49] to on-the-fly compute the states that cannot reach B.

2 In [38], tree-like counterexamples are used for non-linear predicate transformers in IC3.

the strengthened type (ii) obligation must then read: "*the weighted sum over the reachability probabilities of the successors of* s is at most δ". In general, there are infinitely many possible choices of subobligations for the successors of s in order to satisfy the original obligation, because—grossly simplified—there are infinitely many possibilities for a and b to satisfy weighted sums such as 1 3 a + 2 3 b ≤ δ. While we only need one choice of subobligations, picking a *good* one is approximately as hard as solving the entire problem altogether. We hence require a heuristic, which is guided by a *user-provided oracle*. -

*Challenge 4 (Generalization).* "One of the key components of IC3 is [inductive] generalization" [13]. Generalization [36] abstracts single states. It makes IC3 scale, but is *not* essential for correctness. To facilitate generalization, systems should be encoded symbolically, i.e., integer-valued program variables describe states. Frames thus map variables to probabilities. A first aspect is how to effectively present them to an SMT-solver. Conceptually, we use uninterpreted functions and universal quantifiers (encoding program behavior) together with linear real arithmetic to encode the weighted sums occurring when reasoning about probabilities. A second aspect is more fundamental: Abstractly, IC3's generalization guesses an unreachable set of states. We, however, need to guess this set *and* a probability for each state. To be effective, these guesses should moreover eventually yield an inductive frame, which is often highly nonlinear. We propose three SMT-guided interpolation variants for guessing these maps. -

*Structure of this Paper.* We develop PrIC3 gradually: We explain the underlying rationale in Sect. 3. We also describe the core of PrIC3—called PrIC3H—which resembles closely the main loop of standard IC3, but uses adapted frames and termination criteria (Challenge 1). In line with Challenge 3, PrIC3H is parameterized by a heuristic H which is applied whenever we need to select one out of infinitely many probabilities. No requirements on the quality of H are imposed. PrIC3H is *sound* and always terminates: If it returns true, then the maximal reachability probability is bounded by λ. Without additional assumptions about H, PrIC3H is *incomplete*: on returning false, it is unknown whether the returned sub-MDP is indeed a counterexample (Challenge 2). Section 4 details strengthening (Challenge 3). Section 5 presents a sound *and* complete algorithm PrIC3 on top of PrIC3H. Section 6 presents a prototype, discusses our chosen heuristics, and addresses Challenge 4. Section 7 shows some encouraging experiments, but also illustrates need for further progress.

Related Work. Just like IC3 has been a symbiosis of different approaches, PrIC3 has been inspired by several existing techniques from the verification of probabilistic systems.

*BMC.* Adaptions of BMC to Markov chains (MCs) with a dedicated treatment of cycles have been pursued in [57]. The encoding in [24] annotates sub-formulae with probabilities. The integrated SAT solving process implicitly unrolls all paths leading to an exponential blow-up. In [52], this is circumvented by grouping paths, discretizing them, and using an encoding with quantifiers and bit-vectors, but without numerical values. Recently, [56] extends this idea to a PAC algorithm by purely propositional encodings and (approximate) model counting [17]. These approaches focus on MCs and are not mature yet.

*Invariant Synthesis.* Quantitative loop invariants are key in analyzing *probabilistic programs* whose operational semantics are (possibly infinite) MDPs [26]. A quantitative invariant I maps states to probabilities. I is shown to be an invariant by comparing I to the result of applying the MDP's Bellman operator to I. Existing approaches for invariant synthesis are, e.g., based on weakest preexpectations [33,39,40,42,46], template-based constraint solving [25], notions of martingales [3,9,16,55], and solving recurrence relations [10]. All but the last technique require user guidance.

*Abstraction.* To combat state-space explosion, abstraction is often employed. CEGAR for MDPs [37] deals with explicit sets of paths as counterexamples. Game-based abstraction [30,41] and partial exploration [14] exploit that not all paths have to be explored to prove bounds on reachability probabilities.

*Statistical Methods and (deep) Reinforcement Learning.* Finally, an avenue that avoids storing a (complete) model are simulation-based approaches (statistical model checking [2]) and variants of reinforcement learning, possibly with neural networks. For MDPs, these approaches yield weak statistical guarantees [20], but may provide good oracles.

#### 2 Problem Statement

Our aim is to prove that the *maximal probability* of *reaching* a *set* B *of bad states* from the initial state s*I* of a *Markov decision process* M is at most some *threshold* λ.  We refer to [7,50] for a thorough introduction.

Definition 1 (MDPs). *A* Markov decision process *(* MDP*) is a tuple* M = (S, s*I* , Act, P)*, where* S *is a finite set of* states*,* s*I* ∈ S *is the* initial state*,* Act *is a finite set of* actions*, and* P : S × Act × S → [0, 1] *is a* transition probability function*. For state* s*, let* Act (s) = {a ∈ Act | ∃s ∈ S : P(s, a, s ) > 0} *be the* - enabled actions *at* s*. For all states* s ∈ S*, we require* |Act (s)| ≥ 1 *and s*-∈*S* P(s, a, s )=1*.* -

For this paper, we fix an MDP M = (S, s*I* , Act, P), a set of *bad states* B ⊆ S, and a threshold λ ∈ [0, 1]. The *maximal* 3 *(unbounded) reachability probability* to eventually reach a state in B from a state s is denoted by Prmax (s |= ♦B). We characterize Prmax (s |= ♦B) using the so-called *Bellman operator*. Let M*N* denote the set of functions from N to M. Anticipating IC3 terminology, we call a function F ∈ [0, 1]*S* a *frame*. We denote by F[s] the evaluation of frame F for state s.

3 Maximal with respect to all possible resolutions of nondeterminism in the MDP.

Definition 2 (Bellman Operator). *For a set of actions* A ⊆ Act*, we define the* Bellman operator for A *as a frame transformer* Φ*A* : [0, 1]*S* → [0, 1]*S with*

$$\Phi\_A\left(F\right)[s] \ = \begin{cases} 1, & if\ s \in B\\ \max\limits\_{a \in A} \sum\_{s' \in S} P(s, a, s') \cdot F[s']\ , & if\ s \notin B\ . \end{cases}$$

*We write* Φ*a for* Φ{*a*}*,* Φ *for* ΦAct*, and call* Φ *simply* the Bellman operator*.* -

For every state s, the maximal reachability probability Prmax (s |= ♦B) is then given by the least fixed point of the Bellman operator Φ. That is,

$$\forall s: \quad \text{Pr}^{\text{max}}\left(s \mid \vdash \Diamond B\right) \; = \; \left(\text{lfp } \Phi\right)[s] \; ,$$

where the underlying partial order on frames is a complete lattice with ordering

$$\begin{array}{ccccc} F\_1 & \leq & F\_2 & \text{iff} & & \forall s \in S \colon & & F\_1[s] & \leq & F\_2[s] & . \end{array}$$

In terms of the Bellman operator, our formal problem statement reads as follows:

Whenever Prmax (s*I* |= ♦B) ≤ λ indeed holds, we say that the MDP M is *safe* (with respect to the set of bad states B and threshold λ); otherwise, we call it *unsafe*.

Fig. 1. The MDP M serving as a running example.

Recovery Statement 1. *For* λ = 0*, our problem statement is equivalent to the* qualitative reachability *problem solved by (reverse) standard* IC3*, i.e, prove or refute that all bad states in* B *are* unreachable *from the initial state* s*I .*

*Example 1.* The MDP M in Fig. 1 consists of 6 states with initial state s0 and bad states B = {s5}. In s2, actions a and b are enabled; in all other states, one unlabeled action is enabled. We have Prmax (s0 |= ♦B) = 2/3. Hence, M is safe for all thresholds λ ≥ 2/3 and unsafe for λ < 2/3. In particular, M is unsafe for λ = 0 as s5 is *reachable* from s0. -

#### 3 The Core PrIC3 Algorithm

The purpose of PrIC3 is to prove or refute that the maximal probability to reach a bad state in B from the initial state s*I* of the MDP M is at most λ. In this section, we explain the rationale underlying PrIC3. Moreover, we describe the core of PrIC3—called PrIC3H—which bears close resemblance to the main loop of standard IC3 for TSs.

Because of the inherent direction of the Bellman operator, we build PrIC3 on *reverse* IC3 [53], cf. Challenge 3. Reversing constitutes a shift from reasoning along the direction *initial-to-bad* to *bad-to-initial*. While this shift is mostly *inessential* to the fundamentals underlying IC3, the reverse direction is unswayable in the probabilistic setting. Whenever we draw a connection to standard IC3, we thus generally mean *reverse* IC3.

#### 3.1 Inductive Frames

IC3 for TSs operates on (*quali*tative) frames representing sets of states of the TS at hand. A frame F can hence be thought of as a mapping4 from states to {0, 1}. In PrIC3 for MDPs, we need to move from a Boolean to a quantitative regime. Hence, a (*quanti*tative) frame is a mapping from states to probabilities in [0, 1].

For a given TS, consider the frame transformer T that adds to a given input frame F all bad states in B and all predecessors of the states contained in F . The rationale of standard (reverse) IC3 is to find a frame F ∈ {0, 1}*S* such that (i) the initial state s*I* does not belong to F and (ii) applying T takes us down in the partial order on frames, i.e.,

$$\text{(1)}\quad F[s\_I] \quad = \text{ 0} \qquad \text{and} \qquad \text{(ii)}\quad T(F) \le \text{ }F\text{ }\dots$$

Intuitively, (i) postulates the *hypothesis* that s*I* cannot reach B and (ii) expresses that F is closed under adding bad states and taking predecessors, thus affirming the hypothesis.

Analogously, the rationale of PrIC3 is to find a frame F ∈ [0, 1]*S* such that (i) F postulates that the probability of s*I* to reach B is at most the threshold λ and (ii) applying the Bellman operator Φ to F takes us down in the partial order on frames, i.e.,

(i) F[s*I* ] ≤ λ and (ii) Φ(F) ≤ F .

Frames satisfying the above conditions are called *inductive invariants* in IC3. We adopt this terminology. By *Park's Lemma* [48], which in our setting reads

Φ(F) ≤ F implies lfp Φ ≤ F ,

4 In IC3, frames are typically characterized by logical formulae. To understand IC3's fundamental principle, however, we prefer to think of frames as functions in {0, 1}S partially ordered by ≤.

an inductive invariant F would indeed *witness* that Prmax (s*I* |= ♦B) ≤ λ, because

$$\Pr^{\text{max}}\left(s\_I \mid \! \! \! / B\right) \; = \; \left(\text{lfp } \Phi\right)[s\_I] \; \le \; \; F[s\_I] \; \le \; \lambda \; .$$

If no inductive invariant exists, then standard IC3 will find a counterexample: a *path* from the initial state s*I* to a bad state in B, which serves as a witness to refute. Analogously, PrIC3 will find a counterexample, but of a different kind: Since single paths are insufficient as counterexamples in the probabilistic realm (Challenge 2), PrIC3 will instead find a *subsystem* of states of the MDP witnessing Prmax (s*I* |= ♦B) > λ.

#### 3.2 The PrIC3 Invariants

Analogously to standard IC3, PrIC3 aims to find the inductive invariant by maintaining a *sequence of frames* F0 ≤ F1 ≤ F2 ≤ ... such that F*i*[s] overapproximates the maximal probability of reaching B from s within *at most* i *steps*. This i*-step-bounded reachability probability* Prmax s |= ♦≤*i* B can be characterized using the Bellman operator: Φ (**0**) is the 0-step probability; it is 1 for every s ∈ B and 0 otherwise. For any i ≥ 0, we have

$$\Pr^{\text{max}}\left(s \mid \stackrel{<}{\rightharpoonup} \stackrel{<}{B}\right) \\ = \left(\Phi^i\left(\Phi\left(\mathbf{0}\right)\right)\right)[s] \\ = \left(\Phi^{i+1}\left(\mathbf{0}\right)\right)[s],$$

where **0**, the frame that maps every state to 0, is the least frame of the underlying complete lattice. For a finite MDP, the *unbounded* reachability probability is then given by the limit

$$\Pr^{\text{max}}\left(s\vDash\!\!\langle B\rangle\right) = \left(\text{lfp }\Phi\right)[s] \stackrel{(\*)}{=} \left(\lim\_{n\to\infty} \Phi^n\left(\mathbf{0}\right)\right)[s] \\ = \lim\_{n\to\infty} \Pr^{\text{max}}\left(s\vDash\!\!\!\langle \!\*\!\!\langle \!\*\!\!\!\langle \!\*\!\!\langle \!\*\!\!\langle \!\*\!\!\langle \!\*\!\!\langle \!\*\!\langle \!\*\!\langle \!\*\!\langle \!\*\!\langle\!\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\langle\!\/\vert\!\/)\!\/)\!\rangle\!\rangle\!\rangle}\!\rangle\!\rangle]}\right)$$

where (∗) is a consequence of the well-known Kleene fixed point theorem [45].

The sequence F0 ≤ F1 ≤ F2 ≤ ... maintained by PrIC3 should frame-wise overapproximate the increasing sequence Φ (**0**) ≤ Φ2 (**0**) ≤ Φ3 (**0**).... Pictorially:

F0 ≤ F1 ≤ F2 ≤ ... ≤ F*k* ≤ ≤ ≤ ≤

$$\begin{array}{ccccccccc} \mathbf{0} & \leq & \Phi(\mathbf{0}) & \leq & \Phi^2(\mathbf{0}) & \leq & \Phi^3(\mathbf{0}) & \leq & \dots & \leq & \Phi^{k+1}(\mathbf{0}) \end{array}$$

However, the sequence Φ (**0**), Φ2 (**0**), Φ3 (**0**), ... will never explicitly be known to PrIC3. Instead, PrIC3 will ensure the above frame-wise overapproximation property implicitly by enforcing the so-called PrIC3 *invariants* on the frame sequence F0, F1, F2, .... Apart from allowing for a threshold 0 ≤ λ ≤ 1 on the maximal reachability probability, these invariants coincide with the standard IC3 invariants (where λ = 0 is fixed). Formally:

Definition 3 (PrIC3 Invariants). *Frames* F0, ..., F*k, for* k ≥ 0*, satisfy the* PrIC3 invariants*, a fact we will denote by* PrIC3Inv (F0, ..., F*k*)*, if all of the following hold:*

The PrIC3 invariants enforce the above picture: The *chain property* ensures F0 ≤ F1 ≤ ... ≤ F*k*. We have Φ (**0**) = F0 ≤ F0 by *initiality*. Assuming Φ*i*+1 (**0**) ≤ F*i* as induction hypothesis, monotonicity of Φ and *relative inductivity* imply Φ*i*+2 (**0**) ≤ Φ(F*i*) ≤ F*i*+1.

By overapproximating Φ (**0**), Φ2 (**0**), ..., Φ*k*+1 (**0**), the frames F0, ..., F*k* in effect bound the maximal step-bounded reachability probability of every state:

Lemma 1. *Let frames* F0, ..., F*k satisfy the* PrIC3 *invariants. Then*

∀ s ∀ i ≤ k : *Prmax* s |= ♦≤*i* B ≤ F*i*[s].

In particular, Lemma 1 together with *frame-safety* ensures that the maximal step-bounded reachability probability of the *initial state* s*I* to reach B is at most the threshold λ.

As for proving that the *unbounded* reachability probability is also at most λ, it suffices to find two consecutive frames, say F*i* and F*i*+1, that coincide:

Lemma 2. *Let frames* F0, ..., F*k satisfy the* PrIC3 *invariants. Then*

$$\exists \, i < k \colon \quad F\_i \; = \; F\_{i+1} \qquad implies \qquad Pr^{\max}(s\_I \mid \neg \Diamond B) \; \le \; \lambda \; .$$

*Proof.* F*i* = F*i*+1 and *relative inductivity* yield Φ(F*i*) ≤ F*i*+1 = F*i*, rendering F*i inductive*. By Park's lemma (cf. Sect. 3.1), we obtain lfp Φ ≤ F*i* and—by *frame-safety*—conclude

$$\Pr^{\max}(s\_I \mid \mho\_I \Phi) \; = \; \text{(lfp } \Phi) \\
[s\_I] \; \le \; F\_i[s\_I] \; \le \; \lambda \; . \tag{7}$$

#### 3.3 Operationalizing the PrIC3 Invariants for Proving Safety

Lemma 2 gives us a clear angle of attack for *proving* an MDP safe: Repeatedly add and refine frames approximating step-bounded reachability probabilities for more and more steps while enforcing the PrIC3 invariants (cf. Definition 3.2) until two consecutive frames coincide.

Analogously to standard IC3, this approach is taken by the core loop PrIC3H depicted in Algorithm 1; differences to the main loop of IC3 (cf. [23, Fig. 5]) are highlighted in red. A particular difference is that PrIC3H is parameterized by a heuristic H for finding suitable probabilities (see Challenge 3). Since the precise choice of H is irrelevant for the soundness of PrIC3H, we defer a detailed discussion of suitable heuristics to Sect. 4.

```
Data: MDP M, set of bad states B, threshold λ
  Result: true or false and a subset of the states of M
1 F0 ← Φ (0); F1 ← 1; k ← 1; oldSubsystem ← ∅;
2 while true do
3 success, F0,...,Fk, subsystem ← StrengthenH (F0, ..., Fk);
4 if ¬success then returnfalse, subsystem;
5 Fk+1 ← 1;
6 F0,...,Fk+1 ← Propagate (F0,...,Fk+1);
7 if ∃ 1 ≤ i ≤ k : Fi = Fi+1 then returntrue, ;
8 if oldSubsystem = subsystem then returnfalse, subsystem;
9 k ← k + 1; oldSubsystem ← subsystem;
10 end
```
Algorithm 1: PrIC3H (M, B, λ)

As input, PrIC3H takes an MDP M = (S, s*I* , Act, P), a set B ⊆ S of bad states, and a threshold λ ∈ [0, 1]. Since the input is never changed, we assume it to be *globally available*, also to subroutines. As output, PrIC3H returns true if two consecutive frames become equal. We hence say that PrIC3H is *sound* if it only returns true if M is safe.

We will formalize soundness using Hoare triples. For precondition φ, postcondition ψ, and program P, the triple φ P ψ is *valid* (for partial correctness) if, whenever program P starts in a state satisfying precondition φ and terminates in some state s , then s satisfies postcondition ψ. Soundness of PrIC3H then means validity of the triple

$$\{\text{true}\} \text{safe}, \underline{\hspace{0.1cm}} \leftarrow \text{PrlC3}\_{\mathcal{H}}(\mathfrak{M}, B, \lambda) \{\text{safe} \Rightarrow \Pr^{\text{max}} \left(s\_I \mid \!= \lozenge{0}B\right) \leq \lambda\} \;.$$

Let us briefly go through the individual steps of PrIC3H in Algorithm 1 and convince ourselves that it is indeed sound. After that, we discuss why PrIC3H terminates and what happens if it is unable to prove safety by finding two equal consecutive frames.

How **PrIC3***H* works. Recall that PrIC3H maintains a sequence of frames F0, ..., F*k* which is initialized in l. 1 with k = 1, F0 =Φ(**0**), and F1 = **1**, where the frame **1** maps every state to 1. Every time upon entering the while-loop in terms l. 2, the initial segment F0, ..., F*k*−1 satisfies all PrIC3 invariants (cf. Definition 3), whereas the full sequence F0, ..., F*k* potentially violates frame-safety as it is possible that F*k*[s*I* ] > λ.

In l. 3, procedure StrengthenH—detailed in Sect. 4—is called to restore *all* PrIC3 invariants on the *entire* frame sequence: It either returns true if successful or returns false and a counterexample (in our case a subsystem of the MDP) if it was unable to do so. To ensure soundness of PrIC3H, it suffices that StrengthenH restores the PrIC3 invariants whenever it returns true. Formally, StrengthenH must meet the following specification:

Definition 4. *Procedure* StrengthenH *is* sound *if the following Hoare triple is valid:*

$$\begin{aligned} \{ \begin{aligned} \text{\textquotedblleft} \mathsf{PrlC3lnv} \left( F\_0, \ldots, F\_{k-1} \right) & \land \ F\_{k-1} \le F\_k \land \Phi \left( F\_{k-1} \right) \le F\_k \right\} \\ \text{\textit{subclosed}}, \, F\_0, \ldots, F\_k, \underline{\hspace{0.5cm}} & \leftarrow \mathsf{Strengthen}\_{\mathcal{H}} \left( F\_0, \ldots, F\_k \right) \\ \{ \begin{aligned} \text{\textquotedblleft} \text{\textquotedblright} \text{\textquotedblright} \left( F\_0, \ldots, F\_k \right) \end{aligned}} \end{aligned} \} $$

If StrengthenH returns true, then a new frame F*k*+1 = **1** is created in l. 5. After that, the (now initial) segment F0, ..., F*k* again satisfies all PrIC3 invariants, whereas the full sequence F0, ..., F*k*+1 potentially violates frame-safety at F*k*+1. *Propagation* (l. 6) aims to speed up termination by updating F*i*+1[s] by F*i*[s] iff this does not violate relative inductivity. Consequently, the previously mentioned properties remain unchanged.

If StrengthenH returns false, the PrIC3 invariants—premises to Lemma 2 for witnessing safety—cannot be restored and PrIC3H terminates returning false (l. 4). Returning false (also possible in l. 8) has by specification no affect on soundness of PrIC3H.

In l. 7, we check whether there exist two identical consecutive frames. If so, Lemma 2 yields that the MDP is safe; consequently, PrIC3H returns true. Otherwise, we increment k and are in the same setting as upon entering the loop, now with an increased frame sequence; PrIC3H then performs another iteration. In summary, we obtain:

Theorem 1 (Soundness of PrIC3*H*). *If* StrengthenH *is sound and* Propagate *does not affect the* PrIC3 *invariants, then* PrIC3H *is sound, i.e., the following triple is valid:*

 true *safe*, ← PrIC3H (M, B, λ) *safe* =⇒ *Prmax* (s*I* |= ♦B) ≤ λ

**PrIC3***H* Terminates for Unsafe MDPs. If the MDP is unsafe, then there exists a step-bound n, such that Prmax s*I* |= ♦≤*n*B > λ. Furthermore, any sound implementation of StrengthenH (cf. Definition 4) either immediately terminates PrIC3H by returning false or restores the PrIC3 invariants for F0, ..., F*k*. If the former case never arises, then StrengthenH will eventually restore the PrIC3 invariants for a frame sequence of length k = n. By Lemma 1, we have F*n*[s*I* ] ≥ Prmax s*I* |= ♦≤*n*B > λ contradicting frame-safety.

**PrIC3***H* Terminates for Safe MDPs. Standard IC3 terminates on safe finite TSs as there are only finitely many different frames, making every ascending chain of frames eventually stabilize. For us, frames map states to probabilities (Challenge 1), yielding *infinitely many possible frames* even for finite MDPs. Hence, StrengthenH need not ever yield a stabilizing chain of frames. If it continuously fails to stabilize while repeatedly reasoning about the same set of states, we give up. PrIC3H checks this by comparing the subsystem StrengthenH operates on with the one it operated on in the previous loop iteration (l. 8).

Theorem 2. *If* StrengthenH *and* Propagate *terminate, then* PrIC3H *terminates.*

Recovery Statement 2. *For qual. reachability (*λ = 0*),* PrIC3H *never terminates in l. 8.*

**PrIC3***H* is Incomplete. Standard IC3 either proves safety or returns false and a counterexample—a single path from the initial to a bad state. As single paths are insufficient as counterexamples in MDPs (Challenge 2), PrIC3H instead returns a *subsystem* of the MDP M provided by StrengthenH. However, as argued above, we cannot trust StrengthenH to provide a stabilizing chain of frames. Reporting false thus only means that the given MDP *may* be unsafe; the returned subsystem has to be analyzed further.

The full PrIC3 algorithm presented in Sect. 5 addresses this issue. Exploiting the subsystem returned by PrIC3H, PrIC3 returns true if the MDP is safe; otherwise, it returns false and provides a true counterexample witnessing that the MDP is unsafe.

*Example 2.* We conclude this section with two example executions of PrIC3H on a simplified version of the MDP in Fig. 1. Assume that action b has been removed. Then, for every state, exactly one action is enabled, i.e., we consider a Markov chain. Figure 2 depicts the frame sequences computed by PrIC3H (for a reasonable H) on that Markov chain for two thresholds: 5/9 = Prmax (s0 |= ♦B) and 9/10. In particular, notice that *proving the coarser bound of* 9/10 *requires fewer frames than proving the exact bound of* 5/9. -

Fig. 2. Two runs of PrIC3H on the Markov chain induced by selecting action a in Fig. 1. For every iteration, frames are recorded after invocation of StrengthenH.

## 4 Strengthening in **PrIC3***H*

When the main loop of PrIC3H has created a new frame F*k* = **1** in its previous iteration, this frame may violate frame-safety (Definition 3.3) because of F*k*[s*I* ] = 1 ≤ λ. The task of StrengthenH is to restore the PrIC3 invariants on *all* frames F0,...,F*k*. To this end, our first *obligation* is to lower the value in frame i = k for state s = s*I* to δ = λ ∈ [0, 1]. We denote such an obligation by (i, s, δ). Observe that implicitly δ = 0 in the qualitative case, i.e., when proving unreachability. An obligation (i, s, δ) is *resolved* by updating the values assigned to state s in *all frames* F1,...,F*i* to at most δ. That is, for all j ≤ i, we set F*j* [s] to the minimum

```
1 Q ← {(k, sI , λ)} ;
2 while Q not empty do
3 (i, s, δ) ← Q.popMin(); /* pop obligation with minimal frame
      index */
4 if i = 0 ∨ (s ∈ B ∧ δ < 1) then
         /* possible counterexample given by subsystem
            consisting of states popped from Q at some point */
5 return false, , Q.touched();
      /* check whether Fi[s] ← δ violates relative inductivity */
6 if ∃a ∈ Act (s):Φa (Fi−1) [s] > δ then for such an a
7 δ1,...,δn ← H (s, a, δ) ;
8 {s1,...,sn} ← Succs(s, a);
9 Q.push ((i − 1, s1, δ1) ,..., (i − 1, sn, δn) , (i, s, δ));
10 else /* resolve (i, s, δ) without violating relative
      inductivity */
11 F1[s] ← min (F1[s], δ) ; ... ; Fi[s] ← min (Fi[s], δ);
12 end
13 (/* Q empty; all obligations have been resolved */ ) return
  true, F0,...,Fk, Q.touched();
                Algorithm 2: StrengthenH (F0,...,Fk)
```
of δ and the original value F*j* [s]. Such an update affects neither initiality nor the chain property (Definitions 3.1, 3.2). It may, however, violate relative inductivity (Definition 3.4), i.e., Φ (F*i*−1) ≤ F*i*. Before resolving obligation (i, s, δ), we may thus have to further decrease some entries in F*i*−1 as well. Hence, *resolving obligations may spawn additional obligations* which have to be resolved first to maintain relative inductivity. In this section, we present a generic instance of StrengthenH meeting its specification (Definition 4) and discuss its correctness.

**Strengthen***H* by Example. StrengthenH is given by the pseudo code in Algorithm 2; differences to standard IC3 (cf. [23, Fig. 6]) are highlighted in red. Intuitively, StrengthenH attempts to recursively resolve all obligations until either both frame-safety and relative inductivity are restored for *all* frames or it detects a *potential counterexample* justifying why it is unable to do so. We first consider an execution where the latter does not arise:

*Example 3.* We zoom in on Example 2: Prior to the second iteration, we have created the following three frames assigning values to the states s0, s5:

$$F\_0 = \langle 0, 0, 0, 0, 1 \rangle, \qquad F\_1 = \langle 5/9, 1, 1, 1, 1 \rangle, \qquad \text{and} \qquad F\_2 = \mathbf{1}.$$

To keep track of unresolved obligations (i, s, δ), StrengthenH employs a priority queue Q which pops obligations with minimal frame index i first. Our first step is to ensure frame-safety of F2, i.e., alter F2 so that F2[s0] ≤ 5/9; we thus initialize the queue Q with the initial obligation (2, s0, 5/9) (l. 1). To do so, we check whether updating F2[s0] to 5/9 would invalidate relative inductivity (l. 6). This is indeed the case:

$$\Phi(F\_1)\begin{bmatrix} s\_0 \end{bmatrix} = \,^1/2 \cdot F\_1[s\_1] + ^1/2 \cdot F\_1[s\_2] \,^1 = \,^1 \not\le ^5/9.$$

To restore relative inductivity, StrengthenH spawns one new obligation for each relevant successor of s0. These have to be resolved before retrying to resolve the old obligation.5

*In contrast to standard* IC3 *, spawning obligations involves finding suitable probabilities* δ (l. 7). In our example this means we have to spawn two obligations (1, s1, δ1) and (1, s2, δ2) such that 1/2 · δ1 + 1/2 · δ2 ≤ 5/9. There are *infinitely many choices* for δ1 and δ2 satisfying this inequality. Assume some heuristic H chooses δ1 = 11/18 and δ2 = 1/2; we push obligations (1, s1, 11/18), (1, s2, 1/2), and (2, s0, 5/9) (ll. 8, 9). In the next iteration, we first pop obligation (1, s1, 11/18) (l. 3) and find that it can be resolved without violating relative inductivity (l. 6). Hence, we set F1[s1] to 11/18 (l. 11); no new obligation is spawned. Obligation (1, s2, 1/2) is resolved analogously; the updated frame is F1 = (5/9, 11/18, 1/2, 1). Thereafter, our initial obligation (2, s0, 5/9) can be resolved; relative inductivity is restored for F0, F1, F2. Hence, StrengthenH returns true together with the updated frames. -

**Strengthen***H* is Sound. Let us briefly discuss why Algorithm 2 meets the specification of a sound implemenation of StrengthenH (Definition 4): First, we observe that Algorithm 2 alters the frames—and thus potentially invalidates the PrIC3 invariants—only in l. 11 by resolving an obligation (i, s, δ) with Φ (F*i*−1) [s] ≤ δ (due to the check in l. 6).

Let F s → δ denote the frame F in which F[s] is set to δ, i.e.,

$$F\left[s'] = \begin{cases} \delta, & \text{if } s'=s, \\ F[s'], & \text{otherwise}. \end{cases}$$

Indeed, resolving obligation (i, s, δ) in l. 11 lowers the values assigned to state s to at most δ *without* invalidating the PrIC3 invariants:

Lemma 3. *Let* (i, s, δ) *be an obligation and* F0,...,F*i, for* i > 0*, be frames with* Φ (F*i*−1) [s] ≤ δ*. Then* PrIC3Inv (F0,...,F*i*) *implies*

$$\text{Prl\u{C}3\u{\text{inv}}}\left(\left.F\_0\left,\ldots,\left.F\_i\left\right).$$

Crucially, the precondition of Definition 4 guarantees that all PrIC3 invariants except frame safety hold initially. Since these invariants are never invalidated due to Lemma 3, Algorithm 2 is a sound implementation of StrengthenH if it restores frame safety whenever it returns true, i.e., once it leaves the loop with an empty obligation queue Q (ll. 12–13). Now, an obligation (i, s, δ) is only popped from Q in l. 3. As (i, s, δ) is added to Q upon reaching l. 9, the size of Q can only ever be reduced (without returning false) by resolving (i, s, δ) in l. 11. Hence, Algorithm 2 does not return true unless it restored frame safety by resolving, amongst all other obligations, the initial obligation (k, s*I* , λ). Consequently:

5 We assume that the set Succs(s, a) = {s ∈ S | P(s, a, s ) > 0} of *relevant* a*-successors* of state s is returned in some arbitrary, but fixed order.

Lemma 4. *Procedure* StrengthenH *is sound, i.e., it satisfies the specification in Definition 4.*

Theorem 3. *Procedure* PrIC3H *is sound, i.e., satisfies the specification in Theorem 1.*

We remark that, analogously to standard IC3, resolving an obligation in l. 11 may be accompanied by *generalization*. That is, we attempt to update the values of multiple states at once. Generalization is, however, highly non-trivial in a probabilistic setting. We discuss three possible approaches to generalization in Sect. 6.2.

**Strengthen***H* Terminates. We now show that StrengthenH as in Algorithm 2 terminates. The only scenario in which StrengthenH may not terminate is if it keeps spawning obligations in l. 9. Let us thus look closer at how obligations are spawned: Whenever we detect that resolving an obligation (i, s, δ) would violate relative inductivity for some action a (l. 6), we first need to update the values of the successor states s1,...,s*n* ∈ Succs(s, a) in frame i−1, i.e., we push the obligations (i−1, s1, δ1),...,(i−1, s*n*, δ*n*) which have to be resolved first (ll. 7–9). It is noteworthy that, for a TS, a single action leads to a single successor state s1. Algorithm 2 employs a heuristic H to determine the probabilities required for pushing obligations (l. 7). Assume for an obligation (i, s, δ) that the check in l. 6 yields ∃a ∈ Act (s):Φ*a* (F*i*−1) [s] > δ. Then H takes s, a, δ and reports some probability δ*j* for every a-successor s*j* of s. However, an arbitrary heuristic of type H: S × Act × [0, 1] → [0, 1]∗ may lead to non-terminating behavior: If δ1,...,δ*n* = F*i*−1[s1],...F*i*−1[s*n*], then the heuristic has no effect. It is thus natural to require that an *adequate* heuristic H yields probabilities such that the check Φ*a* (F*i*−1) [s] > δ in l. 6 cannot succeed twice for the *same obligation* (i, s, δ) and *same action* a. Formally, this is guaranteed by the following:

Definition 5. *Heuristic* H *is* adequate *if the following triple is valid (for any frame* F*):*

$$\begin{aligned} \left\{ \begin{aligned} &Success(s,a) = s\_1, \dots, s\_n \end{aligned} \right\} \\ \left\{ \begin{aligned} &\delta\_1, \dots, \delta\_n \leftarrow \mathcal{H}(s,a,\delta) \\ &\{ \Phi\_a \left( F \left< s\_1 \mapsto \delta\_1 \right> \dots \left< s\_n \mapsto \delta\_n \right> \right) [s] \le \delta \end{aligned} \right\} \end{aligned}$$

Details regarding our implementation of heuristic H are found in Sect. 6.1.

For an adequate heuristic, attempting to resolve an obligation (i, s, δ) (ll. 3 – 11) either succeeds after spawning it at most |Act(s)| times or StrengthenH returns false. By a similar argument, attempting to resolve an obligation (i > 0, s, ) leads to at most - *a*∈Act(*s*) |{s ∈ S | P(s, a, s ) > 0}| other obligations of the form (i−1, s , ). Consequently, the total number of obligations spawned by Algorithm 2 is bounded. Since Algorithm 2 terminates if all obligations have been resolved (l. 12) and each of its loop iterations either returns false, spawns obligations, or resolves an obligation, we conclude:

Lemma 5. StrengthenH(F0, ..., F*k*) *terminates for every adequate heuristic* H*.*

Recovery Statement 3. *Let* H *be adequate. Then for qualitative reachability (*λ = 0*), all obligations spawned by* StrengthenH *as in Algorithm 2 are of the form* (i, s, 0)*.*

**Strengthen***H* returns **false**. There are two cases in which StrengthenH fails to restore the PrIC3 invariants and returns false. The first case (the left disjunct of l. 4) is that we encounter an obligation for frame F0. Resolving such an obligation would inevitably violate *initiality*; analogously to standard IC3, we thus return false.

The second case (the right disjunct of l. 4) is that we encounter an obligation (i, s, δ) for a bad state s ∈ B with a probability δ < 1 (though, obviously, all s ∈ B have probability =1). Resolving such an obligation would inevitably prevents us from restoring *relative inductivity*: If we updated F*i*[s] to δ, we would have Φ (F*i*−1) [s]=1 > δ = F*i*[s]. Notice that, in contrast to standard IC3, this second case *can* occur in PrIC3:

*Example 4.* Assume we have to resolve an obligation (i, s3, 1/2) for the MDP in Fig. 1. This involves spawning obligations (i−1, s4, δ1) and (i−1, s5, δ2), where s5 is a bad state, such that 1/3 · δ1 + 2/3 · δ2 ≤ 1/2. Even for δ1 = 0, this is only possible if δ2 ≤ 3/4 < 1. -

**Strengthen***H* Cannot Prove Unsafety. If standard IC3 returns false, it proves unsafety by constructing a counterexample, i.e., *a single path from the initial state to a bad state*. If PrIC3 returns false, there are two possible reasons: *Either* the MDP is indeed unsafe, *or* the heuristic H at some point selected probabilities in a way such that StrengthenH is unable to restore the PrIC3 invariants (even though the MDP might in fact be safe). StrengthenH thus only returns a *potential* counterexample which either proves unsafety or indicates that our heuristic was inappropriate.

Counterexamples in our case consist of subsystems rather than a single path (see Challenge 2 and Sect. 5). StrengthenH hence returns the set Q.touched() of all states that eventually appeared in the obligation queue. This set is a conservative approximation, and optimizations as in [1] may be beneficial. Furthermore, in the qualitative case, our potential counterexample subsumes the counterexamples constructed by standard IC3:

Recovery Statement 4. *Let* H0 *be the adequate heuristic mapping every state to* 0*. For qual. reachability (*λ = 0*), if success* = false *is returned by* StrengthenH0 (F0, ..., F*k*)*, then* Q.*touched*() *contains a path from the initial to a bad state.*6

6 Q.touched() might be restricted to only contain this path by some simple adaptions.

Data: global MDP M, set of bad states B, threshold λ Result: true iff Prmax (sI |= ♦B) ≤ λ Ω ← Initialize(); *touched* ← {sI }; 2 do H ← CreateHeuristic(Ω); *safe*, *subsystem* ← PrIC3H(); if *safe* then return true ; if CheckRefutation(*subsystem*) then return false ; *touched* ← Enlarge(*touched*, *subsystem*); Ω ← Refine(Ω, *touched*); while *touched* = S; return Ω(sI ) ≤ λ

Algorithm 3: PrIC3: The outermost loop dealing with possibly imprecise heuristics

#### 5 Dealing with Potential Counterexamples

Recall that our core algorithm PrIC3H is incomplete for a fixed heuristic H: It cannot give a conclusive answer whenever it finds a potential counterexample for two possible reasons: Either the heuristic H turned out to be inappropriate or the MDP is indeed unsafe. The idea to overcome the former is to call PrIC3H finitely often in an outer loop that generates new heuristics until we find an appropriate one: If PrIC3H still does not report safety of the MDP, then it is indeed unsafe. We do not blindly generate new heuristics, but use the potential counterexamples returned by PrIC3H to refine the previous one.

Let consider the procedure PrIC3 in Algorithm 3 which wraps our core algorithm PrIC3H in more detail: First, we create an *oracle* Ω: S → [0, 1] which (roughly) *estimates* the probability of reaching B for every state. A *perfect oracle* would yield *precise* maximal reachability probabilites, i.e., Ω(s) = Prmax (s |= ♦B) for every state s. We construct oracles by user-supplied methods (highlighted in blue). Examples of implementations of all user-supplied methods in Algorithm 3 are discussed in Sect. 7.

Assuming the oracle is good, but not perfect, we construct an adequate heuristic H selecting probabilities based on the oracle7 for all successors of a given state: There are various options. The simplest is to pass-through the oracle values. A version that is more robust against noise in the oracle is discussed in Sect. 6. We then invoke PrIC3H. If PrIC3H reports safety, the MDP is indeed safe by the soundness of PrIC3H.

Check Refutation. If PrIC3H does not report safety, it reports a subsystem that hints to a *potential* counterexample. Formally, this subsystem is a subMDP of states that were 'visited' during the invocation of StrengthenH.

Definition 6 (subMDP). *Let* M = (S, s*I* , Act, P) *be an MDP and let* S ⊆ S *with* s*I* ∈ S *. We call* M*S*- = (S , s*I* , Act, P ) *the* subMDP induced by M and S *, where for all* s, s ∈ S *and all* a ∈ Act*, we have* P (s, a, s ) = P(s, a, s )*.* -

7 We thus assume that heuristic H invokes the oracle whenever it needs to guess some probability.

A subMDP M*S* may be substochastic where missing probability mass never reaches a bad state. Definition 1 is thus relaxed: For all states s ∈ S we require that - *s*-∈*S*- P(s, a, s ) ≤ 1.If the subsystem is unsafe, we can conclude that the original MDP M is also safe.

Lemma 6. *If* M *is a subMDP of* M *and* M *is unsafe, then* M *is also unsafe.* The role of CheckRefutation is to establish whether the subsystem is indeed a true counterexample or a spurious one. Formally, CheckRefutation should ensure:

 true *res* ← CheckRefutation (*subsystem*) *res* = true ⇔ M*subsystem* unsafe .

Again, PrIC3 is backward compatible in the sense that a single fixed heuristic is always sufficient when reasoning about reachability (λ = 0).

Recovery Statement 5. *For qualitative reachability (*λ = 0*) and the heuristic* H0 *from Recovery Statement 4,* PrIC3 *invokes its core* PrIC3H *exactly* once*.*

This statement is true, as PrIC3H returns either *safe* or a subsystem containing a path from the initial state to a bad state. In the latter case, CheckRefutation detects that the subsystem is indeed a counterexample which cannot be spurious in the qualitative setting.

We remark that the procedure CheckRefutation invoked in l. 5 is a classical fallback; it runs an (alternative) model checking algorithm, e.g., solving the set of Bellman equations, for the subsystem. In the worst case, i.e., for S = S, we thus solve exactly our problem statement. Empirically (Table 1) we observe that for reasonable oracles the procedure CheckRefutation is invoked on significantly smaller subMDPs. However, in the worst case the subMDP must include *all* paths of the original MDP, and then thus coincides.

Refine Oracle. Whenever we have neither proven the MDP safe nor unsafe, we refine the oracle to prevent generating the same subsystem in the next invocation of PrIC3H. To ensure termination, oracles should only be refined finitely often. That is, we need some progress measure. The set *touched* overapproximates all counterexamples encountered in some invocation of PrIC3H and we propose to use its size as the progress measure. While there are several possibilities to update *touched* through the user-defined procedure Enlarge (l. 6), every implementation should hence satisfy true *touched* ← Enlarge(*touched*, ) |*touched* | > |*touched*| . Consequently, after finitely many iterations, the oracle is refined with respect to all states. In this case, we may as well rely on solving the characteristic LP problem:

Lemma 7. *The algorithm* PrIC3 *in Algorithm 3 is sound and complete if* Refine(Ω, S) *returns a perfect oracle* Ω *(with* S *is the set of all states).*

Weaker assumptions on Refine are possible, but are beyond the scope of this paper. Moreover, the above lemma does not rely on the abstract concept that heuristic H provides suitable probabilities after finitely many refinements.8

8 One could of course now also create a heuristic that is trivial for a perfect oracle and invoke PrIC3H with the heuristic for the perfect oracle, but there really is no benefit in doing so.

## 6 Practical PrIC3

So far, we gave a conceptual view on PrIC3, but now take a more practical stance. We detail important features of effective implementations of PrIC3 (based on our empirical evaluation). We first describe an implementation without generalization, and then provide a prototypical extension that allows for three variants of generalization.

## 6.1 A Concrete PrIC3 Instance Without Generalization

*Input.* We describe MDPs using the Prism guarded command language9, exemplified in Fig. 3. States are described by valuations to m (integer-valued) program variables vars, and outgoing actions are described by commands of the form

[] guard -> prob1 : update1 & ... & probk : updatek

If a state satisfies guard, then the corresponding action with k branches exists; probabilities are given by probi, the successor states are described by updatei, see Fig. 3b.

Fig. 3. Illustrative Prism-style probabilistic guarded command language example

*Encoding.* We encode frames as logical formulae. Updating frames then corresponds to adding conjuncts, and checking for relative inductivity is a satisfiability call. Our encoding is as follows: States are assignments to the program variables, i.e., States = Z*m*. We use various uninterpreted functions, to whom we give semantics using appropriate constraints. Frames10 are represented by uninterpreted functions Frame : States → R satisfying Frame (s) = d implies F[s] ≥ d. Likewise, the Bellman operator is an uninterpreted function Phi: States → R such that Phi(s) = d implies Φ (F) [s] ≥ d. Finally, we use Bad: States → B with Bad (s) iff s ∈ B.

Among the appropriate constraints, we ensure that variables are within their range, bound the values for the frames, and enforce Phi(s)=1 for s ∈ B. We encode the guarded commands as exemplified by this encoding of the first command in Fig. 3: