Axel Legay

Sean Sedwards

Louis-Marie Traonouez

Scalable Verification of Markov Decision Processes

Markov decision processes (MDP) are useful to model concurrent process optimisation problems, but verifying them with numerical methods is often intractable. Existing approximative approaches do not scale well and are limited to memoryless schedulers. Here we present the basis of scalable verification for MDPSs, using an O(1) memory representation of history-dependent schedulers. We thus facilitate scalable learning techniques and the use of massively parallel verification.

Introduction

Markov decision processes (MDP) describe systems that interleave nondeterministic actions and probabilistic transitions, possibly with rewards or costs assigned to the actions [START_REF] Bellman | Dynamic Programming[END_REF][START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. This model has proved useful in many real optimisation problems and may also be used to represent concurrent probabilistic programs (see, e.g., [START_REF] Bianco | Model checking of probabilistic and nondeterministic systems[END_REF][START_REF] Baier | Principles of model checking[END_REF]). Such models comprise probabilistic subsystems whose transitions depend on the states of the other subsystems, while the order in which concurrently enabled transitions execute is nondeterministic. This order may radically affect the expected reward or the probability that a system will satisfy a given property. It is therefore useful to calculate the upper and lower bounds of these quantities.

Fig. 1 shows a typical fragment of an MDP. Referring in parentheses to the labels in the figure, the execution semantics are as follows. In a given state (s 0 ), an action (a 1 , a 2 , . . .) is chosen nondeterministically to select a distribution of probabilistic transitions (p 1 , p 2 , . . . or p 3 , p 4 , etc.). A probabilistic choice is then made to select the next state (s 1 , s 2 , s 3 , s 4 , . . .). To each of the actions may be associated a reward (r 1 , r 2 , . . .), allowing values to be assigned to sequences of actions.

To calculate the expected total reward or the expected probability of a sequence of states, it is necessary to define how the nondeterminism in the MDP will be resolved. In the literature this is often called a strategy, a policy or an adversary. Here we use the term scheduler and focus on MDPs in the context of model checking concurrent probabilistic systems. Model checking is an automatic technique to verify that a system satisfies a property specified in temporal logic [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]. Probabilistic model checking quantifies the probability that a probabilistic system will satisfy a property [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF]. Classic analysis of MDPs is concerned with finding schedulers that maximise or minimise rewards [START_REF] Bellman | Dynamic Programming[END_REF][START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. The classic verification algorithms for MDPs are concerned with finding schedulers that maximise or minimise the probability of a property, or deciding the existence of a scheduler that ensures the probability of a property is within some bound [START_REF] Bianco | Model checking of probabilistic and nondeterministic systems[END_REF]. Our techniques can be easily extended to include rewards, but in this work we focus on probabilities and leave rewards for future consideration.

Schedulers and State Explosion

The classic algorithms to solve MDPs are policy iteration and value iteration [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. Model checking algorithms for MDPs may use value iteration applied to probabilities [2, Ch. 10] or solve the same problem using linear programming [START_REF] Bianco | Model checking of probabilistic and nondeterministic systems[END_REF].

All consider history-dependent schedulers. Given an MDP with set of actions A, having a set of states S that induces a set of sequences of states Ω = S + , a history-dependent (general) scheduler is a function S : Ω → A. A memoryless scheduler is a function M : S → A. Intuitively, at each state in the course of an execution, a history-dependent scheduler (S) chooses an action based on the sequence of previous states, while a memoryless scheduler (M) chooses an action based only on the current state. History-dependent schedulers therefore include memoryless schedulers. Fig. 2 illustrates a simple MDP for which memoryless and history-dependent schedulers give different optima for logical property X(ψ∧XG t ¬ψ) when p 1 = p 2 and t > 0. The property makes use of the temporal operators next (X) and globally (G). Intuitively, the property states that on the next step ψ will be true and, on the step after that, ¬ψ will be remain true for t+1 time steps. The property is satisfied by the sequence of states s 0 s 1 s 0 s 0 • • •. If p 1 > p 2 , the maximum probability for s 0 s 1 is achieved with action a 2 , while the maximum probability for s 0 s 0 is achieved with action a 1 . Given that both transitions start in the same state, a memoryless scheduler will not achieve the maximum probability achievable with a history-dependent scheduler.

The principal challenge of finding optimal schedulers is what has been described as the 'curse of dimensionality' [START_REF] Bellman | Dynamic Programming[END_REF] and the 'state explosion problem' [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]: the number of states of a system increases exponentially with respect to the number of interacting components and state variables. This phenomenon has led to the design of sampling algorithms that find 'near optimal' schedulers to maximise rewards in discounted MDPs. Probably the best known is the Kearns algorithm [START_REF] Kearns | A sparse sampling algorithm for nearoptimal planning in large Markov decision processes[END_REF], which we briefly review in Section 2.

The state explosion problem of model checking applied to purely probabilistic systems has been well addressed by statistical model checking (SMC) [START_REF] Younes | Probabilistic verification of discrete event systems using acceptance sampling[END_REF]. SMC uses an executable model to approximate the probability that a system satisfies a specified property by the proportion of simulation traces that individually satisfy it. SMC algorithms work by constructing an automaton to accept only traces that satisfy the property. This automaton may then be used to estimate the probability of the property or to decide an hypothesis about the probability. Typically, the probability of property ϕ is estimated by

1 N N i=1 1(ω i |= ϕ),
where ω 1 , . . . , ω N are N independently generated simulation traces and 1(•) is an indicator function that corresponds to the output of the automaton: it returns 1 if the trace is accepted and 0 if it is not. N is chosen a priori to give the required statistical confidence (e.g., using a Chernoff bound [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF], see Section 4.2). Sequential hypothesis tests (e.g., Wald's sequential probability ratio test [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF], see Section 4.1) do not define N a priori, but generate simulation traces until an hypothesis can be accepted or rejected with specified confidence. The state space of the system is not constructed explicitly-states are generated on the fly during simulation-hence SMC is efficient for large, possibly infinite state, systems. Moreover, since the simulations are required to be statistically independent, SMC may be easily and efficiently divided on parallel computing architectures.

SMC cannot be applied to MDPs without first resolving the nondeterminism. Since nondeterministic and probabilistic choices are interleaved in an MDP, schedulers are typically of the same order of complexity as the system as a whole and may be infinite. As a result, existing SMC algorithms for MDPs consider only memoryless schedulers and have other limitations (see Section 2).

Our Approach

We have created memory-efficient techniques to facilitate Monte Carlo verification of nondeterministic systems, without storing schedulers explicitly. In essence, the possibly infinite behaviour of schedulers is fully specified implicitly by the seed of a pseudo-random number generator. Our techniques therefore require almost no additional memory over standard SMC. In doing this, we are the first to provide the basis for a complete lightweight statistical alternative to the standard numerical verification algorithms for MDPs. A further contribution is our derivation of the statistical confidence bounds necessary to test multiple schedulers. These results suggest obvious solutions to problems encountered with existing algorithms that rely on multiple statistical tests (e.g., [START_REF] Henriques | Statistical model checking for Markov decision processes[END_REF]).

In this work we demonstrate the core ideas of our approach with simple SMC algorithms that repeatedly sample from scheduler space. Practical implementations require more sophisticated algorithms that adopt "smart sampling" (optimal use of simulation budget) and lightweight learning techniques. Some of our results make use of these ideas, but a full exposition is not possible here. The Kearns algorithm [START_REF] Kearns | A sparse sampling algorithm for nearoptimal planning in large Markov decision processes[END_REF] is the classic 'sparse sampling algorithm' for large, infinite horizon, discounted MDPs. It constructs a 'near optimal' scheduler piecewise, by approximating the best action from a current state using a stochastic depth-first search. Importantly, optimality is with respect to rewards, not probability (as required by standard model checking tasks). The algorithm can work with large, potentially infinite state MDPs because it explores a probabilistically bounded search space. This, however, is exponential in the discount. To find the action with the greatest expected reward in the current state, the algorithm recursively estimates the rewards of successive states, up to some maximum depth defined by the discount and desired error. Actions are enumerated while probabilistic choices are explored by sampling, with the number of samples set as a parameter. The error is specified as a maximum difference between consecutive estimates, allowing the discount to guarantee that the algorithm will eventually terminate.

There have been several recent attempts to apply SMC to nondeterministc models [START_REF] Bogdoll | Partial order methods for statistical model checking and simulation[END_REF][START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF][START_REF] Henriques | Statistical model checking for Markov decision processes[END_REF][START_REF] Hartmanns | On-the-fly confluence detection for statistical model checking[END_REF]. In [START_REF] Bogdoll | Partial order methods for statistical model checking and simulation[END_REF][START_REF] Hartmanns | On-the-fly confluence detection for statistical model checking[END_REF] the authors present on-the-fly algorithms to remove 'spurious' nondeterminism, so that standard SMC may be used. This approach is limited to the class of models whose nondeterminism does not affect the resulting probability of a property-scheduling makes no difference. The algorithms therefore do not attempt to address the standard MDP model checking problems related to finding optimal schedulers.

In [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF] the authors first find a memoryless scheduler that is near optimal with respect to a reward scheme and discount, using an adaptation of the Kearns algorithm. This induces a Markov chain whose properties may be verified with standard SMC. By storing and re-using information about visited states, the algorithm improves on the performance of the Kearns algorithm, but is thus limited to memoryless schedulers that fit into memory. The near optimality of the induced Markov chain is with respect to rewards, not probability, hence [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF] does not address the standard model checking problems of MDPs.

In [START_REF] Henriques | Statistical model checking for Markov decision processes[END_REF] the authors present an SMC algorithm to decide whether there exists a memoryless scheduler for a given MDP, such that the probability of a property is above a given threshold. The algorithm has an inner loop that generates candidate schedulers by iteratively improving a probabilistic scheduler according to sample traces that satisfy the property. The algorithm is limited to memoryless schedulers because the improvement process counts state-action pairs. The outer loop tests the candidate scheduler against the hypothesis using SMC and is iterated until an example is found or sufficient attempts have been made. The inner loop does not in general converge to the true optimum, but the outer loop randomly explores local maxima. This makes the number of samples used by the inner loop critical: too many may significantly reduce the scope of the random exploration and thus reduce the probability of finding the global optimum. A further problem is that the repeated hypothesis tests of the outer loop will eventually produce erroneous results. We address this phenomenon in Section 4.

We conclude that (i) no previous approach is able to provide a complete set of SMC algorithms for MDPs, (ii) no previous SMC approach considers history-dependent schedulers and (iii) no previous approach facilitates lightweight sampling from scheduler space.

Schedulers as Seeds of Random Number Generators

Storing schedulers as explicit mappings does not scale, so we have devised a way to represent schedulers using uniform pseudo-random number generators (PRNG) that are initialised by a seed and iterated to generate the next pseudorandom value. In general, such PRNGs aim to ensure that arbitrary subsets of sequences of iterates are uniformly distributed and that consecutive iterates are statistically independent. PRNGs are commonly used to implement the uniform probabilistic scheduler, which chooses actions uniformly at random and thus explores all possible combinations of nondeterministic choices. Executing such an implementation twice with the same seed will produce identical traces. Executing the implementation with a different seed will produce an unrelated set of choices. Individual deterministic schedulers cannot be identified, so it is not possible to estimate the probability of a property under a specific scheduler.

An apparently plausible solution is to use independent PRNGs to resolve nondeterministic and probabilistic choices. It is then possible to generate multiple probabilistic simulation traces per scheduler by keeping the seed of the PRNG for nondetermistic choices fixed while choosing random seeds for a separate PRNG for probabilistic choices. Unfortunately, the schedulers generated by this approach do not span the full range of general or even memoryless schedulers. Since the sequence of iterates from the PRNG used for nondeterministic choices will be the same for all instantiations of the PRNG used for probabilistic choices, the i th iterate of the PRNG for nondeterministic choices will always be the same, regardless of the state arrived at by the previous probabilistic choices. The i th chosen action can be neither state nor trace dependent.

General Schedulers Using Hash Functions

Our solution is to construct a per-step PRNG seed that is a hash of the an integer identifying a specific scheduler concatenated with an integer representing the sequence of states up to the present.

We assume that a state of an MDP is an assignment of values to a vector of system variables v i , i ∈ {1, . . . , n}. Each v i is represented by a number of bits b i , typically corresponding to a primitive data type (int, float, double, etc.). The state can thus be represented by the concatenation of the bits of the system variables, such that a sequence of states may be represented by the concatenation of the bits of all the states. Without loss of generality, we interpret such a sequence of states as an integer of n i=1 b i bits, denoted s, and refer to this in general as the trace vector. A scheduler is denoted by an integer σ, which is concatenated to s (denoted σ : s) to uniquely identify a trace and a scheduler.

Our approach is to generate a hash code h = H(σ : s) and to use h as the seed of a PRNG that resolves the next nondeterministic choice.

The hash function H thus maps σ : s to a seed that is deterministically dependent on the trace and the scheduler. The PRNG maps the seed to a value that is uniformly distributed but nevertheless deterministically dependent on the trace and the scheduler. In this way we approximate the scheduler functions S and M described in Section 1.1. Importantly, our technique only relies on the standard properties of hash functions and PRNGs. Algorithm 1 is the basic simulation function of our algorithms.

An Efficient Iterative Hash Function

To implement our approach, we have devised an efficient hash function that constructs seeds incrementally. The function is based on modular division [START_REF] Knuth | The Art of Computer Programming[END_REF]Ch. 6], such that h = (σ : s) mod m, where m is a suitably large prime. Since s is a concatenation of states, it is usually very much larger than the maximum size of integers supported as primitive data types. Hence, to generate h we use Horner's method [START_REF] Horner | A new method of solving numerical equations of all orders, by continuous approximation[END_REF][14, Ch. 4]: we set h 0 = σ and find h ≡ h n (n as given in Section 3.1) by iterating the recurrence relation

h i = (h i-1 2 bi + v i ) mod m. (1)
The size of m defines the maximum number of different hash codes. The precise value of m controls how the hash codes are distributed. To avoid collisions, a simple heuristic is that m should be a large prime not close to a power of 2 [START_REF] Cormen | Introduction to Algorithms[END_REF]Ch. 11]. Practically, it is an advantage to perform calculations using primitive data types that are native to the computational platform, so the sum in (1) should be less than or equal to the maximum permissible value. To achieve this, given x, y, m ∈ N, we note the following congruences:

(x + y) mod m ≡ (x mod m + y mod m) mod m (2) (xy) mod m ≡ ((x mod m)(y mod m)) mod m (3)
The addition in (1) can thus be re-written in the form of ( 2), such that each term has a maximum value of m -1:

h i = ((h i-1 2 bi ) mod m + (v i ) mod m) mod m (4)
To prevent overflow, m must be no greater than half the maximum possible integer. Re-writing the first term of (4) in the form of (3), we see that before taking the modulus it will have a maximum value of (m -1) 2 , which will exceed the maximum possible integer. To avoid this, we take advantage of the fact that h i-1 is multiplied by a power of 2 and that m has been chosen to prevent overflow with addition. We thus apply the following recurrence relation:

(h i-1 2 j ) mod m = (h i-1 2 j-1 ) mod m + (h i-1 2 j-1 ) mod m (5)
Equation ( 5) allows our hash function to be implemented using efficient native arithmetic. Moreover, we infer from (1) that to find the hash code corresponding to the current state in a trace, we need only know the current state and the hash code from the previous step. When considering memoryless schedulers we need only know the current state.

Confidence with Multiple Estimates

The Chernoff bound [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF][START_REF] Chernoff | A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations[END_REF] and Wald sequential probability ratio test [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF] are commonly used to bound errors of SMC algorithms. Their guarantees are probabilistic, such that with specified non-zero probability they produce an incorrect result. If such bounds are used on M schedulers, some of whose true probabilities lie in the interval (0, 1), then as M → ∞ the probability of encountering an error is a.s. 1. In particular, the maximum and minimum estimates will tend to 1 and 0, respectively, regardless of the true values.

To overcome this phenomenon, in Sects. 4.1 and 4.2 we derive new confidence bounds to allow SMC algorithms to test multiple schedulers. We illustrate their use with simple algorithms that sample M schedulers at random, where M is a parameter. These algorithms are the basis of a technique we call "smart sampling", which can exponentially improve convergence. The basic idea is to assign part of the simulation budget to obtain a coarse estimate of the extremal probabilities and to use this information to generate a set of schedulers that contains a "good" scheduler with high probability. The remaining budget is used to refine the set to find the best scheduler. Smart sampling has provided improvements of several orders of magnitude with the illustrated examples and is the subject of ongoing development. Lack of space prevents further discussion.

Sequential Probability Ratio Test for Multiple Schedulers

The sequential probability ratio test (SPRT) of Wald [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF] evaluates hypotheses of the form P(ω |= ϕ) ⊲⊳ p, where ⊲⊳∈ {≤, ≥}. The SPRT distinguishes between two hypotheses, H 0 : P(ω |= ϕ) ≥ p 0 and H 1 : P(ω |= ϕ) ≤ p 1 , where p 0 > p 1 . Hence, to evaluate P(ω |= ϕ) ⊲⊳ p, the SPRT requires a region of indecision (an 'indifference region' [START_REF] Younes | Probabilistic verification of discrete event systems using acceptance sampling[END_REF]) which may be specified by parameter ϑ, such that p 0 = p + ϑ and p 1 = p -ϑ. The SPRT also requires parameters α and β, which specify the maximum acceptable probabilities of errors of the first and second kind, respectively. An error of the first kind is incorrectly rejecting a true H 0 ; an error of the second kind is incorrectly accepting a false H 0 . To choose between H 0 and H 1 , the SPRT defines the probability ratio

ratio = n i=1 (p 1 ) 1(ωi|=ϕ) (1 -p 1 ) 1(ωi |=ϕ) (p 0 ) 1(ωi|=ϕ) (1 -p 0 ) 1(ωi |=ϕ) ,
where n is the number of simulation traces ω i , i ∈ {1, . . ., n}, generated so far. The test proceeds by performing a simulation and calculating ratio until one of two conditions is satisfied:

H 1 is accepted if ratio ≥ (1-β)/α and H 0 is accepted if ratio ≤ β/(1 -α).
To decide whether there exists a scheduler such that P(ω |= ϕ) ⊲⊳ p, we would like to apply the SPRT to multiple (randomly chosen) schedulers. The idea is to test different schedulers, up to some specified number M , until an example is found. Since the probability of error with the SPRT applied to multiple hypotheses is cumulative, we consider the probability of no errors in any of M tests. Hence, in order to ensure overall error probabilities α and β, we adopt

α M = 1 -M √ 1 -α and β M = 1 -M √ 1 -β in our stopping conditions. H 1 is accepted if ratio ≥ (1 -β M )/α M and H 0 is accepted if ratio ≤ β M /(1 -α M ).
Algorithm 2 demonstrates the sequential hypothesis test for multiple schedulers. If the algorithm finds an example, the hypothesis is true with at least the specified confidence.

Chernoff Bound for Multiple Schedulers

Given that a system has true probability p of satisfying a property, the Chernoff bound ensures P(| p -p |≥ ε) ≤ δ, i.e., that the estimate p will be outside the interval [p -ε, p + ε] with probability less than or equal to δ. Parameter δ is related to the number of simulations N by δ = 2e -2N ε 2 [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF], giving N = (ln 2 -ln δ)/(2ε 2 ) .

To demonstrate the scalability of our approach, we consider the choice coordination model of [START_REF] Ndukwu | An expectation transformer approach to predicate abstraction and data independence for probabilistic programs[END_REF] and estimate the minimum probability that a group of six tourists will meet within T steps. The model has a parameter (BOUND ) that limits the state space. We set BOUND = 100, making the state space of ≈ 5 × 10 16 intractable to numerical model checking. For T = 20 and T = 25 the true minimum probabilities are respectively 0.5 and 0.75. Using smart sampling and a Chernoff bound of ε = δ = 0.01, we correctly estimate the probabilities to be 0.496 and 0.745 in a few tens of minutes on a standard laptop computer.

Prospects and Challenges

Our techniques are immediately extensible to continuous time MDPs and other models that use nondeterminism. It is also seems simple to consider MDPs with rewards. Although the presented algorithms are not optimised with respect to simulation budget, in a forthcoming work we introduce the notion of "smart sampling" to maximise the chance of finding good schedulers with a finite budget.

A limitation of our approach is that the algorithms sample from only a subset of possible schedulers. It is easy to construct examples where good schedulers are vanishingly rare and will not be found. Our ongoing focus is therefore to develop memory-efficient learning techniques that construct schedulers piecewise, to improve convergence and consider a much larger set of schedulers.

Figure 1 .Figure 2 .

 12 Figure 1. Fragment of a typical Markov decision process.

Acknowledgement This work was partially supported by the European Union Seventh Framework Programme under grant agreement no. 295261 (MEALS).

Algorithm 2: Hypothesis testing with multiple schedulers

Input:

M, ϕ: the MDP and property of interest H ∈ {H0, H1}: the hypothesis of interest with threshold p ± ϑ α, β: the desired error probabilities of H M : the maximum number of schedulers to test Output: The result of the hypothesis test Let p 0 = p + ϑ and p 1 = p -ϑ be the bounds of H Let αM = 1 -M √ 1 -α and βM = 1 -M √ 1 -β Let A = (1 -βM )/αM and B = βM /(1 -αM ) Let U seed be a uniform PRNG and σ be its sample

The user specifies ε and δ and the SMC algorithm calculates N to guarantee the estimate accordingly. Equation ( 6) is derived from equations

giving N = (ln δ)/(2ε 2 ) to satisfy either inequality. We consider the strategy of sampling M schedulers to estimate the optimum probability. We thus generate M estimates {p 1 , . . . , pM } and take either the maximum (p max ) or minimum (p min ), as required. To overcome the cumulative probability of error with the standard Chernoff bound, we specify that all estimates pi must be within ε of their respective true values p i , ensuring that any pmin , pmax ∈ {p 1 , . . . , pM } are within ε of their true value. Given [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF] and the fact that all estimates pi are statistically independent, the probability that all estimates are less than their upper bound is expressed by P( M i=1 pi -

This leads to the following expression for N , given parameters M , ε and δ:

Since the case for p min is symmetrical, (8) also ensures P(p min -pmin ≥ ε) ≤ δ. Hence, to ensure the more usual conditions that P(|

N scales logarithmically with M (e.g., for ε = δ = 0.01, N ≈ log 1.0002 (M ) + 26472), making it tractable to consider many schedulers. Algorithm 3 is the resulting extremal probability estimation algorithm for multiple schedulers. 1 -δ ))/(2ε 2 ) be the no. of simulations per scheduler Let U seed be a uniform PRNG and σ its sample Initialise pmin ← 1 and pmax ← 0 Set seed of U seed randomly for i ∈ {1, . . . , M } do Iterate U seed to generate σ i Let truecount = 0 be the initial number of traces that satisfy ϕ for j ∈ {1, . . . , N } do

No schedulers were found to satisfy ϕ

Experiments

We implemented Algorithms 2 and 3 in our statistical model checking platform Plasma [1] and performed a number of experiments. Figure 3 shows the empirical cumulative distribution of schedulers generated by Algorithm 3 applied to the MDP of Fig. 2, using p 1 = 0.9, p 2 = 0.5, ϕ = X(ψ∧ XG 4 ¬ψ), ε = 0.01, δ = 0.01 and M = 300. The vertical red and blue lines mark the true probabilities of ϕ under each of the history-dependent and memoryless schedulers, respectively. The grey rectangles show the ±ε error bounds, relative to the true probabilities. There are multiple estimates per scheduler, but all estimates are within their respective confidence bounds. Note that the confidence is specified with respect to estimates, not with respect to optimality. Defining confidence with respect to optimality remains an open problem.

In Fig. 4 we consider a reachability property of the Wireless LAN (WLAN) protocol model of [START_REF] Kwiatkowska | Probabilistic Model Checking of the IEEE 802.11 Wireless Local Area Network Protocol[END_REF]. The protocol aims to minimise "collisions" between devices sharing a communication channel. We estimated the probability of the second collision at time steps {0, 10, . . . , 100}, using Algorithm 3 with M = 4000 schedulers per point. Maximum and minimum estimated probabilities are denoted by blue and red circles, respectively. Maximum probabilities calculated by numerical model checking are denoted by black crosses. The shaded areas indicate the ±ε error of the estimates (Chernoff bound ε = δ = 0.01) and reveal that our estimates are very close to the true values.


Benoît Boyer

Kevin Corre

Axel Legay

Sean Sedwards

PLASMA-lab: a flexible, distributable statistical model checking library

We present PLASMA-lab, a statistical model checking (SMC) library that provides the functionality to create custom statistical model checkers based on arbitrary discrete event modelling languages. PLASMAlab is written in Java for maximum cross-platform compatibility and has already been incorporated in various performance-critical software and embedded hardware platforms. Users need only implement a few simple methods in a simulator class to take advantage of our efficient SMC algorithms. PLASMA-lab may be instantiated from the command line or from within other software. We have constructed a graphical user interface (GUI) that exposes the functionality of PLASMA-lab and facilitates its use as a standalone application with multiple 'drop-in' modelling languages. The GUI adds the notion of projects and experiments, and implements a simple, practical means of distributing simulations using remote clients.

Background and motivation

Statistical model checking (SMC) is a form of probabilistic model checking that employs Monte Carlo methods to avoid the state explosion problem. SMC uses a number of independent simulation traces of a discrete event model to estimate the probability of a property. The traces may be generated on different machines, so SMC can efficiently exploit parallel computation (see Fig. 2). Reachable states are generated on-the-fly and the length of simulations is only weakly related to the size of the state space. Hence SMC tends to scale polynomially with respect to system description (see Fig. 1). Properties may be specified in bounded versions of the same temporal logics used in probabilistic model checking. Since SMC is thus applied to finite traces, it is also possible to use logics and functions that would otherwise be intractable or undecidable.

SMC abstracts the probabilistic model checking problem to one of estimating the parameter of a Bernoulli random variable with well defined confidence (e.g., using a Chernoff bound). The complexity of the estimation problem with respect to confidence is largely independent of the total number of possible traces. Hence SMC may also be applied to stochastic models with continuous states.

Dedicated SMC tools, such as YMER1 , VESPA, APMC2 and COSMOS3 , have been joined by statistical extensions of established tools such as PRISM 4and UPPAAL 5 . In the case of UPPAAL-SMC, this has required the definition of stochastic timed semantics. The tool MRMC6 has both numerical and statistical functionality, but takes as input a low level textual description of a Markov chain. Many other tools are available or under development, with most using a single high level modelling language related to a specific semantics. Our previous tool [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF] suffered the same limitation, prompting us to develop a radically new tool with modular architecture.

PLASMA-lab

PLASMA-lab [3] is an efficient SMC library written in Java, featuring a customisable simulator class. This allows SMC functionality to be added to existing domain-specific modelling platforms, such as DESYRE7 , and allows rapid prototyping of formal verification solutions using, e.g., Scilab8 and MATLAB9 . High performance standalone model checkers can also be constructed with PLASMAlab by including a suitable language parser in the simulator class. PLASMA-lab's integrated development environment facilitates distributed simulation and can work with multiple user-defined language plug-ins.

Properties PLASMA-lab accepts properties described in a form of bounded linear temporal logic (BLTL) extended with custom temporal operators based on concepts such as minimum, maximum and mean of a variable over time.

Model checking modes PLASMA-lab offers three basic modes of model checking: simple Monte Carlo, Monte Carlo using a Chernoff confidence bound and sequential hypothesis testing. There is also a simulation mode for debugging. Rare event model checking modes, such as importance sampling and importance splitting, can be implemented as part of the simulator class when the modelling semantics support them.

-Monte Carlo: the user explicitly specifies the number of simulations that PLASMA-lab must use to estimate the probability of a property. -Chernoff: the user specifies an absolute error ε and a probability δ. PLASMAlab calculates the number of simulations required to ensure that the resulting estimate is within ±ε of the correct value with minimum probability δ. -Sequential: PLASMA-lab adopts the sequential hypothesis ratio test of [START_REF] Younes | Verification and planning for stochastic processes with asynchronous events[END_REF] to verify that the probability of a property is above a user-specified threshold.

The user also specifies a level of indifference and parameters to control errors of Types I and II. The number of simulations is not specified a priori: simulations are performed as necessary. See [START_REF] Younes | Verification and planning for stochastic processes with asynchronous events[END_REF] for details.

Usage PLASMA-lab may be invoked from the command line or embedded in other software as a library. PLASMA-lab is provided as a pre-compiled jar file (plasmalab.jar) and a source template (Simulator.java) to create the simulator class. The minimum requirement is to implement the methods newTrace() and nextState(), that initiate a new simulation and advance the simulation by one step, respectively. Language parsers are typically invoked in the constructor.

Graphical user interface

The GUI provides an integrated development environment (IDE) to facilitate the use of PLASMA-lab as a standalone statistical model checker with multiple 'drop-in' modelling languages. To demonstrate this, we have included a biochemical language and a language based on reactive modules. The website [3] includes other examples. The GUI implements the notion of a project file, that links the description of a model to a specific modelling language simulator and a set of associated properties and experiments. The GUI also provides 2D and 3D graphical output of results and implements a distributed algorithm that will work with any of its associated modelling languages.

Distributed algorithm

The administrative time needed to distribute SMC on parallel computing architectures is often a deterrent. To overcome this, the PLASMA-lab GUI implements a simple and robust client-server architecture, based on Java Remote Method Invocation (RMI) using IPv4/6 protocols. The algorithm will work on dedicated clusters and grids, but can also take advantage of ad hoc networks of heterogeneous computers. The minimum requirement is that the IP address of the GUI is available to the clients. PLASMA-lab implements the SMC distribution algorithm of [START_REF] Younes | Verification and planning for stochastic processes with asynchronous events[END_REF], which avoids the statistical bias that might otherwise occur from load balancing. Distributed performance in illustrated in Fig. 2. The user selects the distributed mode via the GUI and publishes the IP address of the instance of PLASMA-lab GUI that is acting as server. Clients (instances of the PLASMA-lab service application) willing to participate respond by sending a message to the published IP address. The server sends an encapsulated version of the model and property to each of the participating clients, which then wait to be told how many simulations to perform. When sufficient clients are available, the user initiates the analysis by causing the server to broadcast the simulation requirements to each client.

Figure 1 .Figure 2 .

 12 Figure 1. Exponential scaling of numerical model checking vs. linear scaling of PLASMA-lab SMC, considering a fairness property of the probabilistic dining philosophers protocol.

www.tempastic.org/ymer

sylvain.berbiqui.org/apmc

www.lsv.ens-cachan.fr/˜barbot/cosmos/

www.prismmodelchecker.org

www.uppaal.org

www.mrmc-tool.org

www.ales.eu.com

www.scilab.org

www.mathworks.com

www.ict-dali.eu

www.danse-ip.eu

www.ales.eu.com

Applications

PLASMA-lab has been applied to problems from, e.g., systems biology, rare events, performance, reliability, motion planning and systems of systems [3]. PLASMA-lab is the focus of ongoing collaborations with companies Dassault, Thales, IBM, and EADS. PLASMA-lab is also used by several European projects. The following examples relate to the DALi 10 and DANSE 11 projects.

Motion planning PLASMA-lab is used by the DALi project in a novel motion planning application of SMC. DALi aims to develop an autonomous device to help those with impaired ability to negotiate complex crowded environments (e.g. shopping malls). High level constraints and the objectives of the user are expressed in temporal logic, while low level behaviour is predicted by the 'social force model' [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF].

PLASMA-lab was integrated with MATLAB to develop the prototype algorithm. The final version is implemented directly in C on embedded hardware and finds the optimum trajectory in a fraction of a second. PLASMA-lab improves the social force model's ability to avoid collisions by a factor of five.

Systems of systems

The DANSE project is concerned with the design and analysis of 'systems of systems' (SoS). SoS feature a dynamicity of configurations that introduces significant additional complexity (the state and state space of the model are not necessarily known a priori). PLASMA-lab is now an integral part of the DANSE software platform, using a Simulator class that wraps the DESYRE 12 hybrid simulation engine to make dynamicity transparent to SMC.


Roberto Amadini
email: amadini@cs.unibo.it

Maurizio Gabbrielli

Jacopo Mauro
email: jmauro@cs.unibo.it

Portfolio Approaches for Constraint Optimization Problems

Within the Constraints Satisfiability Problems (CSP) context, a methodology that has proved to be particularly performant consists in using a portfolio of different constraint solvers. Nevertheless, comparatively few studies and investigations has been done in the world of Constraint Optimization Problems (COP). In this work, we provide a generalization to COP as well as an empirical evaluation of different state of the art existing CSP portfolio approaches properly adapted to deal with COP. Experimental results confirm the effectiveness of portfolios even in the optimization field, and could give rise to some interesting future research.

Introduction

Constraint Programming (CP) is a declarative paradigm that allows to express relations between different entities in form of constraints that must be satisfied. One of the main goals of CP is to model and solve Constraint Satisfaction Problems (CSP) [START_REF] Mackworth | Consistency in Networks of Relations[END_REF]. Several techniques and constraint solvers were developed starting from the 1960s for solving CSPs and simplified CSPs problems such as the well-known boolean satisfiability problem (SAT), the Satisfiability Modulo Theories (SMT), and Answer Set Programming (ASP). The interest in this research area is still alive and active and, in particular, one of the more recent trends -especially in the SAT field -is trying to solve a given problem by using a portfolio approach [START_REF] Gomes | Algorithm portfolios[END_REF]. An algorithm portfolio is a general methodology that exploits a number of different algorithms in order to get an overall better algorithm. A portfolio of CP solvers can therefore be seen as a particular solver, dubbed portfolio solver, that exploits a collection of m > 1 different constituent solvers s 1 , . . . , s m in order to obtain a globally better CP solver. When a new unseen instance i comes, the portfolio solver tries to predict which are the best constituent solvers s 1 , . . . , s k (k ≤ m) for solving i and then runs such solver(s) on i. This solver selection process is clearly a fundamental part for the success of the approach and it is usually performed by exploiting Machine Learning (ML) techniques.

Exploiting the fact that different solvers are better at solving different problems, portfolios have proved to be particularly effective. For example, the overall winners of international solving competitions like [START_REF]Third International CSP Solver Competition[END_REF][START_REF]SAT Challenge[END_REF] are often portfolio solvers. Despite the proven effectiveness of the portfolio approach in the CSP case, and in particular in the SAT field, a few studies have tried to apply portfolio techniques to Constraint Optimization Problems (COPs). In these problems constraints are used to narrow the space of admissible solutions and then one has to find a solution that minimizes (maximizes) a specific objective function. This is done by using suitable constraint solvers integrated with techniques for comparing different solutions. Clearly a COP is more general than a CSP. Moreover, when considering portfolio approaches, some issues which are obvious for CSPs are less clear for COPs. For example, as we discuss later, defining a suitable metric which allows to compare different solvers is not immediate. These difficulties explain in part the lack of exhaustive studies on portfolios consisting of different COP solvers. Indeed, to the best of our knowledge, a few works deal with portfolios of COP solvers and some of them refer only to a specific problem like the Traveling Salesman Problem, while others use runtime prediction techniques for tuning the parameters of a single solver.

Nevertheless, this area is of particular interest since in many real-life applications we do not want to find just "a" solution for a given problem but "the" optimal solution, or at least a good one. In this work we tackle this problem and we perform a first step toward the definition of COP portfolios. We first formalize a suitable model for adapting the "classical" satisfaction-based portfolios to address COPs, providing also a metric to measure portfolio performances. Then, by using an exhaustive benchmark of 2670 instances, we test the performances of different portfolio approaches using portfolios composed from 2 to 12 different solvers. In particular, we adapt two among the best effective SAT portfolios, namely SATzilla and 3S, to the optimization field. We compare their performances w.r.t. some off-the-shelf approaches -built on top of the widely used ML classifiers -and w.r.t. SUNNY, a recent portfolio approach introduced in [3] that (on the contrary of those mentioned above) does not require an offline training phase.

Empirical results indicate that these approaches always significantly outperform the Single Best Solver available. Moreover, taking as reference the performances in the CSP field, we notice a closure of the performance gap between SATzilla and 3S w.r.t. off-the shelf approaches. Hence, we argue that a study of ad-hoc techniques for COPs may lead to further improvements of portfolio approaches like SATzilla and 3S in order to replicate the excellent results they achieved in the satisfaction field [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF][START_REF] Samulowitz | Snappy: A simple algorithm portfolio[END_REF]. Finally, we observe that the generalization of SUNNY to COPs appears to be particularly effective, since this algorithm has indeed reached the peak performances in our experiments.

Paper structure. In Section 2 we introduce the metrics adopted to evaluate the portfolio approaches for COPs. Section 3 presents the methodology and the portfolio algorithms we used to conduct the tests. The obtained results are detailed in Section 4 while related work is discussed in Section 5. We finally give concluding remarks and discuss feature extensions in Section 6.

Solution quality evaluation

When satisfaction problems are considered, the definition and the evaluation of a portfolio solver is straightforward. Indeed, the outcome of a solver run for a given time on a given instance can be either 'solved' (i.e. a solution is found or the unsatisfiability is proven) or 'not solved' (i.e. the solver does not say anything about the problem). Building and evaluating a CSP portfolio is then conceptually easy: the goal is to maximize the number of solved instances, solving them in less time as possible. Unfortunately, in the COP world the dichotomy solved/not solved is no longer suitable. A COP solver in fact can provide suboptimal solutions or even give the optimal one without being able to prove its optimality. Moreover, on the contrary of what happens in the satisfiability world, in order to speed up the search COP solvers could be executed in a non independent way. Indeed, the knowledge of a sub-optimal solution can be used by a solver to further prune its search space, and therefore to speed up the search process. Thus, the independent (even parallel) execution of a sequence of solvers may differ from a "cooperative" execution where the best solution found by a given solver is used as a lower bound by the solvers that are lunched afterwards.

Although the ideal goal is to prove the optimality in less time as possible, in the real world there is often the need of compromises. For many real life applications it is far better to get a good solution in a relatively short time rather than consume too much time to find the optimal value (or proving its optimality). In order to study the effectiveness of the portfolio approaches we therefore need new and more sophisticated evaluation metrics. In this work we then propose to give to each COP solver (portfolio based or not) a reward proportional to the distance between the best solution it finds and the best known solution. An additional reward is given if the optimality is proven, while a punishment is given if no solution are found without proving unsatisfiability. In particular, given an instance i, we assign to a solver s a score of 1 if it proves optimality for i, 0 if s does not find solutions. Otherwise, we give to s a score corresponding to the value of its best solution scaled into the range [0.25, 0.75], weighting 0.25 and 0.75 respectively the worst and the best known solutions of the known COP solvers.

In order to formally define the scoring function and to evaluate the quality of a solver, we denote with U the universe of the available solvers and with T the solving timeout in seconds that we are willing to wait at most. We use the function sol to define the solver outcomes. In particular we associate to sol(s, i, t) the outcome of the solver s for the instance i at time t. The value sol(s, i, t) can be either unk, if s does not find any solution for i; sat, if s finds at least a solution for i but does not prove the optimality; opt or uns if s proves optimality or unsatisfiability. Similarly, we use the function val to define the values of the objective function. In particular, with val(s, i, t) we indicate the best value of the objective function found by solver s for instance i at time t. If s does not find any solution for i at time t, the value val(s, i, t) is undefined. We assume the solvers behave monotonically, i.e. as time goes the solution quality gradually improves and never degrades.

We are now ready to associate to every instance i and solver s a weight that quantitatively represents how good is s when solving i. We define the scoring value of s (shortly, score) on the instance i at a given time t as a function score such that score(s, i, t) can be either:

(i) 0 if sol(s, i, t) = unk (ii) 1 if sol(s, i, t) ∈ {opt, uns} (iii) 0.75 if sol(s, i, t) = sat and MIN(i) = MAX(i) (iv) 0.75 - (val(s, i, t) -MIN(i)) • 0.5 MAX(i) -MIN(i) if sol(s, i, t) = sat, MIN(i) = MAX(i)
and i is a minimization problem

(v) (val(s, i, t) -MIN(i)) • 0.5 MAX(i) -MIN(i) + 0.25 if sol(s, i, t) = sat, MIN(i) = MAX(i)
and i is a maximization problem where MIN(i) and MAX(i) are the minimal and maximal objective function values found by any solver s at the time limit T1 and x = x if x ≥ 0; 0 otherwise. As an example, consider the scenario in Fig. 1 depicting the performances of three different solvers run on the same minimization problem. By choosing T = 500 as time limit, the score assigned to s 1 is 0.75 because it finds the solution with minimal value (40), the score of s 2 is 0.25 since it finds the solution with maximal value (50), and the score of s 3 is 0 because it does not find a solution. If instead T = 800, the score assigned to s 1 becomes 0.75-(40-10) * 0.5/(50-10) = 0.375 while the score of s 2 is 0.25 and the score of s 3 is 0.75. If instead T = 1000, since s 3 proves the optimality of the value 10 at time 900 (see the point marked with a star in Fig. 1 ) it receives a corresponding reward reaching then the score 1.

The score of a solver is therefore a measure in the range [0, 1] that is linearly dependent on the distance between the best solution it founds and the best solutions found by every other available solver. We decided to scale the values of the objective function in a linear way essentially for the sake of simplicity. Other choices, like for instance using the logarithm of the objective function for scaling or considering the subtended area T 0 val(s, i, t) dt could be equally useful and justifiable in a real scenario. The exploration of the impact of such alternative choices is however outside the scope of this paper, and left as a future work. Moreover in this work we assume the independent execution of the solvers, leaving as a future research the study of portfolio approaches that exploit the collaboration between different solver in order to boost the search speed.

In order to compare different portfolio approaches, we then considered the following evaluation metrics:

-Average Score (AS): the average of the scores achieved by the selected solver(s) on all the instances of the dataset; -Percentage of Optimums Found (POF): the percentage of instances of the dataset for which optimality is proven; -Average Optimization Time (AOT): the average time needed for proving optimality on every instance of the dataset, using a time penalty of T seconds when optimality is nor proven.

Methodology

Taking as baseline the methodology and the results of [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF] in this section we present the main ingredients and the procedure that we have used for conducting our experiments and for evaluating the portfolio approaches.

Solvers, dataset, and features

In order to build our portfolios we used a universe U of 12 solvers retrieved from the MiniZinc Challenge 2012, namely: BProlog, Fzn2smt, CPX, G12/FD, G12/LazyFD, G12/MIP, Gecode, izplus, JaCoP, MinisatID, Mistral and OR-Tools. We used all of them with their default parameters, their global constraint redefinitions when available, and keeping track of each solution found by every solver within a timeout of T = 1800 seconds. To conduct our experiments on a dataset of instances as realistic and large as possible, we have considered all the COPs of the MiniZinc 1.6 benchmark [24]. In addition, we have also added all the instances of the MiniZinc Challenge 2012, thus obtaining an initial dataset of 4977 instances in MiniZinc format.

In order to reproduce the portfolio approaches, we have extracted for each instance a set of 155 features by exploiting the features extractor mzn2feat [START_REF] Amadini | Features for Building CSP Portfolio Solvers[END_REF]. We filtered these features by scaling their values in the range [-1, 1] and by removing all the constant features. In this way, we ended up with a reduced set of 130 features on which we conducted our experiments. We have also filtered the initial dataset by removing, on one hand, the "easiest" instances (i.e. those for which the optimality was proven during the feature extraction) and, on the other, the "hardest" (i.e. those for which the features extraction has required more than T /2 = 900 seconds). These instances were discarded essentially for two reasons. First, if an instance is already optimized during the features extraction, no solver prediction is needed. Second, if the extraction time exceeds half of the timeout it is reasonable to assume that the recompilation of the MiniZinc model into FlatZinc format3 would end up in wasting the time available to solve the instance. The final dataset ∆ on which we conducted our experiments was thus constituted by 2670 instances.

Portfolio composition

After running every solver on each instance of the dataset ∆ keeping track of all the solutions found, we built portfolios of different size m = 2, . . . , 12. While in the case of CSPs the ideal choice is typically to select the portfolio of solvers that maximizes the number of solved instances, in our case such a metric is no longer appropriate since we have to take into account the quality of the solutions. Since maximizing the number of optimality proven instances appeared to be a too restrictive policy, we decided to select for each portfolio size m = 2, . . . , 12 the portfolio P m that maximizes the total score (possible ties have been broken by minimizing the solving time). Formally:

P m = arg max P ∈{S⊆U . |S|=m} i∈∆ max{score(s, i, T ) . s ∈ P }
We then elected a backup solver, that is a solver designated to handle exceptional circumstances like the premature failure of a constituent solver. After simulating different voting scenarios, the choice fell on CPX4 that in the following we refer also as Single Best Solver (SBS) of the portfolio. As a baseline for our experiments, we have also introduced an additional solver called Virtual Best Solver (VBS), i.e. an oracle solver that for every instance always selects and runs the solver of the portfolio having higher score (by using the solving time for breaking ties).

Portfolio Approaches

We tested different portfolio techniques. In particular, we have considered two state of the art SAT approaches (SATzilla and 3S) as well as some relatively simple off-the-shelf ML classifiers used as solver selector. Moreover, we have also implemented a generalization of the recently introduced CSP portfolio solver SUNNY [START_REF] Amadini | SUNNY: a Simple and Dynamic Algorithm Portfolio for Solving CSPs[END_REF] in order to deal with optimization problems. We would like to underline that in the case of 3S and SATzilla approaches we did not use the original methods which are tailored for the SAT domain. As later detailed, we have instead adapted these two approaches for the optimization world trying to modify them as little as possible. A study of alternative adaptations is outside the scope of this paper.

In the following we then provide an overview of these algorithms.

Off the shelf As in the case of satisfiability field [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF], off the shelf approaches were implemented by simulating the execution of a solver predicted by a ML classifier.

We then built 5 different approaches using 5 well-known ML classifiers, viz.: IBk, J48, PART, Random Forest, and SMO, and exploiting their implementation in WEKA [START_REF] Hall | The WEKA data mining software: an update[END_REF] with default parameters. In order to train the models we added for each instance of the dataset a label corresponding to the best constituent solver w.r.t. the score for such instance; for all the instances not solvable by any solver of the portfolio we used a special label no solver. In the cases where the solver predicted by a classifier was labeled no solver, we directly simulated the execution of the backup solver.

3S (SAT Solver Selector) [START_REF] Kadioglu | Algorithm Selection and Scheduling[END_REF] is a SAT portfolio solver that conjugates a fixedtime solver schedule with the dynamic selection of one long-running component solver: first, it executes for 10% of the time limit short runs of solvers; then, if a given instance is not yet solved after such time, a designated solver is selected for the remaining time by using a k-NN algorithm. 3S originally used a portfolio of 21 SAT solvers and was the best-performing dynamic portfolio at the International SAT Competition 2011. The major issue when adapting 3S for optimization problems is to compute the fixed-time schedule since, on the contrary of SAT problems, in this case the schedule should also take into account the quality of the solutions. We then tested different minimal modifications, trying to be as less invasive as possible and mainly changing the objective metric of the original Integer Programming (IP) problem used to compute the schedule. The performances of the different versions we tried were similar. Among those considered, the IP formulation that has achieved the best performance (with a peak AS of 0.78% more than the original one) is the one that: first, tries to maximize the solved instances; then, tries to maximize the sum of the score of the solved instances; finally, tries to minimize the solving time. 5 SATzilla [START_REF] Xu | SATzilla-07: The Design and Analysis of an Algorithm Portfolio for SAT[END_REF] is a SAT solver that relies on runtime prediction models to select the solver that (hopefully) has the fastest running time on a given problem 5 The objective function of the best approach considered was obtained by replacing that of the IP problem defined in [START_REF] Kadioglu | Algorithm Selection and Scheduling[END_REF] instance. Its last version [START_REF] Xu | SATzilla2012: Improved algorithm selection based on cost-sensitive classification models[END_REF] uses a weighted random forest approach provided with an explicit cost-sensitive loss function punishing misclassifications in direct proportion to their impact on portfolio performance. SATzilla won the 2012 SAT Challenge in the Sequential Portfolio Track. Unlike 3S, reproducing this approach turned out to be more straightforward. The only substantial difference concerns the construction of the runtimes matrix that is exploited by SATzilla to constructs its selector based on m(m -1)/2 pairwise cost-sensitive decision forests. 6 Since our goal is to maximize the score rather than to minimize the runtime, instead of using such a matrix we have defined a matrix of "anti-scores" P in which every element P i,j corresponds to the score of solver j on instance i subtracted to 1, that is P i,j = 1-score(j, i, T ).

SUNNY [START_REF] Amadini | SUNNY: a Simple and Dynamic Algorithm Portfolio for Solving CSPs[END_REF] is a lazy algorithm portfolio that, on the contrary of previously mentioned approaches, does not need an offline training phase. For a given instance i and a given portfolio P , SUNNY uses a k-NN algorithm to select from the training set a subset N (i, k) of the k instances closer to i. Then, on-the-fly, it computes a schedule of solvers by considering the smallest sub-portfolio S ⊆ P able to solve the maximum number of instances in the neighborhood N (i, k) and by allocating to each solver of S a time proportional to the number of solved instances in N (i, k).

Even in this case, we faced some design choices to tailor the algorithm for optimization problems. In particular, we decided to select the sub-portfolio S ⊆ P that maximizes the score in the neighborhood and we allocated to each solver a time proportional to its total score in N (i, k). In particular, while in the CSP version SUNNY allocates to the backup solver an amount of time proportional to the number of instances not solved in N (i, k), here we have instead assigned to it a slot of time proportional to k -h where h is the maximum score achieved by the sub-portfolio S.

Validation

In order to validate and test each of the above approaches we used a 5-repeated 5-fold cross validation [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF]. The dataset ∆ was then randomly partitioned in 5 disjoint folds ∆ 1 , . . . , ∆ 5 treating in turn one fold ∆ i , for i = 1, . . . , 5, as the test set and the union of the remaining folds j =i ∆ j as the training set. In order to avoid possible overfitting problems we repeated the random generation of the folds for 5 times, thus obtaining 25 different training sets (consisting of 534 instances each) and 25 different training sets (consisting of 2136 instances). For every instance of every test set we then computed the solving strategy proposed by the particular portfolio approach and we simulated it using a time cap of 1800 seconds. For estimating the solving time we have taken into account both the time needed for converting a MiniZinc model to FlatZinc and the time needed for extracting the features. In order to evaluate the performances, we then computed the metrics introduced in the previous section. (a) Results considering all the approaches and the VBS.

Results

In this section we present the obtained experimental results. 7 . In Fig. 2, 3, 4 we summarize the results obtained by the various techniques considering portfolios of different sizes and by using the Average Score (AS), the Percentage of Optimum Found (POF), and the Average Optimization Time (AOT) metrics introduced in Section 2. For ease of reading, in all the plots we report only the two best approaches among all the off-the-shelf classifiers we evaluated, namely Random Forest (RF) and SMO.

The source code developed to conduct and replicate the experiments is available at http://www.cs.unibo.it/ ~amadini/lion_2014.zip 4.1 Average Score Fig. 2a shows the AS performances of the various approaches, setting as baseline the performances of the Virtual Best Solver (VBS). Fig. 2b for the sake of readability visualizes the same results considering the VBS baseline, the two best approaches only (SUNNY and 3S) and the Single Best Solver (SBS) performance as additional baseline.

All the considered approaches have good performances and they greatly outperform the SBS. As in the case of CSP [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF][START_REF] Amadini | Features for Building CSP Portfolio Solvers[END_REF], it is possible to notice that off-the-shelf approaches have usually lower performances even though the gap between the best approaches and them is not that pronounced. The best portfolio approach is SUNNY that reaches a peak performance of 0.8802 using a portfolio of just 6 solvers. SUNNY is able to close the 91.35% of the gap between the SBS and VBS. 3S however has performances close to those of SUNNY and in particular its best performance (0.8718 with 6 and 12 solvers) is very close to the peak performance of SUNNY. Strangely enough, we can notice that both SUNNY and 3S have non monotonic performances when the portfolio sizes increases. This is particularly evident looking at their performance decrease when a portfolio of size 7 is used instead of one with just 6 solvers. This instability is obviously a bad property for a portfolio approach. We think that in this case it may be due to the noise of the addition of a solver that may disrupt the choices made by the k-NN algorithm on which SUNNY and 3S rely.

SATzilla does not reach the performances of SUNNY or 3S, even though for big portfolio sizes its performances are rather close. Moreover its behavior is monotonic w.r.t. the increase of the size of the portfolio. Hence, it seems that SATzilla is more reliable and scalable and, as also noticed in [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF], it is the only approach that does not present a degradation of performances for portfolios with more than 6 solvers.

Percentage of Optimum Found

Looking at the number of optimum found it is clear from Fig. 3a and 3b that there is a sharp demarcation of SUNNY w.r.t. other approaches. SUNNY appears to find many more optimum w.r.t. the other techniques, reaching a maximum POF of 57.03%. We think that the performances of SUNNY exploit the fact that it schedules more than one solver reducing the risk of making wrong choices. (a) Results considering all the approaches and the VBS. Moreover, it uses this schedule for the entire time window (on the contrary, 3S a static schedule only for 10% of the time window). Another interesting fact is that SUNNY mimics the behavior of the VBS. Thus, SUNNY seems able to properly use the addition of a solver to prove the optimality of instances exploiting the capability of the newly added solver.

Regarding other approaches, it can be observed by the overlapping of their curves in Fig. 3a that they are basically equivalent. What may seem surprising is that the best among them is SMO, which instead turned out to be the worst by considering the AS (Fig. 2a).

Even in this case, as shown in Fig. 3b all the portfolio approaches greatly outperform the SBS. SUNNY in particular is able to close the 85.73% of the gap between the SBS and VBS. Finally, note that there is a significant correlation (the Pearson coefficient is about 0.78) between AS and POF. Hence, maximizing the score is almost equivalent to maximizing the number of optimum to find.

Average Optimization Time

When the AOT metric is considered we can notice that the 3S approach does not perform very well compared to the other approaches. We think that this is due to the fact that 3S is a portfolio that uses more than one solver and it does not employ heuristics to decide which solver has to be executed first. SUNNY instead does not suffer from this problem since it schedules the solvers according to their performances on the already known instances. However, 3S is still able to always outperform the SBS for each portfolio size.

While the performance of SATzilla and the off-the-shelf approaches appear to be very similar, even in this case we can observe the good performances of SUNNY that is able to close the 94.27% of the gap SBS/VBS reaching a peak performance of 786.05 seconds with a portfolio of 10 solvers.

The (anti-)correlation between AOT and AS is lower than the one between POF and AS (the Pearson coefficient is -0.73) but still considerable. It is instead very strong, as can be guessed looking at the definition of these two metrics, the anti-correlation between AOT and POF (the Pearson coefficient is -0.98). This means that trying to maximize the average percentage score is like trying to minimizing the average solving time and to maximize the number of proven optimum.

Finally, we would like to mention that the AOT metric could be too strict and not so significant. In fact, if a solver finds the best value after few seconds and stops its execution without proving optimality it is somewhat over penalized with the timeout value T . In future it may therefore be interesting to study other ways to weight and evaluate the solving time (e.g. a rationale metric could be to consider a properly normalized area under the curve time/value defined by each solver behavior).

Related work

As far as the evaluation of optimization solvers and portfolio approaches is concerned, there exist a variety of metrics used to rank them. Among those used in practice by well known solving competitions worth mentioning are those that rank the solvers by using the number of the solved instances first, considering solving time later in case of ties [START_REF][END_REF][START_REF]SAT Challenge[END_REF]. In [START_REF]MiniZinc Challenge[END_REF] instead the ranking is performed by using a Borda count, i.e. a single-winner election method in which voters rank candidates in order of preference. Differently from the metrics defined in Section 2, these metrics address the quality of the solutions in a less direct way (i.e. by making pairwise comparisons between the score of the different solvers).

In the previous section we have already mentioned SATZilla [START_REF] Xu | SATzilla-07: The Design and Analysis of an Algorithm Portfolio for SAT[END_REF] and 3S [START_REF] Kadioglu | Algorithm Selection and Scheduling[END_REF] as two of the most effectives portfolio approaches in the SAT and CSP domain. For a comprehensive survey on portfolio approaches applied to SAT, planning, and QBF problems we refer the interested reader to the comprehensive survey [START_REF] Kotthoff | Algorithm Selection for Combinatorial Search Problems: A Survey[END_REF] and to [START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF] for CSPs.

As far as optimization problems are concerned, in the 2008 survey on algorithm selection procedures [START_REF] Smith-Miles | Cross-disciplinary perspectives on meta-learning for algorithm selection[END_REF] the authors observe that "there have been surprisingly few attempts to generalize the relevant meta-learning ideas developed by the machine learning community, or even to follow some of the directions of Leyton-Brown et al. in the constraint programming community." To the best of our knowledge, we think that the situation has not improved and are not aware of more recent works addressing explicitly the construction of portfolio solvers for COPs. Indeed, in the literature, we are aware of portfolio approaches developed just for some specific instances of COP. For instance, problems like Mixed Integer Programming, Scheduling, Most Probable Explanation (MPE) and Travel Salesman Problem (TSP) are addressed by means of portfolio techniques exploiting ML methods in [START_REF] Guo | A machine learning approach to algorithm selection for NP-hard optimization problems: a case study on the MPE problem[END_REF][START_REF] Hutter | Algorithm Runtime Prediction: The State of the Art[END_REF].

Other related works target the analysis of the search space of optimization problems by using techniques like landscape analysis [START_REF] Knowles | Towards Landscape Analyses to Inform the Design of Hybrid Local Search for the Multiobjective Quadratic Assignment Problem[END_REF], Kolmogorov complexity [START_REF] Borenstein | Kolmogorov complexity, Optimization and Hardness[END_REF], and basin of attractions [START_REF] Merz | Advanced Fitness Landscape Analysis and the Performance of Memetic Algorithms[END_REF]. Some approaches like [START_REF] Leyton-Brown | Learning the Empirical Hardness of Optimization Problems: The Case of Combinatorial Auctions[END_REF][START_REF] Amanda | Towards insightful algorithm selection for optimisation using meta-learning concepts[END_REF] also use ML techniques to estimate the search space of some algorithms and heuristics on optimization problems. These works look interesting because precise performance evaluations can be exploited in order to built portfolios as done, for instance, by SATzilla [START_REF] Xu | SATzilla-07: The Design and Analysis of an Algorithm Portfolio for SAT[END_REF] in the SAT domain or by [START_REF] Lobjois | Branch and Bound Algorithm Selection by Performance Prediction[END_REF] for optimization problems solved by using branch and bound algorithms.

Another related work is [START_REF] Telelis | Combinatorial Optimization through Statistical Instance-Based Learning[END_REF] where ML algorithms are used to solve the Knapsack and the Set Partitioning problems by a run-time selection of different heuristics during the search. In [START_REF] Kadioglu | ISAC -Instance-Specific Algorithm Configuration[END_REF][START_REF] Xu | Hydra-MIP: Automated Algorithm Configuration and Selection for Mixed Integer Programming[END_REF] automated algorithm configurators based on AI techniques are used to boost the solving process of MIP and optimization problems. In [START_REF] Carchrae | Applying Machine Learning to Low-Knowledge Control of Optimization Algorithms[END_REF] a low-knowledge approach that selects solvers for optimization problems is proposed. In this case, decisions are based only on the improvement of the solutions quality, without relying on complex prediction models or on extensive set of features.

Conclusions and Extensions

In this paper we tackled the problem of developing a portfolio approach for solving COPs. In particular, in order to evaluate the performances of a COP solver we proposed a scoring function which takes into account the solution quality of the solver answers. We then proposed three different metrics to evaluate and compare COP solvers. These criteria were used to compare different portfolio techniques adapted from the satisfiability world with others based on classifiers and with a recently proposed lazy portfolio approach.

The results obtained clearly indicate that exploiting portfolio approaches leads to better performances w.r.t. using a single solver. We conjecture that, especially when trying to prove optimality, the number of times a solver can not give an answer is not negligible and that the solving times have a heavy-tail distribution typical of complete search methods [START_REF] Gomes | Heavy-Tailed Distributions in Combinatorial Search[END_REF]. Hence, as proven in [START_REF] Gomes | Algorithm portfolios[END_REF], a COP setting can be considered an ideal scenario to apply a portfolio approach and obtain statistically better solvers exploiting existing ones.

We noticed that, even though at a first glance it can seem counterintuitive, the best performances were obtained by the portfolio approach which uses a schedule of more than one solver. In these cases the risk of choosing the wrong solver is reduced and, apparently, this is more important than performing part of the computation again, as could happen when two (or more) solvers are lunched on the same instance.

We also noticed that the adaptation of ad-hoc methods (deriving from SAT) does not lead to the same gain of performance that these methods provide in the CSP and SAT field. We believe that the study of new, specific ad-hoc techniques for COPs should be done in order to obtain the same advantages of the SAT field, thus obtaining better portfolios than those which use off-the shelf ML classifiers. This is however left as a future work, as well as adapting and testing other promising portfolio approaches like [START_REF] Kadioglu | ISAC -Instance-Specific Algorithm Configuration[END_REF][START_REF] Malitsky | Algorithm Portfolios Based on Cost-Sensitive Hierarchical Clustering[END_REF][START_REF] Omahony | Using case-based reasoning in an algorithm portfolio for constraint solving[END_REF].

Another direction for further research is the study of how cooperative strategies can be used among the constituent solvers, both in the sequential case and in a parallel setting, where more than one solver of the portfolio is allowed to be executed at the same time.

As previously said, we would also like to study the impact of using other metrics to evaluate the solution quality of the solvers. On the basis of the empirical correlation among the metrics so far considered we are confident that the performance of portfolio approaches should be robust, i.e. the rank of good portfolios approaches does not depend on the specific metric used, provided that the metric is enough "realistic".

The final aim of this work and its future extension would be to exploit our results in order to create a first COP portfolio solver to be submitted to the next MiniZinc Challenge.

2

 2

3 Fig. 1 :

 31 Fig. 1: Solver performances example.

  Results considering SBS, VBS, and the best two approaches.

Fig. 2 :

 2 Fig. 2: AS performances.

Fig. 3 :

 3 Fig. 3: POF performances.

  Results considering SBS, VBS, and the best two approaches.

Fig. 4 :

 4 Fig. 4: AOT performances.

Formally, MIN(i) = min Vi and MAX(i) = max Vi where Vi = {val(s, i, T ) . s ∈ U }.

Avoiding negative scores was necessary because a solver executed by a portfolio solver for t < T seconds could produce a solution that is worse than any solution found by any solver of the universe after T seconds. In this case, the score would be less than 0.25 and therefore could potentially be less than 0. This however is very uncommon: in our experiments we noticed that the 0 score was assigned only to the solvers that did not find any solution.

FlatZinc[START_REF] Becket | Specification of FlatZinc -Version 1[END_REF] is the low level language that each solver uses for solving a given MiniZinc instance. A key feature of FlatZinc is that, starting from a general MiniZinc model, every solver can produce a specialized FlatZinc by redefining the global constraints definitions. We noticed that, especially for huge instances, the time needed for extracting features was strongly dominated by the FlatZinc conversion. However, for the instances of the final dataset this time was in average 10.36 seconds, with a maximum value of 504 seconds and a median value of 3.17

seconds.[START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF] Following[START_REF] Amadini | An Empirical Evaluation of Portfolios Approaches for Solving CSPs[END_REF] methodology, CPX won all the elections we simulated using different criteria, viz.: Borda, Approval, and Plurality.

For more details, we defer the interested reader to[START_REF] Xu | SATzilla2012: Improved algorithm selection based on cost-sensitive classification models[END_REF]

To conduct the experiments we used Intel Dual-Core 2.93GHz computers with 2 GB of RAM and Ubuntu operating system.

(a) Results considering all the approaches and the VBS.