Yves Mocquard
email: yves.mocquard@irisa.fr

Bruno Sericola
email: bruno.sericola@inria.fr

Emmanuelle Anceaume
email: emmanuelle.anceaume@irisa.fr

Population Protocols with Convergence Detection

Keywords: Population protocols, Detection of convergence, Large scale systems, Anonymous agents, Probabilistic analysis

This paper focuses on pairwise interaction-based protocols, and proposes an universal mechanism that allows each agent to locally detect that the system has converged to the sought configuration with high probability. To illustrate our mechanism, we use it to detect the instant at which the proportion problem is solved. Specifically, let nA (resp. nB) be the number of agents that initially started in state A (resp. B) and γA = nA/n, where n is the total number of agents. Our protocol guarantees, with a given precision ε > 0 and any high probability 1 -δ, that after O(n ln(n/δ)) interactions, any queried agent that has set the detection flag will output the correct value of the proportion γA of agents which started in state A, by maintaining no more than O(ln(n)/ε) integers. We are not aware of any such results. Simulation results illustrate our theoretical analysis.

I. INTRODUCTION

The main line of research in the population protocol model has so far been the design of pairwise interaction-based protocols that converge to a given configuration of the system as fast as possible while minimizing the number of states needed to converge to that sought configuration. Actually, since the seminal work of Aspnes [START_REF] Aspnes | An introduction to population protocols[END_REF], a considerable amount of work has been done so far to determine which properties can emerge from pairwise interactions between finite-state nodes, together with the derivation of lower bounds on the time and space needed to reach such properties (e.g., [START_REF] Alistarh | Fast and exact majority in population protocols[END_REF], [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF], [START_REF] Angluin | A simple population protocol for fast robust approximate majority[END_REF], [START_REF] Draief | Convergence speed of binary interval consensus[END_REF], [START_REF] Mertzios | Determining majority in networks with local interactions and very small local memory[END_REF], [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF], [START_REF] Mocquard | Balanced allocations and global clock in population protocols: An accurate analysis[END_REF], [START_REF] Mocquard | Analysis of the propagation time of a rumour in largescale distributed systems[END_REF], [START_REF] Perron | Using three states for binary consensus on complete graphs[END_REF]).

In this paper we go a step further by proposing a mechanism that allows each agent to locally detect that the system has converged to the sought configuration with high probability. As an application, we propose to use this mechanism to detect the instant at which the proportion problem is solved. Specifically, let n A (resp. n B ) be the number of agents that initially started in state A (resp. B) and γ A = n A /n be the majority, where n is the total number of agents. Our protocol guarantees, with a given precision ε > 0 and any high probability 1-δ, that after O(n ln(n/δ)) interactions, any queried agent that has set the detection flag will output the correct value of the proportion γ A of agents which started in state A, by maintaining no more than O(ln(n)/ε) states.

To allow each node to locally detect that its computation of the proportion has converged, we combine three algorithms, each one being run at each node of the system. The first one This work was partially funded by the French ANR project SocioPlug (ANR-13-INFR-0003), and by the DeSceNt project granted by the Labex CominLabs excellence laboratory (ANR-10-LABX-07-01).

is dedicated to the computation of the sought property, that is the computation of the proportion γ A . The second algorithm, run in parallel with the proportion one, aims at constructing the global clock of the system to detect the instant at which convergence is reached at all the nodes. Briefly, when the local clock of at least two nodes have reached a given threshold T max , this means that the number of global interactions in the system is large enough so that all the nodes of the system are able to compute the proportion value with high precision. Thus both nodes can start the propagation of a signal to inform each other node i of the system that i can derive from its local state a good approximation of the proportion γ A . This dissemination is implemented by a randomized propagation algorithm.

We provide in Section V, a new theoretical analysis of the performance of all these three pairwise interaction-based protocols that improve upon existing ones. As will also be proven in Section V, this detection mechanism is universal in the sense that any pairwise interaction-based population protocol can be augmented with this mechanism to safely detect convergence with high probability. The only requirement that must satisfied is that an upper bound on the convergence time of that protocol must be explicitly known.

The remaining of the paper is organized as follows. Section II formalizes the addressed problem. The model of the system together with the different notations adopted in the paper are presented in Section III. The orchestration of the different ingredients of our detection mechanism are presented in Section IV. A deep theoretical analysis of the performance of our detection mechanism is presented in Section V, and a summary of the simulation results is given in Section VI. Finally, Section VII concludes the paper.

II. THE ADDRESSED PROBLEM

We consider a set of n agents, interconnected by a complete graph, that asynchronously start their execution in one of two distinct states A and B. Let n A (resp. n B ) be the number of agents whose initial state is A (resp. B), and let γ A = n A /n be the proportion of the system, with n the total number of agents. We formalize the addressed problem as follows.

Definition 1 (Proportion with Convergence Detection): A population protocol ran by all the nodes of the system solves the proportion with convergence detection problem if with probability at least 1 -δ, for any δ ∈ (0, 1), any node of the system is capable of computing the proportion γ A and detecting the instant at which the computed proportion is an ε-approximation of γ A . The number of interactions and the number of states needed to guarantee the convergence detection must respectively be O(ln(n)) and O(ln(n)/ε), where ε is the precision of the computed proportion, and n is the number of nodes.

III. MODEL AND NOTATIONS

In this work we assume that the collection of nodes communicate through pairwise and asynchronous interactions. Initially, each node starts with an initial symbol A or B represented by ι. The input function of each node initializes its local state according to its initial symbol, and then at each interaction its state is updated using a transition function denoted by f . Interactions between nodes are random: at each discrete time, any two agents are randomly selected to interact. The notion of time in population protocols refers to as the successive steps at which interactions occur, while the parallel time refers to as the total number of interactions averaged by n, see Aspnes et al. [START_REF] Aspnes | An introduction to population protocols[END_REF]. Note that nodes do not maintain nor use identifiers, however for ease of presentation, they are numbered 1, 2, . . . , n.

We will denote by (C

(i) t , T (i) t , S (i)
t ) the state of node i at time t, where C (i) t is used to evaluate the current value of the proportion at node i, T

(i) t
represents the global clock of the system, and S (i) t is a Boolean variable indicating whether proportion convergence has been globally reached or not. Let m and T max be two systems parameters, respectively used to define the initial configuration of the nodes and to determine the global number of interactions after which convergence is reached for all the nodes. Values of both parameters are analyzed in Section V. At any time t, the state set of C The configuration of the system at time t is the state of each node at time t and is denoted by (C t , T t , S t ) where

C t = (C (1) t , . . . , C (n) t ), T t = (T (1) t , . . . , T (n) t ), and S t = (S (1) t , . . . , S (n) t ).
Interactions between nodes are orchestrated by a random scheduler: at each discrete time t ≥ 0, any two indices i and j are randomly chosen to interact with probability p i,j (t). The successive choices of the interacting pair of nodes are supposed to independent and uniformly distributed, which means that we have

P i,j (t) = 1 {i =j} n(n -1)
.

IV. ALGORITHM RUN AT EACH NODE

Each node i maintains, as its current state, a vector made of three components (C (i) , T (i) , S (i) ), initialized according to Algorithm 1.

Pairs of nodes interact randomly (see Algorithm 2) and during their interaction update their state by computing the average of their C values and by incrementing their clock values T 1 Function Init(i):

2 if ι (i) = A then C (i) := m; 3 if ι (i) = B then C (i) := -m; 4 T (i) := 0; 5 S (i) := 0;
Algorithm 1: Initialization of node i's state (input function) as respectively described in Algorithms 3 and 4. The transition function f is given by f (x, y) = ( (x + y)/2 ), (x + y)/2 ). When two interacting nodes have both their clock equal to T max -1 (i.e. the number of global interactions in the system is large enough to allow all the nodes of the system to locally compute the proportion value with high precision), they both set their signal value to 1, which indicates the starting point of the spreading (see Algorithm 5). If during an interaction, a node updates its signal value, i.e. S = 1, then it stops updating both its C and T values, and at any subsequent interactions, this node will only "propagate" signal S.

1 Function UpdateState(i, j):

2 if Spreading (i, j) = 0 then 3 if Clock (i, j) = 0 then 4 Average (i, j); 5 end 6 end
Algorithm 2: Update of the state of two nodes during their interaction 1 Function Average(i,j):

2 (C (i) , C (j) ) := C (i) +C (j) 2
, C (i) +C (j)

2

;

Algorithm 3: Update of C values of two interacting nodes

1 Function Clock(i, j):

2 if T (i) = T (j) = (T max -1) then 3 S (i) := S (j) := 1; 4 return 1; 5 end 6 if T (i) ≤ T (j) then T (i) := T (i) + 1 ; 7 else T (j) := T (j) + 1 ; 8 return 0;
Algorithm 4: Update of T values of two interacting nodes Upon query, a node i returns its estimation ω A of the proportion of initial A according to its current value of C (i) . We have

ω A (C (i) ) = (m + C (i) )/(2m)
Note that in addition to the proportion, node i also returns its signal S (i) (see Algorithm 6). As demonstrated in Section V, the proportion computed by i is an ε-approximation of γ A with any high probability 1 -δ if S (i) = 1. Note that if S (i) = 0 then i does not know how far its computation is from γ A .

1 Function Spreading(i, j):

2 if S (i) t = S (j) t
= 0 then return 0 ; This section is devoted to the analysis of our solution. We split the analysis into five parts, the first one devoted to the analysis of the rumor spreading function (see Section V-A), the second one to the analysis of the average function (see Section V-B) and the third one to the global clock function (see Section V-C). The analysis presented in the fourth part consists in evaluating the global behavior of our protocol by combining the previous evaluations (see Section V-D). We end this section by showing under which hypothesis our convergence detection mechanism can be applied to any pairwise interaction-based protocol (see Section V-E).

A. Analysis of the rumor spreading function

The spreading starts at the first instant t at which two interacting nodes, say i and j, have both their spreading values, S (i) t and S (j) t , equal to 1. This instant occurs when both nodes have their clock equal to T max -1. In order to analyze this spreading time, we first prove a lemma derived from Theorem 4 of [START_REF] Mocquard | Analysis of the propagation time of a rumour in largescale distributed systems[END_REF]. let Y t the number of informed nodes at time t and Θ n the first instant at which all nodes know the rumor. We have

Θ n = inf{t ≥ 0 | Y t = n}.
Note that in our case we have Y 0 = 2. This lemma gives a maximal value of spreading time with a probability less than or equal to any fixed probability δ ∈ (0, 1) when the system initially starts with 2 nodes knowing the rumor. Lemma 2: For all δ, we have

P{Θ n ≥ n (ln(n) -ln(δ)/2) | Y 0 = 2} ≤ δ.
Proof. Applying Theorem 4 of [START_REF] Mocquard | Analysis of the propagation time of a rumour in largescale distributed systems[END_REF] with i = 2 leads, for every k ≥ 0, to

P{Θ n ≥ k | Y 0 = 2} ≤ (n -2) 2 1 -2 n k . Setting k = n (ln(n) -ln(δ)/2) , we obtain 1 - 2 n n(ln(n)-ln(δ)/2) ≤ 1 - 2 n n(ln(n)-ln(δ)/2) = e n(ln(n)-ln(δ)/2) ln(1-2/n) .
Using the fact that ln(1 -x) ≤ -x, for all x ∈ [0, 1), we get

1 - 2 n n(ln(n)-ln(δ)/2) ≤ e -2 ln(n)+ln(δ) = δ n 2
and thus

P{Θ n ≥ n (ln(n) -ln(δ)/2) | Y 0 = 2} ≤ (n -2) 2 δ n 2 ≤ δ,
which completes the proof.

Note that the proof of this lemma does not make use of the Markov inequality. The approximations have all been made with equivalents, which means that the result of this lemma is quite close to the reality. This will be illustrated in section VI-A.

