Sequential Convex Programming for the Efficient Verification of Parametric MDPs
Murat Cubuktepe Nils Jansen Sebastian Junges Joost-Pieter Katoen Ivan Papusha Hasan A. Poonawala Ufuk Topcu


Theory of Scanning Tunneling Spectroscopy: from Kondo Impurities to Heavy Fermion Materials
    Dirk K. Morr
    December 30, 2023
===========================================================================================


Multi-objective verification problems of parametric Markov decision
processes under optimality criteria can be naturally expressed as
nonlinear programs. We observe that many of these computationally
demanding problems belong to the subclass of signomial programs.
This insight allows for a sequential optimization algorithm to
efficiently compute sound but possibly suboptimal solutions. Each stage
of this algorithm solves a geometric programming problem. These geometric
programs are obtained by convexifying the nonconvex constraints of the
original problem. 
Direct applications of the encodings as nonlinear programs are model repair and parameter synthesis.

We demonstrate the scalability and quality of our
approach by well-known benchmarks.





§ INTRODUCTION



We study the applicability of convex optimization to the formal verification of systems that exhibit randomness or stochastic uncertainties. Such systems are formally represented by so-called parametric Markov models. 

In fact, many real-world systems exhibit random behavior and stochastic uncertainties. One major example is in the field of robotics, where the presence of measurement noise or input disturbances requires special controller synthesis techniques <cit.> that achieve robustness of robot actions against uncertainties in the robot model and the environment.
On the other hand, formal verification offers methods for rigorously proving or disproving properties about the system behavior, and synthesizing strategies that satisfy these properties.
In particular, model checking <cit.> is a well-studied technique that provides guarantees on appropriate behavior for all possible events and scenarios.


Model checking can be applied to systems with stochastic uncertainties, including discrete-time Markov chains (MCs), Markov decision processes (MDPs), and their continuous-time counterparts <cit.>.
Probabilistic model checkers are able to verify reachability properties like “the probability of reaching a set of unsafe states is ≤ 10%” and expected costs properties like “the expected cost of reaching a goal state is ≤ 20.”
A rich set of properties, specified by linear- and branching-time logics, reduces to such properties <cit.>.
Tools like PRISM <cit.>, STORM <cit.>, and 
iscasMc <cit.> are probabilistic model checkers capable of handling a wide range of large-scale problems.

Key requirements for applying model checking are a reliable system model and formal specifications of desired or undesired behaviors. 
As a result, most approaches assume that models of the stochastic uncertainties are precisely given. For example, if a system description includes an environmental disturbance, the mean of that disturbance should be known before formal  statements are made about expected system behavior.
However, the desire to treat many applications where uncertainty measures (, faultiness, reliability, reaction rates, packet loss ratio) are not exactly known at design time 


gives rise to parametric probabilistic models <cit.>. Here, transition probabilities are expressed as functions over system parameters, , descriptions of uncertainties.
In this setting, parameter synthesis addresses the problem of computing parameter instantiations leading to satisfaction of system specifications. 
More precisely, parameters are mapped to concrete probabilities inducing the resulting instantiated model to satisfy  specifications.

A direct application is model repair <cit.>, where a concrete model (without parameters) is changed (repaired) such that specifications are satisfied.





Dedicated tools like PARAM <cit.>, PRISM <cit.>, or  <cit.> compute rational functions over parameters that express reachability probabilities or expected costs in a parametric Markov chain (pMC). These optimized tools work with millions of states but are restricted to a few parameters, as the necessary computation of greatest common divisors does not scale well with the number of parameters.
	Moreover, the resulting functions are inherently nonlinear and often of high degree. Evaluation by an SMT solver over nonlinear arithmetic such as Z3 <cit.> suffers from the fact that the solving procedures are exponential in the degree of polynomials and the number of variables. 


This paper takes an alternative perspective. 

We discuss a general nonlinear programming formulation for the verification of parametric Markov decision processes (pMDPs).
 The powerful modeling capabilities of nonlinear programs (NLPs) enable incorporating multi-objective properties and penalties on the parameters of the pMDP.
However, because of their generality, solving NLPs to find a global optimum is difficult. Even feasible solutions (satisfying the constraints) cannot always be computed efficiently <cit.>. 
In contrast, for the class of NLPs called convex optimization problems, efficient methods to compute feasible solutions and global optima even for large-scale problems are available <cit.>. 

We therefore propose a novel automated method of utilizing convex optimization for pMDPs.


