Parameter Synthesis for Markov Models

Markov chain analysis is a key technique in reliability engineering. A practical obstacle is that all probabilities in Markov models need to be known. However, system quantities such as failure rates or packet loss ratios, etc. are often not---or only partially---known. This motivates considering parametric models with transitions labeled with functions over parameters. Whereas traditional Markov chain analysis evaluates a reliability metric for a single, fixed set of probabilities, analysing parametric Markov models focuses on synthesising parameter values that establish a given reliability or performance specification $\varphi$. Examples are: what component failure rates ensure the probability of a system breakdown to be below 0.00000001?, or which failure rates maximise reliability? This paper presents various analysis algorithms for parametric Markov chains and Markov decision processes. We focus on three problems: (a) do all parameter values within a given region satisfy $\varphi$?, (b) which regions satisfy $\varphi$ and which ones do not?, and (c) an approximate version of (b) focusing on covering a large fraction of all possible parameter values. We give a detailed account of the various algorithms, present a software tool realising these techniques, and report on an extensive experimental evaluation on benchmarks that span a wide range of applications.

Markov models play a central role in reliability engineering [1].Markov chain (MC) analysis is employed to compute reliability measures such as the mean time between failures in fault trees [2], [3] and the probability of a system breakdown within a time limit.Numerical as well as simulative approaches are used.In addition, reliability engineers have to make numerous decisions that affect future actions such as when to inspect, repair, or replace faulty components [4].To support decision making over possibly multiple objectives, Markov decision processes (MDPs) are used, as already argued in the original paper on MDPs by Bellman in 1957 [5].
A major practical obstacle is that classical Markov chain analysis requires that all probabilities (or rates) in the Markov model are precisely known a priori.In many cases, this 1 For a complete and formal treatment, we refer to Sect.II.

arXiv:1903.07993v1 [cs.LO] 16 Mar 2019
assumption is too severe.System quantities such as component fault rates, molecule reaction rates, packet loss ratios, etc. are often not, or at best partially, known.Let us give a few examples.The quality of service of a (wireless) communication channel may be modelled by e.g., the popular Gilbert-Elliott model, a two-state Markov chain in which packet loss has an unknown probability depending on the channel's state [6].Other examples include the back-off probability in CSMA/CA protocols determining a node's delay before attempting a transmission [7], the bias of used coins in self-stabilising protocols [8], [9], and the randomised choice of selecting the type of time-slots (sleeping, transmit, or idle) in the birthday protocol, a key mechanism used for neighbour discovery in wireless sensor networks [10] to lower power consumption.
The probabilities in all these systems are deliberately left unspecified.They can later be determined in order to optimise some performance or dependability measure.Likewise, in early stages of reliable system design, the concrete failure rate of components [11] is left unspecified.Optimally, analyses in this stage may even guide the choice of a concrete component from a particular manufacturer.

B. Parametric probabilistic models
What do these examples have in common?The random variables for packet loss, failure rate etc. are not fully defined, but are parametric.Whether a parametric system satisfies a given property or not-"is the probability that the system goes down within k steps below 10 −8 "-depends on these parameters.Relevant questions are then: for which concrete parameter values is such a property satisfied-the (parameter) synthesis problem-and, in case of decision-making models, which parameter values yield optimal designs?That is, for which fixed probabilities do such protocols work in an optimal way, i.e., lead to maximal reliability, maximise the probability for nodes to be discovered, or minimise the time until stabilisation, and so on.These questions are intrinsically hard as parameters can take infinitely many different values that, in addition, can depend on each other.
This paper faces these challenges and presents various algorithmic techniques to treat different variations of the (optimal) parameter synthesis problem.To deal with uncertainties in randomness, parametric probabilistic models are adequate.These models are just like Markov models except that the transition probabilities are specified by arithmetic expressions over real-valued parameters.Transition probabilities are thus functions over a set of parameters.A simple instance is to use intervals over system parameters imposing constant lower and upper bounds on every parameter [12], [13].The general setting as considered here is more liberal as it e.g., includes the possibility to express complex parameter dependencies.This paper considers the analysis of parametric Markov models where probability distributions are functions over system parameters.Specifically, we consider parametric discrete-time Markov chains (pMCs) and parametric discrete-time Markov decision processes (pMDPs).
Example 1.The Knuth-Yao randomised algorithm [14] uses repeated coin flips to model a six-sided die.It uses a fair coin In gray states an unfair coin is flipped with probability 2 /5 for 'heads'; for the unfair coin in the white states this probability equals 7 /10.On the right, the two biased coins have parametric probabilities.
to obtain each possible outcome ('one', 'two', ..., 'six') with probability 1 /6.Fig. 1(a) depicts a Markov chain (MC) of a variant in which two unfair coins are flipped in an alternating fashion.Flipping the unfair coins yields heads with probability 2 /5 (gray states) or 7 /10 (white states), respectively.Accordingly, the probability of tails is 3 /5 and 3 /10, respectively.The event of throwing a 'two' corresponds to reaching the state in the MC.Assume now a specification that requires the probability to obtain 'two' to be larger than 3 /20.Knuth-Yao 's original algorithm accepts this specification as using a fair coin results in 1 /6 as probability to end up in .The biased model however, does not satisfy the specification; in fact, a 'two' is reached with probability 1 /10.

C. Probabilistic model checking
The analysis algorithms presented in this paper are strongly related to (and presented as) techniques from probabilistic model checking.Model checking [15], [16] is a popular approach to verify the correctness of a system by systematically evaluating all possible system runs.It either certifies the absence of undesirable (dangerous) behaviour or delivers a system run witnessing a violating system behaviour.Traditional model checking typically takes two inputs: a finite transition system modelling the system at hand and a temporal logic formula specifying a system requirement.Model checking then amounts to checking whether the transition system satisfies the logical specification.LTL and CTL (linear temporal and computation tree logic, respectively) are popular logics for model checking.Model checking is nowadays a successful analysis technique adopted by mainstream hardware and software industry [17], [18].
To cope with real-world systems exhibiting random behaviour, model checking has been extended to deal with probabilistic, typically Markov, models.Probabilistic model checking [15], [19], [20] takes as input a Markov model of the system at hand together with a quantitative specification specified in some probabilistic extension of LTL or CTL.Example specifications are e.g., "is the probability to reach some bad (or degraded) state below a safety threshold λ?" or "is the expected time until the system recovers from a fault bounded by some threshold κ".Popular logics are extensions of CTL with discrete probabilities [21] and, additionally, real-time constraints [22], [23].Extensions thereof with rewards [24], [25] have been considered too.Efficient probabilistic modelchecking techniques do exist for models such as discrete-time Markov chains (MCs), Markov decision processes (MDPs), and their continuous-time counterparts [19].Probabilistic model checking extends and complements long-standing analysis techniques for Markov models.It has been adopted in the field of performance analysis to analyse stochastic Petri nets [26], [27], in dependability analysis for analysing architectural system descriptions [28], in reliability engineering for fault tree analysis [29], [30], as well as in security [31], distributed computing [9], and systems biology [32].Unremitting algorithmic improvements employing the use of symbolic techniques to deal with large state spaces have led to powerful and popular software tools realising probabilistic model checking techniques such as PRISM [33] and Storm [34].

D. Problem statements
We now give a more detailed description of the parameter synthesis problems considered in this paper.We start off by establishing the connection between parametric Markov models and concrete ones, i.e., ones in which the probabilities are fixed such as MCs and MDPs.Each parameter in a pMC or pMDP (where p stands for parametric) has a given parameter range.The parameter space of the parametric model is the Cartesian product of these parameter ranges.Instantiating the parameters with a concrete value in the parameter space to the parametric model results in an instantiated model.The parameter space defines all possible parameter instantiations, or equivalently, the instantiated models.A parameter instantiation that yields a Markov model, e.g., results in probability distributions, is called well-defined.In general, a parametric Markov model defines an uncountably infinite family of Markov models, where each family member is obtained by a well-defined instantiation.A region R is a fragment of the parameter space; it is well defined if all instantiations in R are well defined.
We are now in a position to describe the three problems considered in this paper.
1) The verification problem is defined as follows: The verification problem.Given a parametric Markov model D, a well-defined region R, and a specification ϕ, the verification problem is to check whether all instantiations of D within R satisfy ϕ.
Consider the following possible outcomes: • If R only contains instantiations of D satisfying ϕ, then the verification problem evaluates to true and the Markov model D on region R accepts specification ϕ.Whenever D and ϕ are clear from the context, we call R accepting.• If R contains an instantiation of D refuting ϕ, then the problem evaluates to false.If R contains only instantiations of D refuting ϕ, then D on R rejects ϕ.Whenever D and ϕ are clear from the context, we call R rejecting.
• If R contains instantiations satisfying ϕ as well as instantiations satisfying ¬ϕ, then D on R is inconclusive w. r. t. ϕ.In this case, we call R inconsistent.In case the verification problem yields false for ϕ, one can only infer that the region R is not accepting, but not conclude whether R is inconsistent or rejecting.To determine whether R is rejecting, we need to consider the verification problem for the negated specification ¬ϕ.Inconsistent regions for ϕ are inconsistent for ¬ϕ too.
Example 3 (Verification problem).Consider the pMC D, the well-defined region R from Example 2, and the specification ϕ := ¬ϕ that constrains the probability to reach to be at most 3 /20.The verification problem is to determine whether all instantiations of D in R satisfy ϕ .As there is no instantiation within R for which the probability to reach is above 3 /20, the verification problem evaluates to true.Thus, R accepts ϕ .
A (simple) region comprising a large range of parameter values may likely be inconsistent, as it contains both instantiations satisfying ϕ, and some satisfying ¬ϕ.Thus, we generalise the problem to synthesise a partition of the parameter space.
2) The exact synthesis problem is described as follows: The synthesis problem.Given a parametric Markov model D and a specification ϕ, the (parameter) synthesis problem is to partition the parameter space of D into an accepting region R a and a rejecting region R r for ϕ.
The aim is to obtain such a partition in an automated manner.
A complete sub-division of the parameter space into accepting and rejecting regions provides deep insight into the effect of parameter values on the system's behaviour.The exact division typically is described by non-linear functions over the parameters, referred to as solution functions.
Example 4. Consider the pMC D, the region R, and the specification ϕ as in Example 3. The solution function: describes the probability to eventually reach .Given that ϕ imposes a lower bound of 3 /20, we obtain The example illustrates that exact symbolic representations of the accepting and rejecting regions may be complex and hard to compute algorithmically.The primary reason is that the boundaries are described by non-linear functions.A viable alternative therefore is to consider an approximative version of the synthesis problem.3) The approximate synthesis problem: The aim is to use simpler and more tractable representations of regions such as (sets of) rectangles, rather than non-linear functions; we refer to such regions as simple.As such shapes ultimately approximate the exact solution function, simple regions become infinitesimally small when getting close to the border between accepting and rejecting areas.For computational tractability, we are thus interested in approximating a partition of the parameter space in accepting and rejecting regions, where we allow also for a (typically small) part to be covered by possibly inconsistent regions.Practically this means that c % of the entire parameter space is covered by simple regions that are either accepting or rejecting.Altogether this results in the following problem description: The approximate synthesis problem.Given a parametric Markov model, a specification ϕ, and a percentage c, the approximate (parameter) synthesis problem is to partition the parameter space of D into a simple accepting region R a and a simple rejecting region R r for ϕ such that R a ∪ R r cover at least c% of the entire parameter space.
Example 5. Consider the pMC D, the region R, and the specification ϕ as in Example 3. The parameter space in Fig. 2 is partitioned into simple regions (rectangles).The green (dotted) area-the union of a number of smaller rectangular accepting regions-indicates the parameter values for which ϕ is satisfied, whereas the red (hatched) area indicates the set of rejecting regions for ϕ.The white area indicates the unknown regions.The indicated partition covers 95% of the parameter space.The sub-division into accepting and rejecting (simple) regions approximates the solution function f ϕ (p, q) given before.

E. Solution approaches
We now outline our approaches to solve the verification problem and the two synthesis problems.For the sake of convenience, we start with the synthesis problem.

1) Synthesis:
The most straightforward description of the sets R a and R r is of the form: The satisfaction relation (denoted |=) can be concisely described by a set of linear equations over the transition probabilities [15].As in the parametric setting the transition probabilities are no longer fixed, but rather defined over a set of parameters, the equations become non-linear.
Example 6 (Non-linear equations for reachability).Take the MC from Fig. 1(a).To compute the probability of eventually reaching, e.g., state , one introduces a variable p s for each transient state s encoding that probability for s.For state s 0 and variable p s0 , the corresponding linear equation reads: where p s1 and p s2 are the variables for s 1 and s 2 , respectively.
The corresponding equation for the pMC from Fig. 1(b) reads: The multiplication of parameters in the model and equation variables leads to a non-linear equation system.
Thus, we can describe the sets R a and R r colloquially as: R a , R r = {u | u satisfies a set of non-linear constraints}.
We provide further details on these constraint systems in Sect.V.
A practical drawback of the resulting equation system is the substantial number of auxiliary variables p s , one for each state in the pMDP.For pMCs, a viable possibility is to simplify the equations by (variants of) state elimination [35].This procedure successively removes states from the pMC until only a start and final state (representing the reachability objective) remain that are connected by a transition whose label is (a mild variant of) the solution function f ϕ that exactly describes the probability to reach a target state: We recapitulate the state-elimination procedure and present several alternatives in Sect.IV.
2) Verification: The basic approach to the verification problem is depicted in Fig. 3.We use a description of the accepting region as computed via the synthesis procedure above.Then, we combine the description of the accepting region with the region R to be verified, as follows: The existence of a rejecting instance in R is thus of relevance; if such a point does not exist, the region is accepting.Using R a and R r as obtained above, the query "is R ∩ R r = ∅?" can be solved via satisfiability modulo theories (SMT) over non-linear arithmetic, checking the conjunction over the corresponding constraints for unsatisfiability.With the help of SMT solvers over this theory like Z3 [36], MathSAT [37], or SMT-RAT [38], this can be solved in a synthesise description of: accepting region R a , and yes for R a → reject, yes for R r → accept, otherwise → unknown Fig. 3. Verification via exact synthesis fully automated manner.This procedure is complete, and is computationally involved.Details of the procedure are discussed in Sect.V.
Parameter lifting [39] is an alternative, approximative solution to the verification problem.Intuitively, this approach over-approximates R r for a given R, by ignoring parameter dependencies.Region R is accepted if the intersection with the over-approximation of R r is empty.This procedure is sound but may yield false negatives as a rejecting point may lie in the over-approximation but not in R r .Tightening the overapproximation makes the approach complete.A major benefit of parameter lifting (details in Sect.VI and Sect.VII) is that the intersection can be investigated by standard probabilistic model-checking procedures.This applicability of mature tools results-as will be shown in Sect.X-in a practically efficient procedure.
3) Approximate synthesis: Here, the central issue is to obtain representations of R a and R r by simple regions such as linear or rectangular regions.Our approach for this parameter space partitioning therefore iteratively obtains partial partitions of the parameter space.The main idea is to compute a sequence R i a i of simple accepting regions that successively extend each other.Similarly, an increasing sequence R i r i of simple rejecting regions is computed.At the i-th iteration, R i a ∪ R i r is the covered fragment of the parameter space.The iterative approach halts when this fragment forms at least c % of the entire parameter space.Termination is guaranteed (under some mild conditions on the order of processing regions) as in the limit a solution to the exact synthesis problem is obtained as The typical approach is to let R i+1 a be the union of R i a , the approximations in the previous iteration, together with some accepting region with a simple representation.Rejecting regions are handled analogously.
Fig. 4 outlines a procedure to address the approximate synthesis problem.As part of our synthesis method, we algorithmically guess a (candidate) region R and guess whether it is accepting or rejecting.We then exploit one of our verification methods to verify whether R is indeed accepting (or rejecting).If it is not accepting (rejecting), we exploit this information together with any additional information obtained during verification to refine the candidate region.This process is repeated until an accepting or rejecting region results.We discuss the method and essential improvements in Sect.VIII.
Example 7. Consider the pMC D and the specification ϕ as in Example 2. The parameter space in Fig. 2 is partitioned into regions.The green (dotted) area-the union of a number of smaller rectangular accepting regions-indicates the parameter values for which ϕ is satisfied, whereas the red (hatched) area indicates the set of rejecting regions for ϕ.Checking whether a region is accepting, rejecting, or inconsistent is done by verification.The small white area consists of regions that are unknown (i.e. , not yet considered) or inconsistent.