B. Analysis of the average function

The average function is modelled by vector C t . This transition function is given, for interacting nodes i and j, by

C (i) t+1 , C (j) t+1 = C (i) t + C (j) t 2 , C (i) t + C (j) t 2 (1)
and

C (r) t+1 = C (r) t for r = i, j.
The following lemma, which states that the sum of the entries of vector C t is constant, is straightforward.

Lemma 3: For every t ≥ 0, we have

n i=1 C (i) t = n i=1 C (i) 0 .
Proof. The proof is immediate since the transformation from C t to C t+1 described in Relation (1) does not change the sum of the entries of C t+1 . Indeed, from Relation (1), we have

C (i) t+1 + C (j) t+1 = C (i) t + C (j) t
and the other entries do not change their values.

We denote by the mean value of the sum of the entries of C t and by L the row vector of R n with all its entries equal to , that is

= 1 n n i=1 C (i) t
and L = ( , . . . , ).

Clearly, from Lemma 3, the value of is independent of the time t. The following theorem shows that, after a given amount of time, the distance between all the C (i) t and is less than 3/2 with any high probability. Recall that the infinite norm is defined for any n-dimensional vector v by v ∞ = max i=1,...,n v i . Theorem 4: For all δ ∈ (0, 1), if there exists a constant K such that C 0 -L ≤ K then, for every t ≥ n (2 ln(K) + 1.78 ln(n) -7.60 ln(δ) + 2.70), we have

P { C t -L ∞ < 3/2} ≥ 1 -δ.
Proof. See Appendix.

We now apply these results to compute the proportion γ A of agents whose initial input was A, which is given by γ

A = n A /(n A + n B ) = n A /n. Recall that the output function ω A is given, for all x ∈ Q C , by ω A (x) = (m + x)/(2m).
The following theorem shows that, after a given amount of time, the value of all the ω A (C (i) t ) is an ε-approximation of the proportion γ A with any probability.

Theorem 5: For all ε, δ ∈ (0, 1), setting m = 3/(4ε) , we have, for all t ≥ n (2.78 ln(n) -2 ln(ε) -7.60 ln(δ) + 2.70),

P {|ω A (C (i) t ) -γ A | < ε for all i = 1, . . . , n} ≥ 1 -δ. Proof. We have C 0 -L ≤ m √ n. Applying Theorem 4, with K = √ n/ε ≥ 3/(4ε) √ n = m √ n,
we obtain for all δ ∈ (0, 1) and t ≥ n (2.78 ln(n) -2 ln(ε) -7.60 ln(δ) + 2.70),

P { C t -L ∞ ≥ 3/2} ≤ δ or equivalently P {|C (i) t -(γ A -γ B )m| < 3/2, for all i = 1, . . . , n} ≥ 1 -δ. Since γ A + γ B = 1 we have |C (i) t -(γ A -γ B )m| = |C (i) t -(2γ A -1)m| = |m + C (i) t -2mγ A | = 2m|ω A (C (i) t ) -γ A |. Then P {|ω A (C (i) t ) -γ A | < 3/(4m), for all i = 1, . . . , n} ≥ 1 -δ. So P {|ω A (C (i) t ) -γ A | < ε, for all i = 1, . . . , n} ≥ 1 -δ,
which completes the proof.

C. Analysis of the clock function

The clock function is modelled by vector T t , and for any two interacting nodes i and j, we have

T (i) t+1 , T (j) t+1 =        T (i) t + 1, T (j) t if T (i) t ≤ T (j) t T (i) t , T (j) t + 1 if T (i) t > T (j) t .
and

T (r) t+1 = T (r) t for r = i, j.
The maximal difference between any two clocks at time t is also called the gap at time t. It is denoted by Gap(t) and is defined by

Gap(t) = max 1≤i≤n T (i) t -min 1≤i≤n T (i) t .
The following theorem gives a maximal value of the gap with any fixed probability. Note that this value is independent of the global time t.

Theorem 6: For all δ ∈ (0, 1), we have

P {Gap(t) ≥ 10 ln(n) -10 ln(δ) + 74} ≤ δ.
Proof. From Relation (1) and Figure 1 of [START_REF] Mocquard | Balanced allocations and global clock in population protocols: An accurate analysis[END_REF] in which we take a = 10 and b = 74, we obtain, for all σ > 0,

P {Gap(t) ≥ 10(1 + σ) ln(n) + 74} ≤ 1/n σ .
Let δ ∈ (0, 1). By taking σ = -ln(δ)/ ln(n), we get σ ln(n) = -ln(δ) and 1/n σ = δ, that is

P {Gap(t) ≥ 10 ln(n) -10 ln(δ) + 74} ≤ δ,
which completes the proof.

The following properties will also be used in the next section. Since at each time only one node has its clock incremented by one we have

n i=1 T t (i) = t.
It follows easily that at each instant t ≥ 0, we have

min 1≤i≤n T (i) t ≤ t n ≤ max 1≤i≤n T (i) t (2)

D. Analysis of the proportion protocol with convergence detection

We now combine all the previous analyses to evaluate the behavior of our proportion protocol with convergence detection. For every n ≥ 2 and for all δ ∈ (0, 1), we introduce the following constants:

• τ 1 = ln(n) -0.5 ln(δ) + 0.55.

• τ 2 = 2.78 ln(n) -2 ln(ε) -7.60 ln(δ) + 11.05.

• τ 3 = 10 ln(n) -10 ln(δ) + 84.99. Constant τ 1 is the constant used in Lemma 2 with δ/3 instead of δ. It is the maximal parallel time for the spreading protocol to converge with probability greater than 1 -δ/3.

Constant τ 2 is the constant used in Theorem 5 with δ/3 instead of δ. It is the maximal parallel time for the proportion protocol to converge with probability greater than 1 -δ/3.

Constant τ 3 is the constant used in Theorem 6 with δ/3 instead of δ. It is the maximal gap obtained with probability greater than 1 -δ/3.

With these constants, we set T max = τ 2 + τ 3 . The following theorem is the main result of the paper. It states that, after time n(T max + τ 1 ), all the nodes have an ε-approximation of γ A and that the spreading is terminated, with probability greater then 1 -δ. More practically, it also states that, if at any instant t a node has its spreading value equal to 1, then all the nodes have an ε-approximation of γ A with probability greater then 1 -δ.

Theorem 7: For every δ ∈ (0, 1) and t ≥ n(T max + τ 1 ), we have

P ω A (C (i) t ) -γ A ≤ ε, S (i) t = 1 ∀i ∈ 1, n ≥ 1 -δ.
Moreover, for every δ ∈ (0, 1) and t ≥ 0, we have

P ω A (C (i) t ) -γ A ≤ ε ∀i ∈ 1, n ∃j such that S (j) t = 1 ≥ 1 -2δ/3.
Proof. The average protocol and the clock protocol both start at time 0 and run independently of each other. We consider first the clock protocol. Let Γ be the first time where two interacting nodes both have their clock value equal to T max -1. Applying Theorem 6 at instant Γ with δ/3 instead of δ, we get P {Gap(Γ) < 10 ln(n) -10 ln(δ/3) + 74} ≥ 1 -δ/3, that is, by definition of τ 3 , P {Gap(Γ) < τ 3 } ≥ 1 -δ/3. By definition of the gap and since that at instant Γ, we have max 1≤i≤n T (i) Γ = T max -1, we obtain

P T max -1 -min 1≤i≤n T (i) Γ < τ 3 ≥ 1 -δ/3, that is, by definition of T max , P min 1≤i≤n T (i) Γ > τ 2 ≥ 1 -δ/3.
From Relation (2) we have min 1≤i≤n T (i) Γ ≤ Γ/n, which leads to

P {Γ > nτ 2 } ≥ 1 -δ/3. (3)
For what concerns the average protocol, to simplify the writing we introduce the events E t defined by

E t = ω A (C (i) t ) -γ A < ε for all i ∈ 1, n .
Applying Theorem 5 with δ/3 instead of δ we get, by definition of τ 2 , for all t ≥ nτ 2 ,

P {E t } ≥ 1 -δ/3.

The random variables C (i)

t and Γ being independent, we have for every t ≥ 0,

P {E Γ+t , Γ > nτ 2 } = ∞ s=nτ2+1 P {E s+t , Γ = s} = ∞ s=nτ2+1 P {E s+t } P {Γ = s} ≥ (1 -δ/3) ∞ s=nτ2+1 P {Γ = s} = (1 -δ/3)P {Γ > nτ 2 } ≥ (1 -δ/3) 2 .
It follows that, for every t ≥ 0,

P {E Γ+t } ≥ P {E Γ+t , Γ > nτ 2 } ≥ (1 -δ/3) 2 . (4)
The starting point of the spreading period is instant Γ + 1. By definition of Γ, instant Γ + 1 is the first time at which exactly two agents have their spreading values equal to 1. More precisely, for every t ≥ 0, we introduce the random variable Y t defined by

Y t = n i=1 S (i) t .
We have Y t = 0 for every t ≤ Γ and Y Γ+1 = 2. The spreading time Θ n is thus the first instant at which the spreading values of all the agents are equal to 1. It is then defined by

Θ n = inf{t ≥ 0 | Y t = n} -(Γ + 1) = inf{t ≥ Γ + 1 | Y t = n} -(Γ + 1).
From Lemma 2 in which we use δ/3 instead of δ, we have, since

Y Γ+1 = 2, P{Θ n < nτ 1 } ≥ 1 -δ/3.
Again, to simplify the writing we introduce the events F t defined by

F t = S (i) t = 1 for all i ∈ 1, n .
By definition of the Boolean S (i) t and of the spreading time Θ n , we have, for all t ≥ 0,

P {F Γ+1+nτ1+t | E Γ+1+nτ1+t } = P {Θ n ≤ nτ 1 + t | E Γ+1+nτ1+t } ≥ 1 -δ/3,
where the last inequality follows from Lemma 2. Unconditioning and using Relation (4), we obtain

P {F Γ+1+nτ1+t , E Γ+1+nτ1+t } = P {F Γ+1+nτ1+t | E Γ+1+nτ1+t } P {E Γ+1+nτ1+t } ≥ (1 -δ/3)(1 -δ/3) 2 = (1 -δ/3) 3 . Recalling that Γ = n i=1 T (i) Γ
and that T (i) t

≤ T max -1, for all t ≥ 0, we get Γ ≤ nT max . Considering instant t = s + nT max -Γ which is positive, finally leads, for all s ≥ 0, to

P {F nTmax+nτ1+s+1 , E nTmax+nτ1+s+1 } ≥ (1 -δ/3) 3 ,
which is equivalent to say that, for all t ≥ n(T max + τ 1 ), we have

P {F t , E t } ≥ (1 -δ/3) 3 ≥ 1 -δ,
which completes the first part of the proof. For the second part of the proof, note that

∃j such that S (j) t = 1 ⇐⇒ Γ ≤ t.
We thus have applying Relation ( 4)

P E t | ∃j such that S (j) t = 1 = P {E t | Γ ≤ t} = P {E Γ+t } ≥ (1 -δ/3) 2 ≥ 1 -2δ/3.
This completes the second part of the proof.

This theorem shows that the convergence is O(ln(n)) and that then number of states needed is equal to

|Q T × Q C × Q S | = 2 (2 3/(4ε) + 1) T max = O (ln(n)/ε).

E. Generalizing the convergence detection mechanism

We now show that our detection mechanism can be applied to any pairwise interaction-based protocol P so that any node of the system can safely detect the instant at which convergence is reached by all the nodes of the system. The only requirement for this mechanism to be applied is that the convergence time of P must be precisely known.

Specifically, let us consider the transition function f of the protocol P such that Relation (1) is replaced by

C (i) t+1 , C (j) t+1 = f C (i) t , C (j) t and C (r) t+1 = C (r) t for r = i, j.
Line 2 of Algorithm 3 is thus replaced by