Many NLP problems for pMDPs belong to the class of signomial programs (SGPs), a certain class of nonconvex optimization problems.
For instance, all benchmarks available at the PARAM–webpage <cit.> belong to this class.




Restricting the general pMDP problem accordingly yields a direct and efficient synthesis method—formulated as an NLP—for a large class of pMDP problems. 

 We list the two main technical results of this paper:

	
  * We relax nonconvex constraints in SGPs and apply a simple transformation to the parameter functions. The resulting programs are geometric programs (GPs) <cit.>, a class of convex programs. We show that a solution to the relaxed GP induces feasibility (satisfaction of all specifications) in the original pMDP problem. Note that solving GPs is polynomial in the number of variables. 
	
  * Given an initial feasible solution, we use a technique called sequential convex programming <cit.> to improve a signomial objective. This local optimization method for nonconvex problems leverages convex optimization by solving a sequence of convex approximations (GPs) of the original SGP. 

Sequential convex programming is known to efficiently find a feasible solution with
good, though not necessarily globally optimal, objective values <cit.>.
We initialize the sequence with a feasible solution (obtained from the GP) of the original problem and compute a trust region. Inside this region, the optimal value of the approximation of the SGP is at least as good as the objective value at the feasible solution of the GP. 
The optimal solution of the approximation is then the initial point of the next iteration with a new trust region.
This procedure is iterated to approximate a local optimum of the original problem. 



Utilizing our results, we discuss the concrete problems of parameter synthesis and model repair for multiple specifications for pMDPs.
Experimental results with a prototype implementation show the applicability of our optimization methods to benchmarks of up to 10^5 states.

As solving GPs is polynomial in the number of variables, our approaches are relatively insensitive to the number of parameters in pMDPs. This is an improvement over state-of-the-art approaches that leverage SMT, which—for our class of problems—scale exponentially in variables and the degree of polynomials. 
This is substantiated by our experiments.




  
Related work.
Several approaches exist for pMCs <cit.> while the number of approaches for pMDPs <cit.> is limited.
Ceska et al. <cit.> synthesize rate parameters in stochastic biochemical networks. 
Multi-objective model checking of non-parametric MDPs <cit.> is a convex problem <cit.>. 
Bortolussi et al. <cit.> developed a Bayesian statistical algorithm for properties on stochastic population models.

Convex uncertainties in MDPs without parameter dependencies are discussed in <cit.>.

Parametric probabilistic models are used to rank patches in the repair of software <cit.> and to compute perturbation bounds <cit.>. 






















































§ PRELIMINARIES
















A probability distribution over a finite or countably infinite set  is a function → with ∑_∈()=1. 

The set of all distributions on  is denoted by ().








  Let V={x_1,…,x_n} be a finite set of strictly positive real-valued variables.
  A monomial over V is an expression of the form
  
    =c· x_1^a_1⋯ x_n^a_n ,

  where c∈_>0 is a positive coefficient, and a_i∈ are exponents for 1≤ i≤ n. 
  A posynomial over V is a sum of one or more monomials:
  
    =∑_k=1^K c_k· x_1^a_1k⋯ x_n^a_nk .

  If c_k is allowed to be a negative real number for any 1≤ k≤ K, then the expression (<ref>) is a signomial. The sets of all monomials, posynomials, and signomials over V are denoted by [V], [V], and [V], respectively.
 

This definition of monomials differs from the standard algebraic definition where exponents are positive integers with no restriction on the coefficient sign. A sum of monomials is then called a polynomial. Our definitions are consistent with <cit.>.




  
  For a set of real-valued variables V, a valuation u over V is a function u V →.
  The set of all valuations over V is ^V.

 Applying valuation u to monomial  over V yields a real number [u]∈ by replacing each occurrence of variables x∈ V in  by u(x); the procedure is analogous for posynomials and signomials using standard arithmetic operations.




 


A parametric Markov decision process (pMDP) is a tuple  with a finite set S of states, an initial state ∈ S, a finite set  of actions, a finite set of real-valued variables , and a transition function S ×× S →[] satisfying for all s∈ S(s) ≠∅,  where (s) = {∈|∃ s'∈ S. (s, , s') ≠ 0}.
If for all s∈ S it holds that |(s)| = 1,  is called a parametric discrete-time Markov chain (pMC)
.





(s) is the set of enabled actions at state s; as (s) ≠∅, there are no deadlock states.




Costs are defined using a state–action cost function S ×→_≥ 0.