F. Overview of the paper
Sect.II introduces the required formalisms and concepts.Sect.III defines the notions of regions and formalises the three problems.The section also shows how to ensure welldefinedness and graph-preservation, two important prerequisites for the verification procedures.The section ends with a bird's eye view of the verification approaches that are later discussed in detail.Sect.IV shows how to do exact synthesis by computing the solution function.Sections V-VII present algorithms for the verification problem.Sect.VIII details the approach to reduce the synthesis problem to a series of verification problems.Sections IX and X contain information about the implementation of the approaches, as well as an extensive experimental evaluation.Sect.XI contains a discussion of the approaches and related work.Sect.XII concludes with an outlook.

G. Contributions of this paper
The paper is loosely based on the conference papers [40] and [39] and extends these works in the following way.It gives a uniform treatment of the solution techniques to the synthesis problem, and treats all techniques uniformly for all different objectives-bounded and unbounded reachability as well as expected reward specifications.The material on SMT-based region verification has been extended in the following way: The paper gives the complete characterisations of the SMT encoding with or without solution function.Furthermore, it is the first to extend this encoding to MDPs under angelic and demonic nondeterminism and includes an explicit and in-depth discussion on exact region checking via SMT checkers.It presents a uniform treatment of the linear equation system for Markov chains and its relation to state elimination and Gaussian elimination.It presents a novel and simplified description of state elimination for expected rewards, and a version of state elimination that is targeted towards MTBDDs.The paper contains a correctness proof of approximate verification for a wider range of pMDPs and contains proofs for expected rewards.It also supports expected-time properties for parametric continuoustime MDPs (via the embedded pMDP).Novel heuristics have been developed to improve the iterative synthesis loop.All presented techniques, models, and specifications are realised in the state-of-the-art tool PROPhESY2 .

A. Basic notations
We denote the set of real numbers by R, the rational numbers by Q, and the natural numbers including 0 by N. Let [0, 1] ⊆ R denote the closed interval of all real numbers between 0 and 1, including the bounds; (0, 1) ⊆ R denotes the open interval of all real numbers between 0 and 1 excluding 0 and 1.
Let X, Y denote arbitrary sets.If X ∩ Y = ∅, we write X Y for the disjoint union of the sets X and Y .We denote the power set of X by 2 X = {X | X ⊆ X}.Let X be a finite or countably infinite set.A probability distribution over X is a function µ :

B. Polynomials, rational functions
Let V denote a finite set of parameters over R and dom(p) ⊆ R denote the domain of parameter p ∈ V .Definition 1 (Polynomial, rational function).For a finite set Let Mon[V ] denote the set of monomials over V .A polynomial g (over V ) with t terms is a weighted sum of monomials: Instantiations replace parameters by constant values in polynomials or rational functions.
Definition 2 (Parameter instantiations).A (parameter) instantiation u of parameters V is a function u : V → R.
We abbreviate the parameter instantiation u with u(p i ) = a i ∈ R by the n-dimensional vector (a 1 , . . ., a n ) ∈ R n for ordered parameters p 1 , . . ., p n .Applying the instantiation u on V to polynomial g ∈ Q[V ] yields g[u] which is obtained by replacing each p ∈ V in g by u(p), with subsequent application of + and •.For rational function

C. Probabilistic models
Let us now introduce the probabilistic models used in this paper.We first define parametric Markov models and present conditions such that their instantiations result in Markov models with constant probabilities.Then, we discuss how to resolve non-determinism in decision processes.
1) Parametric Markov models: The transitions in parametric Markov models are equipped with rational functions over the set of parameters.Although this is the general setting, for some of our algorithmic techniques we will restrict ourselves to linear polynomials 3 .We consider parametric MCs and MDPs as subclasses of a parametric version of classical two-player stochastic games [43].The state space of such games is partitioned into two parts, S and S .At each state, a player chooses an action upon which the successor state is determined according to the (parametric) probabilities.Choices in S and S are made by player and , respectively.pMDPs and pMCs are parametric stochastic one-and zero-player games respectively.Definition 3 (Parametric models).A parametric stochastic game (pSG) is a tuple G = (S, V , s I , Act, P) with a finite set S of states with S = S S , a finite set V of parameters over R, an initial state s I ∈ S, a finite set Act of actions, and a transition function P : S × Act × S → Q(V ) ∪ R ∪ {⊥} with |Act(s)| ≥ 1 for all s ∈ S, where Act(s) = {α ∈ Act | ∃s ∈ S. P(s, α, s ) ≡ 0} is the set of enabled actions at state s.
A parametric state-action reward function rew : S × Act → Q(V ) ∪ R ∪ {⊥} associates rewards with state-action pairs.It is assumed that deadlock states are absent, i.e., Act(s) = ∅ for all s ∈ S. Entries in R ∪ {⊥} in the co-domains of the functions P and rew ensure that the model is closed under instantiations, see Def. 5 below.We silently assume that the input pSGs do not contain symbols from R \ Q or ⊥.A model is called parameter-free if all its transition probabilities are constant.
A pSG intuitively works as follows.In state s ∈ S , player non-deterministically selects an action α ∈ Act(s).With (parametric) probability P(s, α, s ) the play then evolves to state s .On leaving state s via action α, the reward rew(s, α) 3 /5 2 /5 3 /5 is earned.If s ∈ S , the choice is made by player , and as for player , the next state is determined in a probabilistic way.
As by assumption no deadlock states occur, this game goes on forever.A pMDP is a game with one player, whereas a pMC has no players; a pMC thus evolves in a fully probabilistic way.Let D denote a pMC, M a pMDP, and G a pSG.
Example 8. Fig. 5(a)-(c) depict a pSG, a pMDP, and a pMC respectively over parameters V = {p, q}.The states of the players and are drawn as circles and rectangles, respectively.The initial state is indicated by an incoming arrow without source.We omit actions in state s if |Act(s)| = 1.In state s 0 of Fig. 5(a), player can select either action α or β.On selecting α, the game moves to state s 1 with probability p, and to s 2 with probability 1−p.In state s 2 , player can select α or β; in s 1 there is a single choice only.
A transition (s, α, s ) exists if P(s, α, s ) ≡ 0. As pMCs have a single enabled action at each state, we omit this action and just write P (s, s ) for P(s, α, s ) if Act(s) = {α}.A state s is a successor of s, denoted s ∈ succ(s), if P(s, α, s ) ≡ 0 for some α; in this case, s ∈ pred(s ) is a predecessor of s .Given two pSGs G = (S, V , s I , Act, P) and G = (S , V , s I , Act , P ), G is a sub-pSG of G if S ⊆ S, V ⊆ V , s I = s I ∈ S , Act ⊆ Act, and P (s, α, s ) ∈ {P(s, α, s ), 0} for all s, s ∈ S and α ∈ Act .Note that for a given state s ∈ S and action α ∈ Act(s), the sub-pSG might not contain s or α might not be enabled in s, but it is also possible that the sub-pSG omits some but not all successors of α in s.
Remark 1. Parametric stochastic games are the most general model used in this paper.They subsume pMDPs and pMCs.For the sake of readability, we introduce the formal foundations for pSGs and indicate how these apply to subclasses.Several algorithmic approaches later on in this paper are not directly applicable to pSGs, but tailored to either pMDPs or pMCs.This is indicated when introducing these techniques.
A state-action reward function rew : S × Act → R ≥0 associates (non-negative, finite) rewards to outgoing actions.Analogously, Markov chains (MCs) and Markov decision processes (MDPs) are defined as special cases of pMCs and pMDPs, respectively.We use D to denote a MC, M for an MDP and G for an SG.
2) Paths and reachability: An infinite path of a pSG G is an infinite sequence π = s 0 α 0 s 1 α 1 . . . of states s i ∈ S and actions α i ∈ Act(s i ) with P(s i , α i , s i+1 ) ≡ 0 for i ≥ 0. A finite path of a pSG G is a non-empty finite prefix s 0 α 0 . . .s n of an infinite path s 0 α 0 . . .s n α n . . . of G for some n ∈ N. Let Paths G denote the set of all finite or infinite paths of G while Paths G fin ⊆ Paths G denotes the set of all finite paths.For paths in (p)MCs, we omit the actions.The set Paths G (s) contains all paths that start in state s ∈ S. For a finite path π ∈ Paths G fin , last(π) = s n denotes the last state of π.The length |π| of a path π is |π| = n for π ∈ Paths G fin and |π| = ∞ for infinite paths.The accumulated reward along the finite path s 0 α 0 . . .α n−1 s n is given by the sum of the rewards rew(s i , α i ) for 0 ≤ i < n.
The probability measure Pr s over sets of paths can be defined using a standard cylinder construction Pr s (s 0 α 0 . . .s n ) = Π n−1 i=0 P(s i , α i , s i+1 ); details can be found in [15,Ch. 10].
A set of states T ⊆ S is reachable from s ∈ S, written s ∈ ♦T , iff there is a path from s to some s ∈ T .A state s is absorbing iff P(s, α, s) = 1 for all α ∈ Act(s).
3) Model instantiation: Instantiated parametric models are obtained by instantiating the rational functions in all transitions as in Def.Remark 2. The instantiation of a pSG at u is a pSG, but not necessarily an SG.This is due to the fact that an instantiation does not ensure that P(s, α, •) is a probability distribution.In fact, instantiation yields a transition function of the form P : S × Act × S → R ∪ {⊥}.Similarly, there is no guarantee that the rewards rew[u] are non-negative.Therefore, we impose restrictions on the parameter instantiations.Definition 6 (Well-defined instantiation).An instantiation u is well-defined for a pSG G if the pSG G[u] is an SG.
The reward function rew is well-defined if it does only associate non-negative reals to state-action pairs.Example 10.Consider again the pMC in Fig. 5(c).The instantiation u with u(p) = 4 /5 and u(q) = 3 /5 is well-defined and induces the MC D[u] depicted in Fig. 5(d).
From now on, we silently assume that every pSG we consider has at least one well-defined instantiation.This condition can be assured through checking the satisfiability of the conditions in Def. 4, which we discuss in Sect.III-D.
Our methods necessitate instantiations that are not only welldefined, but also preserve the topology of the pSG.
Definition 7 (Graph-preserving).A well-defined instantiation u for pSG G = (S, V , s I , Act, P) is graph-preserving if for all s, s ∈ S and α ∈ Act, Example 11.The well-defined instantiation u with u(p) = 1 and u(q) = 3 /5 for the pMC in Fig. 5(c) is not graph-preserving.

4)
Resolving non-determinism: Strategies 4 resolve the nondeterministic choices in stochastic games with at least one player.For the objectives considered here, it suffices to consider so-called deterministic strategies [44]; more general strategies can be found in [15,Ch. 10].We define strategies for pSGs and assume well-defined instantiations as in Def. 6. Definition 8 (Strategy).A (deterministic) strategy σ i for player i ∈ { , } in a pSG G with state space S = S S is a function Let Str G denote the set of strategies σ = (σ , σ ) for pSG G and Str G i the set of strategies of player i.
A pMDP has only a player-i strategy for the player with S i = ∅; in this case the index i is omitted.A player-i strategy σ i is memoryless if last(π) = last(π ) implies σ i (π) = σ i (π ) for all finite paths π, π .A memoryless strategy can thus be written in the form σ i : S i → Act.A pSG-strategy σ = (σ , σ ) is memoryless if both σ and σ are memoryless.
Remark 3. From now on, we only consider memoryless strategies and refer to them as strategies. 4Also referred to as policies, adversaries, or schedulers.
The notions of strategies for pSGs and pMDPs and of induced pMCs naturally carry over to non-parametric models; e.g., the MC G σ is induced by strategy σ ∈ Str G on SG G.

