{'fraction_non_alphanumeric': 0.035871318761269096, 'fraction_numerical': 0.008730585518615759, 'mean_word_length': 4.347428958051421, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated Model Checking of Parametric Markov Chains

Parametric Markov chains occur quite naturally in various applications: they can be used for a conservative analysis of probabilistic systems (no matter how the parameter is chosen, the system works to specification); they can be used to find optimal settings for a parameter; they can be used to visualise the influence of system parameters; and they can be used to make it easy to adjust the analysis for the case that parameters change. Unfortunately, these advancements come at a cost: parametric model checking is---or rather was---often slow. To make the analysis of parametric Markov models scale, we need three ingredients: clever algorithms, the right data structure, and good engineering. Clever algorithms are often the main (or sole) selling point; and we face the trouble that this paper focuses on -- the latter ingredients to efficient model checking. Consequently, our easiest claim to fame is in the speed-up we have often realised when comparing to the state of the art.


Introduction
The analysis of parametric Markov models is a young and growing field of research. As not only the research direction but also the term 'parametric Markov models' is attractive, it has been used for various generalisations of traditional Markov models. We use Markov chains, where the parameter is used to determine the probabilities and rewards, such that we can reason about the likelihood of obtaining simple temporal properties like safety and reachability as well as standard reward functions, such as long-run average.
What we do not intend to do in this paper is to use parameters to change the size of the system or the shape of the Markov chain. (The latter can, of course, be encoded by using parameters to assign a probability of 0 to an edge, effectively removing it. This would, however, come at the cost of efficiency and is not what we want to use the parameters for.) Using parameters to describe the probabilities of transitions is not quite as easy as it sounds: even when parameters appear in a simple way, like 'p' or '1 − p', the terms that represent the likelihood of obtaining a temporal property or an expected reward can quickly become quite intricate. One ends up with rational functions. We make a virtue of necessity by using this as a motivation to allow for using rational functions of the occurring parameters to represent the probabilities and payoffs.
To allow for an efficient analysis of such complex parametrised systems, we have taken a look at different strategies for the evaluations of-parametric and non-parametric-Markov chains, and considered their suitability for our purposes. We found the stepwise elimination of vertices from a model to be the most attractive approach to port.
Broadly speaking, this approach works like the transformation from finite automata to regular expressions: a vertex is removed, and the new structure has all successors of this state as-potentially new-successors of the predecessors of this vertex. In the transformation from finite automata to regular expressions, one changes the expressions on the edges, while we adjust the probabilities and, if applicable, the rewards on the edges.
When using this approach with explicit probabilities and rewards, one ends up with a Directed Acyclic Graph (DAG) structure in the evaluation. This DAG structure has been exploited to reduce the cost of re-calculating the probabilities for simple temporal properties or expected rewards, and it proves that it also integrates nicely into our framework, where the probabilities and rewards are provided as rational functions. In fact it integrates so naturally that it seems surprising in hindsight that it has not been discovered earlier.
The natural connection occurs when choosing a similar data structure to represent the rational functions that represent the probabilities and rewards. To make full use of the DAG structure that comes with the elimination, we represent these functions in the form of arithmetic circuits-which are essentially DAGs. We have integrated the resulting representation organically in a small extension of ePMC, and tested it on a range of case studies. We have obtained a speed-up of a hefty factor of 20 to 120 when compared to storing functions in terms of coprime numerator and denominator polynomials.
Related work. For (discrete-time) Markov chains (MCs), Daws [6] has devised a language-theoretic approach to solve this problem. In this approach, the transition probabilities are considered as letters of an alphabet. Thus, the model can be viewed as a finite automaton. Then, based on the state elimination method [22], a regular expression that describes the language of such an automaton is calculated. In a post-processing step, this regular expression is recursively evaluated, resulting in a rational function over the parameters of the model. One of the authors has been involved in extending and tuning this method [16] so as to operate with rational functions, which are stored as coprime numerator and denominator polynomials rather than with regular expressions.
The process of computing a function that describes properties (like reachability probabilities or long-run average rewards) that depend on model parameters is often costly. However, once the function has been obtained, it can very efficiently be evaluated for given parameter instantiations. Because of this, parametric model checking of Markov models has also attracted attention in the area of runtime verification, where the acceptable time to obtain values is limited [3,11].
Other works in the area are centred around deciding the validity of boolean formulas depending on the parameter range using SMT solvers or extending these techniques to models that involve nondeterminism [7,14,5,27].  As an example for a parametric model, consider Figure 1. Knuth and Yao [24] have shown how a six-sided dice can be simulated by repeatedly tossing a coin. The idea is to build a Markov chain with transition probabilities of only 0.5 or 1. Borrowing a model from the PRISM website, we have extended this example to a biased dice, simulated by tossing a biased coin. With probability x we see heads, while with probability 1 − x we see tails. This way, we move around in the Markov chain until we obtain a result.
Organisation of the paper. After formalising our setting in Section 2, we describe how we exploit DAGs in the representation of rational functions, and exploit them using synergies with the DAG-style state elimination technique, in Section 3. We then describe how to expand this technique to determine long-run average rewards in Section 4. In Section 5, we evaluate our approach on a range of benchmarks, and discuss the results briefly in Section 6.

Parametric Markov chains with state rewards
Let V = {v 1 , . . . , v n } denote a set of variables over R. A polynomial g over V is a sum of monomials where each i j ∈ N and each a i1 , . . . , in ∈ R. A rational function f over a set of variables V is a fraction f (v 1 , . . . , v n ) = f1(v1,...,vn) f2(v1,...,vn) of two polynomials f 1 , f 2 over V . We denote the set of rational functions from V to R by F V .
where S is a finite set of states, s is the initial state, V = {v 1 , . . . , v n } is a finite set of parameters, and P is the probability matrix P : S × S → F V . A path ω of a PMC D = (S, s, P, V ) is a non-empty finite, or infinite, sequence s 0 , s 1 , s 2 , . . . where s i ∈ S and P(s i , s i+1 ) > 0 for i 0. We let Ω denote the set of infinite paths. With Pr s , we denote the parametric probability measure over Ω assuming that we start in state s, with Pr = Pr s . We use Exp s , Exp to while E = ∅ do 5: se ← choose(E) 6: for all s ∈ post D (se) do 8: P(se, s) ← P(se, s)/(1 − P(se, se)) 9: end for 10: P(se, se) ← 0 11: for all (s1, s2) ∈ pre D (se) × post D (se) do 12: P(s1, s2) ← P(s1, s2) + P(s1, se)P(se, s2) 13: end for 14: if se = s ∧ se / ∈ B ∧ post D (se) = ∅ then 15: Eliminate Given a state s, we denote the set of all immediate predecessors and successors of s in the underlying graph of D by pre D (s) and post D (s), respectively, excluding s itself. We write reach D (s, s ′ ) if s ′ is reachable from s in the underlying graph of D.
Given a PMC D = (S, s, P, V ) we are interested in computing the function that represents the probability of reaching some set of target states B ⊂ S.
Our base algorithm to obtain this value is described in Algorithm 1. A state s e ∈ S is selected, and then eliminated by considering each pair (s 1 , s 2 ) ∈ pre D (s e ) × post D (s e ) and updating the existing probability P(s 1 , s 2 ) by the probability of reaching s 2 from s 1 via s e . Heuristics to determine the order in which states are chosen for elimination by the choose function are discussed in Section 5.5.

Definition 3. A parametric reward function for a PMC
The reward function labels states in D with a rational function over V that corresponds to the reward that is gained if that state is visited. Given a PMC D = (S, s, P, V ) and a reward function r : S → F V , we are interested in the parametric expected accumulated reward defined as r(X(D)) s,i or a variation [25], the parametric expected accumulated reachability reward given B ⊆ S defined as This can, however, be transformed to the former.
Algorithm 1 can be extended to compute the parametric expected accumulated reward. In addition to updating the probability matrix for each predecessor and successor pair, we also update the reward function as follows: The updated value for r(s 1 ) reflects the reward that would be accumulated if a transition would be taken from s 1 to s e , where the expected number of self-loops would be P(se,se) 1−P(se,se) . Upon termination, the algorithm returns the value r(s).

Representing Formulas using Directed Acyclic Graphs
In existing tools for parametric model checking of Markov models, rational functions have traditionally been represented in the form f (v 1 , . . . , v n ) = f1(v1,...,vn) f2(v1,...,vn) , where f 1 (v 1 , . . . , v n ) and f 2 (v 1 , . . . , v n ) [15,8,26] are coprime. As a result, for some cases the representations of such functions are very short. Often, during the state elimination phase, large common factors can be cancelled out, such that one can operate with relatively small functions throughout the whole algorithm. There are, however, many cases without-or with very few-large common factors. The nominator-denominator representations then become larger and larger during the analysis. In this case, the analysis is slowed down severely, mostly by the time taken for the cancellation of common factors. Cancelling out such factors is non-trivial, and indeed a research area in itself. In addition, if formulas become large, this can also lead to out-of-memory problems.
To overcome this issue, we propose the representation of rational functions by arithmetic circuits. These arithmetic circuits are directed acyclic graphs (DAG).
Terminal nodes are labelled with either a variable of the set of parameters V , or with a rational number. Non-final nodes are labelled with a function to be applied on the nodes it has edges to. In our setting, we require two unary functions, additive inverse and multiplicative inverse, and two binary functions, addition and multiplication. All functions used are represented using a single DAG, and a function is represented by a reference to a node of this common DAG.
This representation has two advantages. Firstly, all operations are practically constant time: to apply an operator on two functions, one simply introduces a new node labelled with the according operator, with edges pointing to the two nodes to connect. In particular, we do not have to use expensive methods to cancel out common factors. Secondly, because we are using a DAG and not a tree, common sub-expressions can be shared between different formulas, which is not possible when representing rational functions in terms of two polynomials represented as a list of monomials. For illustration, let us consider the example from Figure 1. We analyse the probability that the final result is . This probability can be described by the function In our DAG-based representation, we would represent the function as in Figure 2; When operating with arithmetic circuits, there are a number of ways to reduce their memory footprint, which will, however, lead to a higher running time. The simplest one is that, while creating a new node to represent a function, it might turn out that there already exists a node with exactly the same operator, and exactly the same operand. In this case, it is better to drop the newly created node and use a reference to the existing node to counter the growth of the DAG. In case we use hash maps for the lookup, we can also still keep the overhead close to constant time. Another optimisation is to use simple algebraic equivalences. This includes computing the values of constant functions. E.g. instead of creating a node representing 2 + 3 we introduce a new terminal node labelled 5, and if we are about the create a new node for y + x but we already have a node for x + y we reuse this node instead. We also take the additive and multiplicate neutral elements into account (rather than creating a new node for 0 + x, we return the one for x, and the like). Another optimisation method is to evaluate functions of the DAG at random points and then to identify functions if the result of this evaluation is the same. Using the Schwartz-Zippel Lemma [29,30,9], we can then bound the probability that we mistakenly identify two nodes although they do not represent the same function. We can again minimise the overheads incurred by this method by using hash maps.
Arithmetic circuits sometimes become very large, consisting of millions of nodes. This way, they cannot serve as a concise, human-readable description of the analysis result. Compared to performing a non-parametric analysis, it is however often still beneficial to obtain a function representation in this form. Even for large dags, obtaining evaluating parameter instantiations is very fast, and linear in the number of DAG nodes. This is useful in particular if a large number of points is required, for instance for plotting a graph. In this case, results can be obtained much faster than using non-parametric model checking, as demonstrated in Section 5. In particular, for any instantiation, values can be obtained in the same, predictable, time. This is quite in contrast to value iteration, where the number of iterations required to obtain a certain precision varies with the concrete values of parameters.
For this reason, parametric model checking is particularly useful for online model checking or runtime verification [3]. Here, one can precompute the DAG before running the actual system, while concrete values can be instantiated at runtime, with a running time that can be precisely calculated offline. Using arithmetic circuits expands the range of systems for which this method is applicable. Evaluation of parameter instantiations can be performed using exact arithmetic or floating-point arithmetic. From our experience, the quality of the floating-point results using DAGs is often better than the one using the representation of rational functions as coprime numerator and denominator, which has been used so far in known implementations. The reason is that, in the latter approach, one often runs into numerical problems such as cancellation, which often forces the use of expensive exact arithmetic to be used for evaluation. The DAG-based method seems to be more robust against such problems.
It has recently come noted by the verification community that the usual way in which value iteration is implemented is not safe, and solutions have already been proposed [1]. While this solves the problem, it requires more complex algorithms and leads to increased model checking time. In case arithmetic circuits are used, it is easy to obtain conservative upper and lower bounds for parameter instantiations. One only has to use interval arithmetic and provide implementations for the basic operations used (addition, multiplication, additive and multiplicative inverse). The increase in the time to evaluate functions is small. In our experiments, the largest interval diameter we have obtained is around 10 −13 . , LRA(D, r u , r l )=LRA(D, r u , r l , s).

Computation of Fractional Long-Run Average Values
In a simple case, r l (·) = 1, which means that we compute the long-run average reward where each step is assumed to take the same amount of time. Solution methods for this property has been implemented (but to the best of our confidence, not been published) for parametric models in PRISM and Storm. The fractional long-run average reward is more general and allows to express values like the average energy usage per task performed more easily. Given a reward structure r : S → F V , we define the recurrence reward as It is known [4] that this value is the same for all states of a BSCC. Furthermore, for r l (·) = 1 we have Return(D, r u , s) Return(D, r l , s) = LRA(D, r u , r l , s), which immediately extends to the general case.
In Section 2, we have discussed how state elimination can be used to obtain values for the expected accumulated reward values. For this, we have repeatedly eliminated states so as to bring the PMC of interest into a form in which reward values can be obtained in a trivial way. It is easy to see that the transformations for the expected accumulated rewards also maintains the recurrence rewards. After having handled each state of our model, we have two possible outcomes.
· · · P = p 1 P = p n P = 1, r u = u s , r l = l s P = 1, r u = u 1 , r l = l 1 P = 1, r u = u n , r l = l n In the simpler case, the remaining model consists of the initial state s with a self-loop with probability one and r u = u s , r l = l s . In this case, we have LRA(D, r u , r l ) = us ls . In the other case, the remaining model consists of the initial state s which has a probability of p i to move to one of the other n remaining states s i , i = 1, . . . , n, which all have a self-loop with probability one and r u (s i ) = u i , r l (s i ) = l i . In this case, we have LRA(D, r u , r l ) = i=1,...,n p i ui li .

Experiments
We now consider four case studies that illustrate the efficiency and scalability of our approach. Three models [21,23,28] are taken from the PRISM benchmark suite 1 , and the last is taken from the authors' work on synchronisation protocols [12,13]. All experiments were conducted on a PC with an Intel Core i7-2600 (tm) processor at 3.4GHz, equipped with 16GB of RAM, and running Ubuntu 16.04. For each case study we compare the performance times obtained for model analysis when using the parametric engine of the model checker ePMC [17] 2 , using either polynomial fractions or DAGs to represent the functions corresponding to transition probabilities and state rewards. Basically, the DAG is implemented as an array of 64-bit integers. Functions are represented as indices to this array. 4 bits describe the type of the node. For terminal nodes, the remaining bits denote the parameter or number used. For non-terminal nodes, 2 × 30 bits are used to refer to the operands within the DAG. We also compare our results to those obtained using the parametric engine of PRISM [26], and the parametric and sampling engines of Storm [8] 3 . Given a parametric model, and a set of valuations for its parameters, we are interested in the total time taken to check some property of interest for every valuation for the parameters. Since our primary concern is the efficiency of multiple evaluations of an existing model, we omit model construction times and restrict our analysis to the total time taken for the evaluation of all parameter valuations. For the parametric engines of ePMC, PRISM, and Storm, we record the total time taken for both state elimination and the evaluation of the resulting function for all parameter valuations. For the sampling engine of Storm, we record the total time taken for value iteration, using default settings to determine convergence. For Storm, we set the precision to 10 −10 rather than the default of 10 −6 . This had a very minor influence on the runtime, and allowed a better comparison to ePMC, the results of which have a precision of < 10 −13 .

Crowds Protocol
The Crowds protocol [28] provides anonymity for a crowd consisting of N Internet users, of whom M are dishonest, by hiding their communication via random routing, where there are R different path reformulates. The model is a PMC parametrised by B = M M+N , the probability that a member of the crowd is untrustworthy, and P , the probability that a member sends a package to a randomly selected receiver. With probability 1−P the packet is directly delivered to the receiver. The property of interest is the probability that the untrustworthy members observe the sender more than they observe others.    values for B, P taken from 0.002, 0.004, . . . , 0.998. There is a substantial increase in the performance of ePMC when using non-simplified DAGs (ePMC(D)), and using DAGs (ePMC(DS)) simplified by evaluating random points (cf. Section 3), instead of polynomial fractions (ePMC) to represent functions. Here, ePMC clearly outperforms the parametric engines of both PRISM and Storm. In some instances, ePMC turns out to be the fastest choice, while the sampling engine of Storm proves to be faster for other instances. Processes that exceeded the time limit of one hour are indicated by -T-, and processes that ran out of memory are indicated by -M-. In Fig. 4 (left) we plot the results for N = 5 and R = 7.

Bounded Retransmission Protocol
The bounded retransmission protocol [21] divides a file, which is to be transmitted, into N chunks. For each chunk, there are at most MAX retransmissions over two lossy channels K and L that send data and acknowledgements, respectively.  Table 2: Performance statistics for bounded retransmission protocol.
The model is a PMC parametrised by pK and pL, the reliability of the channels. We are interested in the probability that the sender reports an unsuccessful transmission after more than 8 chunks have been sent successfully. The performance statistics for different values of N and MAX are shown in Table 2, where each entry shows the total time taken to check all pairwise combinations of values for pK, pL taken from 0.002, 0.004, . . . , 0.998. Here, ePMC with DAGs again has the best performance: the running time remains approximately constant when using this data structure, even for much larger problem instances. In contrast, the running time for both engines of Storm scale linearly. Both the parametric engines of PRISM and ePMC, with polynomial fraction representation, run out of memory for all larger problem instances. Figure 4 (right) plots the results obtained for N = 256 and M AX = 4. As we see the probability of interest first increases with increasing channel reliability, but then decreases again. The reason is that, on the one hand, if the channel reliability is low, then we do not send many chunks successfully. On the other hand, if the channel reliability is high, then it is unlikely that the transmission will fail in the end.

Cyclic Polling Server
This cyclic server polling model [23] is a model of a network, described as a continuous-time Markov chain. There are two parameters, µ and γ. The model consists of one server and N clients. When a client is idle, then a new job arrives at this client with a rate of µ/N . The server 'polls' the clients in a cyclic manner. At each point of time, it observes a single client. If there is a job waiting for a given client, the server servers its job (provided there is one) with a rate of µ. When the client it observes is idle, then the server moves on to observe the next client with a rate of γ. Even though our method targets discrete-time models, we can handle this model by computing the embedded DTMC.
In this case study, we consider the probability that, in the long run, Station 1 is idle. That is, the expected limit average of the time that Station 1, or, due to symmetry, any other station, is idle. We compute this long-run average value using the method described in Section 4. Probabilities are displayed as a   function of the parameters in Figure 5, and Table 3 shows how the various tools perform on this benchmark. With increasing γ the likelihood that Station 1 is idle increases: if we increase γ, then the server will more quickly find stations to be served. As the long-run average idle time only depends on the rate between µ and γ, the likelihood that Station 1 is idle falls with increasing µ. For the current configuration, classic parametric model checking does not seem to be advantageous. Using our DAG-based implementation, however, is much more efficient than classic parametric model checking, but it is space consuming. With the chosen number of parameter instantiations, our method does not quite compete with non-parametric model checking.

Oscillator Synchronisation
The models of [12,13] encode the behaviour of a population of N coupled nodes in a network. Each node has a clock that progresses, cyclically, through a range of discrete values 1, . . . , T . At the end of each clock cycle a node transmits a message to other nodes in the network. Nodes that receive this message adjust their clocks to more closely match those of the firing node. The model is a PMC, parametrised by the likelihood ML that a firing message is lost in the communication medium.   Table 4 shows the results for different values of N and T , where each entry shows the total time taken to check all values of ML taken from 10 −5 , 2 · 10 −5 , . . . , 1 − 10 −5 . Figure 5 (right) plots the results obtained for N = 6 and T = 8. For extremal values of ML, the network is expected to use much more energy to synchronise, because the expected time required for this to occur increases. Very high values of ML result in nearly all firing messages being lost, and hence nodes cannot communicate well enough to coordinate, while very low values of ML lead to perpetually asynchronous states for the network, an artefact of the discreteness of the clock values [12].
In this case study, the DAG-based method, in particular with random points evaluation, performs best, followed by the sampling-based method of Storm. We note that the time required for each value iteration is relatively high, while the cost of evaluating a point for the DAG-based method is quite low. Therefore, the advantage of our approach would have been even more pronounced, if we had evaluated more instantiations in the experiments above. The method of choice thus depends mostly on whether such a high number of instantiations is required.
We have performed value iteration with a (local) precision of 10 −10 for Storm. This does, however, not guarantee any global precision [1]. Obtaining guaranteed results using value iteration is relatively expensive while, as discussed in Section 3, extending our approach to obtain conservative guarantees is relatively simple-and inexpensive-to achieve by using basic interval arithmetic.

Heuristics
An important consideration when performing state elimination is the order, in which different states are eliminated from the graph. Using different elimination orders to evaluate the same model can result in functions, whose representations  (nominator-denominator or DAGs) vary greatly in size, and hence also in the corresponding memory footprint and analysis time. Heuristics for efficient state elimination have been studied in automata theory, to obtain shorter regular expressions from finite-state automata [19,18], and in graph theory, for efficient peeling of a probabilistic network [20]. We employ the following heuristics, consisting of both existing schemes taken from the literature, and novel schemes that prove to be effective for some models.
-NumNew: each state is weighted by the number of new transitions that are introduced to the model when that state is eliminated. That is, we consider each predecessor-successor pair for that state, and add one to the weight if the transition from the predecessor to the successor was not already defined in the underlying graph before state elimination. States with the lowest weight are eliminated first. The aim here is to minimise the total number of transitions as elimination progresses. -MinProd: similarly to NumNew, we consider each predecessor-successor pair.
However, one is added to the weight irrespective of whether that transition already existed in the underlying graph. Again states with the lowest weight are considered first. -TargetBFS: states are eliminated in the order in which they are discovered when conducting a breadth-first search backwards from the target states. -Random: a state is selected uniformly at random for elimination from the set of remaining states. -BFS: states are eliminated in the order in which they are discovered when conducting a breadth-first search from the initial state(s) of the model. -ReverseBFS: similar to BFS, except states are eliminated in reverse order.
In Table 5, we compare the different heuristics described. We have applied each of them for each considered model, and provide the time in seconds required for medium-sized instances. As seen, it turns out that TargetBFS is in general a good choice. In one case, however, NumNew turns out to be faster.

Conclusion and Future Work
We have implemented an approach for the evaluation of parametric Markov chains that exploits the synergies of using DAGs in a state-elimination based analysis and using DAGs in an encoding of rational functions as arithmetic circuits. Our experimental evaluation suggests that these two approaches integrated so seamlessly that they often provide a notable speedup. The nicest observation is that this seems so natural in hindsight that it is almost more surprising that this has not been attempted before than that it works so well. We therefore hope to have discovered one of these simple and natural approaches that will stand the test of time.
The next step in exploiting our approach could be an integration into applications. One of the applications we have in mind is to use it in the context of parameter extraction, which we expect to work similar to Model extraction, for online Model checking. The growing knowledge of the model can be used to refine or adjust the parameters in this application. Our application can help to provide the speed required to make the approach scale, and to keep the analysis and, if required, the visualisation 4 of the effect of the learnt parameters (and the confidence area around them) efficient.
We also note that interval arithmetic could be used to evaluate boxeshyperrectangles [a 1 , b 1 ]×· · · [a n , b n ] of parameter ranges-so as to obtain bounds on the lower and upper values taken by any occurring function value in the box. This approach could be used instead of using SMT solvers (as in [7,14]) to decide PCTL properties. A similar approach to avoid using SMT solvers has been proposed [27], which is however not based on computing a function depending on the parameters but on value iteration. We assume that the DAG-based approach will perform better when a high coverage of the parameter space is required.
It would also be straightforward to parallelise evaluation of points using SIMD approaches such as GPGPU.