Flow Faster: Efficient Decision Algorithms for Probabilistic Simulations

Strong and weak simulation relations have been proposed for Markov chains, while strong simulation and strong probabilistic simulation relations have been proposed for probabilistic automata. However, decision algorithms for strong and weak simulation over Markov chains, and for strong simulation over probabilistic automata are not efficient, which makes it as yet unclear whether they can be used as effectively as their non-probabilistic counterparts. This paper presents drastically improved algorithms to decide whether some (discrete- or continuous-time) Markov chain strongly or weakly simulates another, or whether a probabilistic automaton strongly simulates another. The key innovation is the use of parametric maximum flow techniques to amortize computations. We also present a novel algorithm for deciding strong probabilistic simulation preorders on probabilistic automata, which has polynomial complexity via a reduction to an LP problem. When extending the algorithms for probabilistic automata to their continuous-time counterpart, we retain the same complexity for both strong and strong probabilistic simulations.


Introduction
Many verification methods have been introduced to prove the correctness of systems exploiting rigorous mathematical foundations.As one of the automatic verification techniques, model checking has successfully been applied to automatically find errors in complex systems.The power of model checking is limited by the state space explosion problem.Notably, minimizing the system to the bisimulation [34,35] quotient is a favorable approach.As a more aggressive attack to the problem, simulation relations [33] have been proposed for these models.In particular, they provide the principal ingredients to perform abstractions of the models, while preserving safe CTL properties (formulas with positive universal path-quantifiers only) [16].
Simulation relations are preorders on the state space such that whenever state s ′ simulates state s (written s s ′ ) then s ′ can mimic all stepwise behaviour of s, but s ′ may perform steps that cannot be matched by s.One of the interesting aspects of simulation relations is that they allow a verification by "local" reasoning.Based on this, efficient algorithms for deciding simulation preorders have been proposed in [10,23].
Randomisation has been employed widely for performance and dependability models, and consequently the study of verification techniques of probabilistic systems with and without nondeterminism has drawn a lot of attention in recent years.A variety of equivalence and preorder relations, including strong and weak simulation relations, have been introduced and widely considered for probabilistic models.In this paper we consider discrete-time Markov chains (DTMCs) and discrete-time probabilistic automata (PAs) [39].PAs extend labelled transition systems (LTSs) with probabilistic selection, or, viewed differently, extend DTMCs with nondeterminism.They constitute a natural model of concurrent computation involving random phenomena.In a PA, a labelled transition leads to a probability distribution over the set of states, rather than a single state.The resulting model thus exhibits both non-deterministic choice (as in LTSs) and probabilistic choice (as in Markov chains).
Strong simulation relations have been introduced [26,30] for probabilistic systems.For s s ′ (s ′ strongly simulates s), it is required that every successor distribution of s via action α (called α-successor distribution) has a corresponding α-successor distribution at s ′ .Correspondence of distributions is naturally defined with the concept of weight functions [26].In the context of model checking, strong simulation relations preserve safe PCTL formulas [40].Probabilistic simulation [40] is a relaxation of strong simulation in the sense that it allows for convex combinations of multiple distributions belonging to equally labelled transitions.More concretely, it may happen that for an α-successor distribution µ of s, there is no α-successor distribution of s ′ which can be related to µ, yet there exists a so-called α-combined transition, a convex combination of several α-successor distributions of s ′ .Probabilistic simulation accounts for this and is thus coarser than strong simulation, but still preserves the same class of PCTL-properties as strong simulation does.
Apart from discrete time models, this paper considers continuous-time Markov chains (CTMCs) and continuous-time probabilistic automata (CPAs) [29,42].In CPAs, the transition delays are governed by exponential distributions.CPAs can be considered also as extensions of CTMCs with nondeterminism.CPAs are natural foundational models for various performance and dependability modelling formalisms including stochastic activity networks [37], generalised stochastic Petri nets [32] and interactive Markov chains [24].Strong simulation and probabilistic simulation have been introduced for continuous-time models [9,44].For CTMCs, s s ′ requires that s s ′ holds in the embedded DTMC, and additionally, state s ′ must be "faster" than s which manifests itself by a higher exit rate.Both strong simulation and probabilistic simulation preserve safe CSL formulas [4], which is a continuous stochastic extension of PCTL, tailored to continuous-time models.
Weak simulation is proposed in [9] for Markov chains.In weak simulation, the successor states are split into visible and invisible parts, and the weight function conditions are only imposed on the transitions leading to the visible parts of the successor states.Weak simulation is strictly coarser than the afore-mentioned strong simulation for Markov chains, thus allows further reduction of the state space.It preserves the safe PCTL-and CSLproperties without the next state formulas for DTMCs and CTMCs respectively [9].
Decision algorithms for strong and weak simulations over Markov chains, and for strong simulation over probabilistic automata are not efficient, which makes it as yet unclear whether they can be used as effectively as their non-probabilistic counterparts.In this paper we improve efficient decision algorithms, and devise new algorithms for deciding strong and strong probabilistic simulations for probabilistic automata.Given the simulation preorder, the simulation quotient automaton is in general smaller than the bisimulation quotient automaton.Then, for safety and liveness properties, model checking can be performed on this smaller quotient automata.The study of decision algorithms is also important for specification relations: The model satisfies the specification if the automaton for the specification simulates the automaton for the model.In many applications the specification cannot be easily expressed by logical formulas: it is rather a probabilistic model itself.Examples of this kind include various recent wireless network protocols, such as ZigBee [21], Firewire Zeroconf [11], or the novel IEEE 802.11e,where the central mechanism is selecting among different-sided dies, readily expressible as a probabilistic automaton [31].
The common strategy used by decision algorithms for simulations is as follows.The algorithm starts with a relation R which is guaranteed to be coarser than the simulation preorder .Then, the relation R is successively be refined.In each iteration of the refinement loop, pairs (s, s ′ ) are eliminated from the relation R if the corresponding simulation conditions are violated with respect to the current relation.In the context of labelled transitions systems, this happens if s has a successor state t, but we cannot find a successor state t ′ of s ′ such that (t, t ′ ) is also in the current relation R. For DTMCs, this correspondence is formulated by the existence of a weight function for distributions (P(s, •), P(s ′ , •)) with respect to the current relation R. Checking this weight function condition amounts to checking whether there is a maximum flow over the network constructed out of (P(s, •), P(s ′ , •)) and the current relation R. The complexity for one such check is however rather expensive, it has time complexity O(n 3 / log n).If the iterative algorithm reaches a fix-point, the strong simulation preorder is obtained.The number of iterations of the refinement loop is at most O(n 2 ), and the overall complexity [3] amounts to O(n 7 / log n) in time and O(n 2 ) in space.
Fixing a pair (s, s ′ ), we observe that the networks for this pair across iterations of the refinement loop are very similar: They differ from iteration to iteration only by deletion of some edges induced by the successive clean up of R. We exploit this by adapting a parametric maximum flow algorithm [18] to solve the maximum flow problems for the arising sequences of similar networks, hence arriving at efficient simulation decision algorithms.The basic idea is that all computations concerning the pair (s, s ′ ) can be performed in an incremental way: after each iteration we save the current network together with maximum flow information.Then, in the next iteration, we update the network, and derive the maximum flow while using the previous maximum flow function.The maximum flow problems for the arising sequences of similar networks with respect to the pair (s, s ′ ) can be computed in time O(|V | 3 ) where |V | is the number of nodes of the network.This leads to an overall time complexity O(m 2 n) for deciding the simulation preorder.Because of the storage of the networks, the space complexity is increased to O(m 2 ).Especially in the very common case where the state fanout of a model is bounded by a constant g (and hence m ≤ gn), our strong simulation algorithm has time and space complexity O(n 2 ).The algorithm can be extended easily to handle CTMCs with same time and space complexity.For weak simulation on Markov chains, the parametric maximum flow technique cannot be applied directly.Nevertheless, we manage to incorporate the parametric maximum flow idea into a decision algorithm with time complexity O(m 2 n 3 ) and space complexity O(n 2 ).An earlier algorithm [6] uses LP problems [27,38] as subroutines.The maximum flow problem is a special instance of an LP problem but can be solved much more efficiently [1].
We extend the algorithm to compute strong simulation preorder to also work on PAs.It takes the skeleton of the algorithm for Markov chains: It starts with a relation R which is coarser than , and then refines R until is achieved.In the refinement loop, a pair (s, s ′ ) is eliminated if the corresponding simulation conditions are violated with respect to the current relation.For PAs, this means that there exists an α-successor distribution µ of s, such that for all α-successor distribution µ ′ of s ′ , we cannot find a weight function for (µ, µ ′ ) with respect to the current relation R. Again, as for Markov chains, the existence of such weight functions can be reduced to maximum flow problems.Combining with the parametric maximum flow algorithm [18], we arrive at the same time complexity O(m 2 n) and space complexity O(m 2 ) as for Markov chains.The above maximum flow based procedure cannot be applied to deal with strong probabilistic simulation for PAs.The reason is that an αcombined transition of state s is a convex combination of several α-successor distributions of s, thus induces uncountable many such possible combined transitions.The computational complexity of deciding strong probabilistic simulation has not been investigated before.We show that it can be reduced to solving LP problems.The idea is that we introduce for each α-successor distribution a variable, and then reformulates the requirements concerning the combined transitions by linear constraints over these variables.This allows us to construct a set of LP problem such that whether a pair (s, s ′ ) should be thrown out of the current pair R is equivalent to whether each of the LP problem has a solution.
The algorithms for PAs are then extended to handle their continuous-time analogon, CPAs.In the algorithm, for each pair (s, s ′ ) in the refinement loop, an additional rate condition is ensured by an additional check via comparing the appropriate rates of s and s ′ .The resulting algorithm has the same time and space complexity.Related Works.In the non-probabilistic setting, the most efficient algorithms for deciding simulation preorders have been proposed in [10,23].The complexity is O(mn) where n and m denote the number of states and transitions of the transition system respectively.For Markov chains, Derisavi et al. [17] presented an O(m log n) algorithm for strong bisimulation.Weak bisimulation for DTMCs can be computed in O(n 3 ) time [5].For strong simulation, Baier et al. [3] introduced a polynomial decision algorithm with complexity O(n 7 / log n), by tailoring a network flow algorithm [20] to the problem, embedded into an iterative refinement loop.In [6], Baier et al. proved that weak simulation is decidable in polynomial time by reducing it to linear programming (LP) problems.For a subclass of PAs (reactive systems), Huynh and Tian [25] presented an O(m log n) algorithm for computing strong bisimulation.Cattani and Segala [12] have presented decision algorithms for strong and bi simulation for PAs.They reduced the decision problems to LP problems.To compute the coarsest strong simulation for PAs, Baier et al. [3] presented an algorithm which reduces the query whether a state strongly simulates another to a maximum flow problem.Their algorithm has complexity O((mn 6 + m 2 n 3 )/ log n) 1 .Recently, algorithm for computing simulation and bisimulation metrics for concurrent games [13] has been studied.
Outline of The Paper.The paper proceeds by recalling the definition of the models and simulation relations in Section 2. In Section 3 we give a short interlude on maximum flow problems.In Section 4 we present a combinatorial method to decide strong simulations.In this section we also introduce new decision algorithms for deciding strong probabilistic simulations for PAs and CPAs.In Section 5 we focus on algorithms for weak simulations.Section 6 concludes the paper.