D. Specifications and solution functions
1) Specifications: Specifications constrain the measures of interest for (parametric) probabilistic models.Before considering parameters, let us first consider MCs.Let D = (S, s I , P) be an MC and T ⊆ S a set of target states that (without loss of generality) are assumed to be absorbing.Let ♦T denote the path property to reach T .We overload this notation to also denote the set of states that have a positive probability to reach the target states: ♦T = {s ∈ S | ∃π ∈ Paths D fin (s).last(π) ∈ T }.We consider three kinds of specifications: 1) Unbounded probabilistic reachability The specification ϕ r = P ≤λ (♦ T ) asserts that the probability to reach T from the initial state s I shall be at most λ, where λ ∈ Q ∩ [0, 1].Formally, specification ϕ r is satisfied by MC D, written: where Pr D sI (♦ T ) is the probability mass of all infinite paths that start in s I and visit any state from T .2) Bounded probabilistic reachability In addition to reachability, these specifications impose a bound on the maximal number of steps until reaching a target state.Specification ϕ b = P ≤λ (♦ ≤n T ) asserts that in addition to P ≤λ (♦ T ), states in T should be reached within n ∈ N steps.The satisfaction of P ∼λ (♦ ≤n T ) is defined similar as above.3) Expected reward until a target The specification ϕ e = E ≤κ (♦ T ) asserts that the expected reward until reaching a state in T shall be at most κ ∈ R. Formally, let ER D sI (♦ T ) denote the expected accumulated reward until reaching a state in T ⊆ S from state s I ; if Pr D sI (♦ T ) < 1 then we set ER D sI (♦ T ) := ∞ [15]5 .Then we define We do not treat the accumulated reward to reach a target within n steps, as this is not a very useful measure.In case there is a possibility to not reach the target within n steps, this yields ∞.We omit the superscript D if it is clear from the context.We write ¬ϕ to invert the relation: D |= ¬P ≤λ (♦ T ) is thus equivalent to D |= P >λ (♦ T ).An SG G satisfies specification ϕ under strategy σ if the induced MC G σ |= ϕ.Unbounded reachability and expected rewards are prominent examples of indefinite-horizon properties -they measure behaviour up-to some specified event (the horizon) which may be reached after arbitrarily many steps.
Remark 4. Bounded reachability in MDPs can be reduced to unbounded reachability by a technique commonly referred to as unrolling [45].For performance reasons, it is sometimes better to avoid this unrolling, and present dedicated approaches.
2) Solution functions: Computing (unbounded) reachability probabilities and expected rewards for MCs reduces to solving linear equation systems [15] over the field of reals (or rationals).For parametric MCs, a linear equation system over the field of the rational functions over V results.The solution to this equation system is a rational function.(See Examples 4 and 6 in the introduction.)More details on the the solution function and the equation system follow in Sect.IV and Sect.V, respectively.Definition 10 (Solution functions).For a pMC D = (S, V , s I , P), T ⊆ S and n ∈ N, a solution function for a specification ϕ is a rational function , such that for every well-defined graph-preserving instantiation u: Example 13.Consider the reachability probability to reach s 2 for the pMC in Fig. 6(a).Any instantiation u with u(p), u(q) ∈ (0, 1) is well-defined and graph-preserving.As the only two finite paths to reach s 2 are s 0 s 2 and s 0 s 1 s 2 , we have f r D,{s2} = 1 − p + p•q.
For pSGs (and pMDPs), the solution function depends on the resolution of non-determinism by strategies, i. e., they are defined on the induced pMCs.Formally, a solution function for a pSG G, a reachability specification ϕ r = P ≤λ (♦ T ), and a strategy σ ∈ Str G is a function f r G,σ,T ∈ Q(V ) such that for each well-defined graph-preserving instantiations u it holds: These notions are defined analogously for bounded reachability (denoted f b G,σ,T ) and expected reward (denoted f e G,σ,T ) specifications.

E. Constraints and formulas
We consider (polynomial) constraints of the form g ∼ g with g, g ∈ Q[V ] and ∼∈ {<, ≤, =, ≥, >}.We denote the set of all constraints over V with C[V ].A constraint g ∼ g can be equivalently formulated as g − g ∼ 0. A formula ψ over a set of polynomial constraints is recursively defined: Each polynomial constraint is a formula, and the Boolean combination of formulae is also a formula.
Example 15.Let p, q be variables.
The semantics are standard: i.e., an instantiation u satisfies . An instantiation satisfies ψ ∧ ψ if u satisfies both ψ and ψ .The semantics for other Boolean connectives are defined analogously.Moreover, we will write g = g to denote the formula g < g ∨ g > g .
Checking whether there exists an instantiation that satisfies a formula is equivalent to checking membership of the existential theory of the reals [46].Such a check can be automated using SMT-solvers capable of handling quantifier-free non-linear arithmetic over the reals [36], such as [38], [47].
Statements of the form f ∼ f with f, f ∈ Q(V ) are not necessarily polynomial constraints: however, later we are not interested in instantiations u with f [u] = ⊥, and thus later (in Sect.III-D2) we can transform such constraints into formulae over polynomial constraints.

III. REGION VERIFICATION
This section defines the notion of regions and formalises the verification and synthesis problems.It also shows how to obtain graph-preserving instantiations.Finally, it surveys the verification approaches that are detailed later in the paper.
Instantiated models are amenable to standard probabilistic model checking.However, this sampling is very restrictiveverifying an instantiated model gives results for a single point in the parameter space.A more interesting problem is to determine which parts of the parameter space give rise to a model that complies with the specification.Such sets of parameter values are, inspired by their geometric interpretation, called regions.
We start off by introducing a general satisfaction relation for parametric Markov models for a single given instantiation.We then introduce regions and lift these notions to regions.

Definition 11 (Angelic and demonic satisfaction relations).
For pSG G, well-defined instantiation u, and specification ϕ, the satisfaction relations |= a and |= d are defined by: The angelic relation |= a refers to the existence of a strategy to fulfil the specification ϕ, whereas the demonic counterpart |= d requires all strategies to fulfil ϕ.Observe that G, u |= a ϕ if and only if G, u |= d ¬ϕ.Thus, demonic and angelic can be considered to be dual.For pMCs, the relations |= a and |= d coincide and the subscripts a and d are omitted.
Example 16.Consider the pMDP M in Fig. 6

A. Regions
Instead of considering a single instantiated model, we identify sets of instantiated models by regions, which are solution sets of conjunctions of constraints over V .

Definition 12 (Region
Any region which is a subset of a region R is called a subregion of R.
Example 17.Let the region R over V = {p, q} be described by Regions do not have to describe a contiguous area of the parameter space; e.g., consider the region R described by Regions are semi-algebraic sets [46] which yield the theoretical formalisation of notions such as distance, convexity, etc.It also ensures that regions are well-behaved: Informally, a region in the space R n is given by a finite number of connected cells, and (the boundaries of) each connected cell can be described by a finite set of polynomials.The size R of a region R is given by the Lebesgue measure.All regions are Lebesgue measurable.Two classes of regions are relevant in the current context: linear and rectangular.

Definition 13 (Linear region). A region with representation
Linear regions describe convex polytopes.We refer to the vertices (or angular points) of the polytope as the region vertices.

Definition 14 (Rectangular region). A region R with representation
This definition ensures that all instantiations from graphpreserving regions are well-defined and that the instantiated models have the same topology as the parametric model, cf. 5 below.
Our aim is to consider specifications ϕ that hold for all instantiations represented by a region R of a parametric model G.This is captured by the following satisfaction relation.

Definition 16. (Satisfaction relation for regions)
For pSG G, well-defined region R, and specification ϕ, the relation |= ♣ , ♣ ∈ {a, d}, is defined as: • M, R |= a ϕ, as for strategy σ β = {s 0 → β}, we have Regions can be inconsistent w. r. t. a relation, and consistent w. r. t. its dual relation.The region (0, 1) × (0, 1) is inconsistent for M and |= d , as for both ϕ and ¬ϕ, there is a strategy that is not accepting.For |= a , there is a single strategy which accepts ϕ; other strategies do not affect the relation.

Remark 5. Graph-preserving regions (Def. 15) have the nice property that either
This property can be checked by standard graph analysis [15,Ch. 10].It is thus straightforward to check G, R |= ♣ P =1 (♦T ), an important precondition for computing expected rewards.For the remainder, we assume that for expected rewards, within a region the probability to reach a target is one.
In the remainder of the paper, we repeatedly (and often implicitly) use the following properties for regions.

Lemma 1 (Characterisation for inconsistent regions). For any inconsistent region
The statements follow from the universal quantification over all instantiations in the definition of |= ♣ .Remark 6.Another notion in parameter synthesis is the existence of a robust strategy, that is, The relation differs from G, R |= ϕ in the quantifier order, that is, G, R |= ϕ considers potentially different strategies for different parameter instantiations u ∈ R. The notion of robust strategies leads to a series of quite orthogonal challenges.For instance, the notion is not compositional, that is, if in R 1 and R 2 robust strategies exist, then we cannot conclude the existence of a robust strategy in R 1 ∪R 2 .Moreover, memoryless strategies are not sufficient, see [48].Other than in the Sect.VII, robust strategies are not considered here.

B. Formal problem statements
We are now in a position to formalise the two synthesis problems and the verification problem, see page 3.
The formal synthesis problem.Given pSG G, specification ϕ, and well-defined region R, the synthesis problem is to partition R into R a and R r such that: This problem is the topic of Sect.IV.Remark 7. The solution function for pMCs precisely describes how (graph-preserving) instantiations map to the relevant measure.Therefore, comparing the solution function with the threshold divides the parameter space into an accepting region R a and a rejecting region R r and defines the exact result for the formal synthesis problem.Recall therefore also Ex. 4 on page 3.
The formal verification problem.Given pSG G, specification ϕ, and well-defined region R, the verification problem is to check whether: where |= ♥ denotes the dual satisfaction relation of |= ♣ .This problem is the topic of Sect.V-VII.
The verification procedure allows us to utilise a approximate synthesis problem in which verification procedures are used as a backend.
The formal approximate synthesis problem.Given pSG G, specification ϕ, percentage c, and well-defined region R, the approximate synthesis problem is to partition R into regions R a , R u , and R r such that: where R a ∪ R r cover at least c% of the region R.The regions R a , R u and R r should be finite unions of rectangular regions.This problem is the topic of Sect.VIII.
No requirements are imposed on the (unknown) region R u .

C. A bird's eye view on the verification procedures
In the later sections, we will present several techniques that decide the verification problem for pMCs and pMDPs.(Recall that stochastic games were only used to define the general setting.) The verification problem is used to analyse regions-ofinterest.The assumption that this region contains only welldefined instantiations is therefore natural.It can be checked algorithmically as described in Sect.III-D below.Many verification procedures require that the region is graph-preserving.A decomposition result of well-defined into graph-preserving regions is given in Sect.III-E.
Sect.V presents two verification procedures.The first one directly solves the non-linear equation system, see Example 6 on page 4, as an SMT query.The second procedure reformulates the SMT query using the solution function.While this reformulation drastically reduces the number of variables in the query, it requires an efficient computation of the solution function, as described in Sect.IV.
Sect.VI covers an approximate and more efficient verification procedure, called parameter lifting, which is tailored to multi-linear functions and closed rectangular regions.Under these mild restrictions, the verification problem for pMCs (pMDPs) can be approximated using a sequence of standard verification analyses on non-parametric MDPs (SGs) of similar size, respectively.The key steps here are to relax the parameter dependencies, and consider lower-and upper-bounds of parameters as worst and best cases.

D. Checking whether a region is graph preserving
The verification problem for region R requires R to be well-defined.We first address the problem on how to check this condition.In fact, we present a procedure to check graph preservation which is slightly more general and useful later, see also Remark 5. To falsify that region R is graph preserving, we search for points in R violating the conditions in Def. 7. Using the representation of R, the implication Φ(R) =⇒ graph-preserving needs to be valid since any violating assignment corresponds to a non-graph-preserving instantiation inside R. Technically, we consider satisfiability of the conjunction of: • the inequalities C(R) representing the candidate region, and • a disjunction of (in)equalities describing the violation of the graph-preserving property.This conjunction is satisfiable iff the region is not graph preserving.
1) An equation system for graph preservation: The following equation system is only valid for pSGs with polynomial transition probabilities.We discuss the creation of an equation system for pSGs with rational functions as transition probabilities at the end of the section.The following constraints (1)-( 4) capture the notion of graph preservation: The constraints ensure that (1) all non-zero entries are evaluated to a probability, (2) transition probabilities are probability distributions, (3) rewards are non-negative, and (4) non-zero entries remain non-zero.They can be simplified to: rew(s, α) ≥ 0.
Satisfiability of ( 1)-( 4), or equivalently, deciding whether a region is graph preserving, is as hard as the existential theory of the reals [46], if no assumptions are made about the transition probability and reward functions.This checking can be automated using SMT-solvers capable of handling quantifierfree non-linear arithmetic over the reals [36].The complexity drops to polynomial time once both the region R and all transition probability (and reward) functions are linear and the sums of outgoing transitions always (syntactically) sum to 16 : as linear programming has a polynomial complexity and the formula is then a disjunction over linear programs (with trivial optimisation functions).
2) Handling rational functions: In case the transition probability and reward function are not polynomials, the lefthand side of the statements in ( 1)-( 4) would not be polynomials, and the statements would not be constraints.We therefore perform the following transformations on the statements in (1)-( 4): • Transformation of equalities: with c ∈ Q, and equals < for and ≤ for ≥.
• Transformation of g = g (i.e. of the formula g < g ∨g > g ) is follows the application on both inequalities.
The result is a formula with polynomial constraints that correctly describes graph preservation (or well-definedness).
Example 21.Consider a state with outgoing transition probabilities q and p 1+p .The preservation statements are (after some simplification): Transforming the second item as explained before results in: while transforming the third item yields: Finally, we obtain the following formula (after some further simplifications):

E. Reduction to graph-preserving regions
Our methods only allow regions are graph-preserving.If the region R is well-defined, but not graph-preserving, weas a preprocessing-split the region into subregions.Let us illustrate this.Let us formalise the construction from this example.For a given well-defined region R, and model G, let Z R describe the set of constraints: For X ⊆ Z R , the subregion R X ⊆ R is defined as: It follows that X uniquely characterises which transition probabilities in G are set to zero.In fact, each instance in R X is graph-preserving for the unique sub-pSG G of G obtained from G by removing all zero-transitions in R X .The pSG G is well-defined as R on G is well-defined.By construction, it holds that

IV. EXACT SYNTHESIS BY COMPUTATION OF THE SOLUTION FUNCTION
The solution function for pMCs describes the exact accepting and rejecting regions, as discussed in Sect.III-B 7 .In Sect.V, we will also see that the solution function may be beneficial for the performance of SMT-based (region) verification.
This section discusses how to actually compute the solution function.It starts with some essential observations before recapping the original state elimination approach, albeit slightly rephrased.In the last part, we present alternative, equivalent formulations which sometimes allow for superior performance.

A. Observation
The original approach to compute the solution function of pMCs is via state elimination [35], [49], and is analogous to the computation of regular expressions from nondeterministic finite automata (NFAs) [50].It is suitable for a range of indefinitehorizon properties.The core idea behind state elimination and the related approaches presented here is based on two operations: • Adding short-cuts: Consider the pMC-fragment in Fig. 8(a).The reachability probabilities from any state to t are as in Fig. 8(b), where we replaced the transition from s to s by shortcuts from s to t and all other successors of s , bypassing s .By successive application of shortcuts, any path from the initial state to the target state eventually has length 1.
• Elimination of self-loops: A prerequisite for introducing a short-cut is that the bypassed state is loop-free.Recall that the probability of staying forever in a non-absorbing state is zero, and justifies elimination of self-loops by rescaling all other outgoing transitions, as depicted in the transition from Fig. 8(c) to Fig. 8(d).