C (i) , C (j) := f C (i) , C (j) .
The initial value C 0 of vector C t is given and the set of states

Q C of C (i)
t is supposed to be finite. As convergence indicator, we consider the general function ν from (Q C ) n to {0, 1}. We also suppose that we have a general version of Theorem 5 stating that for all δ ∈ (0, 1) and for all t ≥ τ C (n, δ), we have

P {ν(C t ) = 1} ≥ 1 -δ.
(5) Note that by taking τ 2 = τ C (n, δ) and

ν(C t ) = 1 {|ω A (C (i)
t )-γ A |<ε for all i=1,...,n} we arrive to the previous result of Theorem 5.

Under the previous assumptions, the generalization of Theorem 7 is then the following. We set T max = τ C (n, δ/3) + τ 3 .

Theorem 8: For all δ ∈ (0, 1) and for all t ≥ n(T max + τ 1 ) we have

P {ν(C t ) = 1} ≥ 1 -δ.
Moreover, all δ ∈ (0, 1) and t ≥ 0, we have

P ν(C t ) = 1 ∃j such that S (j) t = 1 ≥ 1 -2δ/3.
Proof. By defining E t = {ν(C t ) = 1}, the proof follows exactly the same lines of the proof of Theorem 7, in which τ 2 is replaced by τ C (n, δ/3).

The number of states is |Q

T × Q C × Q S | = 2T max |Q C |.

VI. SIMULATIONS

In this section we first provide simulation results for the spreading, the average, and the clock functions, and then present simulation results for the full protocol.

A. Spreading rumor

This section shows how tight our bound given in Lemma 2 is to our simulation results.

For our purpose, a simulation consists in the following steps: first, all the n nodes are initialized to 0 except for two nodes which are initialized to 1. Then, at each step of the simulation, two nodes are randomly chosen to interact and update their state, by keeping the maximal value of both ones. The simulation stops when all the nodes have their values equal to 1. We have run N independent simulations and have stored and ordered the N values of the spreading times denoted by θ 1 ≤ . . . ≤ θ N . Recall that the spreading time θ i is the total number of interactions to propagate an information to all the nodes of the system. The estimation of the instant τ such that P(Θ n < nτ } ≥ 1 -δ is thus given by the value θ N (1-δ) .

B. Average

For each value of ε, we take m = 3/(4ε) . Next we choose n A = n/4m+n/2 and n B = n-n A . A simulation consists in the steps described in Algorithm 3 and in Section V-B. The simulation stops when the difference between the minimal and the maximal values of the entries of vector C t is less than or equal to 2. We ran N independent simulations and stored the N values of the number of interactions performed which we ordered as θ 1 ≤ . . . ≤ θ N . The estimation of the instant τ such that, for t ≥ τ ,

P |ω A (C (i) t ) -γ A | < ε for all i = 1, . . . , n ≥ 1 -δ is thus given by the value θ N (1-δ) .

C. The clock

For the clock protocol a simulation consists in the steps described in Algorithm 4 and in Section V-C. We start the evaluation of the gap after the first 50n interactions. We then store the gap every 100 interactions. Wa ran x simulations and for each simulation we stored the gap y times. This means that the duration of a simulation is equal to 100y + 50n. The number N of values of the gap obtained is thus N = xy. These N values are stored and reordered as Gap 1 ≤ . . . ≤ Gap N . The estimation of the instant τ such that P {Gap(t) ≥ τ } ≤ δ is thus given by the value Gap N (1-δ) .

D. Optimized protocol derived from simulation

From Relations ( 6) and ( 7) we derive an intuited value T max close to the expression of T max for the proportion protocol with convergence detection. It is given by

T max = τ 2 (n, δ/3, ε) + τ 3 (n, δ/3) = 1.73 ln(n) -1.23 ln(δ) -2 ln(ε) + 1.05.
For different values of n and ε, we ran N = 1000 independent simulations taking δ = 10 -6 , using the value of T max instead of T max . We stored the convergence times θ 1 , . . . , θ N defined, for i = 1, . . . , N , by Figure 8 provides simulation results for different values of ε. The expected value (θ 1 + • • • + θ N )/N , the minimal value min i=1,...,N θ i and the maximal value max i=1,...,N θ i are shown. The most important lesson drawn from all these simulations is that the convergence time of the proportion protocol with the convergence detection mechanism is of the same order of magnitude as the convergence time of the proportion protocol without any detection mechanism. This is a very nice result.

θ i = inf t ≥ 0 | S (i) t = 1, for all ∈ 1, n .

VII. CONCLUSION

In this paper we have presented how we can augment, in the population model, a proportion protocol with a convergence detection mechanism to allow each node of the system to locally detect the instant at which convergence to the sought property is reached. A deep theoretical analysis of the performance of each ingredient of our solution has been presented, and simulation results show the impressively weak impact of our detection mechanism on the convergence time of the proportion protocol. We have also shown the applicability of our convergence detection mechanism to many other pairwise interaction-based protocols.

which gives using ( 9)

B k=-B kα k,t = nλ. (10)
In the same way, again by definition of α k,t , we have

B k=-B ( -λ + k) 2 α k,t = n i=1 C (i) t 2 = C t 2 .
Observing that C t -L 2 = C t 2 -n 2 and using ( 9) and [START_REF] Mocquard | Balanced allocations and global clock in population protocols: An accurate analysis[END_REF], we obtain

B k=-B k 2 α k,t = C t -L 2 + nλ 2 . ( 11
)
Since C t -L 2 is decreasing, using the hypothesis C 0 -L 2 ≤ n/2, we obtain C t -L 2 ≤ n/2 and thus

B k=-B k 2 α k,t ≤ n 2 + nλ 2 . ( 12
)
Let x be defined by

x = B k=0 α k,t .
We then have

-1 k=-B α k,t = n -x and -1 k=-B kα k,t ≤ -1 k=-B -α k,t = -(n -x).
Using [START_REF] Mocquard | Balanced allocations and global clock in population protocols: An accurate analysis[END_REF], we get

B k=1 kα k,t = nλ - -1 k=-B kα k,t ≥ nλ + n -x.
We also have using the two previous inequalities

B k=-B k 2 α k,t = B k=1 k 2 α k,t + -1 k=-B k 2 α k,t ≥ B k=1 kα k,t - -1 k=-B kα k,t ≥ 2(n -x) + nλ.
Combining this inequality with (12) we obtain

2(n -x) + nλ ≤ B k=-B k 2 α k,t ≤ n 2 + nλ 2 .
These two bounds lead to

x ≥ 3n 4 + nλ(1 -λ) 2 .
Since |λ| < 1/2 we have λ(1 -λ) ≥ -3/4, which gives

x = B k=0 α k,t ≥ 3n 8 .
Using the same reasoning to the sum

Let us now introduce the sequence (Φ t ) t≥0 defined by

Φ t = B k=2 k 2 α k,t + -2 k=-B k 2 α k,t .
From Inequality [START_REF] Perron | Using three states for binary consensus on complete graphs[END_REF], we get

Φ t ≤ n 2 + nλ 2 . ( 14
)
We also introduce the sets H + t and H - t defined by

H + t = i ∈ {1, . . . , n} | C (i) t -+ λ ≥ 2 , H - t = i ∈ {1, . . . , n} | C (i) t -+ λ ≤ -2
and we define

H t = H + t ∪ H - t . It is easily checked that Φ t = i∈Ht C (i) t -+ λ 2 . ( 15
)
Let I + t and I - t be the sets defined by

I + t = i ∈ {1, . . . , n} | C (i) t -+ λ ≥ 0 , I - t = i ∈ {1, . . . , n} | C (i) t -+ λ ≤ 0 .
Relations (13) can be rewritten as

|I + t | ≥ 3n 8 and |I - t | ≥ 3n 8 . (16)
Recall that the random variable X t , which is the pair of agents interacting at time t, is uniformly distributed, i.e., for every i, j ∈ {1, . . . , n} with i = j, we have

Pr{X t = (i, j)} = 1 n(n -1)
.

The main way to decrease Φ t is that an agent of H + t interacts with an agent of I - t or that an agent of H - t interacts with an agent of I + t , at time t. So, we consider the events where an agent of H + t interacts with an agent of I - t or where an agent of H - t interacts with an agent of I + t , at time t.

Let E = H + t × I - t ∪ I - t × H + t ∪ H - t × I + t ∪ I + t × H - t
be the set of these interactions. We introduce the notation

G + t = I + t \ H + t = {i ∈ {1, . . . , n} | C (i) t -+ λ ∈ {0, 1}}, G - t = I - t \ H - t = {i ∈ {1, . . . , n} | C (i) t -+ λ ∈ {-1, 0}} and G t = G + t ∪ G - t .
We now consider the difference Φ t -Φ t+1 according to the different interactions taking place at time t. Suppose that X t = (i, j) with i = j. We have the two following different cases. Case 1)

If (i, j) ∈ H + t × G - t ∪ G - t × H + t ∪ H - t × G + t ∪ G + t × H - t ,
and if we set, to simplify the writing,

a = C (i) t -+ λ 1 {i∈Ht} + C (j) t -+ λ 1 {j∈Ht} and b = C (i) t -+ λ 1 {i∈Gt} + C (j) t -+ λ 1 {j∈Gt} ,
we have

Φ t -Φ t+1 = a 2 - a + b -1 {a+b odd} 2 2 1 {i∈Ht+1} - a + b + 1 {a+b odd} 2 2 1 {j∈Ht+1} , (17)
which gives

Φ t -Φ t+1 ≥ a 2 - a + b -1 {a+b odd} 2 2 - a + b + 1 {a+b odd} 2 2 .
Distinguishing successively the cases where a + b is odd and even, we obtain

Φ t -Φ t+1 ≥ a 2 2 -b a + b 2 - 1 {a+b odd} 2 . ( 18
(i, j) ∈ H + t × G - t ∪ G - t × H + t , which means that a ≥ 2. We thus have -b(a + b/2) = a -1/2 ≥ 3/2 and so Φ t -Φ t+1 ≥ a 2 2 ≥ 12a 2 25 . If b = 1 then we necessarily have (i, j) ∈ H - t × G + t ∪ G + t × H - t , which means that a ≤ -2. We thus have -b(a + b/2) = -a -1/2 ≥ 3/2 and so Φ t -Φ t+1 ≥ a 2 2 ≥ 12a 2
25 .

If b = 0 then we distinguish the cases : a is even, |a| = 3 and |a| ≥ 5.

If a is even then, since b = 0, we have, from Relation (18),

Φ t -Φ t+1 ≥ a 2 2 ≥ 12a 2
25 .

If a = 3 then we have, since b = 0, i / ∈ H t+1 (i ∈ G + t+1 ) and j ∈ H t+1 , which gives using Relation (17)

Φ t -Φ t+1 = 9 - 3 + 1 2 2 = 5 ≥ a 2 2 ≥ 12a 2
25 .

If a = -3 then we have, since b = 0, i ∈ H t+1 and j / ∈ H t+1 (i ∈ G - t+1 ), which gives using Relation (17)

Φ t -Φ t+1 = 9 - -3 -1 2 2 = 5 ≥ a 2 2 ≥ 12a 2 25 .
If a is odd and |a| ≥ 5 then, since b = 0, we have,

a 2 ≥ 5 ⇐⇒ a 2 -1 2 ≥ 12a 2 25 ,
which gives, from Relation (18),

Φ t -Φ t+1 ≥ a 2 -1 2 ≥ 12a 2
25 .

Thus we have shown that if

(i, j) ∈ H + t × G - t ∪ G - t × H + t ∪ H - t × G + t ∪ G + t × H - t then Φ t -Φ t+1 ≥ 12a 2 25 . Case 2) If (i, j) ∈ H + t × H - t ∪ H - t × H + t ,
and if we set, to simplify the writing,