Largely due to algorithmic reasons, the transition probabilities in the literature <cit.> are polynomials or rational functions, , fractions of polynomials.  
 Our restriction to signomials is realistic; all benchmarks from the PARAM–webpage <cit.> contain only signomial transition probabilities. 



A pMDP  is a Markov decision process (MDP) if the transition function is a valid probability distribution, , S ×× S → [0,1] and ∑_s'∈ S(s,,s') = 1 for all s ∈ S ∈(s). 


Analogously, a Markov chain (MC) is a special class of a pMC; a model is parameter-free if all  probabilities are constant. Applying a valuation u to a pMDP, denoted [u], replaces each signomial f in  by f[u]; we call [u] the instantiation of  at u.
The application of u is to replace the transition function f by the probability f[u]. 
A valuation u is well-defined for  if the replacement yields probability distributions at all states; the resulting model [u] is an MDP or an MC.

[pMC]
Consider a variant of the Knuth–Yao model of a die <cit.>, where a six-sided die is simulated by successive coin flips. 
We alternate flipping two biased coins, which result in heads with probabilities defined by the monomials p and q, respectively. Consequently, the probability for tails is given by the signomials 1-p and 1-q, respectively. The corresponding pMC is depicted in Fig. <ref>; and the instantiated MC for p = 0.4 and q = 0.7 is given in Fig. <ref>. Note that we omit actions, as the model is deterministic.






In order to define a probability measure and expected cost on MDPs, nondeterministic choices are resolved by so-called schedulers.

For practical reasons we restrict ourselves to memoryless schedulers; details can be found in <cit.>.




	A (randomized) scheduler for an MDP  is a function S→() such that σ(s)(α) > 0 implies α∈(s).

	The set of all schedulers over  is denoted by ^.



Applying a scheduler to an MDP yields a so-called induced Markov chain.



	Let MDP  and scheduler ∈^. The MC induced by  and  is  ^=(S,,,^) where for all s,s'∈ S,
	
    ^(s,s')=∑_α∈(s)(s)()·(s,α,s').
 




We consider reachability properties and expected cost properties.
For MC  with states S, let sT denote the probability of reaching a set of target states T ⊆ S from  state s∈ S; simply [] denotes the probability for initial state .
We use the standard probability measure as in <cit.>.

For threshold λ∈ [0,1], the reachability property asserting that a target state is to be reached with probability at most λ is denoted =.
The property is satisfied by , written , iff []≤λ.


The cost of a path through MC  until a set of goal states G⊆ S is the sum of action costs visited along the path. The expected cost of a finite path is the product of its probability and its cost.
For G  = 1, the expected cost of reaching G is the sum of expected costs of all paths leading to G.

An expected cost property κG is satisfied if the expected cost of reaching T is bounded by a threshold κ∈.
Formal definitions are given in e.g., <cit.>.

If multiple specifications φ_1,…,φ_q are given, which are either reachability properties or expected cost properties of the aforementioned forms, we write the satisfaction of all specifications φ_1,…,φ_q for an MC  as φ_1…φ_q. 

An MDP  satisfies the specifications φ_1,…,φ_q, iff for all schedulers ∈^ it holds that ^φ_1∧…∧φ_q. The verification of multiple specifications is also referred to as multi-objective model checking <cit.>.
We are also interested in the so-called scheduler synthesis problem, where the aim is to find a scheduler  such that the specifications are satisfied (although other schedulers may not satisfy the specifications).






















































§ NONLINEAR PROGRAMMING FOR PMDPS

In this section we formally state a general pMDP parameter synthesis problem and describe how
it can be formulated using nonlinear programming.



 §.§ Formal problem statement
 






-4

Given a pMDP , specifications
φ_1,…,φ_q that are either probabilistic reachability
properties or expected cost properties, and an objective function
f→ over the variables V, compute a well-defined
valuation u∈^V for , and a (randomized) scheduler ∈^
such that the following conditions hold:

	
  *  Feasibility:
		the Markov chain ^[u] induced by scheduler  and
		instantiated by valuation u satisfies the specifications, , ^[u]φ_1 ∧…∧φ_q.
	
	
  *  Optimality:
		the objective f is minimized.





Intuitively, we wish to compute a parameter valuation and a scheduler such
that all specifications are satisfied, and the objective is globally minimized.