B. State elimination
Let T ⊆ S be a set of target states and assume w. l. o. g. that all states in T are absorbing and that s I ∈ T .
1) Reachability probabilities: We describe the algorithm to compute reachability probabilities based on state-elimination in Alg. 1.In the following, P is the transition matrix.The function eliminate selfloop(P, s) rescales all outgoing probabilities of a non-absorbing state s by eliminating its self-loop.The function eliminate transition(P, s 1 , s 2 ) adds a shortcut from s 1 to the successors of s 2 .Both operations preserve reachability to T .The function eliminate state(P, s) "bypasses" a state s by adding shortcuts from all its predecessors.More precisely, we eliminate the incoming transitions of s, and after all incoming transitions are removed, the state s is unreachable.It is thereby effectively removed from the model.
After removing all non-absorbing, non-initial states S ?, the remaining model contains only self-loops at the absorbing states and transitions emerging from the initial state.Eliminating the self-loop on the initial state (by rescaling) yields a pMC.In this pMC, after a single step, an absorbing state is reached.The absorbing states are either a target state or not.The solution function is then the sum over all transition probabilities to target states in T .
Example 23.Consider again the pMC from Example 8, also depicted in Fig. 9(a).Assume state s 2 is to be eliminated.Applying the function eliminate state(P, s 2 ), we first eliminate the transition s 1 → s 2 , which yields Fig. 9(b), and subsequently eliminate the transition s 0 → s 2 (Fig. 9(c)).State s 2 is now unreachable, so we can eliminate s 2 , reducing computational effort when eliminating state s 1 .For state s 1 , we first eliminate the self-loop (Fig. 9(e)) and then eliminate the transition s 0 → s 1 .The final result, after additionally removing the now unreachable s 1 , is depicted in Fig. 9(f).The result, i.e., the probability to eventually reach s 3 from s 0 in the original model, can now be read from the single transition between these two states.
As for computing of regular expressions from NFAs, the order in which the states are eliminated is essential.Computing an optimal order with respect to minimality of the result, however, is already NP-hard for acyclic NFAs, see [51].For state elimination on pMCs, the analysis is more intricate, as the cost of every operation crucially depends on the size and the structure of the rational functions.We briefly discuss the implemented heuristics in Sect.IX-B1.
Remark 8.The elimination of self-loops yields a rational function.In order to keep these functions as small as possible, it is natural to eliminate common factors of the numerator and the denominator.Such a reduction, however, involves the computation of greatest common divisors (gcds).This operation is expensive for multivariate polynomials.In [52], data structures to avoid their computation are introduced, in [53] a method is presented that mostly avoids introducing common factors.
2) Expected rewards: The state elimination approach can also be adapted to compute expected rewards [49].When eliminating a state s, in addition to adjusting the probabilities of the transitions from all predecessors s 1 of s to all successors s 2 of s, it is also necessary to "summarise" the reward that would have been gained from s 1 to s 2 via s.The presentation in [49] describes these operations on so-called transition rewards.Observe that for the analysis of expected rewards in MCs, we can always reformulate transition rewards in terms of state rewards.We preprocess pMCs to only have rewards at the states: this adjustment simplifies the necessary operations considerably.
The treatment of the expected reward computation is easiest from an adapted (and more performant) implementation of state elimination, as outlined in Alg. 2. Here, we eliminate the probabilities to reach a target state in exactly one step, and collect these probabilities in a vector x which we refer to as one-step-probabilities.Then, we proceed similar as before.However, the elimination of a transition from s 1 to s now has two effects: it updates the probabilities within the non-target states as before, and (potentially) updates the probability x(s 1 ) to reach the target within one step from s 1 (with the probability x(s 1 ) := x(s 1 ) + P(s1,s)•x(s) for each s 2 ∈ succ(s), s = s 2 do P(s 1 , s 2 ) := P(s 1 , s 2 ) + P (s1,s)•P(s,s2) 1−P(s,s) P(s 1 , s) := 0 eliminate state(P, x, s ∈ S) assert P(s, s) = 0 for each s 1 ∈ pred(s) do eliminate transition(P, x, s 1 , s) x(s) := t∈T P(s, t) for each s ∈ S ?P(s, t) := 0 for all s ∈ S, t ∈ T while S ?= ∅ do eliminate state(P, x, s) for some s ∈ S ?S ?:= S ?\ {s} return x(s I ) that the target was reached via s in two steps).Upon termination of the outer loop, the vector x contains the probabilities from all states to reach the target, that is, x(s i ) = x si .
Finally, when considering rewards, the one-step-probabilities contain initially the rewards for the states.Eliminating a transition then moves the (expected) reward to the predecessors by the same sequence of arithmetic operations.
3) Bounded reachability: As discussed in Remark 4, bounded reachability can typically be considered by an unfolding of the Markov model and considering an unbounded reachability property on that (acyclic) unfolding.In combination with state-elimination, that yields the creation of many states that are eliminated afterwards, and does not take into account any problem-specific properties.Rather, and analogous to the parameter-free case [15], it is better to do the adequate matrix-vector multiplication (# number of steps often).The matrix originates from the transition matrix, the vector (after i multiplications) encodes the probability to reach a state within i steps.

C. Linear equation system
The following set of equations is a straightforward adaption of the Bellman linear equation system for MCs found in, e.g., [15], [54] to pMCs.For each state s, a variable x s is used to express the probability Pr s (♦T ) to reach a state in T from the state s.Recall that we overloaded ♦T to also denote the set of states from which T is reachable (with positive probability).Analogously, we use ¬♦T to denote the set of states from which T is not reachable, i. e., ¬♦T = S \ ♦T .We have: This system of equations has a unique solution for every well-defined parameter instantiation.In particular, the set of states satisfying ¬♦T is the same for all well-defined graph-preserving parameter instantiations, as instantiations that maintain the graph of the pMC do not affect the reachability of states in T .
For pMCs, the coefficients are no longer from the field of the real numbers, but rather from the field of rational functions.
Example 24.Consider the equations for the pMC from Fig. 9(a).
Bringing the system in normal form yields: Adding q times the second equation to the third equation (concerning state s 2 ) brings the left-hand side matrix in upper triangular form: x 4 = 0.
The equation system yields the same result as the elimination of the transition from s 2 to s 1 (notice the symmetry between s 1 and s 2 ).
The example illustrates that there is no elementary advantage in doing state elimination over resorting to solving the linear equation sytem by (some variant of) Gaussian elimination.If we are only interested in the probability from the initial state, we do not need to solve the full equation system.The stateelimination algorithm, in which we can remove unreachable states, optimises for this observation, in contrast to (standard) linear equation solving.As in state-elimination, the elimination order of the rows has a significant influence.

D. Set-based transition elimination
To succinctly represent large state spaces, Markov chains are often represented by multi-terminal binary decision diagrams (or variants thereof) [55].Such a symbolic representation handles sets of states instead of single states (and thus also sets of transitions), and thereby exploits symmetries and similarities in the underlying graph of a model.To support efficient elimination, we describe how to eliminate sets of transitions at once.The method is similar to the Floyd-Warshall algorithm for all-pair shortest paths [56].The transition matrix contains one-step probabilities for every pair of source and target states.Starting with a self-loop-free pMC (obtained by eliminating all self-loops from the original pMC), we iterate two operations until convergence.By doing a matrix-matrix multiplication, we effectively eliminate all transitions emanating from all nonabsorbing states simultaneously.As this step may reintroduce (a) After first iteration self-loops, we eliminate them in a second step.As before, eventually only direct transitions to absorbing states remain, which effectively yield the unbounded reachability probabilities.The corresponding pseudo-code is given in Alg. 3.
The approach of this algorithm can conveniently be explained in the equation system representation.Let us therefore conduct one step of the algorithm as an example, where we use the observation that the matrix-matrix multiplication corresponds to replacing the variables x s by their defining equations in all other equations.
Example 25.Reconsider the equations from Example 24: Using the equations for x 0 , x 1 , x 2 to replace their occurrences in all other equations yields: We depict the pMC which corresponds to this equation system in Fig. 10(a).Again, notice the similarity to state-elimination.
For completeness, the result after another iteration is given in Fig. 10(b).
The correctness follows from the following argument: After every iteration, the equations describe a pMC over the same state space as before.As all absorbing states have defining Algorithm 3 Set-based transition elimination for pMCs reachability(pMC D = (S, V , s I , P), T ⊆ S) S ?:= {s ∈ S | s = s I ∧ s ∈ ♦T \ T } for each s ∈ S ?do // can be done in parallel for all s eliminate selfloop(P, s) while ∃s, s ∈ S ? .P(s, s ) = 0 do for each s ∈ S ?, s ∈ S do // can be done in parallel for all s, s P (s, s ) := s P(s, s ) • P(s , s ) for each s ∈ S ?do // can be done in parallel for all s eliminate selfloop(P , s) P := P return t∈T P(s I , t) equations x i ∈ {0, 1}, the equation system is known to have a unique solution [15].Moreover, as the equation system in iteration i implies the equation system in iteration i + 1, they preserve the same (unique) solution.
V. SMT-BASED REGION VERIFICATION In this section, we discuss a complete procedure to verify regions.We first introduce a conjunction of constraints for pMCs, and extend the idea towards a formula over polynomial constraints for pMDPs.We discuss how to perform region verification on this formulation by using an SMT-solver over nonlinear arithmetic, and indicate how to reduce the number of variables by precomputing the solution function.Throughout the section, we focus on unbounded reachability, that is, we assume ϕ = P ≤λ (♦T ).As expected rewards can be described by a similar equation system, lifting the concepts is straightforward.Again, we assume a graph-preserving region R.

A. Satisfiability checking for pMC region checking
Recall from Sect.IV-C the equation system for pMCs, exemplified by a running example.
Example 26.Reconsider the pMC D from Fig. 6(a), repeated in Fig. 11(a) for convenience.The equation system for reaching T = {s 2 }, using x i to denote x si , is given by: The conjunction of the equation system for the pMC is an implicitly existential quantified formula to which we refer by Φ(D) (as R is well-defined).By construction, this formula is satisfiable.
Remark 9.If transitions in the pMC are not polynomial but rational functions, the equations are not polynomial constraints, hence their conjunction is not a formula (Sect.II-E).Instead, each x = P(s, s ) has to be transformed by the rules in Sect.III-D2: then, their conjunction is a formula.This transformation can always be applied, in particular, in the equalities we are never interested in the evaluation of instantiations u ∈ R with P(s, s )[u] = ⊥: Recall that we are interested in analysing this equation system on a well-defined parameter region R: Therefore, for any u ∈ R, P(s, s )[u] = ⊥ for each s, s ∈ S. Thus, when Φ(D) is used in conjunction with Φ(R), we do not need to consider this special case.

Satisfiability of the conjunction of:
• the equation system Φ(D), • a comparison of the initial state s I with the threshold λ, and • a formula Φ(R) describing the parameter region R, means that-for some parameter instantiation within the region-the reachability probability from the initial state s I satisfies the bound.Unlike Φ(D), this conjunction may be unsatisfiable.
For the satisfaction relations |= a and |= d as defined in Def.11, we have to certify that all parameter values within a region yield a reachability probability that satisfies the threshold.That means, we have to quantify over all instantiations u, (roughly) leading to a formula of the form ∀u . . .|= ϕ.By negating this statement, we obtain the proof obligation ¬∃u . . .|= ¬ϕ: no parameter value within the region R satisfies the negated comparison with the initial state.Thus, we check the conjunction of: • the equation system Φ(D), • a comparison of the initial state with the threshold, by inverting the given threshold-relation, and • a formula Φ(R) describing the parameter region.This conjunction is formalised in the following definition.
Definition 18 (Equation system formula).Let D be a pMC, ϕ = P ∼λ (♦T ), and R a region.The equation system formula is given by: If this formula is not satisfiable, then D, R |= ϕ.Otherwise, a satisfying solution is a counterexample.
Observe that the number of variables is |S| + |V |, which quickly becomes too large for SMT-solvers dealing with nonlinear real arithmetic.However, many of the variables are auxiliary variables that encode the probability to reach target states from each individual state.We can get rid of these variables by replacing the full equation system by the solution function (Def.10 on page 9).Definition 19 (Solution function formula).Let D be a pMC, ϕ = P ∼λ (♦T ), and R a region.The solution function formula 8  is given by: f r D,T ∼ λ ∧ Φ(R).
Example 29.We consider the same scenario as in Example 27.
The solution function is given in Example 13.The solution function formula is: By construction, the equation system formula and the solution function formula for pMC D and reachability property ϕ are equisatisfiable.

B. Existentially quantified formula for parametric MDPs
We can also utilise an SMT solver to tackle the verification problem on pMDPs.For parametric MDPs, we distinguish between the angelic and the demonic case, cf.Def.16.We use that optimal strategies for unbounded reachability objectives are memoryless and deterministic.
1) Demonic: The satisfaction relation |= d is defined by two universal quantifiers, ∀u∀σ . . .|= ϕ.We therefore try to refute satisfiability of ∃u∃σ . . .|= ¬ϕ.Put in a game-theoretical sense, the same player can choose both the parameter instantiation u and the strategy σ to resolve the non-determinism.We use the 8 Remark 9 applies also here following generalisation of the set of linear equations, where we define a disjunction over all possible nondeterministic choices: α∈Act(s) We denote the conjunction of ( 12)-( 14) as Φ d (M) for pMDP M9 .Instead of a single equation for the probability to reach the target from state s, we get one equation for each action.The solver can now freely choose which (memoryless deterministic) strategy it uses to refute the property.
Definition 20 (Demonic Equation System Formula).Let M be a pMDP, ϕ = P ≤λ (♦T ), and R a region.The demonic equation system formula is given by: Example 30.Let M be the pMDP from Fig. 11(b).Let R, ϕ be as in Example 27.The demonic equation system formula is with Φ(R) as before, and Similarly, when using the (potentially exponential) set of solution functions, we let the solver choose: Definition 21 (Demonic Solution Function Formula).Let M be a pMDP, ϕ = P ∼λ (♦T ), and R a region.The demonic solution function formula is given by: As the set of solution functions can be exponential, the demonic solution function formula can grow exponentially.
Example 31.The demonic solution function formula for M, ϕ, R as in Example 30, is given by: 2) Angelic: The satisfaction relation |= a has two different quantifiers, ∀u∃σ . . .|= ϕ.Again, we equivalently try to refute the satisfiability of ∃u∀σ . . .|= ¬ϕ.The quantifier alternation can be circumvented by lifting the linear programming (LP) formulation for MDPs [54], where for each nondeterministic choice an upper bound on the probability variables x s is obtained: Intuitively, the conjunction in constraint (17) eliminates the freedom of choosing any strategy from the solver and forces it to use the strategy that minimises the reachability probability.This means that the constraint system is only satisfiable if all strategies violate the probability bound.We denote the conjunction of ( 15)-( 17) as Φ a (M).Notice that, as for parameter-free MDPs, the optimisation objective of the LP formulation can be substituted by the given probability bound.
Definition 22 (Angelic Equation System Formula).Let M be a pMDP, ϕ = P ≤λ (♦T ), and R a region.The angelic equation system formula is given by: Example 32.Let M, ϕ, R as in Example 30.The angelic equation system formula is given by When using the set of solution functions, all strategies have to be considered: Definition 23 (Angelic Solution Function Formula).Let M be a pMDP, ϕ = P ≤λ (♦T ), and R a region.The angelic solution function formula is given by: Example 33.The angelic solution function formula for M, ϕ, R as in Example 30 is given by:

VI. MODEL-CHECKING-BASED REGION VERIFICATION
OF PARAMETRIC MARKOV CHAINS In this section, we discuss an alternative approach to the verification problem for a pMC, a region and a specification.We first treat reachability probabilities, and then extend the approach to the treatment of expected rewards.
In a nutshell, the idea presented in this section is to transform a pMC into an MDP whose minimal (maximal) reachability probability under-approximates (over-approximates) the reachability probability of the pMC.In particular, consider the pMC  From this result, we infer that max u∈R Pr D[u] (♦T ) ≤ 47 /60.Details follow below.

A. Observation
For an instantiation u ∈ R, Pr D[u] (♦T ) can be expressed as a rational function f = g1 /g2 with polynomials g 1 , g 2 due to Def. 10.Recall that we assume region R to be graph-preserving.Therefore, g 2 [u] = 0 for all u ∈ R and f is continuous on any closed region R. Hence, there is an instantiation u ∈ R that induces the maximal (or minimal) reachability probability: To infer that R is accepting (i.e.all instantiations u ∈ R induce probabilities at most λ), it suffices to show that the maximal reachability probability over all instantiations is at most λ: One way to determine the maximum reachability probability is to first determine which u ∈ R induces the maximum, and then compute the probability on the instantiated model D However, constructing an oracle that determines the u that induces the maximum is difficult in general.
Example 35.Consider a three-state pMC where the probability from initial state s I to target state t is a non-linear, nonmonotone transition function, as, e.g., the transition probability from s 0 to s 3 of the pMC in Fig. 9(f).Finding the maximum requires an analysis of the derivative of the solution function, and is (approximately) as hard as the exact verification problem.
Therefore, we assume monotonic transition probabilities, and consider a slightly restricted class of pMCs.
Definition 24 (Locally monotone pMCs).A pMC D = (S, V , s I , P) is locally monotone iff for all s ∈ S there is a multilinear polynomial g s ∈ Q This restriction only constrains the way how a model enters the problem -the resulting reachability probabilities may still be represented by more complicated functions.Moreover, locally monotone pMCs include most pMCs from the literature [39], and also include, e.g., the embedded pMCs of parametric continuous-time Markov chains with multilinear exit rates.Examples of the egligible transition probabilities are p, pq, 1 /p and their complements formed by 1 − p etc.
Thanks to monotonicity, for a locally monotone pMC D = (S, V , s I , P), and a closed rectangular region R we have that for all s, s ∈ S : However, the restriction to local monotonicity does not immediately overcome the challenge of constructing an oracle.
Example 36.Reconsider the locally monotone pMC D in Fig. 5(c)-which is also given in Fig. 13(a)-and the closed rectangular region R = [ 1 /10, 4 /5] × [ 2 /5, 7 /10].We make two observations: s 4 is the only state from which we cannot reach s 3 , furthermore, s 4 is only reachable via s 2 .Hence, it is best to avoid s 2 .From state s 0 , it is thus beneficial if the transition probability to s 2 is as small as possible.Equivalently, it is beneficial if p is as large as possible, as this minimises the probability of reaching s 2 and as p does not occur elsewhere.Now we consider state s 1 : As we want to reach s 3 , the value of q should be preferably low.However, q occurs also at transitions leaving s 2 .From s 2 , q should be assigned a high value as we want to avoid s 4 .In particular, the optimal value for q depends on the probability that we ever visit s 2 , which is directly influenced by the value of p.
As the example indicates, trade-offs in locally monotone pMCs occur due to dependencies where parameters occur at multiple states.These trade-offs make constructing an oracle hard.Summarising, we make the following assumptions throughout the rest of this section: • We restrict the (graph-preserving) region R to be (i) rectangular, and (ii) closed.This restriction makes the bounds of the parameters independent of other parameter instantiations, and ensures that the maximum over the region exists.• We restrict the pMC D to be locally monotone, to exclude difficulties from analysing single transitions.

B. Relaxation
The idea of our approach, inspired by [57], is to drop the aforementioned dependencies between parameters by means of a relaxation of the pMC.Intuitively, the relaxation rel(D) arises from D by equipping each state with its own parameters, thereby eliminating parameter dependencies between different states (if any).This step simplifies finding an optimal instantiation (in the relaxation), but these instantiations might be spurious, i.e., not realisable in the original pMC.The instantiation rel(D)[rel(u)] corresponds to D[u] as depicted in Fig. 5(d).The relaxed region rel(R) contains also instantiations, e.g., ( 4 /5, 1 /2, 3 /5) which are not realisable in R.
For a pMC D and a graph-preserving region R, relaxation increases the set of possible instantiations: Thus, the maximal reachability probability over all instantiations of D within R is bounded by the maximum over the instantiations of rel(D) within rel(R).Lemma 3.For pMC D and region R: Pr rel(D)[u] (♦T ) .
If rel(D) satisfies a reachability property, so does D.
Corollary 1.For pMC D and region R: We now formalise the earlier observation: Without parameter dependencies, finding optimal instantiations in a pMC is simpler.Although rel(D) has (usually) more parameters than D, finding an instantiation u ∈ rel(R) that maximises the reachability probability is simpler than in u ∈ R: For any p s i ∈ rel(V ), we can in state s pick a value in I(p s i ) that maximises the probability to reach T from state s.There is no (negative) effect for the reachability probability at the other states as p s i only occurs at s. Optimal instantiations can thus be determined locally (at the states).
Furthermore, as both D is locally monotone, and there are no parameter dependencies, the maximum reachability probability is relatively easy to find: We only need to consider instantiations u that set the value of each parameter to either the lowest or highest possible value, i. e., u(p s i ) ∈ B(p s i ) for all p s i ∈ rel(V ): Theorem 4. Let D be a locally monotone pMC with states S and T ⊆ S, and a region R.There exists an instantiation u ∈ rel(R) satisfying u(p s i ) ∈ B(p s i ) for all p s i ∈ rel(V ) such that: Pr rel(D) [v] (♦T ).
To prove this statement, we consider the instantiation which assign a value to a parameter strictly between its bounds.Any such instantiation can be modified such that all parameters are assigned to its bound, without decreasing the induced reachability probability.The essential statement is the monotonicity of a parameter without any further dependencies.Lemma 5. Let D be a locally monotone pMC with a single parameter p that only occurs at one state s ∈ S, i.e.P(ŝ, s ) ∈ [0, 1] for all ŝ, s ∈ S with ŝ = s.For region R and T ⊆ S, the probability Pr D (♦T ) is monotonic on R.

Proof. W. l. o. g. let s /
∈ T be the initial state of D and let T be reachable from s. Furthermore, let U denote the standard until-modality and ¬T denote S \ T .Using the characterisation of reachability probabilities as linear equation system (cf.[15]), the reachability probability w. r. t.T (from the initial state) in D is given by: Transposing the equation yields .
The denominator can not be zero as T is reachable from s.
Since D is locally monotone, we have P(s, s ) = f s /gs for s ∈ S and multilinear functions f s , g s ∈ Q[p].We obtain: Hence, Pr D (♦T ) = f1 /f2 is a fraction of two multilinear functions f 1 , f 2 ∈ Q[p] and therefore monotonic on R.
Proof of Theorem 4. By contraposition.Let u ∈ rel(R) with Pr rel(D)[u] (♦T ) = max v∈rel(R) Pr rel(D) [v] (♦T ) .For the contraposition, assume that there exists a parameter p ∈ rel(V ) with u(p) ∈ I R (p) \ B R (p) such that all instantiations u ∈ rel(R) that set p to a value in B R (p) induce a smaller reachability probability, i.e. u (p) ∈ B R (p) and u (q) = u(q) for q = p implies Consider the pMC D = (S, {p}, s, P) with the single parameter p that arises from rel(D) by replacing all parameters q ∈ rel(V ) \ {p} with u(q).We have . Moreover, Pr D (♦T ) is monotonic on I(p) according to Lemma 5. Thus, there is an instantiation u ∈ rel(R) with u (p) ∈ B R (p) and u (q) = u(q) for q = p satisfying This contradicts our assumption for parameter p.