Preliminaries
In Subsection 2.1, we recall the definitions of fully probabilistic systems, discrete-and continuous-time Markov chains [41], and the nondeterministic extensions of these discretetime [40] and continuous-time models [36,7].In Subsection 2.2 we recall the definition of simulation relations.
2.1.Markov Models.Firstly, we introduce some general notations.Let X, Y be finite sets.For f : X → R, let f (A) denote x∈A f (x) for all A ⊆ X.For f : X × Y → R, we let f (x, A) denote y∈A f (x, y) for all x ∈ X and A ⊆ Y , and f (A, y) is defined similarly.Let AP be a fixed, finite set of atomic propositions.
A state s is called stochastic and absorbing if the distribution P(s, •) is stochastic and absorbing respectively.For s ∈ S, let post (s) = Supp(P(s, •)), and let post ⊥ (s) = Supp ⊥ (P(s, •)).Definition 2.2.A labelled discrete-time Markov chain (DTMC) is an FPS M = (S, P, L) where s is either absorbing or stochastic for all s ∈ S.
FPSs and DTMCs are time-abstract, since the duration between triggering transitions is disregarded.We observe the state only at a discrete set of time points 0, 1, 2, . ... We recall the definition of CTMCs which are time-aware: Definition 2.3.A labelled continuous-time Markov chain (CTMC) is a tuple M = (S, R, L) with S and L as before, and R : S × S → R ≥0 is a rate matrix.
For CTMC M, let post (s) = {s ′ ∈ S | R(s, s ′ ) > 0} for all s ∈ S. The rates give the average delay of the corresponding transitions.Starting from state s, the probability that within time t a successor state is chosen is given by 1 − e −R(s,S)t .The probability that a specific successor state s ′ is chosen within time t is thus given by (1 − e −R(s,S)t ) • R(s, s ′ )/R(s, S).A CTMC induces an embedded DTMC, which captures the time-abstract behaviour of it: Definition 2.4.Let M = (S, R, L) be a CTMC.The embedded DTMC of M is defined by emb(M) = (S, P, L) with P(s, s ′ ) = R(s, s ′ )/R(s, S) if R(s, S) > 0 and 0 otherwise.
We will also use P for a CTMC directly, without referring to its embedded DTMC explicitly.If one is interested in time-abstract properties (e. g., the probability to reach a set of states) of a CTMC, it is sufficient to analyse its embedded DTMC.
For a given FPS, DTMC or CTMC, its fanout is defined by max s∈S |post (s)|.The number of states is defined by n = |S|, and the number of transitions is defined by m = s∈S |post (s)|.For s ∈ S, reach(s) denotes the set of states that are reachable from s with positive probability.For a relation R ⊆ S × S and s ∈ S, let R(s) denote the set Markov chains are purely probabilistic.Now we consider extensions of Markov chains with nondeterminism.We first recall the definition of probabilistic automata, which can be considered as the simple probabilistic automata with transitions allowing deadlocks in [39].Definition 2.5.A probabilistic automaton (PA) is a tuple M = (S, Act, P, L) where S and L are defined as before, Act is a finite set of actions, P ⊆ S × Act × Dist(S) is a finite set, called the probabilistic transition matrix.
For (s, α, µ) ∈ P, we use s The fanout of a state s is defined by fan(s) = α∈Act(s) µ∈Stepsα(s) (|µ| + 1).Intuitively, fan(s) denotes the total sum of the sizes of outgoing distributions of state s plus their labelling.The fanout of M is defined by max s∈S fan(s).Summing up over all states, we define the size of the transitions by m = s∈S fan(s).
A Markov decision process (MDP) [36] arises from a PA M if for s ∈ S and α ∈ Act, there is at most one α-successor distribution µ of s, which must be stochastic.
We consider a continuous-time counterpart of PAs where the transitions are described by rates instead of probabilities.A rate function is simply a function r : S → R ≥0 .Let |r| = |{s | r(s) > 0}| denote the size of r.Let Rate(S) denote the set of all rate functions.Definition 2.6.A continuous-time PA (CPA) is a tuple (S, Act, R, L) where S, Act, L as defined for PAs, and R ⊆ S × Act × Rate(S) a finite set, called the rate matrix.
We write s α − → r if (s, α, r) ∈ R, and call r an α-successor rate function of s.For transition s α − → r, the sum r(S) is also called the exit rate of it.Given that the transition s α − → r is chosen from state s, the probability that any successor state is chosen within time t is given by 1 − e −r(S)t , and a specific successor state s ′ is chosen within time t is given by (1 − e −r(S)t ) • r(s ′ ) r(S) .The notion of Act(s), Steps α (s), Steps(s), fanout and size of transitions for PAs can be extended to CPAs by replacing occurrence of distribution µ by rate function r in an obvious way.
The model continuous-time Markov decision processes (CTMDPs) [36,7] can be considered as special CPAs where for s ∈ S and α ∈ Act, there exists at most one rate function r ∈ Rate(S) such that s α − → r.The model CTMDPs considered in paper [42] essentially agree with our CPAs.
Figure 1: An FPS for illustrating the simulation relations2 .

2.2.
Strong and Weak Simulation Relations.We first recall the notion of strong simulation on Markov chains [9], PAs [40], and CPAs [44].Strong probabilistic simulation is defined in Subsection 2.2.2.Weak simulation for Markov chains will be given in Subsection 2.2.3.The notion of simulation up to R is introduced in Subsection 2.2.4.

Strong Simulation.
Strong simulation requires that each successor distribution of one state has a corresponding successor distribution of the other state.The correspondence of distributions is naturally defined with the concept of weight functions [26], adapted to FPSs as in [9].For a relation R ⊆ S × S, we let We write µ ⊑ R µ ′ if there exists a weight function for (µ, µ ′ ) with respect to R.
The first condition requires that only pairs (s, s ′ ) in the relation R ⊥ have a positive weight.In other words, for s, s ′ ∈ S with s ′ ∈ R ⊥ (s), it holds that ∆(s, s ′ ) = 0. Strong simulation requires similar states to be related via weight functions on their distributions [26].Definition 2.8.Let M = (S, P, L) be an FPS, and let R ⊆ S × S. The relation R is a strong simulation on M iff for all s 1 , s 2 with s 1 R s 2 : L(s 1 ) = L(s 2 ) and P(s We say that s 2 strongly simulates s 1 in M, denoted by s 1 M s 2 , iff there exists a strong simulation R on M such that s 1 R s 2 . By definition, it can be shown [9] that M is reflexive and transitive, thus a preorder.Moreover, M is the coarsest strong simulation relation for M. If the model M is clear from the context, the subscript M may be omitted.Assume that s 1 s 2 and let ∆ denote the corresponding weight function.If P(s 2 , ⊥) > 0, we have that P(s 2 , ⊥) = s∈S ⊥ ∆(s, ⊥) = ∆(⊥, ⊥).The second equality follows by the fact that ⊥ can not strongly simulate any real state in S. Another observation is that if s is absorbing, then it can be strongly simulated by any other state s ′ with L(s) = L(s ′ ).
Example 2.9.Consider the FPS depicted in Figure 1.Recall that labelling of states is indicated by colours in the states.Since the yellow (grey) states are absorbing, they strongly simulate each other.The same holds for the green (dark grey) states.We show now that s 1 s 2 but s 2 s 3 .
Since each DTMC is a special case of an FPS, Definition 2.8 applies directly for DTMCs.For CTMCs we say that s 2 strongly simulates s 1 if, in addition to the DTMC conditions, s 2 can move stochastically faster than s 1 [9], which manifests itself by a higher rate.Definition 2.10.Let M = (S, R, L) be a CTMC and let R ⊆ S × S. The relation R is a strong simulation on M iff for all s 1 , s 2 with s 1 R s 2 : L(s 1 ) = L(s 2 ), P(s 1 , •) ⊑ R P(s 2 , •) and R(s 1 , S) ≤ R(s 2 , S).
We say that s 2 strongly simulates s 1 in M, denoted by s 1 M s 2 , iff there exists a strong simulation R on M such that s 1 R s 2 .
Thus, s M s ′ holds if s emb(M) s ′ , and s ′ is faster than s.By definition, it can be shown that M is a preorder, and is the coarsest strong simulation relation for M. For PAs, strong simulation requires that every α-successor distribution of s 1 is related to an α-successor distribution of s 2 via a weight function [40,26]: Definition 2.11.Let M = (S, Act, P, L) be a PA and let R ⊆ S × S. The relation R is a strong simulation on M iff for all s 1 , s 2 with We say that s 2 strongly simulates s 1 in M, denoted s 1 M s 2 , iff there exists a strong simulation R on M such that s 1 R s 2 .
Example 2.12.Consider the PA in Figure 2.Then, it is easy to check s 1 s 2 : each α-successor distribution of s 1 has a corresponding α-successor distribution of s 2 .However, s 1 does not strongly simulate s 2 , as the middle α-successor distribution of s 2 can not be related by any α-successor distribution of s 1 .Now we consider CPAs.For a rate function r, we let µ(r) ∈ Dist(S) denote the induced distribution defined by: if r(S) > 0 then µ(r)(s) equals r(s)/r(S) for all s ∈ S, and if r(S) = 0, then µ(r)(s) = 0 for all s ∈ S. Now we introduce the notion of strong simulation for CPAs [44], which can be considered as an extension of the definition for CTMCs [9]: − → r 1 then there exists a transition s 2 α − → r 2 with µ(r 1 ) ⊑ R µ(r 2 ) and r 1 (S) ≤ r 2 (S).We write s 1 M s 2 iff there exists a strong simulation R on M such that s 1 R s 2 .
Similar to CTMCs, the additional rate condition r 1 (S) ≤ r 2 (S) indicates that the transition s 2 α − → r 2 is faster than s 1 α − → r 1 .As a shorthand notation, we use r 1 ⊑ R r 2 for the condition µ(r 1 ) ⊑ R µ(r 2 ) and r 1 (S) ≤ r 2 (S).For both PAs and CPAs, M is the coarsest strong simulation relation.

Strong Probabilistic Simulations.
We recall the definition of strong probabilistic simulation, which is coarser than strong simulation, but still preserves the same class of PCTLproperties as strong simulation does.We first recall the notion of combined transition [39], a convex combination of several equally labelled transitions: Strong probabilistic simulation is insensitive to combined transitions3 , thus, it is a relaxation of strong simulation.Similar to strong simulation, p M is the coarsest strong probabilistic simulation relation for M. Since MDPs can be considered as special PAs, we obtain the notions of strong simulation and strong probabilistic simulation for MDPs.Moreover, these two relations coincide for MDPs as, by definition, for each state there is at most one successor distribution per action.Example 2.16.We consider again the PA depicted in Figure 2. From Example 2.12 we know that s 2 s 1 .In comparison to state s 1 , state s 2 has one additional α-successor distribution: to states u 4 and v 4 with equal probability 0.5.This successor distribution can be considered as a combined transition of the two successor distributions of s 1 : each with constant 0.5.Hence, we have s 2 p s 1 .
We extend the notion of strong probabilistic simulation for PAs to CPAs.First, we introduce the notion of combined transitions for CPAs.In CPAs the probability that a transition occurs is exponentially distributed.The combined transition should also be exponentially distributed.The following example shows that a straightforward extension of Definition 2.14 does not work.
Example 2.17.For this purpose we consider the CPA in Figure 3. Let r 1 and r 2 denote left and the right α-successor rate functions out of state s 1 .Obviously, they have different exit rates: r 1 (S) = 10, r 2 (S) = 18.Taking each with probability 0.5, we would get the combined transition r = 0.5r 1 + 0.5r 2 : r({u 1 , u 2 }) = 7 and r({v 1 , v 2 }) = 7.However, r is hyper-exponentially distributed: the probability of reaching yellow (grey) states (u 1 or u 2 ) within time t under r is given by: 0.5 ).Similarly, the probability of reaching green (dark grey) states within time t is given by: 0.5 ).From state s 2 , the two α-successor rate functions have the same exit rate 14.Let r ′ 1 and r ′ 2 denote left and the right α-successor rate functions out of state s 2 .In this case the combined transition r ′ = 0.5r ′ 1 + 0.5r ′ 2 is also exponentially distributed with rate 14: the probability to reach yellow (grey) states (u 3 and u 4 ) within time t is 7  14 • (1 − e −14t ), which is the same as the probability of reaching green (dark grey) states within time t.Based on the above example, it is easy to see that to get a combined transition which is still exponentially distributed, we must consider rate functions with the same exit rate: Definition 2.18.Let M = (S, Act, R, L) be a CPA.Let s ∈ S, α ∈ Act(s) and let {r 1 , . . ., r k } ⊆ Steps α (s) where r i (S) = r j (S) for i, j ∈ {1, . . ., k}.The tuple (s, α, r) is a combined transition, denoted by s In the above definition, unlike for the PA case, only α-successor rate functions with the same exit rate are combined together.Similar to PAs, strong probabilistic simulation is insensitive to combined transitions, which is thus a relaxation of strong simulation: We write s 1 p M s 2 iff there exists a strong simulation R on M such that s 1 R s 2 .Recall r 1 ⊑ R r 2 is a shorthand notation for µ(r 1 ) ⊑ R µ(r 2 ) and r 1 (S) ≤ r 2 (S).By definition, the defined strong probabilistic simulation p M is the coarsest strong probabilistic simulation relation for M.
Example 2.20.Reconsider the CPA in Figure 3.As discussed in Example 2.17, the two α-successor rate functions of s 1 cannot be combined together, thus the relation s 0 p s 1 cannot be established.However, s 0 p s 2 holds: denoting the left rate function of s 2 as r 1 and the right rate function as r 2 , we choose as the combined rate function r = 0.5r 1 + 0.5r 2 .Obviously, the conditions in Definition 2.19 are satisfied.