a = C (i) t -+ λ and b = C (j) t -+ λ , we have Φ t -Φ t+1 = a 2 + b 2 - a + b -1 {a+b odd} 2 2 1 {i∈Ht+1} - a + b + 1 {a+b odd} 2 2 1 {j∈Ht+1} , (19)
which gives

Φ t -Φ t+1 ≥ a 2 + b 2 - a + b -1 {a+b odd} 2 2 - a + b + 1 {a+b odd} 2 2 .
Distinguishing successively the cases where a + b is odd and even, we obtain

Φ t -Φ t+1 ≥ a 2 2 + b 2 2 -ab - 1 {a+b odd} 2 . ( 20
)
By definition of H + t and H - t we have -ab ≥ 4, so we obtain

Φ t -Φ t+1 ≥ a 2 2 + b 2 2 ≥ 12a 2 25 + 12b 2 25 .
Putting together the cases 1) and 2), we get

E = H + t × G - t ∪ G - t × H + t ∪ H - t × G + t ∪ G + t × H - t ∪ H + t × H - t ∪ H - t × H + t .
All these six sets are disjoints so we have, using the results obtained in cases 1) and 2) and defining

β t,i = C (i) t -+ λ 2 , (i,j)∈E E(Φ t -Φ t+1 | C t , X t = (i, j)) ≥ 12
(i,j)∈E E(Φ t -Φ t+1 | C t , X t = (i, j)) ≥ 12 25   2|I - t | i∈H + t β t,i + 2|I + t | i∈H - t β t,i   ≥ 9n 25 Φ t . (21)
We have now to show that for the other interactions, at worst Φ t does not increase. More precisely, we show that if, at time t, the pair of agent interacting (i, j) / ∈ E then we have Φ t -Φ t+1 ≥ 0. The condition (i, j) / ∈ E at time t can be splitted into the two following disjoint cases : t . We thus have in these cases, from Relation (15),

• (i, j) ∈ H + t × H + t ∪ H - t × H - t , • C (i)
Φ t -Φ t+1 = C (i) t -+ λ 2 + C (j) t -+ λ 2 - C (i) t + C (j) t 2 -+ λ 2 - C (i) t + C (j) t 2 -+ λ 2 = C (i) t - 2 + C (j) t - 2 - C (i) t + C (j) t 2 - 2 - C (i) t + C (j) t 2 - 2 = C t -L 2 -C t+1 -L 2 ≥ 0.
We consider now the second case where C We go on by distinguishing three subcases.

• If a = 2 we have Φ t -Φ t+1 = 0 because, from [START_REF] Alistarh | Fast and exact majority in population protocols[END_REF] the interaction between values 1 and 2 gives 1 and 2 in one hand, and the interaction between -2 and -1 gives -2 and -1, in the other hand. • If a ≥ 3 and a is odd then Φ t -Φ t+1 = a 2 -2 ((a + 1)/2) 2 = (a 2 -2a -1)/2 ≥ 0 for a ≥ 3 > 1 + √ 2. • If a ≥ 4 and a is even then Φ t -Φ t+1 = a 2 -((a + 2)/2)

2 -(a/2) 2 = (a 2 -2a -2)/2 ≥ 0 for a ≥ 4 ≥ 1 + √ 3. We can then write (i,j) / ∈E,i =j E(Φ t -Φ t+1 | C t , X t = (i, j)) ≥ 0.

(

) 22
All the events {X t = (i, j)} having the same probability, with P{X t = (i, j)} = p i,j = 1 {i =j} /(n(n-1)), so using Inequalities Using the fact that ln(1 -x) ≤ -x, for all x ∈ [0, 1), we get for all t ≥ τ , 1 -9 25(n -1) t = e t ln(1-9/(25(n-1))) ≤ e -9t/(25(n-1)) ≤ e -9τ /(25(n-1)) = 16δ 3n .

Using

P{ C t -L ∞ ≥ 3/2} ≤ P{Φ t = 0} = P{Φ t ≥ 4} ≤ δ,
which completes the proof.

Theorem 4 For all δ ∈ (0, 1), if there exists a constant K such that C 0 -L ≤ K then, for every t ≥ n (2 ln(K) + 1.78 ln(n) -7.60 ln(δ) + 2.70), we have

P { C t -L ∞ < 3/2} ≥ 1 -δ.
Proof. We consider first the case where -= 1/2. Since C 0 -L ∞ ≤ C 0 -L ≤ K and since

  Q C = -m, m , and we have |Q C | = 2m + 1; the state set of T (i) t is defined by Q T = 0, T max -1 , and we have |Q T | = T max ; finally the state set of S (i) t is defined by Q S = {0, 1}, and we have |Q S | = 2. Thus the state set of a node i is equal to |Q C | × |Q T | × |Q S |.

3 S 4 return 1 ; 5 : 6 :

 34156 (i) := S (j) := 1;Algorithm Update of T values of two interacting nodes1 Function ω A (i): 2 return m+C (i)2m , S (i) ; Algorithm Output function of the Algorithm V. ANALYSIS

6 1 Figure 1 .

 611 Figure 1. Parallel convergence time of Rumor Spreading as a function of n, with N = 10 4 .

  Spreading parallel time from two knowings δ τ1 with n = 10 4 θ N (1-δ) with n = 10 4 τ1 with n = 10 3 θ N (1-δ) with n = 10 3 τ1 with n = 10 2 θ N (1-δ) with n = 10 2

Figure 2 .

 2 Figure 2. Parallel convergence time of Rumor Spreading as a function of δ, with N = 10 6 .Recall that the convergence parallel time is equal to the convergence time divided by n. Figures1 and 2depict the convergence parallel times θ N (1-δ) /n and τ 1 for different values of δ for the first one, and for different values n for the second one. Both figures shows that the theoretical results are quite close to the simulation ones.

6 1 Figure 3 . 2

 6132 Figure 3. Parallel convergence time of Proportion Computation as a function of n, with N = 10 4 and ε = 0.01.

Figure 4 .

 4 Figure 4. Parallel convergence time of Proportion Computation as a function of δ, with N = 10 5 and ε = 0.01.

Figures 3 ,

 3 Figures 3, 4 and 5 depict the convergence parallel time θ N (1-δ) /n for different values of δ in the first one, for different values of n for the second one, and for different values of ε for the third one. Note that in both the first and the second figure, we have ε = 0.01, that is m = 75. In each figure the values of θ N (1-δ) /n are compared to an intuited value τ 2 (n, δ, ε) close to the expression of τ 2 whose coefficients have been derived from the simulation results, and given by τ 2 (n, δ, ε) = ln(n) -0.5 ln(δ) -2 ln(ε) -1.80.(6)

2 3 Figure 5 .

 235 Figure 5. Parallel convergence time of Proportion Computation as a function of ε, with N = 10 3 and δ = 0.5.

3 1 Figure 6 .Figures 6

 3166 Figure 6. Gap of the clock as a function of n, with N = 10 6 . Figures 6 and 7 depict the gap Gap N (1-δ) for different values of δ for the first one and for different values of n for the second one. In Figure 6 we chose x = 10 and y = 10000 and in Figure Figures 6 we chose x = 100 and y = 10000. In each figure the values of Gap N (1-δ) /n are compared to an intuited value τ 3 (n, δ) close to the expression of τ 3 whose coefficients have been derived from the simulation results, and given by τ 3 (n, δ) = 0.73 ln(n) -0.73 ln(δ) + 1.5.(7)

3 3 Figure 7 .

 337 Figure 7. Gap of the clock as a function of δ, with N = 10 7 .

Figure 8 .

 8 Figure 8. Parallel convergence time of Proportion Computation with Convergence Detection as a function of n, with N = 10 3 .

)

  We consider the cases b = -1, b = 1 and b = 0 separately. If b = -1 then we necessarily have

t

  -+ λ = 1 and j ∈ H + t or C (i) t -+ λ = -1 and j ∈ H - t or i ∈ H + t and C (j) t -+ λ = 1 or i ∈ H - t and C (j) t -+ λ = -1 .To simplify the notation, we definea = |C (i) t -+ λ|1 {i∈Ht} + |C (j)t -+ λ|1 {j∈Ht} .

( 9 25

 9 21) and (22), we obtainE(Φ t -Φ t+1 | C t ) = j E(Φ t -Φ t+1 | C t , X t = (i, j)) t -Φ t+1 | C t , X t = (i, j)) + (i,j) / ∈E,i =j E(Φ t -Φ t+1 | C t , X t = (i, j)) Since E(Φ t | C t ) = Φ t , this leads to E(Φ t+1 |C t ) ≤ 1n) -ln(δ) -4 ln(2) + ln(3)) .

(n - 1 )

 1 (2 ln(K) + ln(n) -ln δ) ≤ n (2 ln(K) + 1.78 ln(n) -7.60 ln(δ) + 2.70) , Theorem 11 givesP{ C t -L ∞ = 1/2} ≤ δ,for t ≥ n (2 ln(K) + 1.78 ln(n) -7.60 ln(δ) + 2.70).

  25

				i∈H + t j∈G -t	β t,i +	12 25	i∈G -t j∈H + t	β t,j +	12 25	i∈H -t j∈G + t	β t,i +	12 25	i∈G + t j∈H -t	β t,j
	+	12 25	i∈H + t j∈H -t	[β t,i + β t,j ] +	12 25	i∈H -t j∈H + t	[β t,i + β t,j ]
													
	=	12 25	 2|G -t |	i∈H + t	β t,i + 2|G + t |	i∈H -t	β t,i + 2|H + t |	i∈H -t	β t,i + 2|H -t |	t i∈H +	β t,i	 .
	Observing that |G -t | + |H -t | = |I -t |, |G + t | + |H + t | = |I + t | and that |I -t | ≥ 3n/8 and |I + t | ≥ 3n/8, we obtain

  the Markov inequality and Relation (14), which gives, since |λ| < 1/2, Φ 0 ≤ n/2 + nλ 2 ≤ 3n/4, we obtainP{Φ t ≥ 4} ≤ E(Φ t )By definition of Φ t , we have Φ t = 0 ⇐⇒ Φ t ≥ 4. Using moreover the fact that |λ| < 1/2, we haveΦ t = 0 =⇒ α t,k = 0, for every k ∈ H t -≤ 1 -λ, for every i = 1, . . . , n =⇒ -3/2 < C -< 3/2, for every i = 1, . . . , n =⇒ C t -L ∞ < 3/2.

		4	≤	16δ 3n	3n 16	= δ.
	=⇒ -1 ≤ C	(i)
		(i)
		(i)
	This leads to
	C

t -+ λ ≤ 1, for every i = 1, . . . , n =⇒ -1 -λ ≤ C t t t -L ∞ ≥ 3/2 =⇒ Φ t = 0 ⇐⇒ Φ t ≥ 4, that is

APPENDIX

This appendix is dedicated to the proof Theorem 4. To prove it, we first recall some other results obtained in [START_REF] Mocquard | Counting with population protocols[END_REF] and [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF]. Next we prove another theorem, Theorem 13, which will be used for the proof of Theorem 4.

Theorem 9: For every 0 ≤ s ≤ t and y ≥ 0, we have

Proof. See [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF]. Lemma 10: The sequence C t -L 2 is decreasing with t. Proof. See [START_REF] Mocquard | Counting with population protocols[END_REF].

Theorem 11: For all δ ∈ (0, 1), if -= 1/2 and if there exists a constant K such that C 0 -L ∞ ≤ K, then, for every t ≥ (n -1) (2 ln(K) + ln(n) -ln(δ)), we have

Proof. see [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF].

Theorem 12: For all δ ∈ (0, 4/5), if there exists a constant K such that K ≥ max n/2, C 0 -L then, for all t ≥ nθ, we have

Proof. see [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF].

The following theorem is an improvement of a previous result obtained in [START_REF] Mocquard | Optimal Proportion Computation with Population Protocols[END_REF] Theorem 13:

Proof. Let λ be defined by

Note that λ is positive in the first case and negative in the second one. In both cases we have |λ| < 1/2 and -λ is the closest integer to .

Since |λ| ≤ 1/2, this means that, for every i = 1, . . . , n, we have

Let B = 1/2 + n/2 . For k ∈ {-B, -B + 1, . . . , B}, we denote by α k,t the number of agents with the value

where the absolute value of a set is its cardinality. It is easily checked that

Moreover we have, by definition of α k,t ,

are integers and since -= 1/2, we have

Consider now the case where -= 1/2. We apply successively Theorem 12 and Theorem 13 replacing δ by δ/2. We introduce the notation

If C 0 -L < n/2 then we have C 0 -L 2 < n/2 and since C t -L 2 is decreasing (see Lemma 10), we get, for all t ≥ 0,

Let us introduce the instant τ defined by

We have, for all t ≥ τ ,

Using the fact that the Markov chain {C t } is homogeneous and applying Theorem 13, we obtain

Putting together these two results gives, for all t ≥ τ ,

The rest of the proof consists in simplifying the expression of τ . We have


Anicet Bart

Benoît Delahaye

Didier Lime

Éric Monfroy

Charlotte Truchet

Reachability in Parametric Interval Markov Chains using Constraints

Parametric Interval Markov Chains (pIMCs) are a specification formalism that extend Markov Chains (MCs) and Interval Markov Chains (IMCs) by taking into account imprecision in the transition probability values: transitions in pIMCs are labeled with parametric intervals of probabilities. In this work, we study the difference between pIMCs and other Markov Chain abstractions models and investigate the two usual semantics for IMCs: once-and-for-all and at-every-step. In particular, we prove that both semantics agree on the maximal/minimal reachability probabilities of a given IMC. We then investigate solutions to several parameter synthesis problems in the context of pIMCs -consistency, qualitative reachability and quantitative reachability -that rely on constraint encodings. Finally, we propose a prototype implementation of our constraint encodings with promising results.

Introduction

Discrete time Markov chains (MCs for short) are a standard probabilistic modeling formalism that has been extensively used in the litterature to reason about software [START_REF] Courcoubetis | The complexity of probabilistic verification[END_REF] and real-life systems [START_REF] Husmeier | Probabilistic Modeling in Bioinformatics and Medical Informatics[END_REF]. However, when modeling real-life systems, the exact value of transition probabilities may not be known precisely. Several formalisms abstracting MCs have therefore been developed. Parametric Markov chains [START_REF] Alur | Parametric real-time reasoning[END_REF] (pMCs for short) extend MCs by allowing parameters to appear in transition probabilities. In this formalism, parameters are variables and transition probabilities may be expressed as polynomials over these variables. A given pMC therefore represents a potentially infinite set of MCs, obtained by replacing each parameter by a given value. pMCs are particularly useful to represent systems where dependencies between transition probabilities are required. Indeed, a given parameter may appear in several dinstinct transition probabilities, therefore requiring that the same value is given to all its occurences. Interval Markov chains [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF] (IMC for short) extend MCs by allowing precise transition probabilities to be replaced by intervals, but cannot represent dependencies between distinct transitions. IMCs have mainly been studied with two distinct semantics interpretation. Under the once-and-for-all semantics, a given IMC represents a potentially infinite number of MCs where transition probabilities are chosen inside the specified intervals while keeping the same underlying graph structure. The at-every-step semantics, which was the original semantics given to IMCs in [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF], does not require MCs to preserve the underlying graph structure of the original IMC but instead allows an "unfolding" of the original graph structure where different probability values may be chosen (inside the specified interval) at each occurence of the given transition.

Model-checking algorithms and tools have been developed in the context of pMCs [START_REF] Dehnert | PROPhESY: A PRObabilistic ParamEter SYnthesis Tool[END_REF][START_REF] Hahn | PARAM: A model checker for parametric Markov models[END_REF][START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF] and IMCs with the once-and-for-all semantics [START_REF]Model Checking of Open Interval Markov Chains[END_REF][START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF]. State of the art tools [START_REF] Dehnert | PROPhESY: A PRObabilistic ParamEter SYnthesis Tool[END_REF] for pMC verification compute a rational function on the parameters that characterizes the probability of satisfying a given property, and then use external tools such as SMT solving [START_REF] Dehnert | PROPhESY: A PRObabilistic ParamEter SYnthesis Tool[END_REF] for computing the satisfying parameter values. For these methods to be viable in practice, the number of parameters used is quite limited. On the other hand, the model-checking procedure for IMCs presented in [START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF] is adapted from machine learning and builds successive refinements of the original IMCs that optimize the probability of satisfying the given property. This algorithm converges, but not necessarilly to a global optimum. It is worth noticing that existing model checking procedures for pMCs and IMCs strongly rely on their underlying graph structure. As a consequence, to the best of our knowledge, no solutions for model-checking IMCs with the at-every-step semantics have been proposed yet.

In this paper, we focus on Parametric interval Markov chains [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] (pIMCs for short), that generalize both IMCs and pMCs by allowing parameters to appear in the endpoints of the intervals specifying transition probabilities, and we provide four main contributions.

First, we formally compare abstraction formalisms for MCs in terms of succinctness: we show in particular that pIMCs are strictly more succinct than both pMCs and pIMCs when equipped with the right semantics. In other words, everything that can be expressed using pMCs or IMCs can also be expressed using pIMCs while the reverse does not hold. Second, we prove that the once-and-forall and the at-every-step semantics are equivalent w.r.t. rechability properties, both in the IMC and in the pIMC settings. Notably, this result gives theoretical backing to the generalization of existing works on the verification of IMCs to the at-every-step semantics. Third, we study the parametric verification of fundamental properties at the pIMC level: consistency, qualitative reachability, and quantitative reachability. Given the expressivity of the pIMC formalism, the risk of producing a pIMC specification that is incoherent and therefore does not model any concrete MC is high. We therefore propose constraint encodings for deciding whether a given pIMC is consistent and, if so, synthesizing parameter values ensuring consistency. We then extend these encodings to qualitative reachability, i.e., ensuring that given state labels are reachable in all (resp. none) of the MCs modeled by a given pIMC. Finally, we focus on the quantitative reachability problem, i.e., synthesizing parameter values such that the probability of reaching given state labels satisfies fixed bounds in at least one (resp. all) MCs modeled by a given pIMC. While consistency and qualitative reachability for pIMCs have already been studied in [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], the constraint encodings we propose in this paper are significantly smaller (linear instead of exponential). To the best of our knowledge, our results provide the first solution to the quantitative reachability problem for pIMCs. Our last contribution is the implementation of all our verification algorithms in a prototype tool that generates the required constraint encodings and can be plugged to any SMT solver for their resolution. Due to space limitation, proofs and detailed examples are given in [START_REF] Bart | Reachability in Parametric Interval Markov Chains using Constraints[END_REF].

Background

In this section we introduce notions and notations that will be used throughout the paper. Given a finite set of variables X = {x 1 , . . . , x k }, we write D x for the domain of the variable x ∈ X and D X for the set of domains associated to the variables in

X. A valuation v over X is a set v = {(x, d)|x ∈ X, d ∈ D x } of elementary valuations (x, d)
where for each x ∈ X there exists a unique pair of the form (x, d) in v. When clear from the context, we write v(x) = d for the value given to variable x according to valuation v. A rational function f over X is a division of two (multivariate) polynomials g 1 and g 2 over X with rational coefficients, i.e., f = g 1 /g 2 . We write Q the set of rational numbers and Q X the set of rational functions over X. The evaluation v(g) of a polynomial g under the valuation v replaces each variable x ∈ X by its value v(x).

An atomic constraint over X is a Boolean expression of the form f (X) g(X), with ∈ {≤, ≥, <, >, =} and f and g two functions over variables in X and constants. A constraint is linear if the functions f and g are linear. A constraint over X is a Boolean combination of atomic constraints over X.

Given a finite set of states S, we write Dist(S) for the set of probability distributions over S, i.e., the set of functions µ : S → [0, 1] such that s∈S µ(s) = 1. We write I for the set containing all the interval subsets of [0, 1]. In the following, we consider a universal set of symbols A that we use for labelling the states of our structures. We call these symbols atomic propositions. We will use Latin alphabet in state context and Greek alphabet in atomic proposition context. Constraints. Constraints are first order logic predicates used to model and solve combinatorial problems [START_REF] Rossi | Handbook of Constraint Programming (Foundations of Artificial Intelligence)[END_REF]. A problem is described with a list of variables, each in a given domain of possible values, together with a list of constraints over these variables. Such problems are then sent to solvers which decide whether the problem is satisfiable, i.e., if there exists a valuation of the variables satisfying all the constraints, and in this case computes a solution. Checking satisfiability of constraint problems is difficult in general, as the space of all possible valuations has a size exponential in the number of variables.

Formally, a Constraint Satisfaction Problem (CSP) is a tuple Ω = (X, D, C) where X is a finite set of variables, D = D X is the set of all the domains associated to the variables from X, and C is a set of constraints over X. We say that a valuation over X satisfies Ω if and only if it satisfies all the constraints in C. We write v(C) for the satisfaction result of the valuation of the constraints C according to v (i.e., true or false). In the following we call CSP encoding a scheme for formulating a given problem into a CSP. The size of a CSP corresponds to the number of variables and atomic constraints appearing in the problem. Note that, in constraint programming, having less variables or less constraints during the encoding does not necessarily imply faster solving time of the problems. Discrete Time Markov Chains. A Discrete Time Markov Chain (DTMC or MC for short) is a tuple M = (S, s 0 , p, V ), where S is a finite set of states containing the initial state s 0 , V : S → 2 A is a labelling function, and p : S → Dist(S) is a probabilistic transition function. We write MC for the set containing all the discrete time Markov chains.

A Markov Chain can be seen as a directed graph where the nodes correspond to the states of the MC and the edges are labelled with the probabilities given by the transition function of the MC. In this representation, a missing transition between two states represents a transition probability of zero. As usual, given a MC M, we call a path of M a sequence of states obtained from executing M, i.e., a sequence ω = s 1 , s 2 , . . . s.t. the probability of taking the transition from s i to s i+1 is strictly positive, p(s i )(s i + 1) > 0, for all i. A path ω is finite iff it belongs to S * , i.e., it represents a finite sequence of transitions from M.

Example 1. Figure 1 illustrates the Markov chain M 1 = (S, s 0 , p, V ) ∈ MC where S = {s 0 , s 1 , s 2 , s 3 , s 4 }, the atomic proposition are restricted to {α, β}, the initial state is s 0 , and the labelling function V corresponds to {(s 0 , ∅), (s 1 , α), (s 2 , β), (s 3 , {α, β}), (s 4 , α)}. The sequences of states (s 0 , s 1 , s 2 ), (s 0 , s 2 ), and (s 0 , s 2 , s 2 , s 2 ), are three (finite) paths from the initial state s 0 to the state s 2 .

Reachability. A Markov chain M defines a unique probability measure P M over the paths from M. According to this measure, the probability of a finite path ω = s 0 , s 1 , . . . , s n in M is the product of the probabilities of the transitions executed along this path, i.e., P M (ω) = p(s 0 )(s 1 ) • p(s 1 )(s 2 ) • . . . • p(s n-1 )(s n ). This distribution naturally extends to infinite paths (see [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF]) and to sequences of states over S that are not paths of M by giving them a zero probability.