C. Replacing parameters by nondeterminism
In order to determine max u∈rel(R) Pr rel(D)[u] (♦T ), it suffices to make a discrete choice over instantiations u : rel(V ) → R with u(p s i ) ∈ B(p i ).This choice can be made locally at every state, which brings us to the key idea of constructing a (nonparametric) MDP out of the pMC D and the region R, where nondeterministic choices represent all instantiations that have to be considered.In the following, it is convenient to refer to the parameters in a given state s by: Definition 26 (Substitution (pMCs)).For pMC D = (S, V , s I , P) and region R, let the MDP sub R (D) = (S, s I , Act sub , P sub ) with be the (parameter-)substitution of D and R.
Thus, choosing action u in s corresponds to assigning one of the extremal values B(p i ) to the parameters p s i .The number of outgoing actions from state s is therefore 2 |Vs| .
Example 38.Consider pMC D -depicted in Fig. 13(a) -with R = [ 1 /10, 4 /5] × [ 2 /5, 7 /10] as before.The substitution of D and R is shown in Fig. 14(a).In D, each outgoing transition of states s 0 , s 1 , s 2 is replaced by a nondeterministic choice in MDP sub R (D).That is, we either pick the upper or lower bound for the corresponding variable.The solid (dashed) lines depict transitions that belong to the action for the upper (lower) bound.For the states s 3 and s 4 , the choice is unique as their outgoing transitions in D are constant.Fig. 14(b  MC sub R (D) σ which is induced by the strategy σ on MDP sub D (R) that chooses the upper bounds at s 0 and s 2 , and the lower bound at s 1 .Notice that sub R (D) σ coincides with rel(D) [v] for a suitable instantiation v, as depicted in Fig. 13(b).
The substitution encodes the local choices for a relaxed pMC.That is, for an arbitrary pMC, there is a one-to-one correspondence between strategies σ in the MDP sub rel(R) (rel(D)) and instantiations u ∈ rel(R) for rel(D) with u(p s i ) ∈ B(p i ).For better readability, we will omit the superscripts for sets of strategies Str .Combining these observations with Theorem 4, yields the following.
Corollary 2. For a pMC D, a graph-preserving region R, and a set T of target states of D: Furthermore, the nondeterministic choices introduced by the substitution only depend on the values B(p i ) of the parameters p i in R. Since the ranges of the parameters p s i in rel(R) agree with the range of p i in R, we have A direct consequence of these statements yields: Theorem 6.Let D be a pMC, R a graph-preserving region, ϕ a reachability property.Then it holds: Hence, we can deduce via Alg. 4 whether D, R |= ϕ by applying standard techniques for MDP model checking to sub R (D), such as value-and policy iteration, cf.[15], [54].We stress that while the relaxation is key for showing the correctness, equation (18) proves that this step does not actually need to be performed.The approximation error originates from choices where an optimal strategy on sub R (D) chooses actions u 1 and u 2 at states s 1 and s 2 , respectively, with u 1 (p s1 i ) = u 2 (p s2 i ) for some parameter p i , and therefore intuitively disagree on its value.The probability mass that is affected by these choices decreases the smaller the region is.For infinitesimally small regions, the error from the over-approximation vanishes, as the actions for the upper and the lower bound of a parameter become equal up to an infinitesimal.

D. Expected reward properties
The reduction of bounding reachability probabilities on pMCs to off-the-shelf MDP model checking can also be applied to bound expected rewards.To see this, we have to extend the notion of locally monotone parametric Markov chains.

Definition 27 (Locally monotone reward pMC).
A pMC D = (S, V , s I , P) with reward function rew : S → Q(V ) is locally monotone iff for all s ∈ S, there is a multilinear polynomial We now generalise relaxation and substitution to the reward models, and obtain analogous results.
Definition 28 (Substitution for reward pMCs).Let D = (S, V , s I , P) be a pMC, rew : S → Q(V ) a reward function, T ⊆ S a set of target states, and R a region.For s ∈ S, let The MDP sub rew R (D) = (S, s I , Act rew sub , P rew sub ) with reward function rew sub is the (parameter-)substitution of D, rew on R, where • rew sub is given by: The reward approximation of a pMC can be used to identify regions as accepting or rejecting for expected reward properties.Theorem 7. Let D be a pMC with rewards rew, R be a graph-preserving region, and ϕ an expected reward property.
The proof is analogous to the proof of Thm. 6.

VII. MODEL-CHECKING-BASED REGION VERIFICATION OF PARAMETRIC MARKOV DECISION PROCESSES
In the last section, we approximated reachability probabilities in (locally-monotone) pMCs by considering the substitution MDP, see Def. 26.The non-determinism in the MDP encodes the finitely many parameter valuations that approximate the reachability probabilities in the pMC.By letting an adversary player resolve the non-determinism in the MDP, we obtain bounds on the reachability probabilities in the pMC.These bounds can efficiently be computed by standard MDP model checking.
In this section, we generalise approach to pMDPs, which already contain non-determinism.The result naturally leads to a 2-player stochastic game: One player controls the nondeterminism inherent to the MDP, while the other player controls the (abstracted) parameter values.Letting the two players adequately minimise and/or maximise the reachability probabilities in the SG yields bounds on the minimal (and maximal) reachability probabilities in the pMDP.For example, if the player for the original non-determinism maximises and the parameter player minimises, we obtain a lower bound on the maximal probability.These bounds can efficiently be computed by standard SG model checking procedures.
In our presentation below, we discuss the interplay of the two sources of non-determinism.In particular, we show how the generalisation of the method yields an additional source of (over-)approximation.Then, we formalise the construction of the substitution with nondeterminism, analogous to the pMCs from the previous section.In particular, Def.29 is analogous to Def. 26 and Thm. 8 is analogous to Thm. 6. Repeating the concept of relaxation, described in Sect.VI-B, is omitted, as-as also discussed in the last section-it is not a necessary ingredient for the correctness of the approach.

A. Observation
In the following, let M = (S, V , s I , Act, P) be a pMDP and R a graph-preserving, rectangular, closed region.
We analyse R with respect to the demonic relation |= d .We have: The two universal quantifiers can be reordered, and in addition Intuitively, the reformulation states that we have to apply pMC region verification on M σ and R for all σ ∈ Str M .We now want to employ parameter lifting for each strategy.Thus, we want to consider the verification of the substituted pMCs sub R (M σ ).As these substituted pMCs share most of their structure, the set of all such substituted pMCs can be concisely represented as an SG, in which both players cooperate (as witnessed by the same quantifiers).In the scope of this paper, an SG with cooperating players can be concisely represented as an MDP.Consequently, for the demonic relation, pMDP verification can be approximated by MDP model checking.
We now turn our attention to the angelic relation |= a , cf.Def.16.
Here, we cannot simply reorder the quantifiers.However: Now, the left-hand side can be concisely represented as an SG (as in the demonic case).As witnessed by the quantifier elimination, this SG does not reduce to an MDP; the two players have opposing objectives.Nevertheless, we can efficiently analyse this SG (with a variant of value iteration), and thus the left-hand side of the implication above.
Observe that the over-approximation actually computes a robust strategy, as discussed in Remark 6 on page 11.In particular, we now have two sources of approximation: • The approximation that originates from dropping parameter dependencies (as also in the demonic case).• The application of the substitution of parameters with non-determinism on robust strategies rather than of the actual angelic relation.Both over-approximations vanish with declining region size.

B. Replacing parameters by nondeterminism
Example 40.Consider the pMDP M in Fig. 15(a), where state s has two enabled actions α and β.The strategy σ given by {s → α} applied to M yields a pMC, which is subject to substitution, cf.Fig. 15(b).
The parameter substitution of a pMDP (cf.Fig. 15(a)) yields an SG-as in Fig. 15(d).It represents, for all strategies of the pMDP, the parameter-substitution (as in Def.26) of each induced pMC.To ensure that in the SG each state can be assigned to a unique player, we split states in the pMDP which have both (parametric) probabilistic branching and nondeterminism, such that states have either probabilistic branching or non-determinism, but not both.The reformulation is done as follows: After each choice of actions, auxiliary states are introduced, such that the outcome of the action becomes deterministic and the probabilistic choice is delayed to the auxiliary state.This construction is similar to the conversion of Segala's probabilistic automata into Hansson's alternating model [58].More precisely, we • split each state s ∈ S into {s} { s, α | α ∈ Act(s)}, • add a transition with probability one for each s ∈ S and α ∈ Act(s).The transition leads from s to s, α , and • move the probabilistic choice at s w. r. t. α to s, α .Applying this to the pMDP from Fig. 15(a), we obtain the pMDP M in Fig. 15(c), where state s has only nondeterministic choices leading to states of the form s, α with only probabilistic choices.The subsequent substitution on the probabilistic states yields the SG sub R (M ), where one player represents the nondeterminism of the original pMDP M, while the other player decides whether parameters should be set to their lower or upper bound in the region R.For the construction, we generalise V s to state-action pairs: For a pMDP, a state s and action α, let with and, be the (parameter-)substitution of M and R.
We relate the SG sub R (M) under different strategies for player with the substitution in the strategy-induced pMCs of M. We observe that the strategies for player in sub R (M) coincide with strategies in M. Consider the induced MDP (sub R (M)) σ with a strategy σ for player .The MDP (sub R (M)) σ is obtained from sub R (M) by erasing transitions not agreeing with σ.In (sub R (M)) σ player -state have a single enabled action, while player -states have multiple available enabled actions.
Example 41.Continuing Example 40, applying strategy σ to sub R (M) yields (sub R (M)) σ , see Fig. 15(e).The MDP (sub R (M)) σ matches the MDP sub R (M σ ) apart from intermediate states of the form s, α : The outgoing transitions of s in sub R (M σ ) coincide with the outgoing transitions of s, α in (sub R (M)) σ , where s, α is the unique successor of s.
The following corollary formalises that (sub R (M)) σ and sub R (M σ ) induce the same reachability probabilities.
Corollary 3.For pMDP M, graph-preserving region R, target states T ⊆ S, and strategies σ ∈ Str sub R (M) and ρ ∈ Str sub R (M σ ) , it holds that Instead of performing the substitution on the pMC induced by M and σ, we can perform the substitution on M directly and preserve the reachability probability.
Consequently, and analogously to the pMC case (cf.Theorem 6), we can derive whether M, R |= ϕ by analysing a stochastic game.For this, we consider various standard variants of model checking on stochastic games.Definition 30 (Model-relation on SGs).For an SG G, property ϕ, and quantifiers Q 1 , Q 2 , we define G |= Q1,Q2 ϕ as: The order of players, for these games, does not influence the outcome [43], [59].
Theorem 8. Let M be a pMDP, R a region, and ϕ a reachability property.Then: Proof.We only prove the second statement using ϕ = P ≤λ (♦T ), other reachability properties are similar.A proof for the (simpler) first statement can be derived in an analogous manner.We have that M, R |= ¬P ≤λ (♦T ) iff for all u ∈ R there is a strategy σ of M for which the reachability probability in the MC M σ [u] exceeds the threshold λ, i. e., A lower bound for this probability is obtained as follows: The inequality * is due to Corollary 2. The equality * * holds by Corollary 3. Then:

VIII. APPROXIMATE SYNTHESIS BY PARAMETER SPACE PARTITIONING
Parameter space partitioning is our iterative approach to the approximate synthesis problem.It builds on top of region verification, discussed above, and is, conceptually, independent of the methods used for verification discussed later.
Parameter space partitioning is best viewed as a counterexample guided abstraction refinement (CEGAR)-like [60] approach to successively divide the parameter space into accepting and rejecting regions.The main idea is to compute a sequence R i a i of simple accepting regions that successively extend each other.Similarly, an increasing sequence R i r i of simple rejecting regions is computed.At the i-th iteration, is the covered fragment of the parameter space.The iterative approach halts when R i is at least c% of the entire parameter space.Termination is guaranteed: in the limit a solution to the exact synthesis problem is obtained as lim i→∞ R i a = R a and lim i→∞ R i r = R r .Let us describe the synthesis loop for the approximate synthesis as depicted in Fig. 4 on page 5 in detail.In particular, we discuss how to generate candidate regions that can be dispatched to the verifier along with a hypothesis whether the candidate region is accepting or rejecting.We focus on rectangular regions for several reasons: • the automated generation of rectangular regions is easier to generalise to multiple dimensions, • earlier experiments [40] revealed that rectangular regions lead to a more efficient SMT-based verification of regions (described in Sect.V), and • model-checking based region verification (described in Sect.VI) requires rectangular regions.A downside of rectangular regions is that they are neither wellsuited to approximate a region partitioning given by a diagonal, nor to cover well-defined regions that are not rectangular themselves.
Remark 10.In the following, we assume that the parameter space is given by a rectangular well-defined region R.If the parameter space is not rectangular, we over-approximate R by a rectangular region R ⊇ R. If the potential overapproximation of the parameter space R is not well-defined, then we iteratively approximate R by a sequence of welldefined and ill-defined 10 regions.The regions in the sequence 10 A region R is ill-defined if no instantiation in R is well-defined.of well-defined regions are then subject to the synthesis problem.Constructing the sequence of regions is done analogously to the partitioning into accepting and rejecting regions.
Before we present the procedure in full detail, we first outline a naive refinement procedure by means of an example.
Example 42 (Naive refinement loop).Consider the parametric die from Example 5. Suppose we want to synthesise the partitioning as depicted in Fig. 2 on Page 4. We start by verifying the full parameter space R against ϕ.The verifier returns false, as R is not accepting.Since R (based on our knowledge at this point) might be rejecting, we invoke the verifier with R and ¬ϕ, yielding false too.Thus, the full parameter space R is inconsistent.We now split R into four equally-sized regions, all of which are inconsistent.Only after splitting again, we find the first accepting and rejecting regions.After various iterations, the procedure leads to the partitioning in Fig. 16.
Alg. 5 describes this naive region partitioning procedure.It takes a pSG, a region R, a specification ϕ, and a (demonic or angelic) satisfaction relation as input.It first initialises a (priority) queue Q with R. In each iteration, a subregion R of R is taken from the queue, the counter i is incremented, and the sequence of accepted and rejected regions is updated.There are three possibilities.Either R is accepting (or rejecting), and In the latter case, we split R into a finite set of subregions that are inserted into the queue Q. Regions that are not extended are unchanged.
The algorithm only terminates if R a and R r are a finite union of hyper-rectangles.However, the algorithm can be terminated after any iteration yielding a sound approximation.The algorithm ensures lim i→∞ R i = R, if we order Q according to the size of the regions.We omit the technical proof here; the elementary property is that the regions are Lebesguemeasurable (and have a positive measure by construction).
The naive algorithm has a couple of structural weaknesses: • It invokes the verification algorithm twice to determine that the full parameter space is inconsistent.region without samples is obtained: rather than inserting (R , ∅) into Q, we insert the entry (R , sample(R )).
Example 45.After several more iterations, the refinement loop started in Ex. 43 has proceeded to the state in Fig. 17(b).First, we see that the candidate region from Fig. 17(a) was not rejecting.The verification engine gave a counterexample in form of an accepting sample (around p → 0.45, q → 0.52).
Further iterations with smaller regions had some successes, but some additional samples were generated as counterexamples.
The current blue candidate is to be checked next.In Fig. 17(c), we see a further continuation, with even smaller regions being verified.Note the white box on the right border: It has been checked, but the verification timed out without a conclusive answer.Therefore, we do not have a counterexample in this subregion.
It remains to discuss some methods to split a region, and how for some of the constructed regions, verification may be be skipped.We outline more details below.
1) How-to split: Splitting of regions based on the available samples can be done using different strategies.We outline two basic approaches.These approaches can be easily mixed and extended, and their performance heavily depends on the concrete example at hand.a) Equal splitting: This approach splits regions in equallysized regions; the main rationale is that this generates small regions with nice bounds (the bounds are typically powers of two).Splitting in equally sized regions can be done recursively: One projects all samples down to a single dimension, and splits if both accepting and rejecting samples are in the region.The procedure halts if all samples in a region are either accepting or rejecting.The order in which parameters are considered plays a crucial role.Typically, it is a good idea to first split along the larger dimensions.
Example 46.A split in equally-sized regions is depicted in Fig. 18(b), where first the left region candidate is created.The remaining region can be split either horizontally or vertically to immediately generate another region candidate.A horizontal split in the remaining region yields a region without any samples.
The downside of equal splitting is that the position of the splits are not adapted based on the samples.Therefore, the number of splits might be significantly larger than necessary, leading to an increased number of verification calls.
b) Growing rectangles: This approach tries to gradually obtain a large region candidate 11 .The underlying rationale is to quickly cover vast amounts of the parameter space.This is illustrated in Fig. 18(d) (notice that we adapted the samples for a consistent but concise description) where from an initial sampling a large rectangle is obtained as region candidate.
Example 47.Consider the shaded regions in Fig. 18(c).Starting from vertex v = (1, 1), the outer rectangle is maximised to not contain any accepting samples.Taking this outer rectangle as candidate region is very optimistic, it assumes that the accepting samples are on the border.A more pessimistic variant of growing rectangles is given by the inner shaded region.It takes a rejecting sample as vertex v such that the v and v span the largest region.
The growing rectangles algorithm iterates over a subset of the hyper-rectangle's vertices: For each vertex (referred to as anchor), among all possible sub-hyper-rectangles containing the anchor and only accepting or only rejecting samples, the largest is constructed.
Example 48.The growing rectangles approach in its pessimistic fashion takes anchor (0, 0) as anchor and yields the candidate region in Fig. 18(d).
The verification fails more often on large regions (either due to time-outs or due to the over-approximation).Consequently, choosing large candidate regions comes at the risk of failed verification calls, and fragmentation of the parameter space in more subregions.
Furthermore, growing rectangles requires a fall-back splitting strategy: To see why, consider Fig. 16 on page 25.The accepting (green) region does not contain any anchors of the full parameter space, therefore the hypothesis for any created subregion is always rejection.Thus, no subregion containing a (known) accepting sample is ever considered as region candidate.
2) Neighbourhood analysis: Besides considering samples within a region, we like to illustrate that analysis of a region R can and should take information from outside of R into account.First, take Fig. 18(b), and assume that the left region is indeed accepting.The second generated region contains only rejecting samples, but it is only rejecting if all points, including all those on the border to the left region, are rejecting.In other words, the border between the accepting and rejecting regions needs to exactly follow the border between the generated region candidates.The latter case does not occur often, so it is reasonable to shrink or split the second generated region.Secondly, a sensible hypothesis for candidate regions without samples inside is helpful, especially for small regions or in high dimensions.Instead of spawning new samples, we take samples and decided regions outside of the candidate region into account to create a hypothesis.Concretely, we infer the hypothesis for regions without samples via the closest known region or sample.

C. Requirements on verification back-ends
In this section, we have described techniques for iteratively partitioning the parameter space into accepting and rejecting regions.The algorithms rely on verifying regions (and sets of samples) against the specification ϕ.The way in which verification is used in the iterative parameter space partitioning scheme imposes the following requirements on the verification back-end: 1) The verification should work incrementally.That is to say, verification results from previous iterations should be re-used in successive iterations.Verifying different regions share the same model (pMC or pMDP).A simple example of working incrementally is to reuse minimisation techniques for the model over several calls.If a subregion is checked, the problem is even incremental in a more narrow sense: any bounds etc. obtained for the superregion are also valid for the subregion.2) If the verification procedure fails, i.e. if the verifier returns false, obtaining additional diagnostic information in the form of a counterexample is beneficial.A counterexample here is a sample which refutes the verification problem at hand.This wish list is very similar to the typical requirements that theory solvers in lazy SMT frameworks should fulfil [61].Therefore, SMT-based verification approaches naturally match the wish-list.Parameter-lifting can work incrementally: it reuses the graph-structure to avoid rebuilding the MDP, and it may use previous model checking results to improve the time until the model checker converges.due to its approximative nature, does only provide limited diagnostic information.