We refer to a valuation–scheduler pair (u,) that satisfies
condition (<ref>), , only guarantees satisfaction of the
specifications but does not necessarily minimize the objective f, as a
feasible solution to the pMDP synthesis problem. If both
(<ref>) and (<ref>) are satisfied, the pair
is an optimal solution to the pMDP synthesis problem.































 §.§ Nonlinear encoding

We now provide an NLP encoding of Problem <ref>. A general NLP over 
a set of real-valued variables  can be written as

    minimize      f
    subject to   
    ∀ i.  1≤ i≤ m 	      g_i ≤ 0,
    ∀ j.  1≤ i≤ p 	      h_j = 0,

where f, g_i, and h_j are arbitrary functions over , and m and p are the 
number of inequality and equality constraints of the program respectively. Tools like
IPOPT <cit.> solve small instances of such problems.



Consider a
pMDP  with specifications φ_1= and φ_2=κG. We will discuss how additional specifications of either type can be encoded.
The set = ∪ W of variables of the NLP consists of
the variables  that occur in the pMDP as well as a set W of additional variables:

	
  * {^s,α| s ∈ S, ∈(s) },
		which define the randomized scheduler  by (s)()=^s,.
		
	
  * { p_s | s ∈ S }, 
		where p_s is the probability of reaching the target set 
		T⊆ S from state s under scheduler , and
	
  * { c_s | s ∈ S }, where c_s is the expected cost to reach G⊆ S from s under .



 A valuation over  consists of a valuation u∈^V over the
pMDP variables and a valuation w∈^W over the additional variables.

    f 
    
          p_≤λ,				
          c_≤κ,				
    ∀ s∈ S.	     ∑_∈(s)^s,=1, 
    ∀ s∈ S  ∀∈(s).       0 ≤^s,≤ 1,				
    ∀ s∈ S  ∀∈(s).	      ∑_s'∈ S(s,,s')=1,	
    ∀ s, s'∈ S  ∀∈(s).	       0 ≤(s,,s') ≤ 1,	
    ∀ s∈ T.	      p_s=1,															
    ∀ s∈ S∖ T.       p_s=∑_∈(s)σ^s,·∑_s'∈ S(s,,s')· p_s', 
    ∀ s∈ G.	       c_s=0,															
    ∀ s∈ S∖ G.	      c_s= ∑_∈(s)σ^s,·(c(s,) + ∑_s'∈ S (s,,s') · c_s').

The NLP (<ref>)–(<ref>) encodes Problem <ref> in the following way.
The objective function f in (<ref>) is any real-valued function over the variables . 
The constraints (<ref>) and (<ref>) encode the
specifications φ_1 and φ_2, respectively.
The constraints (<ref>)–(<ref>)
ensure that the scheduler obtained is well-defined by requiring that the
scheduler variables at each state sum to unity. 
Similarly, the constraints
(<ref>)–(<ref>) ensure that
for all states, parameters from  are instantiated such that
probabilities sum up to one. 
(These constraints are included if not all probabilities at a state are constant.)
The probability of reaching the target for all states in the target set is
set to one using (<ref>).
The reachability probabilities in each state 
depend on the reachability of the successor states and the transition
probabilities to those states through (<ref>).
Analogously to the reachability probabilities, the cost for each goal state G⊆ S
must be zero, thereby precluding the collection of infinite cost at
absorbing states, as enforced by (<ref>).
Finally, the expected cost for all states except target states is given by
the equation (<ref>), where according to the
strategy  the cost of each action is added to the expected cost of
the successors. 


	
	

	
	
We can readily extend the NLP to include more specifications. If
another reachability property φ'=λ'T' is given, we add the set of probability variables { p'_s |
s ∈ S} to W, and duplicate the 
constraints (<ref>)–(<ref>) accordingly.




To ensure satisfaction of φ', we also add the constraint
p'_≤λ'.
The procedure is similar for additional expected cost properties. 
By construction, we have the following result relating the NLP encoding and Problem <ref>.

	
	The NLP (<ref>)–(<ref>) is sound and complete with respect to Problem <ref>.


We refer to soundness in the sense that each variable assignment that satisfies
the constraints induces a scheduler and a valuation of parameters such that a
feasible solution of the problem is induced. Moreover, any optimal solution to
the NLP induces an optimal solution of the problem. Completeness means that all
possible solutions of the problem can be encoded by this NLP; while
unsatisfiability means that no such solution exists, making the problem
infeasible.




  
Signomial programs. By Def. <ref> and <ref>, all constraints in the NLP consist of signomial functions.
A special class of NLPs known as signomial programs (SGPs) is of the form (<ref>)–(<ref>) where f, g_i and h_j are signomials over , see Def. <ref>. Therefore, we observe that the NLP (<ref>)–(<ref>) is an SGP. We will refer to the NLP as an SGP in what follows.