Given a MC M, the overall probability of reaching a given state s from the initial state s 0 is called the reachability probability and written P M s0 (♦s) or P M (♦s) when clear from the context. This probability is computed as the sum of the probabilities of all finite paths starting in the initial state and reaching this state for the first time. Formally, let reach s0 (s) = {ω ∈ S * | ω = s 0 , . . . s n with s n = s and s i = s ∀0 ≤ i < n} be the set of such paths. We then define P M (♦s) = ω∈reachs 0 (s) P M (ω) if s = s 0 and 1 otherwise. This notation naturally extends to the reachability probability of a state s from a state t that is not s 0 , written P M t (♦s) and to the probability of reaching a label α ⊆ A written P M s0 (♦α). In the following, we say that a state s (resp. a label α ⊆ A) is reachable in M iff the reachability probability of this state (resp. label) from the initial state is strictly positive. In Figure 1 the probability of the path (s 0 , s 2 , s 1 , s 1 , s 3 ) is 0.3 • 0.5 • 0.5 • 0.5 = 0.0375 and the probability of reaching the state s

1 is P M1 (♦s 1 ) = p(s 0 )(s 1 ) + Σ +∞ i=0 p(s 0 )(s 2 )•p(s 2 )(s 2 ) i •p(s 2 )(s 1 ) = p(s 0 )(s 1 ) + p(s 0 )(s 2 )•p(s 2 )(s 1 )•(1/(1 -p(s 2 )(s 2 ))) = 1.
Furthermore, the probability of reaching β corresponds to the probability of reaching the state s 2 .

Markov Chains Abstractions

Modelling an application as a Markov Chain requires knowing the exact probability for each possible transition of the system. However, this can be difficult to compute or to measure in the case of a real-life application (e.g., precision errors, limited knowledge). In this section, we start with a generic definition of Markov chain abstraction models. Then we recall three abstraction models from the literature, respectively pMC, IMC, and pIMC, and finally we present a comparison of these existing models in terms of succinctness. A Markov chain Abstraction Model is a specification theory for MCs. It consists in a set of abstract objects, called specifications, each of which representing a (potentially infinite) set of MCs -implementations -together with a satisfaction relation defining the link between implementations and specifications. As an example, consider the powerset of MC (i.e., the set containing all the possible sets of Markov chains). Clearly, (2 MC , ∈) is a Markov chain abstraction model, which we call the canonical abstraction model. This abstraction model has the advantage of representing all the possible sets of Markov chains but it also has the disadvantage that some Markov chain abstractions are only representable by an infinite extension representation. Indeed, recall that there exists subsets of [0, 1] ⊆ R which cannot be represented in a finite space (e.g., the Cantor set [START_REF] Cantor | ber unendliche, lineare punktmannigfaltigkeiten v (on infinite, linear point-manifolds)[END_REF]). We now present existing MC abstraction models from the literature.

Existing MC Abstraction Models

Parametric Markov Chain is a MC abstraction model from [START_REF] Alur | Parametric real-time reasoning[END_REF] where a transition can be annotated by a rational function over parameters. We write pMC for the set containing all the parametric Markov chains.

Definition 2 (Parametric Markov Chain). A Parametric Markov Chain ( pMC for short) is a tuple I = (S, s 0 , P, V, Y ) where S, s 0 , and V are defined as for MCs, Y is a set of variables (parameters), and P : S × S → Q Y associates with each potential transition a parameterized probability.

Let M = (S, s 0 , p, V ) be a MC and I = (S, s 0 , P, V, Y ) be a pMC. The satisfaction relation |= p between MC and pMC is defined by M |= p I iff there exists a valuation v of Y s.t. p(s)(s ) equals v(P (s, s )) for all s, s in S.

Example 3. Figure 2 shows a pMC I = (S, s 0 , P, V, Y ) where S, s 0 , and V are similar to the same entities in the MC M from Figure 1, the set of variable Y contains only one variable p, and the parametric transitions in P are given by the edge labelling (e.g., P (s 0 , s 1 ) = 0.7, P (s 1 , s 3 ) = p, and P (s 2 , s 2 ) = 1 -p). Note that the pMC I is a specification containing the MC M from Figure 1.

Interval Markov Chains extend MCs by allowing to label transitions with intervals of possible probabilities instead of precise probabilities. We write IMC for the set containing all the interval Markov chains. [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF]). An Interval Markov Chain ( IMC for short) is a tuple I = (S, s 0 , P, V ), where S, s 0 , and V are defined as for MCs, and P : S × S → I associates with each potential transition an interval of probabilities. Example 4. Figure 3 illustrates IMC I = (S, s 0 , P, V ) where S, s 0 , and V are similar to the MC given in Figure 1. By observing the edge labelling we see that P (s 0 , s 1 ) = [0, 1], P (s 1 , s 1 ) = [0.5, 1], and P (s 3 , s 3 ) = [START_REF] Alur | Parametric real-time reasoning[END_REF][START_REF] Alur | Parametric real-time reasoning[END_REF]. On the other hand, the intervals of probability for missing transitions are reduced to [0, 0], e.g., P (s 0 , s 0 ) = [0, 0], P (s 0 , s 3 ) = [0, 0], P (s 1 , s 4 ) = [0, 0].

Definition 3 (Interval Markov Chain

In the literature, IMCs have been mainly used with two distinct semantics: at-every-step and once-and-for-all. Both semantics are associated with distinct satisfaction relations which we now introduce.

The once-and-for-all IMC semantics ( [START_REF] Dehnert | PROPhESY: A PRObabilistic ParamEter SYnthesis Tool[END_REF][START_REF] Wongpiromsarn | Tulip: A software toolbox for receding horizon temporal logic planning[END_REF][START_REF] Puggelli | Polynomial-time verification of pctl properties of mdps with convex uncertainties[END_REF]) is alike to the semantics for pMC, as introduced above. The associated satisfaction relation |= o I is defined as follows: A MC M = (T, t 0 , p, V M ) satisfies an IMC I = (S, s 0 , P, V I ) iff (T, t 0 , V M ) = (S, s 0 , V I ) and for all reachable state s and all state s ∈ S, p(s)(s ) ∈ P (s, s ). In this sense, we say that MC implementations using the onceand-for-all semantics need to have the same structure as the IMC specification.

On the other hand, the at-every-step IMC semantics, first introduced in [START_REF] Jonsson | Specification and refinement of probabilistic processes[END_REF], operates as a simulation relation based on the transition probabilities and state labels, and therefore allows MC implementations to have a different structure I is defined as follows: A MC M = (T, t 0 , p, V M ) satisfies an IMC I = (S, s 0 , P, V I ) iff there exists a relation R ⊆ T × S such that (t 0 , s) ∈ R and whenever (t, s) ∈ R, we have 1. the labels of s and t correspond:

s 0 s 1 α s 2 β s 3 α, β s 4 α [0, 1] [0, 1] [q, 1] [0.3, q] [0, p] [0.2, p] [0, 0.5] 1 [0, 0.5] [0.5, p]
V M (t) = V I (s), 2. there exists a correspondence function δ : T → (S → [0, 1]) s.t. (a) ∀t ∈ T if p(t)(t ) > 0 then δ(t ) is a distribution on S (b) ∀s ∈ S : (Σ t ∈T p(t)(t ) • δ(t )(s )) ∈ P (s, s ), and (c) ∀(t , s ) ∈ T × S, if δ(t )(s ) > 0, then (t , s ) ∈ R. By construction, it is clear that |= a I is more general than |= o I , i.e.
, that whenever M |= o I I, we also have M |= a I I. The reverse is obviously not true in general, even when the underlying graphs of M and I are isomorphic (see [START_REF] Bart | Reachability in Parametric Interval Markov Chains using Constraints[END_REF] for details).

Example 5 (Example 4 continued).

Consider the MC M 1 with state space S from Figure 1 and the MC M 2 with state space T from Figure 4. They both satisfy the IMC I with state space S given in Figure 3. Furthermore, M 1 satisfies I with the same structure. On the other hand, for the MC M 2 given in Figure 4, the state s 3 from I has been "split" into two states t 3 and t 3 in M 2 and the state t 1 from M 2 "aggregates" states s 1 and s 4 in I. The relation R ⊆ T × S containing the pairs (t 0 , s 0 ), (t 1 , s 1 ), (t 1 , s 4 ), (t 2 , s 2 ), (t 3 , s 3 ), and (t 3 , s 3 ) is a satisfaction relation between M 2 and I.

Parametric Interval Markov Chains, as introduced in [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF], abstract IMCs by allowing (combinations of) parameters to be used as interval endpoints in IMCs. Under a given parameter valuation the pIMC yields an IMC as introduced above. pIMCs therefore allow the representation, in a compact way and with a finite structure, of a potentially infinite number of IMCs. Note that one parameter can appear in several transitions at once, requiring the associated transition probabilities to depend on one another. Let Y be a finite set of parameters and v be a valuation over Y . By combining notations used for IMCs and pMCs the set I(Q Y ) contains all parametrized intervals over [0, 1], and for all

I = [f 1 , f 2 ] ∈ I(Q Y ), v(I) denotes the interval [v(f 1 ), v(f 2 )] if 0 ≤ v(f 1 ) ≤ v(f 2 ) ≤ 1
and the empty set otherwise4 . We write pIMC for the set containing all the parametric interval Markov chains.

Definition 4 (Parametric Interval Markov Chain [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]). A Parametric Interval Markov Chain ( pIMC for short) is a tuple P = (S, s 0 , P, V, Y ), where S, s 0 , V and Y are defined as for pMCs, and P : S × S → I(Q Y ) associates with each potential transition a (parametric) interval.

In [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] the authors introduced pIMCs where parametric interval endpoints are limited to linear combination of parameters. In this paper we extend the pIMC model by allowing rational functions over parameters as endpoints of parametric intervals. Given a pIMC P = (S, s 0 , P, V, Y ) and a valuation v, we write v(P) for the IMC (S, s 0 , P v , V ) obtained by replacing the transition function P from P with the function P v : S × S → I defined by P v (s, s ) = v(P (s, s )) for all s, s ∈ S. The IMC v(P) is called an instance of pIMC P. Finally, depending on the semantics chosen for IMCs, two satisfaction relations can be defined between MCs and pIMCs. ). Example 6. Consider the pIMC P = (S,0, P, V, Y ) given in Figure 5. The set of states S and the labelling function are the same as in the MC and the IMC presented in Figures 1 and3 respectively. The set of parameters Y has two elements p and q. Finally, the parametric intervals from the transition function P are given by the edge labelling (e.g., P (s 1 , s 3 ) = [0.3, q], P (s 2 , s 4 ) = [0, 0.5], and P (s 3 , s 3 ) = [START_REF] Alur | Parametric real-time reasoning[END_REF][START_REF] Alur | Parametric real-time reasoning[END_REF]). Note that the IMC I from Figure 3 is an instance of P (by assigning the value 0.6 to the parameter p and 0.5 to q). Furthermore, as said in Example 5, the Markov Chains M 1 and M 2 (from Figures 1 and4 respectively) satisfy I, therefore M 1 and M 2 satisfy P.

In the following, we consider that the size of a pMC, IMC, or pIMC corresponds to its number of states plus its number of transitions not reduced to 0, [0, 0] or ∅. We will also often need to consider the predecessors (Pred), and the successors (Succ) of some given states. Given a pIMC with a set of states S, a state s in S, and a subset S of S, we write:

-Pred(s) = {s ∈ S | P (s , s) / ∈ {∅, [0, 0]}} -Succ(s) = {s ∈ S | P (s, s ) / ∈ {∅, [0, 0]}} -Pred(S ) = s ∈S Pred(s ) -Succ(S ) = s ∈S Succ(s )

Abstraction Model Comparisons