IX. IMPLEMENTATION
All the algorithms and constructions in this paper have been implemented, and are publicly available via PROPhESY12 .In particular, PROPhESY supports algorithms for: • the exact synthesis problem: via computing the solution function, using either of the three variants of stateelimination, discussed in Sect.IV. • the verification problem: via an encoding to an SMTsolver as in Sect.V or by employing the parameter lifting method as in Sect.VI and VII.• the approximate synthesis problem: via parameter space partitioning, that iteratively generates verification calls as described in Sect.VIII.
PROPhESY is implemented in python, and designed as a flexible toolbox for parameter synthesis.PROPhESY internally heavily relies on high-performance routines of probabilistic model checkers and SMT-solvers.In particular, the computation of the solution function and the parameter lifting have been implemented and tightly integrated in the probabilistic model checker Storm [34].
PROPhESY can be divided in three parts: 1) First and foremost, it presents a library consisting of: a) data structures for parameter spaces and instantiations, solution functions, specifications, etc., built around the python bindings of the library carl 13 (featuring computations with polynomials and rational functions), b) algorithms such as guided sampling, various candidate region generation procedures, decomposition of regions, etc., methods that require tight integration with the model are realised via the python bindings of Storm14 , c) abstract interfaces to backend tools, in particular probabilistic model checkers, and SMT-checkers, together with some concrete adapters for the different solvers, see Fig. 19.Storm [34] Storm (python) z3 [47] SMT-RAT [38] Fig. 19.High-level architecture of PROPhESY and its backends 2) An extensive command-line interface which provides simple access to the different core functionalities of the library, ranging from sampling to full parameter synthesis.3) A prototypical web-service running on top of the library, which allows users to interact with the parameter synthesis via a web-interface.PROPhESY is constructed in a modular fashion: besides the python bindings for carl, all non-standard packages and tools (in particular model checkers and SMT solvers) are optional.Naturally, the full power of PROPhESY can only be used if these packages are available.
Besides the methods presented in this paper, PROPhESY contains two further mature parameter synthesis methods: 1) particle-swarm optimisation inspired by [62], and 2) convex optimisation from [63].
The information in the remainder details the implementation and the possibilities provided by PROPhESY.The section contains some notions from probabilistic model checking [15], [19], [20].We refrain from providing detailed descriptions of these notions, as it would go beyond the scope of this paper.

A. Model construction and preprocessing (Realised in Storm)
The model checker Storm supports the creation of pMCs and pMDPs from both PRISM-language model descriptions [33] and JANI-specifications [64].The latter can be used as intermediate format to support, e.g., digital-clock PTAs with parameters written in Modest [65], or to support expected time properties of generalised stochastic Petri nets [66] with parametric rates and/or weights.Parametric models can be built using the matrix-based, explicit representation, as well as the symbolic, decision diagram (dd)-based engine built on top of sylvan [67].Both engines support the computation of qualitative properties, an essential preprocessing step, and bisimulation minimisation on parametric models, as described in [49].We advocate the use of the Storm-python API adapter: Its interactive nature avoids the repetition of expensive steps.In particular, it allows for the incremental usage of parameter lifting and sampling.
The support for rational functions is realised via the library carl 15 .The rational function is stored as a tuple consisting of multivariate polynomials.These polynomials are by default stored in a partially factorised fashion, cf.[52].Each factor (a polynomial) is stored as an ordered sparse sum of terms, each term consists of the coefficient and a sparse representation of variables with their non-zero exponents.For manipulating the (rational) coefficients, we exploit gmp 16 or cln 17 .The former is thread-safe, while the latter performs slightly better with single-thread usage.Computation of GCDs in multivariate polynomials is done either via ginac [68] or cocoa [69].

B. Solution function computation (Realised in Storm)
The computation of solution functions for pMCs as discussed in Sect.IV is implemented for a variety of specifications: • reachability and reach-avoid probabilities, • expected rewards, including expected time of continuoustime Markov chains, • step-bounded reachability probabilities, and • long-run average probabilities and rewards.The computation is realised either via state elimination, or via Gaussian elimination.An implementation of set-based transition elimination is available for symbolic representations of the pMC.
1) State elimination: As the standard sparse matrix representation used by Storm is not suitable for fast removal and insertion of entries, a flexible sparse matrix with faster delete and insert operations is used.
The order in which states are eliminated has a severe impact on the performance [40].Storm supports a variety of static (pre-computed) and dynamic orderings for the elimination: • several static orders (forward (reversed), backward (reversed)) based on the order of state-generation by the model construction algorithms.This latter order is typically determined by a depth-first search through the highlevel model description 18 , • orders based on the topology of the pMC, e.g., based on the decomposition in strongly connected components, • orders (Regex) which take into account the in-degree (the number of incoming transitions at a state), inspired by [51], [70], • orders (SPen, DPen) which take into account the complexity of the rational function corresponding to the transition probability.The complexity is defined by the degree and number of terms of the occurring polynomials.The orders are computed as penalties for states, and the order prefers states with a low penalty.For dynamic orderings (Regex, DPen), the penalties are recomputed as the in-degree of states and complexity of transition probabilities change during stateelimination.
2) Gaussian elimination: Storm supports Eigen [71] as a linear equation system solver over the field of rational functions.It uses the "supernodal" (supernodes) LU factorisation.The matrix is permuted by the column approximate minimum degree permutation (COLAMD) algorithm to reorder the matrix.One advantage is that this solver is based on sparse model-checking algorithm for parameter-free models.The solver therefore, in addition to the properties supported by state-elimination, supports the construction in [72] for conditional probabilities and rewards.
3) Set-based transition elimination: This elimination method is targeted for symbolic representations of the Markov chain.Set-based transition elimination is implemented via matrixmatrix multiplications.In every multiplication, a copy of the dd-representation of a matrix over variables ( s, t) is made.The copy uses renamed dd-variables ( t, t ).Then, a multiplication of the original matrix with the copy can be done on the dd level yielding a matrix ( s, t ).Renaming t to t yields a matrix on the original dd-variables.

C. Parameter lifting (Realised in Storm)
For parameter lifting (Sect.VI and VII), the major effort beyond calling standard model-checking procedures is the construction of the substituted (lifted) model.As parameter lifting for different regions does not change the topology of the lifted model, it is beneficial to create a template of the lifted model once, and to substitute the values according to the region at hand.The substitution operation can be sped up by exploiting the following observation: Typically, transition probability functions coincide for many transitions.Thus, we evaluate each occurring function once and substitute the outcome directly at all occurrences.Moreover, for a growing number of regions to be checked, any one-time preprocessing of the lifted model eventually pays off.In particular, we apply minimisation techniques before construction of the lifted model.We use both bisimulation minimisation as well as stateelimination of parameter-free transitions.These minimisations drastically reduce the run-time of checking a single region.We use numerical methods first: for regions that we want to classify as accepting (or rejecting) we resort to the analysis of MDPs using policy iteration with rational numbers.For that, we initialise the policy iteration with a guess based on the earlier numerical results.

D. SMT-based region verification (Realised in PROPhESY)
This complete region checking procedure is realised by constructing SMT queries, as elaborated in Sect.V. When invoking the SMT solver, we use some features of the SMTlib standard [73].First of all, when checking several regions, we use backtrack-points to only partly reset the solver: More precisely, the problem description is given by a conjunction of subformulae, where the conjunction is represented by a stack.We first push the constraints for the problem to the stack, save a backtrack point, and then store the region.Once we have checked a particular region, we backtrack to the back-track point, that is, we remove the constraints for the particular region from the problem description.This way, we reuse simplifications and data structures the solver constructed for the problem description covering the model (and not the region).To support both verifying the property and its negation, the problem description is slightly extended.We add two Boolean variables (accepting and rejecting).The following gives an example of the encoding together with checking whether a region R 1 is accepting, and a region R 2 is rejecting, using the notation of Sect.V.
We accelerate the selection of regions by getting a rough picture through sampling, as discussed in Sect.VIII.We support two engines for computing the samples: Either via model checking, or by instantiating the solution function.Sampling on the solution function should always be done exactly, as the evaluation of the typically highly-nonlinear solution functions is (again typically) numerically unstable.In each iteration, based on the current set of samples, a new set of sampling candidates is computed.The choice of the new samples can be modified in several ways.The standard used here is via linear interpolation between accepting and rejecting samples.

F. Partitioning (Realised in PROPhESY)
For the construction of region candidates, we split the initial regions according to our heuristic (quads or growing rectangles, cf.Sect.VIII-B) until none of the regions is inconsistent.We sort the candidate regions based on their size in descending order.Furthermore, we prefer regions where we deem verification to be less costly: Candidate regions that are supposed to be accepting and are further away from samples or regions that are rejecting are preferred over those regions which have rejecting samples or regions in their neighbourhood.

X. EMPIRICAL EVALUATION
In this section, we show the applicability of the presented approaches based on a selection of benchmarks.
A. Set-up 1) Benchmarks: We consider five case studies from the literature.The selection represents various application domains.
a) NAND multiplexing: With integrated circuits being build at ever smaller scale, they are more prone to defects and/or to exhibit transient failures [74].One way to overcome these deficiencies is the implementation of redundancy at gatelevel.In particular, one aims to construct reliable devices from unreliable components.NAND multiplexing is such a technique, originally due to von Neumann [75].Automated analysis of NAND multiplexing via Markov chain model checking was considered first in [76].They also studied the influence of gate failures in either of the stages of the multiplexing by sampling various values.We take the pMC from [40], that replaced these probabilities with parameters.We analyse the effect of changing failure probabilities of the gates on the reliability of the multiplexed NAND.b) Herman's self-stabilising protocol: In distributed systems, tokens are used to grant privileges (e.g., access to shared memory) to processes.Randomisation is an essential technique to break the symmetry among several processes [77].Herman's probabilistic algorithm [8] is a token circulation algorithm for ring structures.In each step, every process possessing a token passes the token along with probability p and keeps the token with probability 1−p.The algorithm is self-stabilising, i.e., started from any illegal configuration with more than one token the algorithm recovers to a legal configuration with a unique token.The recovery time crucially depends on the probability of passing the token, and an optimal value for p depends on the size of the system [9].We investigate the expected recovery time by parameter synthesis, inspired by [78].
c) Mean-time-to-failure of a computer system: In reliability engineering, fault trees are a prominent model to describe how a system may fail based on faults of its various components [1], [2].Dynamic fault trees (DFTs, [79]) extend these fault trees with a notion of a state, and allow to model spare management and temporal dependencies in the failure behaviour.State-of-the-art approaches for dynamic fault trees translate such fault trees into Markov chains [29], [30], [80]; evaluation of the mean-time-to failure boils down to the analysis of the underlying Markov chain.Probabilities and rewards originate from the failure rate of the components in the described system.Such failure rates are often not known (precisely), especially during design time.Therefore, they may be represented by parameters.We take the HECS DFT [81] benchmark describing the failure of a computer system with an unknown failure rate for the software interface and the spare processor, as first described in [82].We analyse how this failure rate affects the expected time until the failure (meantime-to-failure) of the complete computer system.d) Network scheduling: This benchmark [83] concerns the wireless downlink scheduling of traffic to different users, with hard deadlines and prioritised packets.The system is time-slotted: time is divided into periods and each period is divided into an equal number of slots.At the start of each time period, a new packet is generated for each user with a randomly assigned priority.The goal of scheduling is to, in each period, deliver the packets to each user before the period ends.Packets not delivered by the end of a period are dropped.Scheduling is non-trivial, as successful transmissions are not stochastically independent, i.e., channels have a (hidden) internal state.The system is described as a partially observable Markov decision process [84], a prominent formalism in the AI community.We take the Network model from [85], and consider the pMC that describes randomised finite memory controllers that solve this scheduling problem, based on a translation from [86].Concretely, the parameters represent how the finite memory controller randomises.We evaluate the effect of the randomisation in the scheduling on the expected packet loss.
e) Bounded retransmission protocol: The bounded retransmission protocol (BRP, [87], [88]) is a variant of the alternating bit protocol.It can be used as part of an OSI data link layer, to implement retransmitting corrupted file chunks between a sender and a receiver.The system contains two channels; from  [89].We consider the parametric version from [49].
The parameters naturally reflect the channel qualities.The model contains non-determinism as the arrival of files on the link layer cannot be influenced.This non-determinism hampers a manual analysis.The combination of parametric probabilities and non-determinism naturally yields a pMDP.We analyse the maximum probability that a sender eventually does not report a successful transmission.
Remark 11.Other benchmarks and a thorough performance evaluation have been presented before in [40] (for stateelimination and parameter space partitioning) and [39] (for parameter lifting).
2) Benchmark statistics: Table I summarises relevant information about the concrete instances that we took from the benchmarks.The id is used for reference.The benchmark refers to the name of the benchmark-set, while the instance describes the particular instance from this benchmark set.We give the 3) Evaluation: We conducted the empirical evaluation on an HP BL685C G7 with Debian 9.6.Each evaluation run could use 8 cores with 2.1GHz each.However, most runs only made use of one core.We set the timeout to 1 hour and the memory limit to 16GB.We used PROPhESY version 2.0, together with the Storm-python bindings version 1.3.1,z3 version 4.8.4.All benchmark files are made available via PROPhESY.