SGPs with equality constraints consisting of functions that are not affine are not convex in general. 
In particular, the SGP (<ref>)–(<ref>) is not necessarily convex. Consider a simple pMC only having transition probabilities of the form p and 1-p, as in Example <ref>. The function in the equality
constraint (<ref>) of the corresponding SGP encoding is not affine in
parameter p and the probability variable p_s for some state s∈ S.
More generally, the equality constraints
(<ref>),
(<ref>), and
(<ref>)
involving  are not necessarily affine, and thus the SGP may not be a convex program <cit.>.


Whereas for convex programs global optimal solutions can be
found efficiently <cit.>, such guarantees are
not given for SGPs. 
However, we can efficiently obtain local optimal solutions for SGPs in our setting, as shown in the following sections.








§ CONVEXIFICATION


We investigate how to transform the SGP (<ref>)–(<ref>) into a convex program by relaxing equality constraints and a lifting of variables of the SGP. 


A certain subclass of SGPs called geometric programs (GPs) can be transformed into convex programs <cit.> and solved efficiently.
A GP is an SGP of the form (<ref>)–(<ref>) where f,g_i∈[] and h_j∈[].
We will refer to a constraint with posynomial or monomial function as a posynomial or monomial constraint, respectively.

















 §.§ Transformation and relaxation of equality constraints

As discussed before, the SGP (<ref>)–(<ref>) is not convex because of the presence of non-affine equality constraints. 
First observe the following transformation <cit.>:

    f≤ h⟺f/h≤ 1,

for f∈[] and h∈ []. Note that monomials are strictly positive (Def. <ref>). This (division-)transformation of f≤ h yields a posynomial inequality constraint.

We relax all equality constraints of SGP (<ref>)–(<ref>) that are not monomials to inequalities, then we apply the division-transformation wherever possible. Constraints (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) are transformed to

    p_/λ≤ 1,
         c_/κ≤ 1,
    ∀ s∈ S.	     ∑_∈(s)^s,≤ 1,
    ∀ s∈ S ∀∈(s).	      ∑_s'∈ S(s,,s')≤ 1,
    ∀ s∈ S∖ T.	     ∑_∈(s)σ^s,·∑_s'∈ S(s,,s')· p_s'/p_s≤ 1,
    ∀ s∈ S∖ G.	     ∑_∈(s)σ^s,·(c(s,) + ∑_s'∈ S (s,,s') · c_s')/c_s≤ 1 .


These constraints are not necessarily posynomial inequality constraints because (as in Def. <ref>) we allow signomial expressions in the transition probability function . Therefore, replacing (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) in the SGP with (<ref>)–(<ref>) does not by itself convert the SGP to a GP.



 §.§ Convexification by lifting





The relaxed equality constraints (<ref>)–(<ref>) involving  are signomial, rather than posynomial, because the parameters enter Problem <ref> in signomial form. Specifically, consider the relaxed equality constraint (<ref>) at s_0 in Example <ref>,

    p· p_s_1 + (1-p)· p_s_2/p_s_0≤ 1.

The term (1-p)· p_s_2 is signomial in p and p_s_2. We lift by introducing a new variable p̅ = 1-p, and rewrite (<ref>) as a posynomial inequality constraint and an equality constraint in the lifted variables:

    p· p_s_1 + p̅· p_s_2/p_s_0≤ 1,
    	  p̅ = 1-p.



We relax the (non-monomial) equality constraint to p + p̅≤ 1. 

More generally, we restrict the way parameters occur in  as follows. Refer to Fig. <ref>. For every state s∈ S and every action ∈(s) we require that there exists at most one state s̅∈ S such that (s,,s̅)∈[] and (s,,s')∈[] for all s'∈ S∖{s̅}. In particular, we require that
	
    (s,,s̅)_∈[] = 1 - ∑_s'∈ S∖{s̅}(s,,s')_∈[] .

This requirement is met by all benchmarks available at the PARAM–webpage <cit.>.
In general, we lift by introducing a new variable p̅_s,,s̅=(s,,s̅) for each such state s∈ S; refer to Fig. <ref>.
We denote this set of lifting variables by L. Lifting as explained above then creates a new transition probability function  where for every s,s'∈ S and ∈ we have (s,,s')∈[∪ L]. 