IMC, pMC, and pIMC are three Markov chain Abstraction Models. In order to compare their expressiveness and compactness, we introduce the comparison operators and ≡. Let (L 1 , |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models containing respectively L 1 and L 2 . We say that L 1 is entailed by L 2 , written L 1 L 2 , iff all the MCs satisfying L 1 satisfy L 2 modulo bisimilarity. (i.e., ∀M |= 1 L 1 , ∃M |= 2 L 2 s.t. M is bisimilar to M ). We say that L 1 is (semantically) equivalent to L 2 , written L 1 ≡ L 2 , iff L 1 L 2 and L 2 L 1 . Definition 5 introduces succinctness based on the sizes of the abstractions. Definition 5 (Succinctness). Let (L 1 , |= 1 ) and (L 2 , |= 2 ) be two Markov chain abstraction models. L 1 is at least as succinct as L 2 , written L 1 ≤ L 2 , iff there exists a polynomial p such that for every L 2 ∈ L 2 , there exists 5 Thus, L 1 is strictly more succinct than L 2 , written

L 1 ∈ L 1 s.t. L 1 ≡ L 2 and |L 1 | ≤ p(|L 2 |).
L 1 < L 2 , iff L 1 ≤ L 2 and L 2 ≤ L 1 .
We start with a comparison of the succinctness of the pMC and IMC abstractions. Since pMCs allow the expression of dependencies between the probabilities assigned to distinct transitions while IMCs allow all transitions to be independant, it is clear that there are pMCs without any equivalent IMCs (regardless of the IMC semantics used), therefore (IMC, |= o I ) ≤ pMC and (IMC, |= a I ) ≤ pMC (see [START_REF] Bart | Reachability in Parametric Interval Markov Chains using Constraints[END_REF] for details). On the other hand, IMCs imply that transition probabilities need to satisfy linear inequalities in order to fit given intervals. However, these types of constraints are not allowed in pMCs. It is therefore easy to exhibit IMCs that, regardless of the semantics considered, do not have any equivalent pMC specification. As a consequence, pMC ≤ (IMC, |= o I ) and pMC ≤ (IMC, |= a I ). We now compare pMCs and IMCs to pIMCs. Recall that the pIMC model is a Markov chain abstraction model allowing to declare parametric interval transitions, while the pMC model allows only parametric transitions (without intervals), and the IMC model allows interval transitions without parameters. Clearly, any pMC and any IMC can be translated into a pIMC with the right semantics (once-and-for-all for pMCs and the chosen IMC semantics for IMCs). This means that (pIMC, |= o pI ) is more succinct than pMC and pIMC is more succinct than IMC for both semantics. Furthermore, since pMC and IMC are not comparable due to the above results, we have that the pIMC abstraction model is strictly more succinct than the pMC abstraction model and than the IMC abstraction model with the right semantics. Our comparison results are presented in Proposition 1. Further explanations and examples are given in [START_REF] Bart | Reachability in Parametric Interval Markov Chains using Constraints[END_REF]. Note that (pMC, |= p ) ≤ (IMC, |= o I ) could be achieved by adding unary constraints on the parameters of a pMC, which is not allowed here. However, this would not have any impact on our other results.

Qualitative Properties

As seen above, pIMCs are a succinct abstraction formalism for MCs. The aim of this section is to investigate qualitative properties for pIMCs, i.e., properties that can be evaluated at the specification (pIMC) level, but that entail properties on its MC implementations. pIMC specifications are very expressive as they allow the abstraction of transition probabilities using both intervals and parameters. Unfortunately, as it is the case for IMCs, this allows the expression of incorrect specifications. In the IMC setting, this is the case either when some intervals are ill-formed or when there is no probability distribution matching the interval constraints of the outgoing transitions of some reachable state. In this case, no MC implementation exists that satisfies the IMC specification. Deciding whether an implementation that satisfies a given specification exists is called the consistency problem. In the pIMC setting, the consistency problem is made more complex because of the parameters which can also induce inconsistencies in some cases. One could also be interested in verifying whether there exists an implementation that reaches some target states/labels, and if so, propose a parameter valuation ensuring this property. Both the consistency and the consistent reachability problems have already been investigated in the IMC and pIMC setting [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF][START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]. In this section, we briefly recall these problems and propose new solutions based on CSP encodings. Our encodings are linear in the size of the original pIMCs whereas the algorithms from [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF][START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] are exponential.

ρ 0 ρ 1 ρ 2 ρ 3 ρ 4 πp ∈ [0, 1] πq ∈ [0, 1] θ 1 0 θ 2 0 θ 1 1 θ 3 1 θ 1 2 θ 2 2

Existential Consistency

A pIMC P is existential consistent iff there exists a MC M satisfying P (i.e., there exists a MC M satisfying an IMC I instance of P). As seen in Section 2, pIMCs are equipped with two semantics: once-and-for-all (|= o pI ) and at-everystep (|= a pI ). Recall that |= o pI imposes that the underlying graph structure of implementations needs to be isomorphic to the graph structure of the corresponding specification. In contrast, |= a pI allows implementations to have a different graph structure. It therefore seems that some pIMCs could be inconsistent w.r.t |= o pI while being consistent w.r.t |= a pI . On the other hand, checking the consistency w.r.t |= o pI seems easier because of the fixed graph structure. In [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF], the author firstly proved that both semantics are equivalent w.r.t. existential consistency, and proposed a CSP encoding for verifying this property which is exponential in the size of the pIMC. Based on this result of semantics equivalence w.r.t. existential consistency from [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF] we propose a new CSP encoding, written C ∃c , for verifying the existential consistency property for pIMCs.

Let P = (S,s 0 ,P,V,Y ) be a pIMC, we write C ∃c (P) for the CSP produced by C ∃c according to P. Any solution of C ∃c (P) will correspond to a MC satisfying P. In C ∃c (P), we use one variable π p with domain [0, 1] per parameter p in Y ; one variable θ s s with domain [0, 1] per transition (s, s ) in {{s}×Succ(s) | s ∈ S}; and one Boolean variable ρ s per state s in S. These Boolean variables will indicate for each state whether it appears in the MC solution of the CSP (i.e., in the MC satisfying the pIMC P). For each state s ∈ S, Constraints are as follows:

(1) ρ s , if s = s 0 (3) ¬ρ s ⇔ Σ s ∈Pred(s)\{s} θ s s = 0, if s = s 0 (5) ρ s ⇒ θ s s ∈ P (s, s ), for all s ∈ Succ(s)

(2) ρ s ⇔ Σ s ∈Succ(s) θ s s = 1 (4) ¬ρ s ⇔ Σ s ∈Succ(s) θ s s = 0
Recall that given a pIMC P the objective of the CSP C ∃c (P) is to construct a MC M satisfying P. Constraint (1) states that the initial state s 0 appears in M. Constraint (3) ensures that for each non-initial state s, variable ρ s is set to false iff s is not reachable from its predecessors. Constraint (2) ensures that if a state s appears in M, then its outgoing transitions form a probability distribution. On the contrary, Constraint (4) propagates non-appearing states (i.e., if a state s does not appear in M then all its outgoing transitions are set to zero). Finally, Constraint [START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF] states that, for all appearing states, the outgoing transition probabilities must be selected inside the specified intervals.

Example 7. Consider the pIMC P given in Figure 5. Figure 6 describes the variables in C ∃c (P): one variable per transition (e.g., θ 1 0 , θ 2 0 , θ 1 1 ), one Boolean variable per state (e.g., ρ 0 , ρ 1 ), and one variable per parameter (π p and π q ). The following constraints correspond to the Constraints (2), ( 3), (4), and (5) generated by our encoding C ∃c for the state 2 of P:

¬ρ 2 ⇔ θ 2 0 = 0 ¬ρ 2 ⇔ θ 1 2 + θ 2 2 + θ 4 2 = 0 ρ 2 ⇔ θ 1 2 + θ 2 2 + θ 4 2 = 1 ρ 2 ⇒ 0 ≤ θ 1 2 ≤ π p ρ 2 ⇒ 0.2 ≤ θ 2 2 ≤ π p ρ 2 ⇒ 0 ≤ θ 4 2 ≤ 0.5
Finally, Figure 7 describes a solution for the CSP C ∃c (P). Note that given a solution of a pIMC encoded by C ∃c , one can construct a MC satisfying the given pIMC by keeping all the states s s.t. ρ s is equal to true and considering the transition function given by the probabilities in the θ s s variables. We now show that our encoding works as expected. Proposition 2. A pIMC P is existential consistent iff C ∃c (P) is satisfiable.

Our existential consistency encoding is linear in the size of the pIMC instead of exponential for the encoding from [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] which enumerates the powerset of the states in the pIMC resulting in deep nesting of conjunctions and disjunctions.

Qualitative Reachability

Let P = (S, s 0 , P, V, Y ) be a pIMC and α ⊆ A be a state label. We say that α is existential reachable in P iff there exists an implementation M of P where α is reachable (i.e., P M (♦α) > 0). In a dual way, we say that α is universal reachable in P iff α is reachable in any implementation M of P. As for existential consistency, we use a result from [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF] that states that both pIMC semantics are equivalent w.r.t. existential (and universal) reachability. We therefore propose a new CSP encoding, written C ∃r , that extends C ∃c , for verifying these properties. Formally, CSP C ∃r (P) = (X ∪ X , D ∪ D , C ∪ C ) is such that (X, D, C) = C ∃c (P), X contains one integer variable ω s with domain [0, |S|] per state s in S, D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

(6) ω s = 1, if s = s 0 (7) ω s = 1, if s = s 0 (8) ρ s ⇔ (ω s = 0) (9) ω s > 1 ⇒ s ∈Pred(s)\{s} (ω s = ω s + 1) ∧ (θ s s > 0), if s = s 0 (10) ω s = 0 ⇔ s ∈Pred(s)\{s} (ω s = 0) ∨ (θ s s = 0), if s = s 0 Recall first that CSP C ∃c (P ) constructs a Markov chain M satisfying P.
Informally, for each state s in M the Constraints ( 6), ( 7), ( 9) and [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF] in C ∃r ensure that ω s = k iff there exists in M a path from the initial state to s of length k -1 with non zero probability; and state s is not reachable in M from the initial state s 0 iff ω s equals to 0. Finally, Constraint (8) enforces the Boolean reachability indicator variable ρ s to bet set to true iff there exists a path with non zero probability in M from the initial state s 0 to s (i.e., ω s = 0).

Let S α be the set of states from P labeled with α. C ∃r (P) therefore produces a Markov chain satisfying P where reachable states s are such that ρ s = true. As a consequence, α is existential reachable in P iff C ∃r (P) admits a solution such that s∈Sα ρ s ; and α is universal reachable in P iff C ∃r (P) admits no solution such that s∈Sα ¬ρ s . This is formalised in the following proposition.

Proposition 3. Let P = (S, s 0 , P, V, Y ) be a pIMC, α ⊆ A be a state label, S α = {s | V (s) = α}, and (X, D, C) be the CSP C ∃r (P).

-CSP (X, D, C ∪ s∈Sα ρ s ) is satisfiable iff α is existential reachable in P -CSP (X, D, C ∪ s∈Sα ¬ρ s ) is unsatisfiable iff α is universal reachable in P As for the existential consistency problem, we have an exponential gain in terms of size of the encoding compared to [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]: the number of constraints and variables in C ∃r is linear in terms of the size of the encoded pIMC. Remark. In C ∃r Constraints (3) inherited from C ∃c are entailed by Constraints [START_REF] De Moura | Z3: An efficient SMT solver[END_REF] and [START_REF] Delahaye | Consistency for parametric interval Markov chains[END_REF] added to C ∃r . Thus, in a practical approach one may ignore Constraints (3) from C ∃c if they do not improve the solver performances.

Quantitative Properties

We now move to the verification of quantitative reachability properties in pIMCs. Quantitative reachability has already been investigated in the context of pMCs and IMCs with the once-and-for-all semantics. Due to the complexity of allowing implementation structures to differ from the structure of the specifications, quantitative reachability in IMCs with the at-every-step semantics has, to the best of our knowledge, never been studied. In this section, we propose our main theoretical contribution: a theorem showing that both IMC semantics are equivalent with respect to quantitative reachability, which allows the extension of all results from [START_REF] Wongpiromsarn | Tulip: A software toolbox for receding horizon temporal logic planning[END_REF][START_REF] Benedikt | LTL model checking of interval Markov chains[END_REF] to the at-every-step semantics. Based on this result, we also extend the CSP encodings introduced in Section 4 in order to solve quantitative reachability properties on pIMCs regardless of their semantics.

Equivalence of |= o

I and |= a I w.r.t quantitative reachability Given an IMC I = (S, s 0 , P, V ) and a state label α ⊆ A, a quantitative reachability property on I is a property of the type P I (♦α)∼p, where 0 < p < 1 and ∼ ∈ {≤, <, >, ≥}. Such a property is verified iff there exists an MC M satisfying I (with the chosen semantics) such that P M (♦α)∼p.

As explained above, all existing techniques and tools for verifying quantitative reachability properties on IMCs only focus on the once-and-for-all semantics. Indeed, in this setting, quantitative reachability properties are easier to compute because the underlying graph structure of all implementations is known. However, to the best of our knowledge, there are no works addressing the same problem with the at-every-step semantics or showing that addressing the problem in the once-and-for-all setting is sufficiently general. The following theorem fills this theoretical gap by proving that both semantics are equivalent w.r.t quantitative reachability. In other words, for all MC M such that M |= a I I and all state label α, there exist MCs M ≤ and M ≥ such that M ≤ |= o I I, M ≥ |= o I I and P M ≤ (♦α) ≤ P M (♦α) ≤ P M ≥ (♦α). This is formalized in the following theorem.

Theorem 1. Let I = (S, s 0 , P, V ) be an IMC, α ⊆ A be a state label, ∼ ∈ {≤, < , >, ≥} and 0 < p < 1. I satisfies P I (♦α)∼p with the once-and-for-all semantics iff I satisfies P I (♦α)∼p with the at-every-step semantics.

The proof is constructive (see [START_REF] Bart | Reachability in Parametric Interval Markov Chains using Constraints[END_REF]): we use the structure of the relation R from the definition of |= a I in order to build the MCs M ≤ and M ≥ .

Constraint Encodings

Note that the result from Theorem 1 naturally extends to pIMCs. We therefore exploit this result to construct a CSP encoding for verifying quantitative reachability properties in pIMCs. As in Section 4, we extend the CSP C ∃c , that produces a correct MC implementation for the given pIMC, by imposing that this MC implementation satisfies the given quantitative reachability property. In order to compute the probability of reaching state label α at the MC level, we use standard techniques from [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF] that require the partitioning of the state space into three sets S , S ⊥ , and S ? that correspond to states reaching α with probability 1, states from which α cannot be reached, and the remaining states, respectively. Once this partition is chosen, the reachability probabilities of all states in S ? are computed as the unique solution of a linear equation system (see [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF], Theorem 10.19, p.766). We now explain how we identify states from S ⊥ , S and S ? and how we encode the linear equation system, which leads to the resolution of quantitative reachability. Let P = (S, s 0 , P, V, Y ) be a pIMC and α ⊆ A be a state label. We start by setting S = {s | V (s) = α}. We then extend C ∃r (P) in order to identify the set S ⊥ . Let C ∃r (P, α) = (X ∪ X , D ∪ D , C ∪ C ) be such that (X, D, C) = C ∃r (P), X contains one Boolean variable λ s and one integer variable α s with domain [0, |S|] per state s in S, D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

(11) α s = 1, if α = V (s) (12) α s = 1, if α = V (s) (13) λ s ⇔ (ρ s ∧ (α s = 0)) (14) α s > 1 ⇒ s ∈Succ(s)\{s} (α s = α s + 1) ∧ (θ s s > 0), if α = V (s) (15) α s = 0 ⇔ s ∈Succ(s)\{s} (α s = 0) ∨ (θ s s = 0), if α = V (
s) Note that variables α s play a symmetric role to variables ω s from C ∃r : instead of indicating the existence of a path from s 0 to s, they characterize the existence of a path from s to a state labeled with α. In addition, due to Constraint (13), variables λ s are set to true iff there exists a path with non zero probability from the initial state s 0 to a state labeled with α passing by s. Thus, α cannot be reached from states s.t. λ s = false. Therefore, S ⊥ = {s | λ s = false}.

Finally, we encode the equation system from [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF] in a last CSP encoding that extends C ∃r . Let C ∃r (P, α) = (X ∪ X , D ∪ D , C ∪ C ) be such that (X, D, C) = C ∃r (P, α), X contains one variable π s per state s in S with domain [0, 1], D contains the domains of these variables, and C is composed of the following constraints for each state s ∈ S:

(

) ¬λ s ⇒ π s = 0 (17) λ s ⇒ π s = 1, if α = V (s) (18) λ s ⇒ π s = Σ s ∈Succ(s) π s θ s s , if α = V 16
s) As a consequence, variables π s encode the probability with which state s eventually reaches α when s is reachable from the initial state and 0 otherwise.