Weak Simulations.
We now recall the notion of weak simulation [9] on Markov chains4 .Intuitively, s 2 weakly simulates s 1 if they have the same labelling, and if their successor states can be grouped into sets U i and V i for i = 1, 2, satisfying certain conditions.Consider Figure 4. We can view steps to V i as stutter steps while steps to U i are visible steps.With respect to the visible steps, it is then required that there exists a weight function for the conditional distributions: P(s 1 , •)/K 1 and P(s 2 , •)/K 2 where K i intuitively is the probability to perform a visible step from s i .The stutter steps must respect the weak simulation relations: thus states in V 2 should weakly simulate s 1 , and state s 2 should weakly simulate states in V 1 .This is depicted by dashed arrows in the figure.For reasons we will explain later in Example 2.22, the definition needs to account for states which partially belong to U i and partially to V i .Technically, this is achieved by functions δ i that distribute s i over U i and V i in the definition below.For a given pair (s 1 , s 2 ) and functions If (s 1 , s 2 ) and δ i are clear from the context, we write U i , V i instead.
Definition 2.21.Let M = (S, P, L) be a DTMC and let R ⊆ S × S. The relation R is a weak simulation on M iff for all s 1 , s 2 with s 1 R s 2 : L(s 1 ) = L(s 2 ) and there exist functions Figure 5: A DTMC where splitting states is necessary to establish the weak simulation.In the model some states are drawn more than once.
(b) if K 1 > 0 and K 2 > 0 then for all states w ∈ S: where (3) for u 1 ∈ U 1 there exists a path fragment s 2 , w 1 , . . ., w n , u 2 with positive probability such that n ≥ 0, s 1 R w j for 0 < j ≤ n, and u 1 R u 2 .We say that s 2 weakly simulates s 1 in M, denoted s 1 M s 2 , iff there exists a weak simulation R on M such that s 1 R s 2 .
Note again that the sets U i , V i in the above definition are defined according to Equation 2.1 with respect to the pair (s 1 , s 2 ) and the functions δ i .The functions δ i can be considered as a generalisation of the characteristic function of U i in the sense that we may split the membership of a state to U i and V i into fragments which sum up to 1.For example, if δ 1 (s) = 1  3 , we say that 1 3 fragment of the state s belongs to U 1 , and 2 3 fragment of s belongs to V 1 .Hence, U i and V i are not necessarily disjoint.Observe that Condition 3 will in the sequel be called the reachability condition.If K 1 > 0 and K 2 = 0, which implies that U 2 = ∅ and U 1 = ∅, the reachability condition guarantees that for any visible step s 1 → u 1 with u 1 ∈ U 1 , s 2 can reach a state u 2 which simulates u 1 while passing only through states simulating s 1 .Assume that we have S = {s 1 , s 2 , u 1 } where L(s 1 ) = L(s 2 ) and u 1 has a different labelling.There is only one transition P(s 1 , u 1 ) = 1.Obviously s 1 s 2 .Dropping Condition 3 would mean that s 1 s 2 .We illustrate the use of fragments of states in the following example: Example 2.22.Consider the DTMC depicted in Figure 5.For states u 1 , u 2 , v 1 , v 2 , obviously the following pairs (u 1 , u 2 ), (u 1 , v 2 ), (v 1 , v 2 ) are in the weak simulation relation.The state u 2 cannot weakly simulate v 1 .Since v 2 weakly simulates v 1 , it holds that s 2 s 5 .Similarly, from u 1 v 1 we can easily show that s 1 s 4 .We observe also that s 2 s 3 : K 1 > 0 and K 2 > 0 since both s 2 and s 3 have yellow (grey) successor states, but the required function ∆ cannot be established since u 2 cannot weakly simulate any successor state of s 2 (which is v 1 ).Thus s 2 s 3 .
Without considering fragments of states, we show that a weak simulation between s 1 and s 3 cannot be established.Since s 2 s 3 , we must have We consider the following two cases: • The case δ 2 (s 4 ) = 1.In this case we have that K 2 = 1.A function ∆ must be defined satisfying Condition 2b in Definition 2.21.Taking w = s 4 , the following must hold: The state s 2 is the only successor of s 1 that can be weakly simulated by s 4 , so ∆(s 2 , s 4 ) = 0.75 must hold.However, the equation does not hold any more, as on the left side we have 0.75 but on the right side we have 0.5 instead.• The case δ 2 (s 4 ) = 0.In this case we have still K 2 > 0. Similar to the previous case it is easy to see that the required function ∆ cannot be defined: the equation does not hold since the left side is 0 (no states in U 2 can weakly simulate s 2 ) but the right side equals 0.5.
Thus without splitting, s 3 does not weakly simulate s 1 .We show it holds that ), (q 1 , q 1 ), (s 2 , s 5 ), (s 2 , s 4 )} is a weak simulation relation.By the discussions above, it is easy to verify that every pair except (s 1 , s 3 ) satisfies the conditions in Definition 2.21.
Weak simulation for CTMCs is defined as follows.
Definition 2.23 ([9, 8]).Let M = (S, R, L) be a CTMC and let R ⊆ S × S. The relation R is a weak simulation on M iff for s 1 R s 2 : L(s 1 ) = L(s 2 ) and there exist functions δ i : S → [0, 1] (for i = 1, 2) satisfying Equation 2.1 and Conditions 1 and 2 of Definition 2.21 and the rate condition: In this definition, the rate condition 3 ′ strengthens the reachability condition of the preceding definition.If U 1 = ∅, we have that K 1 > 0; the rate condition then requires that Figure 6: A simple FPS for illustrating the simulation up to R. K 2 > 0, which implies U 2 = ∅.For both DTMCs and CTMCs, the defined weak simulation is a preorder [9], and is the coarsest weak simulation relation for M.  6.Let R = {(s 1 , s 2 ), (w 1 , w 2 )}.Since L(q 1 ) = L(q 2 ) we have that w 1 w 2 .Thus, R is not a strong simulation.However, s 1 R s 2 , as the weight function is given by ∆(w 1 , w

Maximum Flow Problems
Before introducing algorithms to decide the simulation preorder, we briefly recall the preflow algorithm [20] for finding the maximum flow over the network N = (V, E, c) where V is a finite set of vertices, E ⊆ V × V is a set of edges, and c : E → R >0 ∪ {∞} is the capacity function.V contains a distinguished source vertex and a distinguished sink vertex .We extend the capacity function to all vertex pairs: conservation rule The value of a flow function f is given by f ( , V ).A maximum flow is a flow of maximum value.A preflow is a function f : V × V → R satisfying Conditions 1 and 2 above, and the relaxation of Condition 3: Related to maximum flows are minimum cuts.A cut of a network N = (V, E, c) is a partition of V into two disjoint sets (X, X ′ ) such that ∈ X and ∈ X ′ .The capacity of (X, X ′ ) is the sum of all capacities of edges from X to X ′ , i. e., v∈X,w∈X ′ c(v, w).A minimum cut is a cut with minimal capacity.The Maximum Flow Minimum Cut Theorem [1] states that the capacity of a minimum cut is equal to the value of a maximum flow.
The Preflow Algorithm.We initialise the preflow f by: and 0 otherwise.The preflow algorithm preserves the validity of the preflow f and the distance function d.If there is an active vertex v such that the residual edge (v, w) is admissible, we push δ := min{e(v), c f (v, w)} amount of flow from v toward the sink along the admissible edge (v, w) by increasing f (v, w) (and decreasing f (w, v)) by δ.The excesses of v and w are then modified accordingly by: e(v) = e(v) − δ and e(w) = e(w) + δ.If v is active but there are no admissible edges leaving it, one may relabel v by letting d(v Pushing and relabelling are repeated until there are no active vertices left.The algorithm terminates if no such operations apply.The resulting final preflow f is a maximum flow. Feasible Flow Problem.Let A ⊆ E be a subset of edges of the network N = (V, E, c), and define the lower bound function l : A → R >0 which satisfies l(e) ≤ c(e) for all e ∈ A. We address the feasible flow problem which consists of finding a flow function f satisfying the condition: f (e) ≥ l(e) for all e ∈ A. We briefly show that this problem can be reduced to the maximum flow problem [1].
We can replace a minimum flow requirement on edge v → w by turning v into a demanding vertex (i.e., a vertex that consumes part of its inflow) and turning w into a supplying vertex (i.e., a vertex that creates some outflow ex nihilo).The capacity of edge v → w is then reduced accordingly.Now, we are going to look for a flow-like function for the updated network.The function should satisfy the capacity constraints, and the difference between outflow and inflow in each vertex corresponds to its supply or demand, except for and .To remove that last exception, we add an edge from to with capacity ∞.
We then apply another transformation to the updated network so that we can apply the maximum flow algorithm.We add new source and target vertices ′ and ′ .For each supplying vertex s, we add an edge ′ → s with the same capacity as the supply of the vertex.For each demanding vertex d, we add an edge d → ′ with the same capacity as the demand of the vertex.In [1] it is shown that the original network has a feasible flow if and only if the transformed network has a flow h that saturates all edges from ′ and all edges to ′ .The flow h necessarily is a maximum flow, and if there is an h, each maximum flow satisfies the requirement; therefore it can be found by the maximum flow algorithm.An example will be given in Example 5.4 in Section 5.

Algorithms for Deciding Strong Simulations
We first recall the basic algorithm to compute the largest strong simulation relation in Subsection 4.1.Then, we refine this algorithm to deal with strong simulation on Markov chains in Subsection 4.2, and extend it to deal with probabilistic automata in Subsection 4.3.
SimRel s (M) In Subsection 4.4 we present an algorithm for deciding strong probabilistic simulation for probabilistic automata.
4.1.Basic Algorithm to Decide Strong Simulation.The algorithm in [3], copied as SimRel s in Algorithm 1, takes as a parameter a model, which, for now, is an FPS M. The subscript 's' stands for strong simulation; a very similar algorithm, namely SimRel w , will be used for weak simulation later.To calculate the strong simulation relation for M, the algorithm starts with the initial relation In iteration i, it generates R i+1 from R i by deleting each pair (s 1 , s 2 ) from R i if s 2 cannot strongly simulate s 1 up to R i , i. e., s 1 R i s 2 .This proceeds until there is no such pair left, i. e., R i+1 = R i .Invariantly throughout the loop it holds that R i is coarser than M (i.e., M is a sub-relation of R i ).We obtain the strong simulation preorder M = R i , once the algorithm terminates.
The decisive part of the algorithm is the check in Line 1.6, i. e., whether s 1 R i s 2 .This can be answered via solving a maximum flow problem on a particular network N (P(s 1 , •), P(s 2 , •), R i ) constructed from P(s 1 , •), P(s 2 , •) and R i .This network is the relevant part of a graph containing two copies t ∈ S ⊥ and t ∈ S ⊥ of each state where S ⊥ = {t | t ∈ S ⊥ } as follows: Let (the source) and (the sink) be two additional vertices not contained in S ⊥ ∪ S ⊥ .For µ, µ ′ ∈ Dist(S), and a relation R ⊆ S × S we define the network The capacity function c is defined as follows: c( , s) = µ(s) for all s ∈ Supp ⊥ (µ), c(t, ) = µ ′ (t) for all t ∈ Supp ⊥ (µ ′ ), and c(s, t) = ∞ for all other (s, t) ∈ E. This network is a bipartite network, since the vertices can be partitioned into two subsets V 1 := Supp ⊥ (µ) ∪ { } and V 2 := Supp ⊥ (µ ′ ) ∪ { } such that all edges have one endpoint in V 1 and another in V 2 .Later, we will use two variations of this network: For γ ∈ R >0 , we let N (µ, γµ ′ , R) denote the network obtained from N (µ, µ ′ , R) by setting the capacities to the sink to: c(t, ) = γµ ′ (t).For two states s 1 , s 2 of an FPS or a CTMC, we let N (s 1 , s 2 , R) denote the network N (P(s 1 , •), P(s 2 , •), R).
The following lemma expresses the crucial relationship between maximum flows and weight functions on which the algorithm is based.It is a direct extension of [3, Lemma 5.1]: Lemma 4.1.Let S be a finite set of states and R be a relation on S. Let µ, µ ′ ∈ Dist(S).Then, µ ⊑ R µ ′ iff the maximum flow of the network N (µ, µ ′ , R) has value 1.
Proof.As we introduced the auxiliary state ⊥, µ and µ ′ are stochastic distributions in Dist(S ⊥ ).The rest of the proof follows directly from [3, Lemma 5.1].
Thus we can decide s 1 R i s 2 by computing the maximum flow in N (s 1 , s 2 , R i ) and then check whether it has value 1.We recall the correctness and complexity of SimRel s which will also be used later.Proof.First we show that after the last iteration (say iteration k), it holds that is coarser than Now we show by induction that the loop of the algorithm invariantly ensures that R i is coarser than .Assume i = 1.By definition of strong simulation, s 1 s 2 implies L(s 1 ) = L(s 2 ).Thus, the initial relation R 1 is coarser than the simulation relation .Now assume that R i is coarser than for some 1 ≤ i < k; we will show that also R i+1 is coarser than .Pick a pair (s 1 , s 2 ) ∈ arbitrarily.By Definition 2.8, P(s 1 , •) ⊑ P(s 2 , •), so there exists a weight function for (P(s 1 , •), P(s 2 , •)) with respect to .Inspection of Definition 2.7 shows that the same function is also a weight function with respect to any set coarser than .As R i is coarser than by induction hypothesis, we conclude that P(s 1 , •) ⊑ R i P(s 2 , •), and from Subsection 2.2.4,s 1 R i s 2 .This implies that (s 1 , s 2 ) ∈ R i+1 by line 1.6 for all s 1 s 2 .Therefore, R i+1 is coarser than for all i = 1 . . ., k.
Now we show the complexity.For one network N (s 1 , s 2 , R i ) = (V, E, c), the sizes of the vertices |V | and edges |E| are bounded by 2n + 4 and (n + 1) 2 + 2n, respectively.The number of edges meets the worst case bound O(n 2 ).To the best of our knowledge, the best complexity of the flow computation for the network [14,19].In the algorithm SimRel s , only one pair, in the worst case, is removed from R i in iteration i, which indicates that the test whether times and so on.Altogether it is bounded by ).Hence, the overall time complexity amounts to O(n 7 / log n).The space complexity is O(n 2 ) because of the representation of the transitions in N (s 1 , s 2 , R i ).4.2.An Improved Algorithm for FPSs.We first analyse the behaviour of SimRel s in more detail.For this, we consider an arbitrary pair (s 1 , s 2 ), and assume that (s 1 , s 2 ) stays in relation R 1 , . . ., R k throughout the iterations i = 1, . . ., k, until the pair is either found not to satisfy s 1 R k s 2 or the algorithm terminates with a fix-point after iteration k.Then altogether the maximum flow algorithms are run k-times for this pair.However, the networks N (s 1 , s 2 , R i ) constructed in successive iterations are very similar, and may often be identical across iterations: They differ from iteration to iteration only by deletion of some edges induced by the successive cleanup of R i .For our particular pair (s 1 , s 2 ) the network might not change at all in some iterations, because the deletions from R i do not affect their direct successors.We are going to exploit this observation by an algorithm that reuses the already computed maximum flow, in a way that whatever happens is good: If no Apply the preflow algorithm to calculate the maximum flow for N (s 1 , s 2 , R i ), but initialise the preflow to f i and the distance function to d i .2.6: return 2.12: Apply the preflow algorithm to calculate the maximum flow for N (s 1 , s 2 , R i ).2.13: return Algorithm 2: Algorithm for a sequence of maximum flows.changes occur from N (s 1 , s 2 , R i−1 ) to N (s 1 , s 2 , R i ), then the maximum flow is the same as the one in the previous iteration.If changes do occur, the preflow algorithm can be applied to get the new maximum flow very fast, using the maximum flow and distance function constructed in the previous iteration as a starting point.
To understand the algorithm, we look at the network N (s 1 , s 2 , R 1 ).Let D 1 , . . ., D k be pairwise disjoint subsets of R 1 , which correspond to the pairs deleted from R 1 in iteration i, so denote the maximum flow of the network N (s 1 , s 2 , R i ) for 1 ≤ i ≤ k.We sometimes omit the superscript (s 1 , s 2 ) in the parameters if the pair (s 1 , s 2 ) is clear from the context.We address the problem of checking |f i | = 1 for all i = 1, . . ., k.Our algorithm sequence of maximum flows ) is shown as Algorithm 2. It executes iteration i of a parametric flow algorithm, where N (s 1 , s 2 , R i−1 ) is the network for (s 1 , s 2 ) and f i−1 and d i−1 are the flow and the distance function resulting from the previous iteration i − 1; and D i−1 is a set of edges that have to be deleted from N (s 1 , s 2 , R i−1 ) to get the current network.The algorithm returns a tuple, in which the first component is a boolean that tells whether |f i | = 1; it also returns the new network N (s 1 , s 2 , R i ), flow f i and distance function d i to be reused in the next iteration.Smf is inspired by the parametric maximum algorithm in [18].A variant of Smf is used in the first iteration, shown in lines 2.11-2.13.
This algorithm for sequence of maximum flow problems is called in an improved version of SimRel s shown as Algorithm 3. Lines 3.2-3.7 contain the first iteration, very similar to the first iteration of Algorithm 1 (lines 1.4-1.7).At line 3.4 we prepare for later iterations the set where pre(s) = {t ∈ S | P(t, s) > 0}.This set contains all pairs (u 1 , u 2 ) such that the network N (u 1 , u 2 , R 1 ) contains the edge (s 1 , s 2 ).Iteration i (for i > 1) of the loop (lines 3.10-3.18)calculates R i+1 from R i .In lines 3.11-3.14,we collect edges that should be removed from N (u 1 , u 2 , R i−1 ) in the sets . At line 3.16, the algorithm Smf constructs the maximum flow for parameters using information from iteration i − 1.It uses the set then it constructs the maximum flow f i for the network N (s 1 , s 2 , R i ).If Smf returns true, (s 1 , s 2 ) is inserted into R i+1 and survives this iteration (line 3.18).
Consider the algorithm Smf and assume that i > 1.At lines 2.1-2.4,we remove the edges D i−1 from the network N (s 1 , s 2 , R i−1 ) and generate the preflow f i based on the flow f i−1 , which is the maximum flow of the network N (s 1 , s 2 , R i−1 ), by • setting f i (u 1 , u 2 ) = 0 for all deleted edges (u 1 , u 2 ) ∈ D i−1 , and • reducing f i (u 2 , ) such that the preflow f i becomes consistent with the (relaxed) flow conservation rule.The excess e(v) is increased if there exists (v, w) ∈ D i−1 such that f i−1 (v, w) > 0, and unchanged otherwise.Hence, f i after line 2.4 is a preflow.The distance function d i−1 = d i is still valid for this preflow, since removing the set of edges D i−1 does not introduce new residual edges.This guarantees that, at line 2.5, the preflow algorithm finds a maximum flow over the network N (s 1 , s 2 , R i ).In line 2.6, Smf returns whether the flow has value 1 together with information to be reused in the next iteration.(If |f k | < 1 at some iteration k, then |f j | < 1 for all iterations j ≥ k because deleting edges does not increase the maximum flow.In that case, it would be sufficient to return false.)We prove the correctness and complexity of the algorithm Smf:

and d i−1 be as returned by some earlier call to
) be the set of edges that will be removed from the network N (s 1 , s 2 , R i−1 ) during the (i − 1)th call of Smf.Then, the (i − 1)th call of Smf returns true iff s 1 R i s 2 .
Proof.By Lemma 4.1, Smf init returns true iff |f 1 | = 1, which is equivalent to s 1 R 1 s 2 .Let i > 1.As discussed, at the beginning of line 2.5, the function f i−1 is a flow (thus a preflow) with value 1, and the distance function d i−1 is a valid distance function.It follows directly from the correctness of the preflow algorithm [2] that after line 2.5, f i is a maximum flow for N (s 1 , s 2 , R i ).Thus, Smf returns true (i.e.
Lemma 4.4.Consider the pair (s 1 , s 2 ) and assume that |post (s Proof.In the bipartite network N (s 1 , s 2 , R 1 ), the set of vertices are partitioned into subsets V 1 = post (s 1 )∪{ } and V 2 = post (s 2 )∪{ } as described in Section 4.1.Generating the initial network (line 2.11) takes time in O(|V 1 | |V 2 |).In our sequence of maximum flow problems, the number of (nontrivial) iterations, denoted by k, is bounded by the number of edges, i. e., k ≤ |E| ≤ |V 1 | |V 2 | − 1.We split the work being done by all calls to Smf(i, N (s 1 , s 2 , •), . ..) together with the initial call to the preflow algorithm (line 2.12 and line 2.5) into edge deletions, relabels, non-saturating pushes, saturating pushes.(A non-saturating push along an edge (u, v) moves all excess at u to v; by such a push, the number of active nodes never increases.) All edge deletions take time proportional to k i=1 |D i |, which is less than the number of edges in the network.Therefore, edge deletions take time O(|V 1 ||V 2 |).For all v ∈ V , it holds that d i+1 (v) = d i (v), i.e., the labelling function at the beginning of iteration i + 1 is the same as the labelling function at the end of iteration i.
We discuss the time for relabelling and saturating pushes [2].For a bipartite network, the distance of the source can be initialised to d( ) = 2 |V 1 | instead of |V |, and d(v) never grows above 4 |V 1 | for all v ∈ V .For v ∈ V , let I(v) denote the set of nodes containing w such that either (v, w) ∈ E or (w, v) ∈ E. Intuitively, it represents edges which could be admissible leaving v.The time for relabel operations with respect to node v is thus (4 Altogether, this gives the time for all relabel operations: v∈V ((4 Between two consecutive saturating pushes on (v, w), the distances d(v) and d(w) must increase by 2. Thus, the number of saturating pushes on edge (v, w) is bounded by 4 |V 1 |.Summing over all edges, the work for saturating pushes is bounded by O(|V 1 | |E|).Now we discuss the analysis of the number of non-saturating pushes, which is very similar to the proof of Theorem 2.2 in [22] where Max-d version of the algorithm is used.Assume that in iteration l ≤ k of Smf, the last relabelling action occurs.In the Max-d version [22], always the active node with the highest label is selected, and once an active node is selected, the excess of this node is pushed until it becomes 0. This implies that, between any two relabel operations, there are at most n active nodes processed (otherwise the algorithm terminates and we get the maximum flow).Also observe that at each time an active node is selected, at most one non-saturating push can occur, which implies that there are at most n non-saturating pushes between node label increases.Since d i (v) is bounded by 4|V 1 |, the number of relabels altogether is bounded by O(|V 1 ||V |).Thus, the number of non-saturating pushes before the iteration l is bounded by O(|V 1 ||V | 2 ).Since the distance function does not change after iteration l any more, inside any of the iterations l ′ ≥ l, there are again at most n − 1 non-saturating pushes.Hence, the number of nonsaturating pushes is bounded by Proof.By Lemma 4.3, Smf init (i, s 1 , s 2 , R 1 ) returns true in iteration i = 1 iff s 1 R 1 s 2 ; Smf(i, N (s 1 , s 2 , R i−1 ), . ..) returns true in iteration i > 1 iff s 1 R i s 2 .The rest of the correctness proof is the same as the proof of Theorem 4.2.Proof.We first show the space complexity.In most cases, it is enough to store information from the previous iteration until the corresponding structure for the current iteration is calculated.The size of the set Listener (s Summing over all (s 1 , s 2 ) yields a memory consumption in O(m 2 ) again.Hence, the overall space complexity is O(m 2 ).Now we show the time complexity.We observe that a pair (s 1 , s 2 ) belongs to D i in at most one iteration.Therefore, the time needed in lines 3.11-3.14 in all iterations together is bounded by the size of all sets Listener (s 1 ,s 2 ) , which is O(m 2 ).We analyse the time needed for all calls to the algorithm Smf.Recall that the fanout g equals max s∈S |post (s)|, and therefore |post (s i )| ≤ g for i = 1, 2. By Lemma 4.4, the complexity attributed to the pair (s 1 , s 2 ) is bounded by O(g |post (s 1 )| |post (s 2 )|).Taking the sum over all possible pairs, we get the bound gm 2 ∈ O(m 2 n).If g is bounded by a constant, we have m ≤ gn, and the time complexity is gm 2 ≤ g 3 n 2 ∈ O(n 2 ).In this case the space complexity is also O(n 2 ).
Strong Simulation for Markov Chains.We now consider DTMCs and CTMCs.Since each DTMC is a special case of an FPS the algorithm SimRel FPS s applies directly.Let M = (S, R, L) be a CTMC.Recall that s M s ′ holds if s emb(M) s ′ in the embedded DTMC, and s ′ is faster than s.We can ensure the additional rate condition by incorporating it into the initial relation R.More precisely, initially R contains only those pair (s, s ′ ) such that L(s) = L(s ′ ), and that the state s ′ is faster than s, i. e., we replace line 3.1 of the algorithm by to ensure the additional rate condition of Definition 2.10.In the refinement steps afterwards, only the weight function conditions need to be checked with respect to the current relation in the embedded DTMC.Thus, we arrive at an algorithm for CTMCs with the same time and space complexity as for FPSs.
Assume that we get the maximum flow f 1 which sends 1 2 amount of flow along the path , u 2 , u 4 , and 1  2 amount of flow along , u 1 , u 3 , .Hence, the check for (s 1 , s 2 ) is successful in the first iteration.The checks for the pairs (u 1 , u 3 ), (u 1 , u 4 ) and (u 2 , u 3 ) are also successful in the first iteration.However, the check for the pair (u 2 , u 4 ) fails, as the probability to go from u 4 to q 3 in the embedded DTMC is 2  5 , while the probability to go from u 2 to q 1 in the embedded DTMC is 1.
In the second iteration, the network N (s 1 , s 2 , R 2 ) is obtained from N (s 1 , s 2 , R 1 ) by deleting the edge (u 2 , u 4 ).In N (s 1 , s 2 , R 2 ), the flows on (u 2 , u 4 ) and on (u 4 , ) are set to 0, and the vertex u 2 has a positive excess 1  2 .Applying the preflow algorithm, we push the excess from u 2 , along u 3 , u 1 , u 4 to .We get a maximum flow f 2 for N (s 1 , s 2 , R 2 ) which sends 1  2 amount of flow along the path , u 2 , u 3 , and 1 2 amount of flow along , u 1 , u 4 , .Hence, the check for (s 1 , s 2 ) is also successful in the second iteration.Once the fix-point is reached, R still contains (s 1 , s 2 ).

Strong Simulation for Probabilistic Automata.
In this subsection we present algorithms for deciding strong simulations for PAs and CPAs.It takes the skeleton of the algorithm for FPSs: it starts with a relation R which is coarser than , and then refines R until is achieved.In the refinement loop, a pair (s, s ′ ) is eliminated from the relation R if the corresponding strong simulation conditions are violated with respect to the current relation.For PAs, this means that there exists an α-successor distribution µ of s, such that for all α-successor distribution µ ′ of s ′ , we cannot find a weight function for (µ, µ ′ ) with respect to the current relation R.
Let M = (S, Act, P, L) be a PA.We aim to extend Algorithm 3 to determine the strong simulation on PAs.For a pair (s 1 , s 2 ), assume that L(s 1 ) = L(s 2 ) and that Act(s 1 ) ⊆ Act(s 2 ), which is guaranteed by the initialisation.We consider line 3.17, which checks the condition P(s 1 , •) ⊑ R i P(s 2 , •) using Smf.By Definition 2.11 of strong simulation for PAs, we should instead check the condition Recall the condition µ 1 ⊑ R i µ 2 is true iff the maximum flow of the network N (µ 1 , µ 2 , R i ) has value one.Sometimes, we write N (s 1 , α, µ 1 , s 2 , µ 2 , R i ) to denote the network N (µ 1 , µ 2 , R i ) associated with the pair (s 1 , s 2 ) with respect to action α.Our first goal is to extend Smf to check Condition 4.1 for a fixed action α and αsuccessor distribution µ 1 of s 1 .To this end, we introduce a list Sim (s 1 ,α,µ 1 ,s 2 ) that contains all potential candidates of α-successor distributions of s 2 which could be used to establish if match then 4.10: Sim i ← tail(Sim i ) 4.12: return (false, ∅, NIL, NIL, NIL) Algorithm 4: Subroutine to calculate whether s 1 R i s 2 , as far as s 1 α − → µ 1 is concerned.The parameter Sim denotes the subsets of α-successor distributions of s 2 serving as candidates for possible µ 2 .the condition µ 1 ⊑ R µ 2 for the relation R considered.The set Sim (s 1 ,α,µ 1 ,s 2 ) is represented as a list.This and some subsequent notations are similar to those used by Baier et al. in [3].We use the function head(•) to read the first element of a list; tail(•) to read all but the first element of a list; and empty(•) to check whether a list is empty.As long as the network for a fixed candidate µ 2 = head(Sim (s 1 ,α,µ 1 ,s 2 ) ) allows a flow of value 1 over the iterations, we stick to it, and we can reuse the flow and distance function from previous iterations.If by deleting some edges from N (µ 1 , µ 2 , R), its flow value falls below 1, we delete µ 2 from Sim (s 1 ,α,µ 1 ,s 2 ) and pick the next candidate.
The algorithm ActSmf, shown as Algorithm 4, implements this.It has to be called for each pair (s 1 , s 2 ) and each successor distribution s 1 α − → µ 1 of s 1 .It takes as input the list of remaining candidates Sim , the information from the previous iteration (the network N (µ 1 , µ 2 , R i−1 ), flow f i−1 , and distance function d i−1 ), and the set of edges that have to be deleted from the old network ) be the set of edges that will be removed from the network during the (i − 1)th call of ActSmf.Then, the (i − 1)th call of algorithm ActSmf returns true iff: Proof.Once Smf returns false because the maximum flow for the current candidate µ 2 has value < 1, it will never become a candidate again, as edge deletions cannot lead to increased flow.The correctness proof is then the same as the proof of Lemma 4.3.Proof.We first consider the space complexity.We save the sets D

Proof. The proof follows the same lines as the proof of the correctness of SimRel
, the networks N (s 1 , α, µ 1 , s 2 , µ 2 , R i ) which are updated in every iteration, Listener (s 1 ,s 2 ) and the sets Summing over all (s 1 , α, µ 1 , s 2 , µ 2 ), we get: Similarly, the memory needed for saving the networks has the same bound O(m 2 ).Now we consider the set Listener (s 1 ,s 2 ) for the pair (s . Then, it holds that s 1 ∈ Supp(µ 1 ) and s 2 ∈ Supp(µ 2 ).Hence, the tuple (u 1 , α, µ 1 , u 2 , µ 2 ) can be an element of Listener (s 1 ,s 2 ) of some arbitrary pair (s 1 , s 2 ) at most |µ 1 | |µ 2 | times, which corresponds to the maximal number of edges between the set of nodes Supp(µ 1 ) and Supp(µ 2 ) in N (s 1 , α, µ 1 , s 2 , µ 2 , R 1 ).Summing over all (s 1 , α, µ 1 , s 2 , µ 2 ), we get that memory needed for the set Listener is also bounded by O(m 2 ).For each pair (s 1 , s 2 ) and s has size |Steps α (s 2 )|.Summing up, this is smaller than or equal to m 2 according to Inequality 4.2.Hence, the overall space complexity amounts to O(m 2 ).Now we consider the time complexity.All initialisations (lines 5.1-5.6 of SimRel PA s and the initialisations in ActSmf init , which calls Smf init ) take O(m 2 ) time.We observe that a pair (s 1 , s 2 ) belongs to D i during at most one iteration.Because of the Inequality 4.2, the time needed in lines 5.13-5.17 is in O(m 2 ).The rest of the algorithm is dominated by the time needed for calling Smf in line 4.2 of ActSmf.By Lemma 4.4, the time complexity for successful and unsuccessful checks concerning the tuple (s Taking the sum over all possible tuples (s 1 , α, µ 1 , s 2 , µ 2 ) we get the bound gm 2 according to Inequality 4.2.Hence, the complexity is O(m 2 n).If the fanout g is bounded by a constant, we have m ≤ gn.Thus, the time complexity is in the order of O(n 2 ).In this case the space complexity is also O(n 2 ).The decision algorithm for strong simulation for CPAs can be adapted from SimRel PA s in Algorithm 5 easily: Notations are extended with respect to rate functions instead of distributions in an obvious way.To guarantee the additional rate condition, we rule out successor rate functions of s 2 that violate it by replacing line 5.6 by: For each pair (s 1 , s 2 ), and successor rate functions r i ∈ Steps α (s i ) (i = 1, 2), the subroutine for checking whether r 1 ⊑ R i r 2 is then performed in the network N (µ(r 1 ), µ(r 2 ), R i ).Obviously, the so obtained algorithm for CPAs has the same complexity O(m 2 n).

4.4.
Strong Probabilistic Simulation.The problem of deciding strong probabilistic simulation for PAs has not been tackled yet.We show that it can be computed by solving LP problems which are decidable in polynomial time [27].In Subsection 4.4.1, we first present an algorithm for PAs.We extend the algorithm to deal with CPAs in Subsection 4.4.2.4.4.1.Probabilistic Automata.Recall that strong probabilistic simulation is a relaxation of strong simulation in the sense that it allows combined transitions, which are convex combinations of multiple distributions belonging to equally labelled transitions.Again, the most important part is to check whether s 1 p R s 2 where R is the current relation.By Definition 2.15, it suffices to check L(s 1 ) = L(s 2 ) and the condition: However, since the combined transition involves the quantification of the constants c i ∈ [0, 1], there are possibly infinitely many such µ 2 .Thus, one cannot check µ 1 ⊑ R µ 2 for each possible candidate µ 2 .The following lemma shows that this condition can be checked by solving LP problems which are decidable in polynomial time [27,38].
Proof.First assume that s 1  it suffices to check L(s 1 ) = L(s 2 ) and the condition: Recall that for CPAs only successor rate functions with the same exit rate can be combined together.For state s ∈ S, we let E(s) := {r(S) | s α − → r} denote the set of all possible exit rates of α-successor rate functions of s.For E ∈ E(s) and α ∈ Act(s), we let Steps E α (s) = {r ∈ Steps α (s) | r(S) = E} denote the set of α-successor rate functions of s with the same exit rate E. As for PAs, to check the condition s 1 p R s 2 we resort to a reduction to LP problems.Lemma 4.13.Let M = (S, Act, R, L) be a given CPA, and let R ⊆ S × S. Let (s 1 , s 2 ) ∈ R with L(s 1 ) = L(s 2 ) and that Act(s 1 ) ⊆ Act(s 2 ).Then, s 1 p R s 2 iff for each transition s 1 α − → r either r(S) = 0, or there exists E ∈ E(s 2 ) with E ≥ r(S) such that the following LP has a feasible solution, which consists of Constraints 4.4,4.5,4.6 of Lemma 4.12, and additionally: where k = Steps E α (s) with Steps E α (s) = {r 1 , . . ., r k }.Proof.The proof follows the same strategy as the proof of Lemma 4.12, in which the induced distribution of the corresponding rate function should be used.Now we are able to check Condition 4.9 by solving LP problems.For a CPA M = (S, Act, R, L), and a relation R ⊆ S × S, let (s 1 , s 2 ) ∈ R with L(s 1 ) = L(s 2 ) and Act(s 1 ) ⊆ Act(s 2 ).For s 1 α − → r 1 , and E ∈ E(s 2 ), we introduce the predicate LP ′ (s 1 , α, r 1 , s 2 , E) which is true iff E ≥ r 1 (S) and the corresponding LP problem has a solution.Then, s 1 p R s 2 iff the conjunction α∈Act(s 1 ) r 1 ∈Steps α (s 2 ) E∈E(s 2 ) LP ′ (s 1 , α, r 1 , s 2 , E) is true.The decision algorithm is given in Algorithm 7. As complexity we have to solve O(n 2 m) LP problems and each of them has at most O(n 2 ) constraints.

Algorithms for Deciding Weak Simulations
We now turn our attention to weak simulations.Similar to strong simulations, the core of the algorithm is to check whether s 1 R s 2 , i.e., s 2 weakly simulates s 1 up to the current relation R. As for strong simulation up to R, s 1 R s 2 does not imply s 1 M s 2 , since no conditions are imposed on pairs in R different from (s 1 , s 2 ).By the definition of weak simulation, for fixed characteristic functions δ i (i = 1, 2), the weight function conditions can be checked by applying maximum flow algorithms.Unfortunately, δ i -functions are not known a priori.Inspired by the parametric maximum flow algorithm, in this chapter, we show that one can determine whether such characteristic functions δ i exist with the help of breakpoints, which can be computed by analysing a parametric network constructed out of P(s 1 , •), P(s 2 , •) and R. We present dedicated algorithms for DTMCs in Subsection 5.1 and CTMCs in Subsection 5.2.

5.
1.An Algorithm for DTMCs.Let M = (S, P, L) be a DTMC.Let R ⊆ S × S be a relation and s 1 R s 2 .Whether s 2 weakly simulates s 1 up to R is equivalent to whether there exist functions δ i : S → [0, 1] such that the conditions in Definition 2.21 are satisfied.Assume that we are given the U i -characterising functions δ i .In this case, s 1 R s 2 can be checked as follows: • Concerning Condition 1a we check whether for all v ∈ S with δ 1 (v) < 1 it holds that v R s 2 .Similarly, for Condition 1b, we check whether for all v ∈ S with δ 2 (v) < 1 it holds that s 1 R v. • The reachability condition can be checked by using standard graph algorithms.In more detail, for each u with δ 1 (u) > 0, the condition holds if a state in R(u) is reachable from s 2 via R(s 1 ) states.• Finally consider Condition 2. From the given δ i functions we can compute K i .In case of that K 1 > 0 and K 2 > 0, we need to check whether there exists a weight function for the conditional distributions P(s 1 , •)/K 1 and P(s 2 , •)/K 2 with respect to the current relation R. From Lemma 4.1, this is equivalent to check whether the maximum flow for the network constructed from (P(s 1 , •)/K 1 , P(s 2 , •)/K 2 ) and R has value 1.
To check s 1 R s 2 , we want to check whether such δ i functions exist.The difficulty is that there exist uncountably many possible δ i functions.In this section, we first show that whether such δ i exists can be characterised by analysing a parametric network in Subsection 5.1.1.Then, in Subsection 5.1.2,we recall the notion of breakpoints, and show that the breakpoints play a central role in the parametric networks considered: only these points need to be considered.Based on this, we present the algorithm for DTMCs in Subsection 5.1.3.An improvement of the algorithm for certain cases is reported in Subsection 5.1.4.We introduce some notations.We focus on a particular pair (s 1 , s 2 ) ∈ R, where R is the current relation.We partition the set post (s i ) into M U i (for: must be in U i ) and P V i (for: potentially in V i ).The set P V 1 consists of those successors of s 1 which can be either put into U 1 or V 1 or both:  consists of the successor states which can only be placed in U 1 .The sets P V 2 and M U 2 are defined similarly by: P Example 5.1.Consider the DTMC in Figure 8 and the relation We say a flow function f of N (γ) is valid for N (γ) iff f saturates all edges ( , u 1 ) with u 1 ∈ M U 1 and all edges (u 2 , ) with u 2 ∈ M U 2 .If there exists a valid flow f for N (γ), we say that γ is valid for N (γ).The following lemma considers the case in which both s 1 and s 2 have visible steps: ).Then, s 1 R s 2 iff there exists a valid γ for N (γ).
Proof.By assumption, we have that s ′ i ∈ M U i for i = 1, 2, thus M U i = ∅, and it holds that δ i (s ′ i ) = 1 for i = 1, 2.
We first show the only if direction.Assume s 1 R s 2 , and let δ i , U i , V i , K i , ∆ (for i = 1, 2) as described in Definition 2.21.Since M U i = ∅ for i = 1, 2, both K 1 and K 2 are greater than 0. We let γ = K 1 /K 2 .For s, s ′ ∈ S, we define the function f for N (γ): Since δ i (s) ≤ 1 (i = 1, 2) for s ∈ S, f ( , s) ≤ P(s 1 , s) and f (s, ) ≤ γP(s 2 , s).Therefore, f satisfies the capacity constraints.f also satisfies the conservation rule: Hence, f is a flow function for N (γ).For u 1 ∈ M U 1 , we have δ 1 (u 1 ) = 1, therefore, f ( , u 1 ) = P(s 1 , u 1 ).Analogously, f (u 2 , ) = γP(s 2 , u 2 ) for u 2 ∈ M U 2 .Hence, f is valid for N (γ), implying that γ is valid for N (γ).Now we show the if direction.Assume that there exists γ > 0 and a valid flow f for N (γ).The function δ 1 is defined by: δ 1 (s) equals f ( , s)/P(s 1 , s) if s ∈ post (s 1 ) and 0 otherwise.The function δ 2 is defined similarly: δ 2 (s) equals f (s, )/γP(s 2 , s) if s ∈ post (s 2 ) and 0 otherwise.Let the sets U i and V i be defined as required by Definition 2.21.It follows that Since the amount of flow out of is the same as the amount of flow into , we have 2, both of K 1 and K 2 are greater than 0. We show that the Conditions 1a and 1b of Definition 2.21 are satisfied.For v 1 ∈ V 1 , we have that δ 1 (v 1 ) < 1 which implies that f ( , v 1 ) < P(s 1 , v 1 ).Since f is valid for N (γ), and since the edge ( , v 1 ) is not saturated by f , it must hold that v We define ∆(w, w ′ ) = f (w, w ′ )/K 1 for w, w ′ ∈ S. Assume that ∆(w, w ′ ) > 0.Then, f (w, w ′ ) > 0, which implies that (w, w ′ ) is an edge of N (γ), therefore, (w, w ′ ) ∈ R. By the flow conservation rule, f ( , w) ≥ f (w, w ′ ) > 0, implying that δ 1 (w) > 0. By the definition of U 1 , we obtain that w ∈ U 1 .Similarly, we can show that w ′ ∈ U 2 .Hence, the Condition 2a is satisfied.To prove Condition 2b: where equality ( * ) follows from the flow conservation rule.Therefore, for w ∈ S we have that K 1 ∆(w, U 2 ) = P(s 1 , w)δ 1 (w).Similarly, we can show K 2 ∆(U 1 , w) = P(s 2 , w)δ 2 (w).Condition 2b is also satisfied.As K 1 > 0 and K 2 > 0, the reachability condition holds trivially, hence, s 1 R s 2 .Observe that for γ = 6 7 , the edges to v 1 and o 1 are saturated, i.e., we have used full capacities of the edge ( , v 1 ) and ( , o 1 ).Thus, by a greater value of γ, although the capacities c({v 2 , o 2 , o 3 }, ) increase (become greater than 3  4 ), no additional flow can be sent through {v 1 , o 1 }.For the other breakpoint 2, we observe that for a value of γ ≤ 2, we can still send γ 8 through the path , u 1 , u 2 , , but for γ > 2, the edge to u 1 keeps saturated, thus the amount of flow sent through this path does not increase any more.Thus, for γ ∈ [ 6  7 , 2], the maximum value, or the value of the minimum cut, is 3  4 + γ 8 .The first term 3 4 corresponds to the amount of flow through v 1 and o 1 .The breakpoint 6  7 is not valid since the edge to u 1 can not be saturated.As discussed in Example 5.4, the breakpoint 2 is valid.The curve for κ(γ) is depicted in Figure 11.
In the following lemma we show that if there is any valid γ, then at least one breakpoint is valid.
Proof.Consider the network N (γ * ).Assume that the maximum flow f γ * is a valid maximum flow for N (γ * ).
Assume first γ ′ ∈ (γ * , γ 2 ].We use the augmenting path algorithm [1] to obtain a maximum flow f * in the residual network N f γ * (γ ′ ), requiring that the augmenting path contains no cycles, which is a harmless restriction.Then, f γ ′ := f γ * + f * is a maximum flow in N (γ ′ ).Since f γ * saturates edges from to M U 1 , f γ ′ saturates edges from to M U 1 as well , as flow along an augmenting path without cycles does not un-saturate edges to M U 1 .We choose the minimum cut (X, X ′ ) for N (γ * ) with respect to f γ * such that M U 2 ∩ X ′ = ∅, or equivalently M U 2 ⊆ X.This is possible since f γ * saturates all edges (u 2 , ) with u 2 ∈ M U 2 .The minimum cut for f γ ′ , then, can also be chosen as (X, X ′ ), as (γ ′ , κ(γ ′ )) lies on the same line segment as (γ * , κ(γ * )).Hence, f γ ′ saturates the edges from M U 2 to , which indicates that For the valid maximum flow f γ * we select the minimal cut (X, X ′ ) for N (γ * ) such that M U 1 ∩ X = ∅.Let d denote a valid distance function corresponding to f γ * .We replace f γ * (v, ) by min{f γ * (v, ), c γ ′ (v, )} where c γ ′ is the capacity function of N (γ ′ ).The modified flow is a preflow for the network N (γ ′ ).Moreover, d stays a valid Figure 12: A network in which more than one breakpoint is valid.
distance function as no new residual edges are introduced.Then, we apply the preflow algorithm to get a maximum flow f γ ′ for the network N (γ ′ ).Since no flow is pushed back from the sink, edges from M U 2 to are kept saturated.Since (γ * , κ(γ * )) and (γ ′ , κ(γ ′ )) are on the same line segment, the minimal cut for f γ ′ can also be chosen as (X, X ′ ), which indicates that f γ ′ saturates all edges to M U 1 .This implies that γ ′ is valid for In Example 5.5, only one breakpoint is valid.In the following example we show that it is in general possible that more than one breakpoint is valid.
Example 5.7.Consider the network depicted in Figure 12.By a similar analysis as Example 5.5, we can compute that there are three breakpoints Proof.If there exists a valid γ for N (γ), Lemma 5.6 guarantees that one of the breakpoints of κ(γ) is valid.The other direction is trivial.
For a given breakpoint, we need to solve one feasible flow problem to check whether it is valid.In the network N (γ) the capacities of the edges leading to the sink are increasing functions of a real-valued parameter γ.If we reverse N (γ), we get a parametric network that satisfies the conditions in [18]: The capacities emanating from are non-decreasing functions of γ.So we can apply the breakpoint algorithm [18] to obtain all of the breakpoints of N (γ).5.1.3.The Algorithm for DTMCs.Let M be a DTMC and let SimRel w (M) denote the weak simulation algorithm, which is obtained by replacing line 1.6 of SimRel s (M) in Algorithm 1 by: if s 1 R s 2 .The condition s 1 R s 2 is checked in Ws(M, s 1 , s 2 , R), shown as Algorithm 8.
The first part of the algorithm is the preprocessing part.Line 8.1 tests for the case that s 1 could perform only stutter steps with respect to the current relation R. If line 8.5 Proof.The proof follows exactly the lines of the proof of Theorem 4.2.Let iteration k be the last iteration of Ws.Then, by Lemma 5.9, for each pair (s 1 , s 2 ) ∈ R k , it holds that s 2 weakly simulates s 1 up to R k , so R k is a weak simulation.On the other hand, one can prove by induction that each R i is coarser than .
Complexity.For (s 1 , s 2 ) ∈ R, we have shown that to check whether s 1 R s 2 we could first compute the breakpoints of N (γ), then solve O(|V |) many maximum flow problems.To achieve a better bound, we first prove that applying a binary search method over the breakpoints, we only need to consider O(log |V |) breakpoints, and thus solve O(log |V |) maximum flow problems.
Assume that the sets M U i and P V i for i = 1, 2 are constructed as before for N (γ).Recall that a flow function f of N (γ) is valid for N (γ) iff f saturates all edges ( , u 1 ) with u 1 ∈ M U 1 and all edges (u 2 , ) with u 2 ∈ M U 2 .If f is also a maximum flow, we say that f is a valid maximum flow of N (γ).We first reformulate Lemma 5.2 using maximum flow.Lemma 5.11.There exists a valid flow f for N (γ) iff there exists a valid maximum flow f m for N (γ).
Proof.Assume there exists a valid flow f for N (γ).We let N f (γ) denote the residual network.We use the augmenting path algorithm to get a maximum flow f ′ in the residual network N f (γ).Assume that the augmenting path contains no cycles, which is a harmless restriction.Let f m = f + f ′ .Obviously, f m is a maximum flow for N (γ), and it saturates all of the edges saturated by f .Hence, f m is valid for N (γ).The other direction is simple, since a valid maximum flow is also a valid flow for N (γ).
We first discuss how to get a valid maximum flow provided that γ is valid.Observe that even if γ is valid for N (γ), not all maximum flows for N (γ) are necessarily valid.Consider the DTMC on the left part of Figure 13.Assume that the relation R is given by {(s 1 , s 2 ), (s 1 , v 2 ), (v 1 , s 2 ), (u 1 , u 2 ), (u 1 , v 2 )} and consider the pair (s 1 , s 2 ).Thus, we have that P The network N (γ) is depicted on the right part of Figure 13.The maximum flow f for N (1) has value 0.5.If f contains positive sub-flow along the path , u 1 , v 2 , , it does not saturate the edge (u 2 , ).On the contrary, the flow along the single path , u 1 , u 2 , with value 0.5 would be a valid maximum flow.This example gives us the intuition to use the augmenting path through edges between M U 1 and M U 2 as much as possible.For this purpose we define a cost function cost from edges in N (γ) to real numbers as follows: cost cost (s, s ′ ) = 0 otherwise.The costs of edges starting from source, or ending at sink, or in P V 1 × P V 2 are 0. The cost of a flow f is defined by cost (f ) = e∈E f (e)cost (e).By definition of the cost function, we have the property cost (f ) = f ( , M U 1 ) + f (M U 2 , ), i.e., the cost equals the sum of the amount of flow from into M U 1 and from M U 2 into .Lemma 5.12.Assume that γ > 0 is valid for N (γ).Let f γ denote a maximum flow over N (γ) with maximum cost.Then, f γ is valid for N (γ).
Proof.By Lemma 5.11, provided γ is valid for N (γ), there exists a valid maximum flow function g for N (γ).Since g saturates edges to M U 1 and from M U 2 , obviously, cost (g) = P(s 1 , M U 1 ) + γP(s 2 , M U 2 ).Assume that f γ is not valid, which indicates that f γ does not saturate an edge ( , u 1 ) with u This contradicts the assumption that f γ has maximum cost.
Let N U (γ) denote the subnetwork of N (γ) where the set of vertices is restricted to M U 1 , M U 2 and { , }.The following lemma discusses how to construct a maximum flow with maximum cost.Lemma 5.13.Assume that f * is an arbitrary maximum flow of N U (γ) and f is an arbitrary maximum flow in the residual network N f * (γ) with the residual edges from M U 1 back to removed, as well as the residual edges from back to M U 2 .Then f γ = f * + f is a maximum flow over N (γ) with maximum cost.
Proof.Recall that the cost of f γ is equal to cost(f γ ) = f γ ( , M U 1 ) + f γ (M U 2 , ).Assume that cost(f γ ) is not maximal for the sake of contradiction.Let f be a maximum flow such that cost (f γ ) < cost (f ).Without loss of generality, we assume that f γ ( , M U 1 ) < f ( , M U 1 ).It holds that f γ = f * + f where f * is a maximum flow of N U (γ), and f is a maximum flow in the residual network N f * (γ) with the residual edges from M U 1 back to removed, as well as the residual edges from back to M U 2 .On the one hand, f * sends as much flow as possible along M U 1 in N U (γ).Since in the residual network N f * (γ) edges from M U 1 back to are removed, this guarantees that no flow can be sent back to from M U 1 .On the other hand, f sends as much flow as possible from M U 1 to P V 2 (and also from P V 1 to M U 2 ) in N f * (γ).Thus, f γ ( , M U 1 ) must be maximal which contradicts the assumption f γ ( , M U 1 ) < f ( , M U 1 ).
5.1.4.An Improvement.The algorithm Ws(M, s 1 , s 2 , R) is dominated by the part in which all breakpoints (O(n) many) must be computed, and a binary search is applied to the breakpoints, with O(log n) many feasible flow problems.In this section we discuss how to achieve an improved algorithm if the network N (γ) can be partitioned into sub-networks.Let H denote the sub-relation R ∩ [(post (s 1 ) ∪ {s 1 }) × (post (s 2 ) ∪ {s 2 })], which is the local fragment of the relation R. Now let A 1 , A 2 , . . .A h enumerate the classes of the equivalence relation (H ∪ H −1 ) * generated by H, where h denotes the number of classes.W. l.o. g., we assume in the following that A h is the equivalence class containing s 1 and s 2 , i. e., s 1 , s 2 ∈ A h .The following lemma gives some properties of the sets A i provided that s 1 R s 2 : Lemma 5.16.For (s 1 , s 2 ) ∈ R, assume that there exists a state s ′ 1 ∈ post (s 1 ) such that s ′ 1 ∈ R −1 (s 2 ), and s ′ 2 ∈ post (s 2 ) such that s ′ 2 ∈ R(s 1 ).A 1 , . . ., A h be the sets constructed for (s 1 , s 2 ) as above.If s 1 R s 2 , the following hold: (1) P(s 1 , A i ) > 0 and P(s 2 , A i ) > 0 for all i < h, (2) γ i = K 1 /K 2 for all i < h where γ i = P(s 1 , A i )/P(s 2 , A i ) Proof.Since s 1 R s 2 , we let δ i , U i , V i , ∆ (for I = 1, 2) as described in Definition 2.21.Because of states s ′ 1 and s ′ 2 , we have K 1 > 0, K 2 > 0. Let post i (s j ) = A i ∩ post (s j ) for i = 1, . . ., h and j = 1, 2. We prove the first part.For i < h, the set A i does not contain s 1 nor s 2 , but only (some of) their successors, so it is impossible that both P(s 1 , A i ) = 0 and P(s 2 , A i ) = 0. W. l.o. g., assume that P(s 1 , A i ) > 0. There exists t ∈ post i (s 1 ) such that P(s 1 , t) > 0. Obviously δ 1 (t) = 1, thus: 0 < P(s 1 , t) = (K 1 ∆(t, U 2 ))/δ 1 (t) = K 1 ∆(t, U 2 ) which implies that ∃u 2 ∈ A i with ∆(t, u 2 ) > 0. Hence, P(s 2 , u 2 ) = K 2 ∆(U 1 , u 2 ) ≥ K 2 ∆(t, u 2 ) > 0. Now we prove the second part.It holds that: ∆(a i , U 2 ) where ( * ) follows from Condition 2b of Definition 2.21, (!) follows from the equation δ 1 (a i ) = 1 for all a i ∈ post i (a 1 ) with i < n, and ( †) follows from the fact that if a ∈ post i (s 1 ), then ∆(a, b) = 0 for b ∈ U 2 \post i (s 2 ).In the same way, we get P(s 2 , A i ) = K 2 • a i ∈A i ∆(A i , a i ).Therefore, For the case h > 1, the above lemma allows to check whether s 1 R s 2 efficiently.For this case we replace lines 8.6-8.7 of Ws by the sub-algorithm WsImproved in Algorithm 9.The partition A 1 , . . ., A h is constructed in line 9.1.Lines 9.2-9.10follow directly from Lemma 5.16: if γ i = γ j for some i, j < h, we conclude from Lemma 5.16 that s 1 R s 2 .Line 9.11 follows from the following lemma, which is the counterpart of Lemma 5.2: Lemma 5.17.For (s 1 , s 2 ) ∈ R, assume that there exists a state s ′ 1 ∈ post (s 1 ) such that s ′ 1 ∈ R −1 (s 2 ), and s ′ 2 ∈ post (s 2 ) such that s ′ 2 ∈ R(s 1 ).Assume that h > 1, and assume WsImproved(M, s 1 , s 2 , R) reaches line 9.11.Then, s 1 R s 2 iff γ 1 is valid for N (γ 1 ).
Proof.First, assume that s 1 R s 2 .According to Lemma 5.2, there exists a valid γ * for N (γ * ).As in the proof of Lemma 5.2, γ * = K 1 /K 2 is valid for N (γ * ).If Ws reaches line 9.11, by Lemma 5.16, we have γ 1 = K 1 /K 2 , hence, γ 1 is valid for N (γ 1 ).The other direction follows directly from Lemma 5.2. each such pair (s 1 , s 2 ) we hope to get h i ′ > 1 because some other elements of R i have been eliminated from it.We only perform the expensive computation if an incomplete iteration does no longer refine the relation.

5.
2. An Algorithm for CTMCs.Let M = (S, R, L) be a CTMC.We now discuss how to handle CTMCs.Recall that in Definition 2.23, we have the rate condition 3 ′ : K 1 R(s 1 , S) ≤ K 2 R(s 2 , S).To determine M , we simplify the algorithm for DTMCs.If K 1 > 0 and K 2 = 0, s 1 R s 2 because of the rate condition.Hence, we do not need check the reachability condition, and lines 8.3-8.5 of the algorithm Ws(M, s 1 , s 2 , R) can be skipped.For states s 1 , s 2 and relation R, we use N (γ) to denote the network defined in the embedded DTMC emb(M).To check the additional rate condition we use the following lemma: Lemma 5.19.Let s 1 R s 2 .Assume that there exists s ′ 1 ∈ post(s 1 ) such that s ′ 1 ∈ R −1 (s 2 ).Then, s 1 R s 2 in M iff there exists γ ≤ R(s 2 , S)/R(s 1 , S) such that γ is valid for N (γ).
Proof.Assume first s 1 R s 2 in M. Let δ i , U i , V i , K i , ∆ (for i = 1, 2) as described in Definition 2.23.Obviously, s ′ 1 must be in U 1 , implying that K 1 > 0. Because of the rate condition it holds that K 1 R(s 1 , S) ≤ K 2 R(s 2 , S), which implies that K 2 > 0. It is sufficient to show that γ := K 1 /K 2 is valid for N (γ).Exactly as in the proof of Lemma 5.2 (the only if direction), we can construct a valid flow f for N (γ).Thus, γ is valid for N (γ) and γ ≤ R(s 2 , S)/R(s 1 , S).Now we show the other direction.By assumption, γ is valid for N (γ).We may assume that there exists a valid flow function f for N (γ).We define δ i , V i , U i , K i , ∆ (for i = 1, 2) as in the proof (the if direction) of Lemma 5.2.Recall that s ′ 1 must be in U 1 , implying that f ( , s ′ 1 ) > 0. Thus, there must be a node s in N (γ) with f (s, ) > 0, which implies that s ∈ U 2 .Thus we have K 2 > 0. Using the proof (the if direction) of Lemma 5.2, it holds that s 1 R s 2 in emb(M), moreover, it holds that γ = K 1 /K 2 .By assumption it holds that γ ≤ R(s 2 , S)/R(s 1 , S) which is exactly the rate condition.
To check the rate condition for the case h > 1, we replace line 9.11 of the algorithm WsImproved by: return (γ 1 ≤ γ * ∧ γ 1 is valid for N (γ 1 )) where γ * = R(s 2 , S)/R(s 1 , S) can be computed directly.In case h = 1, we replace line 8.7 of Ws by: return (∃i ∈ {1, . . ., j}. b i ≤ γ * ∧ b i is valid for N (b i )) to check the rate condition.Or, equivalently, we can check whether the minimal valid breakpoint γ m is smaller than or equal to γ * .The binary search algorithm introduced for DTMCs can also be modified slightly to find the minimal valid breakpoint.The idea is that, if we find a valid breakpoint, we first save it, and then continue the binary search on the left side.If another breakpoint is valid, we save the smaller one.As the check for the reachability condition disappears for CTMCs, we get better bound for sparse CTMCs: which is bounded by 2kgm 2 .Since k is bounded by n 2 , the time complexity is bounded by 4gm 2 n 2 .If g is a constant, we have m ≤ gn, hence, the time complexity is 4g 3 n 4 ∈ O(n 4 ).

Conclusion and Future Work
In this paper we have proposed efficient algorithms deciding simulation on Markov models with complexity O(m 2 n).For sparse models where the fanout is bounded by a constant, we achieve for strong simulation the complexity O(n 2 ) and for simulation O(n 4 ) on CTMCs and O(n 5 ) for DTMCs.We extended the algorithms for computing simulation preorders to handle PAs with the complexity O(m 2 n) for strong simulation.For strong probabilistic simulation, we have shown that the preorder can be determined by solving LP problems.We also considered their continuous-time analogon, CPAs, and arrived at an algorithm with same complexities as for PAs.

Figure 4 :
Figure 4: Splitting of successor states for weak simulations.

Algorithm 1 :
Basic algorithm to decide strong simulation.

Algorithm 3 :
Improved algorithm for deciding strong simulation for FPSs.

Theorem 4 . 6 .
The algorithm SimRel FPS s (M) runs in time O(m 2 n) and in space O(m 2 ).If the fanout is bounded by a constant, it has complexity O(n 2 ), both in time and space.

Example 4 . 7 .
Consider the CTMC in the left part of Figure7(it has 10 states).Consider the pair (s 1 , s 2 ) ∈ R 1 .The network N (s 1 , s 2 , R 1 ) is depicted on the right of the figure. s

FPS s in Theorem 4 . 5 .
The only new element is that we now have to quantify over the actions and successor distributions as prescribed by Definition 2.11.This translates to a conjunction in lines 5.8 and 5.21 of the algorithm.Exploiting Lemma 4.8 we get the correctness.Now we give the complexity of the algorithm: Theorem 4.10.The algorithm SimRel PA s (M) runs in time O(m 2 n) and in space O(m 2 ).If the fanout of M is bounded by a constant, it has complexity O(n 2 ), both in time and space.

Remark 4 . 11 .
Let M = (S, Act, P, L) be a PA, and let k = s∈S α∈Act(s) |Steps α (s)|, called the number of transitions in [3], denote the number of all distributions in M. The algorithm for deciding strong simulation introduced by Baier et al. has time complexity O((kn 6 + k 2 n 3 )/ log n), and space complexity O(k 2 ).The number of distributions k and the size of transitions m are related by k ≤ m ≤ nk.The left equality is established if |µ| = 1 for all distributions, and the right equality is established if |µ| = n for all distributions.

Algorithm 7 :
Algorithm for deciding strong probabilistic simulation for CPAs.

Example 5 . 3 . 1 4 1 8Figure 11 :
Figure 11: The value of the maximum flow, or equivalently the value of the minimum cut, as a function of γ for the network in Figure 10.

Figure 13 :
Figure 13: A simple DTMC for illustrating that not all maximum flows are valid.

Theorem 5 . 20 .
If the fanout g of M is bounded by a constant, the time complexity for CTMC is O(n 4 ).Proof.In the proof of Theorem 5.15 we have shown that Ws has complexity O(|V 1 | |V 2 | 2 ).As we do not need to check the reachability condition, the overall complexity of the algorithm SimRel w (M) (see Inequality 5.1) is s 1 ∈S s 2 ∈S l i=1 |post (s 1 )| + |V 1 | |V 2 | 2 2.2.4.Simulation up to R. For an arbitrary relation R on the state space S of an FPS with s 1 R s 2 , we say that s 2 simulates s 1 strongly up to R, denoteds 1 R s 2 , if L(s 1 ) = L(s 2 ) and P(s 1 , •) ⊑ R P(s 2 , •).Otherwise we write s 1 R s 2 .Since only the first step is considered for R , s 1 R s 2 does not imply s 1 M s 2 unless R is a strong simulation.By definition, R is a strong simulation if and only if for all s 1 R s 2 it holds that s 1 R s 2 .Likewise, we say that s 2 simulates s 1 weakly up to R, denoted by s 1 R s 2 , if there are functions δ i and U i , V i , ∆ as required by Definition 2.21 for this pair of states.Otherwise, we write s 1 R s 2 .Similar to strong simulation up to R, s 1 R s 2 does not imply s 1 M s 2 , since no conditions are imposed on pairs in R different from (s 1 , s 2 ).Again, R is a weak simulation if and only if for all s 1 R s 2 it holds that s 1 R s 2 .These conventions extend to DTMCs, CTMCs, PAs and CPAs in an obvious way.For PAs and CPAs, strong probabilistic simulation up to R, denoted by p R , is also defined analogously.Example 2.24.Consider the FPS in Figure 1 ,s 2 ) is bounded by |pre(s 1 )| |pre(s 2 )| where pre(s) = {t ∈ S | P(t, s) > 0}.Summing over all (s 1 , s 2 ), we get s 1 ∈S s 2 ∈S |pre(s 1 )| |pre(s 2 )| = m 2 .Assume we run iteration i.For every pair (s 1 , s 2 ), we generate the set D (s 1 , s 2 , R i ) together with f i and d i .Obviously, the size of D )| |post (s 2 )|.Summing over all (s 1 , s 2 ), we get the bound O(m 2 ).The number of edges of the network N (s 1 , s 2 , R 1 ) (together with f i and d