We call the set of constraints obtained through transformation, relaxation, and lifting of every constraint of the SGP (<ref>)–(<ref>) as shown above the convexified constraints. 
Any posynomial objective subject to the convexified constraints forms by construction a GP over the pMDP parameters , the SGP additional variables W, and the lifting variables L. 



 §.§ Tightening the constraints

A solution of the GP as obtained in the previous section does not have a direct relation to the original SGP (<ref>)–(<ref>).
In particular, a solution to the GP may not have the relaxed constraints satisfied with equality. 
For (<ref>) and (<ref>), the induced parameter valuation and the scheduler are not well-defined, , the probabilities may not sum to one. In particular, this makes it trivial to fulfill the GP, just assign 0 to all schedulers
We need to relate the relaxed and lifted GP to Problem <ref>. By defining a regularization function F over all parameter and scheduler variables, we ensure that the constraints are satisfied with equality; enforcing well-defined probability distributions.













    F = ∑_p∈1/p + ∑_p̅∈ L1/p̅ + ∑_s∈ S,∈(s)1/σ_s, .

The function F is monotone in all its variables. We discard the original objective f in (<ref>) and form a GP with the regularization objective F (<ref>):

    minimize       F
    subject to
         p_/λ≤ 1,
         c_/κ≤ 1,
    ∀ s∈ S.	     ∑_∈(s)^s,≤ 1,
    ∀ s∈ S ∀∈(s).      ^s,≤ 1			,
    ∀ s∈ S ∀∈(s).	      ∑_s'∈ S(s,,s')≤ 1,
    ∀ s, s'∈ S ∀∈(s).	      (s,,s') ≤ 1,	
    ∀ s∈ T.	      p_s=1,															
    ∀ s∈ S∖ T.	     ∑_∈(s)σ^s,·∑_s'∈ S(s,,s')· p_s'/p_s≤ 1,
    ∀ s∈ S∖ G.	     ∑_∈(s)σ^s,·(c(s,) + ∑_s'∈ S(s,,s') · c_s')/c_s≤ 1 .



Since the objective F (<ref>) and the inequality constraints (<ref>) and (<ref>) are monotone in , L, and the scheduler variables, each optimal solution for a feasible problem satisfies them with equality. We obtain a well-defined scheduler  and a valuation u as in Problem 1. Note that variables from (<ref>) are explicitly excluded from the GP by treating them as constants.

The reachability probability constraints (<ref>) and cost constraints (<ref>) need not be satisfied with equality. However, (<ref>) is equivalent to 

    p_s≥∑_∈(s)σ^s,·∑_s'∈ S(s,,s')· p_s'

for all s∈ S∖ T and ∈.
The probability variables p_s are assigned upper bounds on the actual probability to reach the target states T under scheduler  and valuation u. Put differently, the p_s variables cannot be assigned values that are lower than the actual probability; ensuring that  and u induce satisfaction of the specification given by (<ref>) if the problem is feasible and  and u are well-defined. An analogous reasoning applies to the expected cost computation (<ref>). 
A solution consisting of a scheduler or valuation that are not well-defined occurs only if Problem <ref> itself is infeasible. Identifying that such a solution has been obtained is easy.

These facts allow us to state the main result of this section.

	
	A solution of the GP (<ref>)–(<ref>) inducing well-defined scheduler  and valuation u is a feasible solution to Problem 1.

Note that the actual probabilities induced by  and u for the given pMDP  are given by the MC ^[u] induced by  and instantiated by u.
Since all variables are implicitly positive in a GP, no transition probability function will be instantiated to probability zero. 

The case of a scheduler variable being zero to induce the optimum can be excluded by a previous graph analysis.











§ SEQUENTIAL GEOMETRIC PROGRAMMING


We showed how to efficiently obtain a feasible solution

 for Problem <ref> by solving GP (<ref>)–(<ref>). 









We propose a sequential convex programming trust-region method to compute a local optimum of the SGP (<ref>)–(<ref>), following <cit.>, solving a sequence of GPs. We obtain each GP by replacing signomial functions in equality constraints of the SGP (<ref>)–(<ref>) with monomial approximations of the functions. 

























Given a posynomial f ∈[], variables = {x_1,…,x_n}, and a valuation u ∈^, a monomial approximation f̂∈[] for f near u is 

    ∀ i. 1 ≤ i ≤ n
    	  f̂ = f[u]
    	∏_i=1^n(x_iu(x_i))^a_i,
    	  where 
    	a_i=u(x_i)f[u]∂ f∂ x_i[u].
	