B. Exact synthesis via the solution function
To evaluate the exact synthesis approach, we use state elimination with 7 different heuristics, set-based transition elimination, and Gaussian elimination.All configurations are evaluated with and without strong bisimulation.
First, we show the sizes of the solution function: The results are summarised in Table II.The id references the corresponding benchmark instance in Table I.The BRP model is not included: The set of all strategies is significantly too large.The next four columns state properties of the resulting rational function.We give the degree of both the numerator (degree num) and denominator (degree denom), as well as the number of terms in both polynomials (# terms num, # terms denom).The next column gives the number of configurations (out of the 18) which successfully finished within the time limit.The last two columns indicate timings.We give the times (in seconds) to compute the solution function (time mc) and the total time including model building, (optional) bisimulation minimisation and computing the solution function.For these timings we give two numbers per benchmark instance: The upper row describes the median value over all successful configurations and the lower row describes the best result obtained.Thus, while functions often grow prohibitively large, medium-sized functions can still be computed.Contrary to model checking for parameter-free models, model building is typically not the bottleneck.
Furthermore, we see that the selected heuristic is indeed crucial.Consider instance 11: 11 heuristics successfully compute the solution function (and most of them within a second).However, 7 others yield a timeout.That leads us to compare some heuristics in Fig. 20.The plot depicts the cumulative solving times for selected configurations over all 18 benchmark instances (excluding BRP).Gaussian and set-based refer to these approaches, respectively, all other configurations are variants of state-elimination, cf.Sect.IX-B1, (bisim) denotes that bisimulation minimisation is used.The x-axis represents the number of solved instances and the (logarithmic) y-axis represents the time in seconds.A point (x, y) in the plot represents the x fastest instances which could be solved within a total time of y seconds.For 15 instances, one of the depicted configurations was the fastest overall.Regex based configurations were the fastest eight times, DPen based ones four times and three times configurations based on FwRev were fastest.From these numbers, we conclude that the selection of the heuristic is essential, and depending on the model to be analysed.From the graph, we further observe that although using a Gaussian elimination yields good performance, stateelimination based approaches can (significantly) outperform the Gaussian elimination on some benchmarks.The DPen solves all instances (the only configuration to do so), but Regex is overall (slightly) faster.The uninformed FwRev with bisimulation works surprisingly well for these benchmarks (but that is mostly coincidence).The set-based elimination is clearly inferior on the benchmarks considered here, but allows to analyse some models with a very regular structure and a  gigantic state space, e.g., a parametric Markov chain for the analysis of the bluetooth protocol [90].

C. Three types of region verification
We evaluate region verification using two SMT-based approaches (SF: based on first computing the Solution Function, or ETR: encoding the equations into Existential Theory of the Reals), and PLA.In particular, we present some results for the Herman benchmark: it features a single parameter, and therefore is well-suited for the illustration of some concepts.We visualised the results for instance 11 in Fig. 21.The xaxis represents the probability p and the y-axis the expected recovery time.We indicate the solution function in blue.The threshold in the following is set to λ = 5 and indicated by the orange horizontal line.The black columns depict six different regions 19 that are evaluated with region checking.For each region we want to verify whether the expected recovery time is at least 5.The results are summarised in (the upper part of) Table III.The first column id references the benchmark instance and the second column gives the threshold λ.The next columns indicate the considered region and the technique.The last columns give the result of the region verification and the time (in seconds) needed for the computation.The timeout (TO) was set to 120 seconds.
For benchmark instance 11, PLA computes a result within milliseconds and the computation time is independent of the considered region.The SMT-based techniques take longer and the SF technique in particular does not terminate within two minutes.However, the ETR technique could yield a result for region [0.28, 0.35] whereas PLA could not give a conclusive answer due to its inherent over-approximation.
We now consider the region verification on the NAND model with two parameters.We visualised the solution function for instance 13 in Fig. 22.The considered threshold is λ = 0.3.Green coloured parts indicate parameter instantiations leading to probabilities above λ and red parts lie below λ.The results of the verification for different regions are given in (the lower part of) Table III.PLA is again the fastest technique, but for larger regions close to the threshold PLA can often not provide a conclusive answer.Contrary to before, SF is superior to ETR.The performance of the SMT-based techniques (again) greatly depends on the considered region.It is only natural that the size of the region, and the difference to the threshold have a significant influence on the performance of region verification.These observations are general and do hold on all other benchmarks.Furthermore, parameter lifting seems broadly applicable, and in the setting evaluated here, clearly faster than SMT-based approaches.Parameter lifting overapproximates and therefore might only give a decisive result in an refinement loop such as parameter space partitioning.The SMT-based approaches are a valuable fallback.When relying on the SMT techniques, it is heavily model-dependent which performs better.Table IV   gives some additional results, indicating the performance of the different verification techniques.

D. Approximative synthesis via parameter space partitioning
We now evaluate the parameter space partitioning.We use the implementation in PROPhESY with the three verification procedures evaluated above.Therefore, we focus here on the actual parameter space partitioning.
First, consider again Herman for illustration purposes.Region verification is not applicable for instance 10 (with threshold 5), as neither all instantiations accept nor all reject the specification.Instead, parameter space partitioning delivers which of these instantiations accept, and which reject the specification.The resulting parameter space partitioning is visualised in Fig. 23.
Next, we compare the three verification techniques-each with two different methods for selecting candidate regions-in Fig. 24.Fig. 24(a) depicts the computation on the Herman model with 5 processes and threshold λ = 5.The plot depicts the covered area for all three techniques with both quads (straight lines) and rectangles (dashed lines) as regions.The x-axis represents the computation time (in seconds) on a logarithmic scale and the y-axis represents the percentage of covered area.A point (x, y) in the plot represents y percent of the parameter space which could be covered within x seconds.
For Herman, SMT-based techniques perform better than PLA.PLA was able to cover 64% of the parameter space within milliseconds.However, in the remaining hour only 2% more space was covered.The SMT-based techniques were able to cover at least 99% of the parameter space within 15 seconds.Moreover, the rectangles cover the parameter space faster than quads.We also perform the parameter space partitioning on the NAND model with two different thresholds: We compare the parameter space partitioning techniques for threshold λ = 0.1 in Fig. 24(b), and for threshold λ = 0.3 in Fig. 24(c).For NAND, the PLA technique performs better than the SMT-based techniques.For threshold λ = 0.1, PLA could cover at least 99% of the parameter space within 1 second.The main reason is that the border is in a corner of the parameter space.Additionally, the SMT-based techniques with rectangles are significantly faster than the quads for this threshold.For threshold λ = 0.3, more region verification steps were necessary.PLA still outperforms ETR and SF.However, the use of rectangles over quads does not lead to a better performance for this threshold.At any point in time, there can be very significant differences between the heuristics for candidate generation, especially in settings where single region verification calls become expensive.Finally, we summarise an overview of the performance in Table IV.For brevity, we pruned some rows, especially if the present approaches already struggle with smaller instances.The id is a reference to the benchmark instance.The technique is given in the next column.In the next three columns we give for each technique the time (in seconds) needed to cover at least 50%, 90% and 98% of the complete parameter space.The next two columns give the complete covered area-i.e. the sum of the sizes of all accepting or rejecting regionswhen terminating the parameter space partitioning after 1h, together with the safe area, i.e. the sum of the sizes of all accepting regions.The last two columns indicate the percentage of the total time spent in generating the regions (time reg gen) and verifying the regions (time analysis).PLA is almost always superior, but not on all benchmarks (and not on all

XI. RELATED WORK AND DISCUSSION
We discuss related work with respect to various relevant topics.a) Computing a solution function: This approach was pioneered by [35] and significantly improved by [49].Both PRISM [33] and PARAM [42] support the computation of a solution function based on the latter method.It has been adapted in [52] to an elimination of SCCs and a more clever representation of rational functions.This representation has been adapted by Storm [34].In [91], computing a solution function via a computer algebra system was considered.That method targets small, randomly generated pMCs with many parameters.Recently, [53] explored the use of one-step fraction-free Gaussian elimination to reduce the number of GCD computations.For pMDPs, [92] experimented with the introduction of discrete parameters to reflect strategy choicesthis method, however, scales poorly.In [93] and [94], variants of value iteration with a dd-based representation of the solution function are presented.b) Equation system formulation: Regarding pMDPs, instead of introducing a Boolean structure, one can lift the linear program formulation for MDPs to a nonlinear program (NLP).This lifting has been explored in [95], and shown to be not feasible in general.Although the general NLP does not lie in the class of convex problems, a variety of verification related problems can be expressed by a sequence of geometric programs, which is exploited in [96].Alternatively, finding satisfying parameter instantiations in pMDPs under demonic non-determinism and with affine transition probabilities can be approached by iteratively solving a convex-concave program that approximates the original NLP [63].Alternatively, more efficient solvers can be used [97] for subclasses of pMDPs.An alternative parametric model with a finite set of parameter instantiations, but without the assumption that these instantiations are graph-preserving is considered in [98].
c) Model repair: The problem of model repair is related to parameter synthesis.In particular, for a Markov model and a refuted specification the problem is to transform the model such that the specification is satisfied.If repair consists of changing transition probabilities, the underlying model is parametric, where parameters are added to probabilities.The problem was first defined and solved either by a nonlinear program or parameter synthesis in [95].A greedy approach was given in [99] and efficient simulation-based methods are presented in [62].In addition, parametric models are used to rank patches in the repair of software [100].
d) Interval Markov chains: Instead of parametric transitions, interval MCs or MDPs feature intervals at their transitions [101], [102].These models do not allow for parameter dependencies, but verification is necessarily "robust" against all probabilities within the intervals, see for instance [103], where convex optimization is utilised, and [104], where efficient verification of multiple-objectives is introduced.In [105], [106], these models are extended to so-called parametric interval MCs, where interval bounds themselves are parametric.
e) Sensitivity analysis: Besides analysing in which regions the system behaves correctly w. r. t. the specification, it is often desirable to perform a sensitivity analysis [107], [108], i. e., to determine in which regions of the parameter space a small perturbation of the system leads to a relatively large change in the considered measure.In our setting, such an analysis can be conducted with little additional effort.Given a rational function for a measure of interest, its derivations w. r. t. all parameters can be easily computed.Passing the derivations with user-specified thresholds to the SMT solver then allows for finding parameter regions in which the system behaves robustly.Adding the safety constraints described earlier, the SMT solver can find regions that are both safe and robust.
f) Parameters with distributions: Rather than a model in which the parameter values are chosen from a set, they can be equipped with a distribution.The verification outcome are then confidence intervals rather than absolute guarantees.In [109], simulation based methods are used, whereas [110], [111] use statistical methods on a solution function.pMDPs with a distribution over the parameters are considered in [48].
g) Ensuring graph preservation: Checking graphpreservation is closely related to checking whether a welldefined point instantiation exists, which has an exponential runtime in the number of parameters [112].For parametric interval Markov chains, the question whether there exists a well-defined instantiation is referred to as consistency and received attention in [105], [113].
h) Robust strategies: Robust strategies for pMDPs, as mentioned in Remark 6, are considered in, among others, [114], [115].These and other variants of synthesis problems on pMDPs were compared in [116].A variant where parameters are not non-deterministically chosen, but governed by a prior over these parameters, has recently been considered [48].
i) Continuous time: Parametric CTMCs were first considered by [117].A method with similarities to parameter lifting has been proposed in [118].The method was improved in [119] and implemented in PRISM-PSY [120].A combination with sampling-based algorithms to find good parameter instantiations is explored in [121].Parameter synthesis with statistical guarantees has been explored in [122], [123].In [124], finding good parameter instantiations is considered by identifying subsets of parameters which have a strictly positive or negative influence on the property at hand.j) Complexity: For graph-preserving pMCs, many complexity results are collected in [53] and [125].In particular, the verification problem considered in this paper is known to be square-root-sum-hard and in the existential theory of the reals.Furthermore, [86] establishes connections to the computation of strategies in partially observable MDPs [84], a prominent model in AI.For pMDPs, so far only lower bounds (from pMCs) are known.This paper establishes membership in the class ETR (existential theory of the reals) via the encodings in Sect.V.

XII. CONCLUSION AND FUTURE WORK
This paper gives an extensive account of parameter synthesis for discrete-time Markov chain models.In particular, we considered three different variants of parameter synthesis questions.For each problem variant, we give an account of the available algorithms from the literature, together with several extension from our side.All algorithms are available in the open-source tool PROPhESY.For future work, we would like to develop methods which identify and exploit structural properties of many Markov chains and Markov decision processes, and to develop methods that handle pMDPs on regions that are not graph-preserving.

Fig. 1 .
Fig.1.A (a) biased and (b) parametric variant of Knuth-Yao's algorithm.In gray states an unfair coin is flipped with probability 2 /5 for 'heads'; for the unfair coin in the white states this probability equals 7 /10.On the right, the two biased coins have parametric probabilities.

2 . 5 (
Definition Instantiated pSG).For a pSG G = (S, V , s I , Act, P) and instantiation u of V , the instantiated pSG at u is given by G[u] = (S, s I , Act, P[u]) with P[u](s, α, s ) = P(s, α, s )[u] for all s, s ∈ S and α ∈ Act.The instantiation of the parametric reward function rew at u is rew[u] with rew[u](s, α) = rew(s, α)[u] for all s ∈ S, α ∈ Act.Instantiating pMDP M and pMC D at u is denoted by M[u] and D[u], respectively.
are non-strict.Rectangular regions are hyper-rectangles and a subclass of linear regions.A closed rectangular region R can be represented as R = × p∈V [a p , b p ] with parameter intervals [a p , b p ] described by the bounds a p and b p for all p ∈ V .For a region R, we refer to the bounds of parameter p by B R (p) = {a p , b p } and to the interval of parameter p by I R (p) = [a p , b p ].We may omit the subscript R, if it is clear from the context.For a rectangular region R, the size R equals p∈V b p − a p .Regions represent sets of instantiations G[u] of a pSG G.The notions of well-definedness and graph-preservation from Def. 7 trivially lift to regions: By the duality of |= a and |= d , a region is thus rejecting iff ∀u ∈ R. G, u |= ♣ ϕ.Note that this differs from G, R |= ♣ ϕ.Example 19.Reconsider the pMDP in Fig. 6(b), with R = [ 2 /5, 1 /2] × [ 2 /5, 1 /2] and ϕ = P > 4 /5 (♦{s 2 }).The corresponding solution functions are given in Example 14.

Example 22 .
Reconsider the pMC D from Fig. 5(c), and let R = [0, 1] × [0, 1], which is well defined but not graph preserving.Region R can be partitioned into 9 regions, see Fig. 7(a) where each dot, line segment, and the inner region are subregions of R. All subregions are graph-preserving on some sub-pMC of D. Consider, e.g., the line-region R = {u ∈ R | p[u] = 0}.The pMC D is not graph-preserving on R , as the transition s 0 p − → s 1 disappears when p = 0.However, R is graph-preserving on the sub-pMC D in Fig. 7(b), which is obtained from D by removing the transitions on the line-region p=0.

Fig. 12 .
Fig. 12.A pMC D and the substitution sub R (D).

in Fig. 12 (
a) and a region R = [ 1 /10, 4 /5] × [ 2 /5, 7 /10].The method creates the MDP in Fig. 12(b), where different types of arrows reflect different actions.The MDP is created by adding in each state two actions: One reflecting the lower bound of the parameter range, one reflecting the upper bound.Model checking on this MDP yields an maximal probability of 47 /60.

Example 39 . 4
Reconsider Example 38.From sub R (D) in Fig. 14(a), we can derive max σ∈Str Pr sub R (D) σ (♦T ) = 47 /60 and, by Theorem 6, D, R |= P ≤ 4 /5 (♦T ) follows.Despite the large region R, we establish a non-trivial upper bound on the reachability probability over all instantiations in R. If the over-approximation by region R is too coarse for a conclusive answer, region R can be refined, meaning that Algorithm Parameter lifting reachability(pMC D, T ⊆ S, region R, P ≤λ (♦T ) ) Construct sub R (D) if ∀σ ∈ Str sub R (D) |= P ≤λ (♦T ) then // via standard MDP model checking procedures return true else if ∀σ ∈ Str sub R (D) |= P >λ (♦T ) then // via standard MDP model checking procedures return false else return unknown smaller regions are considered.Intuitively, as more potential parameter values are excluded by reducing the region size, the actual choice of the parameter value has less impact on reachability probabilities.The smaller the region gets, the smaller the over-approximation: The optimal instantiation on the pMC D is over-approximated by some strategy on sub R (D).