Let p ∈ [0, 1] ⊆ R be a probability bound. Adding the constraint π s0 ≤ p (resp. π s0 ≥ p) to the previous C ∃r encoding allows to determine if there exists a MC M |= a pI P such that P M (♦α) ≤ p (resp ≥ p). Formally, let ∼ ∈ {≤, <, ≥, >} be a comparison operator, we write ∼ for its negation (e.g., ≤ is >). This leads to the following theorem.

Theorem 2. Let P = (S, s 0 , P, V, Y ) be a pIMC, α ⊆ A be a label, p ∈ [0, 1], ∼ ∈ {≤, <, ≥, >} be a comparison operator, and (X, D, C) be C ∃r (P, α):

-CSP (X, D, C ∪ (π s0 ∼ p)) is satisfiable iff ∃M |= a pI P s.t. P M (♦α) ∼ p -CSP (X, D, C ∪ (π s0 ∼ p)) is unsatisfiable iff ∀M |= a pI P: P M (♦α) ∼ p 6

Prototype Implementation

Our results have been implemented in a prototype tool6 which generates the above CSP encodings, and CSP encodings from [START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] as well. Given a pIMC in a text format inspired from [START_REF] Wongpiromsarn | Tulip: A software toolbox for receding horizon temporal logic planning[END_REF], our tool produces the desired CSP as a SMT instance with the QF NRA logic (Quantifier Free Non linear Real-number Arithmetic). This instance can then be fed to any solver accepting the SMT-LIB format with QF NRA logic [START_REF] Barrett | The Satisfiability Modulo Theories Library (SMT-LIB)[END_REF]. For our benchmarks, we chose Z3 [START_REF] De Moura | Z3: An efficient SMT solver[END_REF] (latest version: 4.5.0). QF NRA does not deal with integer variables. In practice, logics mixing integers and reals are harder than those over reals only. Thus we obtained better results by encoding integer variables into real ones. In our implementations each integer variable x is declared as a real variable whose real domain bounds are its original integer domain bounds; we also add the constraint x < 1 ⇒ x = 0. Since we only perform incrementation of x this preserves the same set of solutions.

In order to evaluate our prototype, we extend the nand model from [START_REF] Norman | Evaluating the reliability of NAND multiplexing with PRISM[END_REF] 7 . The original MC nand model has already been extended as a pMC in [START_REF] Dehnert | PROPhESY: A PRObabilistic ParamEter SYnthesis Tool[END_REF], where the authors consider a single parameter p for the probability that each of the N nand gates fails during the multiplexing. We extend this model to pIMC by considering intervals for the probability that the initial inputs are stimulated and we have one parameter per nand gate to represent the probability that it fails. pIMCs in text format are automatically generated from the PRISM model.

Table 1 summarizes the size of the considered instances of the model (in terms of states, transitions, and parameters) and of the corresponding CSP problems (in terms of number of variables and constraints). In addition, we also present the resolution time of the given CSPs using the Z3 solver. Our experiments were performed on a 2.4 GHz Intel Core i5 processor with time out set to 10 minutes and memory out set to 2Gb.

Conclusion and future work

In this paper, we have compared several Markov Chain abstractions in terms of succinctness and we have shown that Parametric Interval Markov Chain is a strictly more succinct abstraction formalism than other existing formalisms such as Parametric Markov Chains and Interval Markov Chains. In addition, we have proposed constraint encodings for checking several properties over pIMC. In the context of qualitative properties such as existencial consistency or consistent reachability, the size of our encodings is significantly smaller than other existing solutions. In the quantitative setting, we have compared the two usual semantics for IMCs and pIMCs and showed that both semantics are equivalent with respect to quantitative reachability properties. As a side effect, this result ensures that all existing tools and algorithms solving reachability problems in IMCs under the once-and-for-all semantics can safely be extended to the at-every-step semantics with no changes. Based on this result, we have then proposed CSP encodings addressing quantitative reachability in the context of pIMCs regardless of the chosen semantics. Finally, we have developed a prototype tool that automatically generates our CSP encodings and that can be plugged to any constraint solver accepting the SMT-LIB format as input.

We plan to develop our tool for pIMC verification in order to manage other, more complex, properties (e.g., supporting the LTL-language in the spirit of what

Fig. 3 :

 3 Fig. 1: MC M1

Definition 1 (

 1 Markov chain Abstraction Model). A Markov chain abstraction model (an abstraction model for short) is a pair (L, |=) where L is a nonempty set and |= is a relation between MC and L. Let P be in L and M be in MC we say that M implements P iff (M, P) belongs to |= (i.e., M |= P). When the context is clear, we do not mention the satisfaction relation |= and only use L to refer to the abstraction model (L, |=).

Fig. 4 :

 4 Fig. 4: MC M2 satisfying the IMC I from Figure 3 with a different structure

Fig. 5 :

 5 Fig. 5: pIMC P

Proposition 1 .

 1 The Markov chain abstraction models can be ordered as follows w.r.t. succinctness: (pIMC, |= o pI ) < (pMC, |= p ), (pIMC, |= o pI ) < (IMC, |= o I ) and (pIMC, |= a pI ) < (IMC, |= a I ).

Fig. 6 : 5 ⊥Fig. 7 :

 657 Fig.6: Variables in the CSP produced by C ∃c for the pIMC P from Fig.5

  They are written |= a pI and |= o pI and defined as follows: M |= a pI P (resp. |= o pI ) iff there exists an IMC I instance of P s.t. M |= a I I (resp. |= o I

Table 1 :

 1 trans. #par. #var. #cstr. time #var. #cstr. time #var. #cstr. time Benchmarks

		pIMC			C∃c	C∃r	C∃r
	Benchmark #states #nand K=1; N=2 104	147	4	255	1,526 0.17s 170 1,497 0.19s 296	2,457 69.57s
	nand K=1; N=3	252	364	5	621	3,727 0.24s 406 3,557 0.30s 703	5,828 31.69s
	nand K=1; N=5	930	1,371	7	2,308 13,859 0.57s 1,378 12,305 0.51s 2,404 20,165 T.O.
	nand K=1; N=10 7,392 11,207	12 18,611 111,366 9.46s 9,978 89,705 13.44s 17,454 147,015 T.O.

Indeed, when 0 ≤ v(f1) ≤ v(f2) ≤ 1 is not respected, the interval is inconsistent and therefore empty.

|L1| and |L2| are the sizes of L1 and L2, respectively.

All resources, benchmarks, and source code are available online as a Python library at https://github.com/anicet-bart/pimc_pylib

Available online at http://www.prismmodelchecker.com

This work is partially supported by the ANR national research programs PACS (ANR-14-CE28-0002) and Coverif (ANR-15-CE25-0002), and the regional programme Atlanstic2020 funded by the French Region Pays de la Loire and the European Regional Development Fund.

Tulip [START_REF] Wongpiromsarn | Tulip: A software toolbox for receding horizon temporal logic planning[END_REF] does). We also plan on investigating a practical way of computing and representing the set of all solutions to the parameter synthesis problem.