Intuitively, we compute a linearization f̂ of f∈[] around a fixed valuation u.



We enforce the fidelity of monomial approximation f̂ of f ∈[] by restricting valuations to remain within a set known as trust region. We define the following constraints on the variables  with t > 1 determining the size of the trust region:

    ∀ i. 1 ≤ i ≤ n
    	   (1/t)· u(x_i) ≤ x_i ≤ t· u(x_i)










For a given valuation u, we approximate the SGP (<ref>)–(<ref>) to obtain a local GP as follows. 
First, we apply a lifting procedure (Section <ref>) to the SGP ensuring that all constraints consist of posynomial functions. The thus obtained posynomial inequality constraints are included in the local GP.
After replacing posynomials in every equality constraint by their monomial approximations near u, the resulting monomial equality constraints are also included.
Finally, we add trust region constraints (<ref>) for scheduler and parameter variables. The objective function is the same as for the SGP.
The optimal solution of the local GP is not necessarily a feasible solution to the SGP. 
Therefore, we first normalize the scheduler and parameter values to obtain well-defined probability distributions. These normalized values are used to compute precise probabilities and expected cost using PRISM. The steps above provide a feasible solution of the SGP.

We use such approximations to obtain a sequence of feasible solutions to the SGP approaching a local optimum of the SGP. 
First, we compute a feasible solution u^(0) for Problem <ref> (Section <ref>), 


forming the initial point of a sequence of solutions u^(0),…, u^(N), N ∈. 
The solution u^(k) for 0≤ k≤ N is obtained from a local GP defined using u^(k-1) as explained above.







The parameter t for each iteration k is determined based on its value for the previous iteration, and the ratio of f[u^(k-1)] to f[u^(k-2)], where f is the objective function in (<ref>). The iterations are stopped when | f[u^(k)] - f[u^(k-1)] | <  ϵ. Intuitively, ϵ defines the required improvement on the objective value for each iteration; once there is not enough improvement the process terminates.


































































§ APPLICATIONS

We discuss two applications and their restrictions for the general SGP (<ref>)–(<ref>).




  
Model repair.

For MC  and specification φ with φ, the model repair problem <cit.> is to transform  to ' such that 'φ. The transformation involves a change of transition probabilities. Additionally, a cost function measures the change of probabilities. The natural underlying model is a pMC where parameters are added to probabilities. The cost function is minimized subject to constraints that induce satisfaction of φ. In <cit.>, the problem is given as NLP. Heuristic <cit.> and simulation-based methods <cit.> (for MDPs) were presented.

Leveraging our results, one can readily encode model repair problems for MDPs, multiple objectives, and restrictions on probability or cost changes directly as NLPs. The encoding as in <cit.> is  handled by our method in Section <ref> as it involves signomial constraints. We now propose a more efficient approach, which encodes the change of probabilities using monomial functions.

Consider an MDP  and specifications φ_1,…,φ_q with φ_1…φ_q. For each probability (s,,s')=a∈ that may be changed, introduce a parameter p, forming the parameter set . We define a parametric transition probability function by '(s,,s')=p· a∈[V]. The quadratic cost function is for instance  f=∑_p∈ V p^2∈[V]. 

By minimizing the sum of squares of the parameters (with some regularization), the change of probabilities is minimized.

By incorporating these modifications into SGP (<ref>)–(<ref>), our approach is directly applicable. Either we restrict the cost function f to an upper bound, and efficiently solve a feasibility problem (Section <ref>), or we compute a local minimum of the cost function (Section <ref>).

In contrast to <cit.>, our approach works for MDPs and has an efficient solution. While <cit.> uses fast simulation techniques, we can directly incorporate multiple objectives and restrictions on the results while offering an efficient numerical solution of the problem.



  
Parameter space partitioning.
For pMDPs, 

tools like PRISM <cit.> or PROPhESY <cit.> aim at partitioning the parameter space into regions with respect to a specification.
A parameter region is given by a convex polytope defined by linear inequalities over the parameters, restricting valuations to a region. Now, for pMDP  a region is safe regarding a specification φ, if no valuation u inside this region and no scheduler  induce ^[u]φ. Vice versa, a region is unsafe, if there is no valuation and scheduler such that the specification is satisfied. 
In <cit.>, this certification is performed using SMT solving.
More efficiency is achieved by using an approximation method <cit.>.

Certifying regions to be unsafe is directly possible using our approach.

Assume pMDP , specifications φ_1,…,φ_q, and a region candidate defined by a set of linear inequalities.

We incorporate the inequalities in the NLP (<ref>)–(<ref>). If the feasibility problem (Section <ref>) has no solution, the region is unsafe. This yields the first efficient numerical method for this problem of which we are aware.
Proving that a region is safe is more involved. Given one specification φ=, we maximize the probability to reach T. If this probability is at most λ, the region is safe. For using our method from Section <ref>, one needs domain specific knowledge to show that a local optimum is a global optimum.









§ EXPERIMENTS


We implemented a prototype using the Python interfaces of the probabilistic model checker STORM <cit.> and the optimization solver MOSEK <cit.>.
 All experiments were run on a 2.6 GHz machine with 32 GB RAM. 
 We used PRISM <cit.> to correct approximation errors as explained before. 

We evaluated our approaches using mainly examples from the PARAM–webpage <cit.> and from PRISM <cit.>.
We considered several parametric instances of the Bounded Retransmission Protocol (BRP) <cit.>,  NAND Multiplexing <cit.>, and the Consensus protocol (CONS) <cit.>. For BRP, we have a pMC and a pMDP version, NAND is a pMC, and CONS is a pMDP.
For obtaining feasibility solutions, we compare to the SMT solver Z3 <cit.>. For additional optimality criteria, there is no comparison to another tool possible as IPOPT <cit.> already fails for the smallest instances we consider.

Fig. <ref> states for each benchmark instance the number of states (#states) and the number of parameters (#par). We defined two specifications consisting of a expected cost property () and a reachability property (). For some benchmarks, we also maximized the probability to reach a set of “good states” (*). We list the times taken by MOSEK; for optimality problems we also list the times PRISM took to compute precise probabilities or costs (Section <ref>). For feasibility problems we list the times of Z3. The timeout (TO) is 90 minutes.










We observe that both for feasibility with optimality criteria we can handle most benchmarks of up to 10^5 states within the timeout, while we ran into a timeout for CONS. The number of iterations N in the sequential convex programming is less than 12 for all benchmarks with ϵ=10^-3.
As expected, simply solving feasibility problems is faster by at least one order of magnitude. Raising the number of parameters from 2 to 4 for BRP does not cause a major performance hit, contrary to existing tools. For all benchmarks except NAND, Z3 only delivered results for the smallest instances within the timeout. 


To demonstrate the insensitivity of our approach to the number of parameters, we considered a pMC of rolling multiple  Knuth–Yao dice with 156 states, 522 transitions and considered instances with up to 8 different parameters. The timeout is 100 seconds.
In Fig. <ref> we compare our encoding in MOSEK for this benchmark to the mere computation of a rational function using PROPhESY <cit.> and again to Z3. PROPhESY already runs into a timeout for 4 parameters[Due to the costly computation of greatest common divisors employed in .].

Z3 needs around 15 seconds for most of the tests. Using GPs with MOSEK proves far more efficient  as it needs less than one second for all instances.

In addition, we test model repair (Section <ref>) on a BRP instance with 17415 states for φ=0.9T. The initial parameter instantiation violates φ. We performed model repair towards satisfaction of φ. The probability of reaching T results in 0.79 and the associated cost is 0.013. The computation time is 21.93 seconds. We compare our result to an implementation of <cit.>, where the probability of reaching T is 0.58 and the associated cost is 0.064. However, the time for the simulation-based method is only 2.4 seconds, highlighting the expected trade-off between optimality and computation times for the two methods.
	
Finally, we encode model repair for the small pMC from Example <ref> in IPOPT, see <cit.>. For ψ=0.125T where T represents the outcome of the die being 2, the initial instantiation induces probability 1/6. With our method, the probability of satisfying ψ is 0.1248 and the cost is 0.0128. With IPOPT, the probability is 0.125 with cost 0.1025, showing that our result is nearly optimal.



















§ CONCLUSION AND FUTURE WORK

We presented a way to use convex optimization in the field of parameter synthesis for parametric Markov models. Using our results, many NLP encodings of related problems now have a direct and efficient solution.

Future work will concern the integration of these methods into mature tools like PRISM or PROPhESY to enable large-scale benchmarking by state space reduction techniques and advanced data structures. Moreover, we will explore extensions to richer models like continuous-time Markov chains <cit.>.


plainyr











