of the same job will have the same duration; and (ii) randomly and uniformly sampled from a given probability distribution law D, whose PDF is f and cumulative distribution function (CDF) is F . The probability distribution is assumed to be nonnegative, since we model execution times, and is defined either on a finite support [a, b], where 0 ≤ a < b, or on an infinite support [a, ∞) where 0 ≤ a. Hence, the execution time of a job is a random variable X, and P(X ≤ T ) = F (T ) = T a f (t)dt. For notational convenience, we sometimes extend the domain of f outside the support of D by letting f (t) = 0 for t ∈ [0, a].

B. Cost model

), where α > 0, β ≥ 0 and γ ≥ 0. If t > t 1 , another reservation should be paid for. Hence, the user needs to make a (possibly infinite) sequence of reservations S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ), where: 1) t i < t i+1 for all i ≥ 1. Indeed, because jobs are deterministic, it is redundant to have a duration in the sequence that is not strictly larger than the previous one, hence that duration can be removed from the sequence; 2) all possible execution times of the job are indeed smaller than or equal to some t i in the sequence. This simply means that the sequence must tend to infinity if job execution times are not upper-bounded.

Throughout the paper, we assume that both properties hold when speaking of a reservation sequence. For notational convenience, we define t 0 = 0, in order to simplify summations. Now, for a sequence S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ), and for a job execution time t, the cost is

C(k, t) = k-1 i=1 (αt i + βt i + γ) + αt k + βt + γ ( 2
)
where k is the smallest index in the sequence such that t ≤ t k (or equivalently, t k-1 < t ≤ t k ; recall that t 0 = 0).

C. Objective

The goal is to find a scheduling strategy, i.e., a sequence of increasing reservation durations, that minimizes the cost in expectation. Formally, the expected cost for a sequence S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ) can be written as:

E(S) = ∞ k=1 t k t k-1 C(k, t)f (t)dt (3)
Indeed, when t k-1 < t ≤ t k , the cost is C(k, t), weighted with the corresponding probability. Here are two examples:

• UNIFORM(a, b): for a uniform distribution over the interval [a, b] where 0 < a < b, we have f (t) = 1 b-a if a ≤ t ≤ b, and f (t) = 0 otherwise. Given a finite sequence S = ( a+b 2 , b), the expected cost is

E(S) = a+b 2 a (α a+b 2 + βt + γ) 1 b-a dt + b a+b 2 (α a+b 2 + β a+b 2 + γ) + (αb + βt + γ) 1 b-a dt
The first term is for values of t that are in [a, a+b 2 ] and the second term is for larger values of t in [ a+b 2 , b]. For the latter term, we pay a constant cost α a+b 2 + β a+b 2 + γ for the first unsuccessful reservation, and then a cost that depends upon the value of t for the second reservation if β = 0.

• EXP(λ): for an exponential distribution with rate λ and support in [0, ∞), we have f (t) = λe -λt for all t ≥ 0. Given an infinite and unbounded sequence S = ( 1 λ , 2 λ , . . . , i λ , i+1 λ , . . . ), the expected cost is

E(S) = ∞ k=1 k λ k-1 λ k-1 i=1 (α i λ +β i λ +γ)+α k λ +βt+γ λe -λt dt
Again, when t ∈ [ k-1 λ , k λ ], we pay a fixed cost for the k -1 first reservations, and a possibly variable cost for the kth reservation. Looking at the expression of E(S) above, we easily see that the given sequence S has a finite expected cost E(S). In fact, there are many sequences with finite expected cost, such as those defined by t i = ui + v for i ≥ 1, where u and v are positive constants.

We are now ready to state the optimization problem:

Definition 1 (STOCHASTIC). Given a probability distribution (with CDF F ) for the execution times of stochastic jobs, and given a cost function given by Equation (1) (with parameters α, β and γ), find a reservation sequence S with minimal expected cost E(S) as given in Equation [START_REF]AWS pricing information[END_REF].

We further define RESERVATIONONLY to be the instance of STOCHASTIC where the cost is a linear function of the reservation length only, i.e., when β = γ = 0. For RESER-VATIONONLY, we can further consider α = 1 without loss of generality. For instance, such costs are incurred when making reservations of resources to schedule jobs on some cloud platforms, with hourly or daily rates. Throughout the paper, we focus on the usual probability distributions, hence we assume that the density function f and the CDF F of D are smooth (infinitely differentiable), and that D has finite expectation.

III. CHARACTERIZING THE OPTIMAL SOLUTION

In this section, we establish key properties of an optimal solution in the general setting.

A. Cost function

We start by establishing a simpler expression for the cost function of STOCHASTIC.

Theorem 1. Given a sequence S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ), the cost function given by Equation (3) (with parameters α, β and γ) can be rewritten as (with t 0 = 0):

E(S) = β • E[X]+ ∞ i=0 (αt i+1 +βt i +γ) P(X ≥ t i ) (4)
Proof. We first expand Equation (3) as follows:

E(S) = ∞ k=1 t k t k-1 k i=1 (αt i +γ)+ k-1 i=1 βt i +βt f (t)dt (5)
We compute the three terms on the right-hand side separately. By defining t 0 = 0, the first term can be expressed as:

∞ k=1 t k t k-1 k i=1 (αti + γ) f (t)dt = ∞ k=1 k i=1 (αti + γ) t k t k-1 f (t)dt = ∞ k=1 k i=1 (αti + γ) P(X ∈ [t k-1 , t k ]) = ∞ i=1 ∞ k=i (αti + γ) P(X ∈ [t k-1 , t k ]) = ∞ i=1 (αti + γ) P(X ≥ ti-1)
Similarly, we obtain the second term:

∞ k=1 t k t k-1 k-1 i=1 βti f (t)dt = ∞ i=1 βtiP(X ≥ ti)
and the third term:

∞ k=1 t k t k-1 βtf (t)dt = β • E[X]
Plugging these three terms back into Equation ( 5), we get the desired expression for the cost function as given by Equation (4).

B. Upper bound on t o

1 and finite expected cost In this section, we extract an upper bound for the first request t o 1 of an optimal sequence S o to STOCHASTIC, which allows us to show that the expected cost E(S o ) is upper bounded too, and hence finite. This result holds in a general setting, namely, for any distribution D such that E(X 2 ) < ∞.

Obviously, if the distribution's support is upper bounded, such as for UNIFORM(a, b), a solution is to choose that upper bound for t o 1 (e.g., t o 1 ≤ b for UNIFORM(a, b)). Hence, we focus on distributions with infinite support [a, ∞) and aim at restricting the search for an optimal t o 1 to a bounded interval [a, A 1 ] for some A 1 . We derive the following result.

Theorem 2. For any distribution

D with infinite support [a, ∞) such that E[X 2 ] < ∞, the value t o 1 of an optimal sequence S o = (t o 1 , t o 2 , . . . , t o i , t o i+1 , . . . ) satisfies t o 1 ≤ A 1 , and E(S o ) ≤ A 2 , where A 1 = E[X]+1+ α+β 2α (E[X 2 ]-a 2 )+ α+β +γ α (E[X]-a) (6)
A 2 = β • E(X) + αA 1 + γ (7)
Proof. We consider the sequence S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ) with t i = a + i for i ≥ 1 (and t 0 = 0), and compute

E(S) -β • E[X] = ∞ i=0 (αti+1 + βti + γ)P(X ≥ ti) = ∞ i=0 (α(a + i + 1) + β(a + i) + γ)P(X ≥ a + i) = α(a + 1) + γ + ∞ i=1 (α + β)(a + i)P(X ≥ a + i) + (α + γ) ∞ i=1 P(X ≥ a + i) = α(a + 1) + γ + (α + β) ∞ i=1 a+i a+i-1 (a + i)P(X ≥ a + i)dt + (α + γ) ∞ i=1 a+i a+i-1 P(X ≥ a + i)dt
Note that for all t ∈ [a+i-1, a+i], we have both a+i ≤ t+1 and P(X ≥ a + i) ≤ P(X ≥ t), thus (a + i)P(X ≥ a + i) ≤ (t + 1)P(X ≥ t) Hence, we can write:

E(S) -β • E[X] ≤ α(a + 1) + γ + (α + β) ∞ i=1 a+i a+i-1 (t + 1)P(X ≥ t)dt + (α + γ) ∞ i=1 a+i a+i-1 P(X ≥ t)dt = α(a + 1) + γ + (α + β) ∞ a (t + 1)P(X ≥ t)dt + (α + γ) ∞ a P(X ≥ t)dt ≤ α(a + 1) + γ + (α + β) ∞ a t • P(X ≥ t)dt + (2α + β + γ) ∞ a P(X ≥ t)dt
For the last inequality, we have split

∞ a (t + 1)P(X ≥ t)dt into ∞ a tP(X ≥ t)dt and ∞ a P(X ≥ t)dt.
Extending the support of D to [0, ∞) by letting f (t) = 0 for 0 ≤ t ≤ a, and hence P(X ≥ t) = 1 for 0 ≤ t ≤ a, we have the following property for any integer p ≥ 1:

∞ 0 t p-1 • P(X ≥ t)dt = ∞ t=0 t p-1 ∞ x=t f (x)dxdt = ∞ x=0 f (x) x t=0 t p-1 dtdx = ∞ 0 x p p f (x)dx = E[X p ]
p Hence, using p = 1, we have:

∞ a P(X ≥ t)dt = ∞ 0 P(X ≥ t)dt- a 0 P(X ≥ t)dt = E[X]-a
and using p = 2, we get

∞ a t • P(X ≥ t)dt = ∞ 0 t • P(X ≥ t)dt - a 0 t • P(X ≥ t)dt = E[X 2 ]-a 2 2
Altogether, we derive that

E(S) ≤ β • E[X] + αA 1 + γ (8)
where A 1 is given by Equation [START_REF] Canon | Evaluation and optimization of the robustness of dag schedules in heterogeneous environments[END_REF]. From Equation (4), the expected cost of any sequence S satisfies E(S) ≥ β • E[X] + αt 1 + γ (cost of expected execution time and cost of first request). Hence, necessarily in an optimal sequence, the first reservation t o 1 satisfies to t o 1 ≤ A 1 . Thus, Equation ( 6) gives the desired bound on t o 1 .

C. Properties of optimal sequences

We now derive a recurrence relation between the successive requests in the optimal sequence for STOCHASTIC. Theorem 3. Let S o = (t o i ) i≥1 denote an optimal sequence for STOCHASTIC. For all i ≥ 1, if F (t o i ) = 1, we have the following property:

αt o i+1 + βt o i + γ = α 1 -F (t o i-1 ) f (t o i ) + β 1 -F (t o i ) f (t o i ) (9)
Proof. We fix an index j ≥ 1 such that F (t o j ) = 1 and consider the expected cost when we replace t o j by an arbitrary value t ∈ [t o j-1 , t o j+1 ]. This amounts to using the sequence

S o j (t) = (t o 1 , t o 2 , • • • , t o j-1 , t, t o j+1 , • • • )
whose expected cost, according to Equation (4), is the following:

E(S o j (t)) = β • E[X] + i =j-1,j (αt o i+1 + βt o i + γ)P(X ≥ t o i ) + (αt + βt o j-1 + γ)P(X ≥ t o j-1 ) + (αt o j+1 + βt + γ)P(X ≥ t)
which we can rewrite as:

E(S o j (t)) = C j +αt(1-F (t o j-1 ))+(αt o j+1 +βt+γ)(1-F (t))
where C j is some constant independent of t. By definition, the minimum of

E(S o j (t)) on [t o j-1 , t o j+1 ] is achieved at t = t o j
(and potentially at other values). Because E(S o j (t)) is smooth, we have that its derivative at t o j , which is not an extremity of the interval [t o j-1 , t o i+1 ], must be equal to zero, i.e., ∂E(S o j (t)) ∂t = 0. This gives:

α(1 -F (t o j-1 ))+β(1 -F (t o j )) -(αt o j+1 +βt o j +γ)f (t o j ) = 0 (10)
To get the final result, it remains to show that f (t o j ) = 0. Otherwise, we would get from Equation (10) that α(1 -

F (t o j-1 ))+β(1-F (t o j )) = 0, which implies that F (t o j-1 ) = 1 because α > 0 (and β(1 -F (t o j )) ≥ 0). But then, F (t o j ) ≥ F (t o j-1 ) = 1
, which contradicts the initial assumption. Hence, f (t o j ) = 0, and rewriting Equation ( 10) directly leads to Equation [START_REF] Dong | Resource load based stochastic DAGs scheduling mechanism for Grid environment[END_REF].

Note that the condition F (t o i ) = 1 in Theorem 3 applies to distributions with finite support, such as UNIFORM(a, b), where F (b) = 1. For the usual distributions with infinite support, such as EXP(λ), we have F (t) < 1 for all t ∈ [0, ∞) and an optimal sequence must be infinite. In essence, Theorem 3 suggests that an optimal sequence is characterized solely by its first value t o 1 : Proposition 1. For a smooth distribution with infinite support, solving STOCHASTIC reduces to finding

t o 1 that minimizes ∞ i=0 (αt i+1 +βt i +γ) P(X ≥ t i )
where t o 0 = 0, and for all i ≥ 2,

t o i = 1 -F (t o i-2 ) f (t o i-1 ) + β α 1 -F (t o i-1 ) f (t o i-1 ) -t o i-1 - γ α ( 11
)
For a smooth distribution with finite support, the recurrence in Equation [START_REF] Friedman | Report on a multicenter fMRI quality assurance protocol[END_REF] still holds but the optimal sequence stops as soon as it reaches t o i with F (t o i ) = 1. Proposition 1 provides an optimal algorithm for general smooth distributions, up to the determination of t o 1 . However, computing the optimal t o 1 , remains a difficult problem, except for simple distributions such as UNIFORM(a, b) (see Section III-D).

D. Uniform distributions

In this section, we discuss the optimal strategy for a uniform distribution UNIFORM(a, b), where 0 < a < b. Intuitively, one could try and make a first reservation of duration, say, t 1 = a+b 2 , and then a second reservation of duration t 2 = b. However, we show that the best approach is to make a single reservation of duration t 1 = b, for any value of the parameters α, β and γ: Proof. We study here only the sequences of length smaller than or equal to 2, i.e., S = (t 1 , b), where t 1 ≤ b. The full proof is available in the companion report [START_REF] Aupy | Reservation Strategies for Stochastic Jobs[END_REF].

Note that necessarily, in a sequence of length 2,

t 2 = b otherwise t 2 < b and E((t 1 , t 2 )) = ∞ because the interval [t 2 , b] has non-zero measure.
Using Equation ( 4) we obtain:

E((t 1 , b)) = βE(X) + (αt 1 + γ)P(X ≥ 0) + (αb + βt 1 + γ)P(X ≥ t 1 ) = βE(X) + (αt 1 + γ) + b-t1 b-a (αb + βt 1 + γ) = βE(X) + α b 2 -t1a b-a + b-t1 b-a (βt 1 + γ) + γ
We can verify easily that this is minimized when t 1 = b (and that E((t 1 , b)) = E((t 1 )). Hence S o = (b) is optimal amongst the sequences of length smaller than or equal to 2.

E. Exponential distributions

In this section, we provide partial results for the RESER-VATIONONLY problem (β = γ = 0 and α = 1) with an exponential distribution EXP(λ). From Theorem 2 (and the example in Section II-C), we know that there exist sequences of finite expected cost. We further characterize the optimal solution as follows:

Proposition 2. Let S 1 = (s 1 , s 2 , . . . , s i , s i+1 , . . . ) denote the optimal sequence for RESERVATIONONLY with an EXP(1) distribution. Then, s 2 = e s1 , and for i ≥ 3,

s i = e si-1-si-2 = e si-1 Π i-1 j=2 s j (12)
The expected cost of

S 1 is E(S 1 ) = s 1 + 1 + ∞ i=1 e -si
. Furthermore, the optimal sequence for RESERVATIONONLY with EXP(λ) distribution is the infinite sequence S λ = (t 1 , t 2 , . . . , t i , t i+1 , . . . ) such that t i = si λ for i ≥ 1. Its expected cost is E(S λ ) = 1 λ E(S 1 ). Proof. Consider an EXP(λ) distribution. From Equation (4), the expected cost of the optimal sequence S λ is E(S) = ∞ i=0 t i+1 e -λti where t 0 = 0, t 1 is unknown, and the value of t i for i ≥ 2 is given by Equation [START_REF] Friedman | Report on a multicenter fMRI quality assurance protocol[END_REF] as t i = e λ(t i-1 -t i-2 ) λ for i ≥ 2. Introducing s i = λt i for all i ≥ 0, we derive that

E(S λ ) = 1 λ ∞ i=0 s i+1 e -si
with s i = e si-1-si-2 for all i ≥ 2. We have s i+1 e -si = e -si-1 for i ≥ 1, which gives the desired value for E(S λ ). Now, we prove Equation ( 12) by induction. It holds for i = 3, because s 3 = e s2-s1 = e s 2 e s 1 = e s 2 s2 . Assume that it holds for any j ≤ i. Then s i+1 = e si-si-1 = e s i e s i-1 , and by induction e si-1 = s i s 2 . . . s i-1 , hence the result.

Again, the optimal sequence is fully characterized by the value of t 1 or s 1 . Here, s 1 is independent of λ. In other words, the solution for EXP(1) is generic, and the solution for EXP(λ) for an arbitrary λ can be directly derived from it. Unfortunately, we do not know how to compute s 1 analytically. However, a brute-force search provides the value s 1 ≈ 0.742, which means that the first reservation for EXP(λ) should be approximately three quarters of the mean value 1 λ of the distribution, for any λ > 0.

IV. HEURISTICS FOR ARBITRARY DISTRIBUTIONS

The results of the preceding section provide a strategy to compute the optimal sequence up to the determination of t o 1 , since Theorem 3 and Proposition 1 allow us to compute the subsequent t o i 's. However, while we have derived an upper bound on t o 1 , we do not know how to compute its exact value for an arbitrary distribution. In this section, we introduce several heuristics for the STOCHASTIC problem under arbitrary probability distributions.

A. Brute-force procedure

We first present a procedure called BRUTE-FORCE that simply tries different values for the first reservation length t 1 in a sequence S, and then computes the subsequent values according to Equation [START_REF] Friedman | Report on a multicenter fMRI quality assurance protocol[END_REF]. Specifically, we try M different values of t 1 on the interval [a, b], where a is the lower bound of the distribution and b is the upper bound if the distribution is finite. Otherwise, we set b = A 1 , which is an upper bound on the optimal t o 1 as given in Equation [START_REF] Canon | Evaluation and optimization of the robustness of dag schedules in heterogeneous environments[END_REF]. As an example, when α = 1, β = γ = 0, we have

A 1 ≤ 1 + 2E(X) + E(X 2 )
2 . For each m = 1, . . . , M , we generate a sequence that starts with

t 1 = a + m • b-a M .
Given a sequence S, its expected cost is evaluated via a Monte-Carlo process, as described in Section V-A: we randomly draw N execution times from the distribution, and compute the expected cost incurred by the sequence over the N samples. We finally return the minimum expected cost found over all the M values of t 1 . Note that some values of t 1 may not lead to any result, because the sequence computed based on it and using Equation ( 11) may not be strictly increasing. In this case, we simply ignore the sequence. The complexity of this heuristic is O(M N ).

We point out that the actual optimal value for the first request t o 1 would possibly lie in between two successive values of t 1 that we try. However, because we deal with smooth probability distributions, we expect to return a t 1 and an associated expected cost that are close to the optimal when M and N are sufficiently large. In the performance evaluation, we set M = 5000 and N = 1000.

B. Discretization-based dynamic programming

We now present a heuristic that approximates the optimal solution for STOCHASTIC by first discretizing the continuous distribution and then computing an optimal sequence for the discrete problem via dynamic programming.

1) Truncating and discretizing continuous distributions: If a continuous distribution has finite support [a, b], where 0 ≤ a < b, then we can directly discretize it. Otherwise, for a distribution with infinite support [a, ∞), where 0 ≤ a, we need to first truncate it in order to operate on a bounded interval. In the latter case, we define b = Q(1 -ǫ), where Q(x) = inf{t|F (t) ≥ x} is the quantile function. That is, we discard the final ǫ ∈ (0, 1) quantile of the distribution, which for usual distributions ensures that b is finite. In either case, the discretization will then be performed on the interval [a, b]. Let n denote the number of discrete values we will sample from the continuous distribution. The result will be a set of n pairs (v i , f i ) i=1...n , where the v i 's represent the possible execution times of the jobs, and the f i 's represent the corresponding probabilities. We envision two schemes for the discretization:

• EQUAL-PROBABILITY: This scheme ensures that all the discrete execution times have the same probability. Thus, for all i = 1, 2, . . . , n, we can compute

v i = Q i • F (b) n and f i = F (b)
n .

• EQUAL-TIME: This scheme makes the discrete execution times equally spaced in the interval [a, b]. Thus, for all i = 1, 2, . . . , n, the execution times and their probabilities are computed as

v i = a + i • b-a n and f i = F (v i ) -F (v i-1
). Note that when the continuous distribution has infinite support, the probabilities for the n discrete execution times do not sum up to 1, i.e.,

n i=1 f i = F (b) = 1 -ǫ.
A smaller value of ǫ and a larger number n will provide a better sampling of the continuous distribution in either discretization scheme. In the performance evaluation, we set ǫ = 10 -7 and n = 1000.

2) Dynamic programming for discrete distributions: We now present a dynamic programming algorithm to compute the optimal sequence for any discrete probability distribution. It will be used with the discretization schemes to approximate the optimal solution for an arbitrary continuous distribution.

Theorem 5 (Discrete distribution). If X ∼ (v i , f i ) i=1...n , then STOCHASTIC can be solved optimally in polynomial time.

Proof. Let E *

i denote the optimal expected cost given that X ≥ v i . In this case, to compute the optimal expected cost, the probability distribution of X needs to be first updated as

f ′ k = f k n j=i fj , ∀k = i, . . . , n, which guarantees that n k=i f ′ k = 1.

We can then express E *

i based on the following dynamic programming formulation:

E * i = min i≤j≤n αv j +γ + j k=i f ′ k • βv k + n k=j+1 f ′ k βv j +E * j+1
In particular, to compute E * i , we make a first reservation of all possible discrete values (v j ) j=i...n and select the one that incurs the minimum total expected cost. For each v j considered, if the job's actual execution time is greater than v j (with probability n k=j+1 f ′ k ), the total cost also includes the optimal cost E * j+1 for making subsequent reservations. The dynamic program is initialized with E * n = αv n + βv n + γ, and the optimal total expected cost is given by E *

1 . The complexity is O(n 2 ), since each E * i depends on n -i other expected costs, with associated probability updates and summations that can be computed in O(n -i) time. The optimal sequence of reservations can be obtained by backtracking the decisions made at each step.

Note that the sequence obtained by dynamic programming always ends with the largest value v n = b. When applying it back to a continuous distribution with infinite support, more values will be needed, because the sequence must tend to infinity as explained in Section II-B. In this case, additional values can be appended to the sequence by using other heuristics, such as the ones presented next in Section IV-C.

C. Other heuristics

We finally present some simple heuristics that are inspired by common resource allocation strategies in the literature. These heuristics do not explore the structure of the optimal solution nor the probability distribution, but rely on simple incremental methods to generate reservation sequences.

In the following, we will use µ = E(X) to denote the mean of a given distribution, σ2 = E(X 2 ) -µ 2 to denote its variance, and m = Q( 12 ) to denote its median, where Q(x) = inf{t|F (t) ≥ x} represents the quantile function. The different heuristics are defined as follows:

• MEAN-BY-MEAN: start with the mean (i.e., t 1 = µ) and then make each subsequent reservation request by computing the conditional expectation of the distribution in the remaining interval, i.e.,

t i = E(X|X > t i-1 ) = ∞ t i-1 tf (t)dt 1-F (ti-1) for all i ≥ 2.
• MEAN-STDEV: start with the mean and then increment the reservation length by one standard deviation (σ) for each subsequent request, i.e., t i = µ + (i -1)σ for all i ≥ 2.

• MEAN-DOUBLING: start with the mean and then double the reservation length for each subsequent request, i.e., t i = 2 i-1 µ for all i ≥ 2.

• MEDIAN-BY-MEDIAN: each request is the median of the distribution in the remaining interval, i.e.,

t i = Q(1 -1 2 i ) for all i ≥ 2.
Note that deriving the sequence for MEAN-BY-MEAN is straightforward for some distributions (e.g., exponential, uniform), but more involved for others. Recursive formulas are provided in [START_REF] Aupy | Reservation Strategies for Stochastic Jobs[END_REF] to compute the sequence using this heuristic for the considered distributions, along with key parameters (e.g., mean, variance, quantile) to facilitate the sequence computations for the other heuristics.

V. PERFORMANCE EVALUATION

In this section, we evaluate the different heuristics presented in Section IV, and compare their performance. The code and setup of the experiments presented in this section are publicly available on https://gitlab.inria.fr/vhonore/ipdps 2019 stochastic-scheduling.

A. Evaluation methodology

For each heuristic that generates a reservation sequence S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ) under a particular probability distribution D, we approximate its expected cost via a Monte-Carlo process 2 : we randomly sample N possible execution times from the distribution, and then average over the cost of all the N samples, i.e.,

Ẽ(S) = 1 N N i=1 C(k, t)| t←D ( 13
)
where C(k, t) is the cost for a specific execution time t drawn from the distribution, computed using Equation [START_REF] Afanasyev | Service provider competition: Delay cost structure, segmentation, and cost advantage[END_REF]. For the presented evaluation results, we set N = 1000.

-λt λ = 1.0 t ∈ [0, ∞) Weibull(λ, κ) κ λ t λ κ-1 e -( t λ ) κ λ = 1.0 κ = 0.5 t ∈ [0, ∞) Gamma(α, β) β α Γ(α) t α-1 e -βt α = 2.0 β = 2.0 t ∈ [0, ∞) LogNormal ν, κ 2 1 tκ √ 2π e - (ln t-ν) 2 2κ 2 ν = 3.0 κ = 0.5 t ∈ (0, ∞) TruncatedNormal(ν, κ 2 , a) 1 κ 2 π • e -1 2 t-ν κ 2 1-erf a-ν κ √ 2 ν = 8.0 κ 2 = 2.0 a = 0.0 t ∈ [a, ∞)
Pareto(ν, α)

αν α t α+1 ν = 1.5 α = 3.0 t ∈ [ν, ∞)
Distributions with finite support Uniform(a, b)

1 b-a a = 10.0 b = 20.0 t ∈ [a, b] Beta(α, β) t α-1 •(1-t) β-1 B(α,β) α = 2.0 β = 2.0 t ∈ [0, 1]
BoundedPareto(L, H, α)

αL α t -α-1 1-L H α L = 1.0 H = 20.0 α = 2.1 t ∈ [L, H]
To get uniform results, we normalize the expected cost of each heuristic by the expected cost of an omniscient scheduler, which knows the job execution time t a priori, and thus would make a single request of length t 1 = t. Averaging over all possible values of t from the distribution D, the omniscient scheduler has an expected cost:

E o = ∞ 0 (αt + βt + γ)f (t)dt = (α + β)E[X] + γ
Hence, the normalized ratio will always be greater than or equal to 1, and a smaller ratio means a better result.

We perform the evaluation of the heuristics under two different reservation-based scenarios.

• RESERVATIONONLY (Section V-B): This scenario is based on the Reserved Instance pricing scheme available in AWS [START_REF]AWS pricing information[END_REF], where the user pays exactly what is requested. Hence, we set α = 1, β = γ = 0. We consider nine probability distributions in this case, six of which have infinite support and the remaining three have finite support. Table I lists these distributions with instantiations of their parameters used in the evaluation. • NEUROHPC (Section V-C): This scenario is based on executing large jobs on HPC platforms, where the cost, as represented by the total turnaround time of a job, is the sum of its waiting time in the queue and its actual execution time. We set β = 1 for the execution time and instantiate the waiting time function (α, γ) by curvefitting the data from Fig. 2b. The probability distribution is derived from the execution traces of the neuroscience application shown in Fig. 1b.

B. Results for RESERVATIONONLY scenario

Table II presents, for each heuristic, the normalized expected cost, i.e., Ẽ(S)/E o , under different probability distributions. The BRUTE-FORCE heuristic tries M = 5000 values of t 1 , and both discretization heuristics set the truncation parameter to be ǫ = 10 -7 and use n = 1000 samples. First, the normalized costs allow us to compare the performance of these heuristics generally with long-running jobs that require a large amount of resources having higher priorities. Jobs that are kept in the waiting queue for a long period of time could also be upgraded and moved up in the queue. Slurm schedules the jobs from the top of the high-priority queue and moves down. Even though larger jobs (in term of time and space) have higher priorities, generally the lack of resource availability in the system leads to longer wait times. On the other hand, smaller jobs, despite having lower priorities, are usually scheduled quickly thanks to the backfilling algorithms that place them in the unused time slots between successive large jobs.

Some studies (e.g., [START_REF] Mu | Utilization, predictability, workloads, and user runtime estimates in scheduling the IBM SP2 with backfilling[END_REF], [START_REF] Rodrigo Álvarez | HPC system lifetime story: Workload characterization and evolutionary analyses on NERSC systems[END_REF], [START_REF] Yamamoto | The K computer operations: Experiences and statistics[END_REF]) have analyzed the impact of scheduling strategies on the performance of applications in large HPC centers. Some of these studies show that the penalty for jobs with longer requested walltimes and/or larger numbers of nodes is higher than that for jobs with shorter elapsed times and smaller numbers of nodes. This is observed, for example, in [START_REF] Yamamoto | The K computer operations: Experiences and statistics[END_REF] for the K computer from Riken Advanced Institute for Computational Science. The study shows that, for applications requesting similar computing resources, the wait time generally increases with larger requested processing times and can cause delays of hours for large scientific applications, although it is also dependent on other workloads submitted to the system. Some HPC centers divide the resources into seasons for users to utilize the reserved resources. Users tend to submit more jobs toward the end of a season causing contention at the scheduler level which results in even longer waiting times. The study in [START_REF] Rodrigo Álvarez | HPC system lifetime story: Workload characterization and evolutionary analyses on NERSC systems[END_REF] presents a trend of the evolution of the workload of HPC systems and the corresponding scheduling policies as we move from monolithic MPI applications to high-throughput and data-intensive jobs. The paper shows that the cost paid in terms of the wait time of applications in the queue has generally increased over the years with less uniform workloads. The study in [START_REF] Mu | Utilization, predictability, workloads, and user runtime estimates in scheduling the IBM SP2 with backfilling[END_REF] shows that systems that give each job a partition of the resources for exclusive use and allocate such partitions in the order of job arrivals could suffer from severe fragmentation, leading to low utilization. The authors propose an aggressive backfilling algorithm for dealing with such fragmentation. However, users are still expected to provide accurate runtime estimates. The study shows that over-estimation may lead to a long wait time and possibly excessive CPU quota loss, while under-estimations can lead to job terminations before completion. Some recent schedulers [START_REF] Park | 3sigma: distribution-based cluster scheduling for runtime uncertainty[END_REF] consider the distribution of execution time of the submitted jobs to take their scheduling decision in order to increase their overall utility.

b) Stochastic job scheduling: Many works deal with stochastic job scheduling (e.g., [START_REF] Canon | Evaluation and optimization of the robustness of dag schedules in heterogeneous environments[END_REF], [START_REF] Dong | Resource load based stochastic DAGs scheduling mechanism for Grid environment[END_REF], [START_REF] Tang | A stochastic scheduling algorithm for precedence constrained tasks on grid[END_REF]- [START_REF] Weiss | Turnpike optimality of smith's rule in parallel machines stochastic scheduling[END_REF]). Various models [START_REF] Canon | Scheduling stochastic tasks under deadline and budget constraints[END_REF] have been proposed to model the performance of executing stochastic jobs on computing platforms. For instance, in [START_REF] Li | Scheduling precedence constrained stochastic tasks on heterogeneous cluster systems[END_REF], stochastic jobs are modeled as a DAG of tasks whose execution times and communication times are stochastically independent. In this paper, we model jobs by an execution time following a probability distribution. The authors in [START_REF] Tang | A stochastic scheduling algorithm for precedence constrained tasks on grid[END_REF] propose a model based on resource load in grid systems. Several refinements can be envisioned, such that improving scheduler performances by including distribution features in order to optimize final performance. Also, dealing with heterogeneous nodes increases problem complexity [START_REF] Topcuoglu | Performance-effective and lowcomplexity task scheduling for heterogeneous computing[END_REF]. We refer the reader to the book by Pinedo [START_REF] Pinedo | Scheduling: Theory, Algorithms, and Systems[END_REF] which contains a comprehensive survey of stochastic scheduling problems, and to the book chapter [START_REF] Gaujal | Comparisons of stochastic task-resource systems[END_REF] for a detailed comparison of stochastic task-resource systems.

c) Pricing and reservation schemes for cloud computing: Cloud computing platforms have emerged as another option for executing HPC applications. Job scheduling in the cloud has an even bigger challenge [START_REF] Gupta | Evaluating and improving the performance and scheduling of HPC applications in cloud[END_REF], since it needs to deal with highly heterogeneous resources with a wide range of processor configurations, interconnects, virtualization environments, etc.

Different pricing and reservation schemes are also available for users who submit jobs to a cloud service. Several works have been conducted to study these schemes in the cloud, and from a computer science perspective, many of these studies focus on the pricing strategies and service management of platform providers [START_REF] Afanasyev | Service provider competition: Delay cost structure, segmentation, and cost advantage[END_REF], [START_REF] Chen | Pricing schemes in cloud computing: Utilization-based versus reservation-based[END_REF], [START_REF] Dierks | Cloud pricing: the spot market strikes back[END_REF], [START_REF] Xu | Dynamic cloud pricing for revenue maximization[END_REF]. Some works consider modeling the delays for users [START_REF] Afanasyev | Service provider competition: Delay cost structure, segmentation, and cost advantage[END_REF] and how providers manage the idle resources [START_REF] Dierks | Cloud pricing: the spot market strikes back[END_REF]. The work in [START_REF] Xu | Dynamic cloud pricing for revenue maximization[END_REF] studies the pricing practices of Amazon AWS [START_REF]AWS pricing information[END_REF] when the price is dynamically adapted to real-time demand and idle resources. In [START_REF] Chen | Pricing schemes in cloud computing: Utilization-based versus reservation-based[END_REF], authors provide an analytical model of pricing for reservation-based scheme (used by Amazon AWS) and utilization-based scheme (used by Google GCP [START_REF]GCP pricing information[END_REF]). They show that the effective price mainly depends on the variation of platform usage and the competition for customers. Some tools are also provided for users to perform cost evaluation in order to select which type of platform to use. They show that users with high-volatility demand should consider using AWS offers while one should use GCP in the other case. The experimental results in this paper suggest that, compared with on-demand or utilizationbased services, reservation strategies can provide cost-effective options for executing stochastic jobs when there is significant difference in the offered price.

VII. CONCLUSION

In this paper, we have studied the problem of scheduling stochastic jobs on a reservation-based platform. We have shown the existence of an optimal reservation sequence when the job execution time follows a set of classical distributions, and we have characterized the optimal solution up to the duration of the first reservation. We do not know how to compute this duration analytically, but we have provided an upper bound and a brute-force procedure to generate a solution that is close to the optimal. We have also introduced several heuristics, one based upon discretizing the continuous distribution, and some relying on standard measures, such as the mean, variance and quantiles of the distribution. We have demonstrated the effectiveness of these heuristics via comprehensive simulations conducted using both classical distributions and execution traces of a real neuroscience application.

Future work will include allowing requests with variable amount of resources, hence offering a combination of a reservation time and a number of processors. Another interesting direction is to include checkpoint snapshots at the end of some reservations. We expect the solutions such as the one introduced in this work not to work because of the difficulty of choosing which reservations to checkpoint. Indeed we do not expect the strategy "checkpoint all reservations" to be optimal. Hence the checkpointing approach calls for a complicated trade-off between doing useful work through the reservations and sacrificing some time/budget in order to avoid restarting the job whenever its execution time exceeds the length of the current reservation. The cost will then depend both on the length of the reservation and on a conditional probability based on previous checkpointing decisions.

Theorem 4 .

 4 For a uniform distribution UNIFORM(a, b), the optimal sequence for STOCHASTIC is S o = (b).

Table I .

 I Probability distributions and parameter instantiations

	Distribution	PDF f (t)	Instantiation	Support
		Distributions with infinite support
	Exponential (λ)	λe

Table II .

 II Normalized expected costs of different heuristics in the RESERVATIONONLY scenario under different distributions. The values in the brackets show the expected costs normalized by those of the BRUTE-FORCE heuristic (best online heuristic).

	Distribution	BRUTE-FORCE	MEAN-BY-MEAN	MEAN-STDEV	MEAN-DOUB.	MED-BY-MED	EQUAL-TIME	EQUAL-PROB.
	Exponential	2.15	2.36 (1.10)	2.39 (1.11)	2.42 (1.13)	2.83 (1.32)	2.31 (1.07)	2.36 (1.10)
	Weibull	2.12	2.76 (1.30)	3.58 (1.69)	3.03 (1.43)	3.05 (1.44)	2.40 (1.13)	2.22 (1.05)
	Gamma	2.02	2.26 (1.12)	2.18 (1.08)	2.24 (1.11)	2.51 (1.24)	2.20 (1.09)	2.13 (1.05)
	Lognormal	1.85	2.19 (1.19)	2.09 (1.13)	1.95 (1.06)	2.30 (1.24)	1.87 (1.01)	1.93 (1.04)
	TruncatedNormal	1.36	1.98 (1.46)	1.83 (1.35)	1.98 (1.46)	2.16 (1.60)	1.38 (1.02)	1.36 (1.00)
	Pareto	1.62	1.82 (1.12)	2.18 (1.34)	1.75 (1.08)	2.26 (1.39)	1.71 (1.05)	1.66 (1.03)
	Uniform	1.33	2.21 (1.66)	1.90 (1.43)	1.67 (1.26)	2.21 (1.66)	1.33 (1.00)	1.33 (1.00)
	Beta	1.75	2.02 (1.15)	2.11 (1.20)	1.98 (1.13)	2.45 (1.40)	1.79 (1.02)	1.80 (1.02)
	BoundedPareto	1.80	1.84 (1.02)	2.09 (1.16)	1.83 (1.01)	2.81 (1.56)	2.00 (1.11)	1.91 (1.06)

The possibly infinite sequence prevents us from analytically evaluating its expected cost.

Acknowledgments: We thank Bennett Landman and his MASI Lab at Vanderbilt for sharing the medical imaging database used to extract the execution time distributions. This research was supported in part by the National Science Foundation grant CCF1719674, Vanderbilt Institutional Fund, and Inria-Vanderbilt associated team Keystone. Part of this work was done while Valentin Honoré was visiting Vanderbilt University.


Miriam Pescador-Rojas

ESCOM Denis Pallez
email: denis.pallez@unice.fr

Carlos Ignacio
email: carlos@mutuofinanciera.comccoello@cs.cinvestav.mx

Hernández Castellanos

Carlos A Coello Coello

Studying the Effect of Robustness Measures in Offline Parameter Tuning for Estimating the Performance of MOEA/D

Keywords: Offline parameter tuning, robustness measures, multiobjective evolutionary algorithms, scalarizing functions

Offline parameter tuning (OPT) of multi-objective evolutionary algorithms (MOEAs) has the goal of finding an appropriate set of parameters for solving a large number of problems. According to the no free lunch theorem (NFL), no algorithm can have the best performance in all classes of optimization problems. However, it is possible to find an appropriate set of parameters of an algorithm for solving a particular class of problems. For that sake, we need to study how to estimate the aggregation quality function for an algorithmic configuration assessed on a set of optimization problems.

In this paper, we study robustness measures for dealing with the parameter settings of stochastic algorithms. We focus on decomposition-based MOEAs and we propose to tune scalarizing functions for solving some classes of problems based on the Pareto front shapes using up to 7 objective functions. Based on our experimental results, we were able to derive interesting guidelines to evaluate the quality of algorithmic configurations using a combination of descriptive statistics.

I. INTRODUCTION

The challenge for offline parameter tuning (OPT) techniques is to find the best parameter configuration of an evolutionary algorithm that successfully solves a large number of multiobjective optimization problems (MOPs). In multiobjective evolutionary algorithms (MOEAs), each parameter has a significant influence on their performance such as its rate of convergence and the quality of the solutions obtained. OPT is subject to uncertainty which originates from two main sources: 1) the stochastic procedures involved in MOEAs, 2) the optimal configuration for an MOEA optimizing a set of problem instances at the same time.

The goal of this paper is to study robustness measures used in the area of uncertainty [START_REF] Beyer | Robust optimization a comprehensive survey[END_REF], [START_REF] Gaspar-Cunha | Robustness in multiobjective optimization using evolutionary algorithms[END_REF] to compute the performance of an algorithmic configuration used in the offline tuning methods. Our main contribution is to provide a set of guidelines on the use of robustness measures for solving multi-objective optimization problems (MOPs). We adopt a particular case study to improve the performance of a MOEA based on decomposition (MOEA/D) and we employ an evolutionary OPT tool which has been adopted in other studies [START_REF] Montero | Towards a Method for Automatic Algorithm Configuration: A Design Evaluation Using Tuners[END_REF]- [START_REF] Pescador-Rojas | An Overview of Weighted and Unconstrained Scalarizing Functions[END_REF]. We focus on an important component of MOEA/D: the scalarizing function (SF), which is a strategy to transform a MOP into several single objective problems. Some SFs involve additional parameters in R that require to be tuned according to particular characteristics on the MOP such as the Pareto front shape or the number of objective functions.

The remainder of this paper is organized as follows. Previous works on parameter settings are briefly discussed in Section II with a particular focus on the MOPs detailed in Section IV. Then, uncertainty and robustness measures are introduced in Section III. Decomposition-based MOEAs based on scalarizing functions are described in Section V. Section VI presents our methodology for studying robustness. Section VII presents our results and their discussion. We finally conclude the paper and present some possible paths for future work in Section VIII.

II. PARAMETERS SETTINGS

The performance of EAs is strongly related to the definition of appropriate parameter values. Typically, before running an EA, the user empirically determines the proper encoding, selection mechanism, and evolutionary operators, as well as their numerical parameters such as the population size, mutation and crossover rates. There exist two main possibilities for dealing with the parameter setting problem: the use of offline parameter tuning and the use of online parameter control strategies [START_REF] Eiben | Parameter tuning for configuring and analyzing evolutionary algorithms[END_REF]. The first one refers to selecting a set of parameter values which can be established by hand according to the user's prior experience or by applying experimental design methods. In these cases, the same set of parameters is used in all the iterations of an EA. The second strategy includes adaptive mechanisms where the parameter values are modified using information gathered during the evolutionary search process. Adaptive and self-adaptive strategies are feedbackbased approaches. The former uses indicators to monitor the performance of the evolutionary process. The second strategy encodes the parameter values in each individual of the population and solves at the same time the optimization and the configuration problems.

In this work, we focus on the OPT problem which can be formalized as a 6-tuple T = (A, P, C, O, ψ, b max ) as follows:

• A is a target algorithm to be tuned. In this work, we concentrate on MOEA/D, which is the most popular decomposition-based MOEA; • P is the set of parameters of A to be tuned. For instance, if A is an EA, then population size, maximum number of generations, crossover rate, etc. belong to P . • C is the set of values for each parameter of P . An instance c ∈ C is called a configuration of A; • O is called a scenario and is a set of several problem instances (for example, {DTLZ1, DTLZ3, DTLZ5, DTLZ7} from the Deb-Thiele-Laumanns-Zitzler (DTLZ) test suite [START_REF] Deb | Scalable Test Problems for Evolutionary Multiobjective Optimization[END_REF]); • ψ is a scalar fitness measure used by T to assess the quality of a configuration e.g., hypervolume (HV) [START_REF] Zitzler | Multiobjective Optimization Using Evolutionary Algorithms-A Comparative Study[END_REF], inverted generational distance (IGD) [START_REF] Coello Coello | Solving Multiobjective Optimization Problems using an Artificial Immune System[END_REF] or the R2 [START_REF] Hansen | Evaluating the quality of approximations to the non-dominated set[END_REF] indicator. • b max is the maximum budget of function evaluations available to tune A on O;

Here, the goal is to find the 'best' configuration c ∈ C, based on ψ, for tuning A in an scenario O with b max function evaluations.

OPT is computationally expensive, but it is useful for deriving knowledge about the relationships among the parameters involved in an EA. This procedure can tackle the problem from two perspectives: specialization or generalization of algorithms. From the first perspective, given an algorithm's configuration, the aim is to find a subset of optimization problems on which this algorithm obtains good results. In the case of generalization of algorithms, the aim is to find an algorithm configuration which solves the major quantity of problems with different features [START_REF] Smit | Parameter tuning of evolutionary algorithms: Generalist vs. specialist[END_REF]. Both perspectives need to deal with the so-called no free lunch (NFL) theorem which proves that, under certain assumptions, no optimization algorithm is superior to any other on all possible optimization problems [START_REF] Wolpert | No free lunch theorems for optimization[END_REF].

In recent years, several tuning methods have been designed to search automatically the most appropriate configuration in the parameter values of the stochastic search based algorithms. One of the first was sequential parameter optimization (SPO) [START_REF]Sequential parameter optimization[END_REF] which is based on statistical techniques such as design and analysis of computer experiments. In [START_REF] Frank | Automatic algorithm configuration based on local search[END_REF]- [START_REF] Hoos | Automated Algorithm Configuration and Parameter Tuning[END_REF], the so-called ParamILS method is introduced. This approach consists of a steady state algorithm that uses an iterated local search to improve only categorical parameters using an initial configuration. The relevance estimation and value calibration (REVAC) method was proposed in [START_REF] Nannen | Relevance estimation and value calibration of evolutionary algorithm parameters[END_REF]. This approach works with a population of parameter configurations for estimating the distribution of the target algorithm. Here, the relevance of each parameter is determined according to the entropy measurement. Examples of tuning tools that work on categorical and numerical parameters at the same time are the iterated racing for configuration (I-RACE) framework [START_REF] López-Ibáez | The irace package: Iterated racing for automatic algorithm configuration[END_REF] and the evolutionary algorithm called EVOCA [START_REF] Riff | A New Algorithm for Reducing Metaheuristic Design Effort[END_REF]. I-RACE samples a set of configurations according to a particular distribution and incorporates an iterative racing procedure for selecting the best configurations. EVOCA is a steady state EA that, at each iteration generates two new individuals via wheel-crossover and a local search procedure. It can find good parameter values without requiring an in-depth knowledge of parameter tuning methods.

III. UNCERTAINTY

The typical goal in optimization is to identify optimal solutions. In the case of parameters settings, we want to optimize the configuration c ∈ C of an MOEA in an scenario O. As defined in Section II, O contains several problems o 1..k to be optimized, one at a time. However, the optimal configuration for solving one problem o i can be poor for another problem o j of one scenario O. Thus, the additional challenge is to have a set of configurations that have a 'good' performance in a set of scenarios. The main issue is how to define this performance. For that sake, we propose to use robustness measures, some of the concepts of robust optimization [START_REF] Gorissen | A practical guide to robust optimization[END_REF], to solve the OPT problem represented by Equation [START_REF] Beyer | Robust optimization a comprehensive survey[END_REF].

max c∈C,o∈O ψ A (c, o) (1)
In order to do it, we substitute the nominal objective for a robust measure that aggregates on the results from the different problems measured by the scalar fitness function ψ. In the following, we present the robust measures that we use in this work.

• Mean: it is one of the most commonly used measures in the literature. With this measure, we are looking for the configuration that works best in the mean of the cases. However, this measure is not appropriate when the results contain outliers since an outlier value has a big impact on the mean value itself. • Median: it is quite useful since it removes outliers.

However, removing outliers would mean that we are not interested in all cases but instead, in most of them. Thus, using this measure would mean that we are interested in the configuration that works well for most scenarios but it could fail completely in the others. • Worst case: in this case, the fitness would be represented by the worst result in the given scenario. This would optimize over worst cases and thus ensure that the configuration would work at least with that quality. However, this measure can be over-conservative, since the worst case could never happen in practice. Each of those measures has some advantages and drawbacks. The decision to use them should include the preference of the decision maker as well as the aim of the algorithm that we are tuning.

Moreover, since the MOEA has uncertainty itself (the same configuration can give different quality), it is also required to treat it in the same form as for the scenario. We first approximate the fitness for the MOEA and afterwards, for the scenario. Namely, we combine: mean-mean, median-median, best-worst, mean-worst, median-worst and worst-worst. We emphasize the worst-case since we are aiming for a configuration that works in all problem instances.

IV. MULTI-OBJECTIVE OPTIMIZATION A MOP can be represented as follows:

Minimize f (x) := f 1 (x), f 2 (x), . . . , f m (x) (2) subject to x ∈ S, (3)
where S ⊂ IR n is the feasible region and each decision variables vector x ∈ S is related to an objective vector f (x) ∈ IR m . The solution to a MOP, called the Pareto optimal set (POS), is given by:

P OS := {x ∈ S : ∃y ∈ S, y ≺ x}, (4)
and its image in IR m , is called the Pareto optimal front (POF). Some MOPs have complicated characteristics that cause difficulties to MOEAs for converging to the P OS such as a large number of decision variables (large-scale MOPs) or a large number of objective functions (many-objective problems). Furthermore, there are several Pareto front shapes such as linear, convex, concave, mixed, disconnected and degenerated geometries (see Fig. 1). Additionally, multi-frontal problems can cause a premature convergence for MOEAs, which results in the generation of false Pareto fronts.

There are three main classes of MOEAs: (1) based on Pareto dominance, (2) based on a performance indicator, and (3) based on decomposition. The last class has an efficient search ability to deal with many-objective optimization problems in a low computational time. However, the performance of decomposition-based MOEAs strongly depends on a scalarizing function and a set of target directions. The most popular decomposition-based MOEA is MOEA/D [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], which have been widely used and modified with the aim of improving its performance when solving a wide variety of MOPs.

MOEA/D transforms a MOP into several single-objective subproblems which are solved in a collaborative manner. Each solution from the population is associated with a search direction (weight vector) to optimize a transformation method called the scalarizing function (SF). Two goals should be accomplished during the search process. The first goal is to minimize the distance between a candidate solution and the reference point (typically, the ideal vector) in order to achieve convergence towards the POF. The second goal is the use of different search directions defined by a uniform distribution around all regions in the objective space. It aims to promote diversity for covering the entire POF. MOEA/D establishes neighborhoods at each target direction and applies iteratively a mating selection mechanism and evolutionary operators to enhance individuals. After that, an update process modifies the reference points and the current population according to the SF values. The first version of MOEA/D used the simulated binary crossover (SBX) and polynomial-based mutation (PM).

V. SCALARIZING FUNCTIONS A scalarizing function (SF) is a mathematical programming technique that transforms a MOP into a single-objective one. Its goal is to combine a vector of objective functions y = [f 1 , . . . , f m ] T with a reference direction λ = {λ 1 , . . . , λ m } to obtain a scalar value g(y) : IR m → IR.

This paper is focused on SFs that can be stated as:

minimize g(f (x); λ) (5) subject to x ∈ S, (6)
where λ is a weight vector that must satisfy m i=1 λ i = 1 and λ i ≥ 0 for all i ∈ {1, . . . , m}. Each component value represents the relative importance assigned to each objective function.

f (x) is a transformation function (see equation ( 7)) to handle negative or incommensurable objective functions.

f (x) := f (x) -z * , (7)
where

z * = [z 1 , • • • , z m ]
T is the ideal point that denotes the optimum value of the i th objective functions of a MOP. An estimation of the ideal point can be obtained by minimizing each of the objective functions individually.

Next, we describe a set of SFs commonly used in decomposition-based MOEAs, and we provide some information that aims to provide a better understanding of the impact of their model parameters for solving MOPs with Scalarizing -Function (ASF) [START_REF] Wierzbicki | The Use of Reference Objectives in Multiobjective Optimization[END_REF] Augmented max

y i λ i + α i y i λ i
small ASF (AASF) [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF] α > 0 Penalty d1 + θd2, θ ∈ (0, ∞) Boundary where d1 := y • λ λ suggested Intersection (PBI) [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] and d2 := y -d1 λ λ θ = 5

Exponential i e p λ i -1 e p y i p ∈ [1, ∞) Weighted Criteria (EWC) [START_REF] Athan | A Note on Weighted Criteria Methods for Compromise Solutions in Multi-objective Optimization[END_REF] Vector Angle The parameter α in AASF should take small values. In [START_REF] Tutum | A Multimodal Approach for Evolutionary Multi-objective Optimization (MEMO): Proof-of-Principle Results[END_REF], it was recommended to set α ≈ 10 -4 . MOEAs that adopt this scalarizing function are [START_REF] Tutum | A Multimodal Approach for Evolutionary Multi-objective Optimization (MEMO): Proof-of-Principle Results[END_REF], [START_REF] Ruiz | A Preference-based Evolutionary Algorithm for Multiobjective Optimization: The Weighting Achievement Scalarizing Function Genetic Algorithm[END_REF]. The EWC function requires a large value of p to achieve Pareto optimality, but this can lead to numerical overflow [START_REF] Marler | Survey of Multi-objective Optimization Methods for Engineering[END_REF]. In [START_REF] Carpinelli | Exponential Weighted Method and a Compromise Programming Method for Multiobjective Operation of plug-in Vehicle Aggregators in Microgrids[END_REF], EWC was used to solve a problem related to a voltage distribution network. For the VADS function, t = 100 was recommended by its authors. The PBI function can produce uniformly distributed solutions in objective space by setting appropriate values for θ. There are several studies [START_REF] Sato | Inverted PBI in MOEA/D and its Impact on the Search Performance on Multi and Many-Objective Optimization[END_REF]- [START_REF] Ishibuchi | Relation between weight vectors and solutions in moea/d[END_REF] that provide a sensitivity analysis of PBI, varying θ ∈ [0.1, 100].

Figure 2 illustrates some boxplots that represent the effect of the parameter sensitivity in the AASF and PBI functions. Here, we use the MOEA/D framework with the same parameters settings, only varying the SF. We compute the normalized hypervolume on 30 independent runs when solving the DTLZ1, DTLZ2 and DTLZ2 -1 test problems. We can see that one parameter value can be appropriate for a particular Pareto front shape and can work poorly for others. This effect occurs in the same way if we vary the MOP's dimensionality. Such sensitivity to the parameters settings justifies the need to find and optimize the 'best' configuration when dealing with several MOPs using the same MOEA configuration.

VI. EXPERIMENTAL METHODOLOGY

The goal of our experiments was to study the effect that robustness measures cause in the OPT methods. To achieve this goal, we used the descriptive statistics mentioned in Section III and the tuning tool called evolutionary calibrator (EVOCA) was adopted to find the suitable scalarizing function (SF) and its model parameters to reach the maximum performance within the MOEA/D framework. Our two particular goals were: 1) to verify which measures are the most useful in three scenarios characterized according to the Pareto front geometry, and 2) to identify the most robust configurations when scaling up the number of objective functions in a MOP.

According to the tuple presented in Section II, we use the next elements in our experimental OPT: Laumanns-Zitzler (DTLZ) test suite [START_REF] Huband | A Review of Multiobjective Test Problems and a Scalable Test Problem Toolkit[END_REF] and the Lamé superspheres (LS) test problems [START_REF] Emmerich | Test Problems Based on Lamé Superspheres[END_REF] varying the parameter γ = {0.3, 0.5, 1.0, 2.0, 4.0} with the goal of achieving diverse Pareto front shapes. We classified our test MOPs into three different scenarios according to their Pareto front geometry as follows.

1

) convex = {DTLZ3 -1 , LS (γ = 0.3), LS (γ = 0.5)} 2) linear = {DTLZ1, LS (γ = 1.0)} 3) concave = {DTLZ3, LS (γ = 2.0), LS (γ = 4.0)}
We tested this approach for m = 3, 5, 7 objectives. The number of decision variables was set to n = m + 4.

• ψ := the normalized hypervolume (NHV), defined as:

N HV (A) := HV i r i , (8)
where ) , r] is the hypervolume indicator [START_REF] Zitzler | Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications[END_REF]. L measures convergence towards the PF and maximum spread through the union of hypercubes formed by all non-dominated elements in A and a reference point r := (r 1 , . . . , r m ) T . It is set to (2, . . . , 2) T in all our test problems. A high N HV value is better.

HV (A, r) = L ∪ µ i=1 [A (i
• b max := The maximum budget established in the EVOCA algorithm is 10, 000 function evaluations and ten seeds per test problem at each generation.

Table II shows the parameter values used by MOEA/D, where m is the number of objectives, n is the number of decision variables, H is a parameter used by the simplex lattice design (SLD) method to generate uniform weight vectors. popsize and N F E are the population size and number of function evaluations adopted. B is the neighborhood size per weight vector. Algorithm 1 shows the general structure of EVOCA while Algorithm 2 shows how we introduce the different robustness measures. Mean, median or worst are aggregation functions that could be used instead of σ or ω. In EVOCA, the chromosome is represented by a string where each element corresponds to a parameter, and its value. Thus, the string length is the number of parameters to be tuned. In our experiments, we adopted 6 SFs and 4 model parameters. The population size is computed considering the number of parameters to be tuned and their initial domain sizes. The population size adopted in EVOCA is 25.

Σ ← ∅ 3 foreach o ∈ O do 4 Ω ← ∅ 5 foreach s ∈ S do 6 apf ← A(c, o, s) // MOEA execution 7 Ω ← Ω ∪ ψ(apf ) 8 Σ ← Σ ∪ ω(Ω) 9 return σ(Σ)
EVOCA uses two operators: a wheel-crossover operator that constructs one child from the whole population. It uses a roulette wheel procedure [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] to select the value of the gene of each offspring. The child generated replaces the worst individual on the current population. The crossover procedure is performed at each iteration. The mutation operator is a hill climbing first improvement procedure, which takes a copy of the child generated by the crossover operator and tries to improve it by modifying one of its parameter values. The mutation operator is always applied. When a numerical parameter is selected, it tries to randomly take a new value within a continuous range that represents the parameter domain. The mutated child generated replaces the second worst individual on the current population, in case of producing a better performance.

VII. RESULTS

Table III shows the scalarizing functions obtained by EVOCA for each of the approaches in the different scenarios related to the Pareto front shapes. It is interesting to see that the approaches select not only different parameters for the scalarizing function but also different functions in several cases. This is an interesting result since it shows the impact of changing the robustness measure that is adopted.

For all scenarios, the experimental results are shown in Table IV. In gray, we show the best result among the different robustness measures. An arrow upwards indicates that the approach is outperformed in a statistically significant way by the baseline algorithm (EVOCA with mean-mean). An arrow downwards means that the baseline algorithm outperforms the robustness measure. Furthermore, Figures 3, 4 and5 show the box plots for the different scenarios using the hypervolume indicator without applying a normalization process. Marks a, b, c, d, e, f correspond to each robustness measures presented in Table IV. We can see in Figure 3 that EVOCA's recommendation for convex scenarios outperforms the baseline version only in the Lamé supersphere with α = 0.3 and α = 0.5 using median-median and best-worst measures. But in multifrontal MOPs, the median-median fails and other measures obtain

DTLZ3X, 07D

HV Fig. 3: Scenario with convex geometry and 3, 5, 7 objectives similar results. This is because the median statistics discard the outlier results. Figure 4 shows the linear scenario, where there is an evident trend: the best-worst measure outperforms the baseline version in MOPs with more than 3 objectives. Finally, we obtain similar results in Figure 5, where EVOCA's recommendation improved the baseline version only in unimodal MOPs.

From these results, it is possible to observe that the bestworst approach, is capable of outperforming the baseline algorithm in most cases when the problems are unimodal. A similar case occurs with the median. However, for the multimodal problems we can observe a deterioration of the quality of the results. This suggests that such problems should be in a different scenario and have their own configuration to find good solutions.

Also, we can notice that the configurations found work well when increasing the number of objectives. Furthermore, we can observe that the different robustness measures outperform the baseline algorithm in most of the problems (except for the multi-modal problems). This suggests that alternative measures to treat uncertainty can have a positive impact while searching for configurations for a set of scenarios.

VIII. CONCLUSIONS AND FUTURE WORK

In this work, we have presented a comparative study to analyze the effect of different robustness measures in the offline parameter tuning procedure applied to MOEA/D. We found that AASF showed the best results for the convex case, VADS for the concave case and PBI for the linear case. This shows evidence that the shape of the Pareto front has an effect on the scalarizing function to be used. Furthermore, the values for the scalarizing functions are different among the robustness measures. We have observed that the median and best-worst cases have interesting results for the unimodal problems. Thus, they are worth studying in more depth.

As part of our future work, we are interested in scaling this study regarding the number of objectives in order to find robust configurations that scale up with a higher number of objectives. Furthermore, we intend to increase the number of test problems and their characteristics. Our preliminary results suggest that multi-modal problems deserve a special scenario. Finally, we plan to use more sophisticated robustness measures and to adopt them as a starting point for the design of selfadaptive models.

Fig. 1 :

 1 Fig. 1: Different Pareto front shapes

Fig. 2 :

 2 Fig. 2: Examples of parameter sensitivity for the AASF and PBI functions. Figures a) and b) show the normalized hypervolume indicator for linear (DTLZ1), concave (DTLZ3) and convex (DTLZ3 -1 ) Pareto fronts.

• A :

 : The multi-objective evolutionary algorithm based on decomposition (MOEA/D). • P := {CHE, ASF, AASF, EW C, V ADS, P BI} with their corresponding model parameter values defined in the ranges α ∈ [0, 10.0], p ∈ [0.1, 10.0], t ∈ [1, 100], and θ ∈ [0.1, 50.0]. • O := DTLZ1, DTLZ3, DTLZ3 -1 from the Deb-Thiele-

Algorithm 1 : 6 Evaluate

 16 EVOCA(A,P ) Input : Metaheuristic A, Parameters P Output: Population 1 Pop = Generate Initial Population; 2 Evaluate each configuration ∈ Pop using Algorithm 2; 3 while not termination criterion met do 4 New-Pop = Pop; 5 Child = Wheel-crossover(Pop); Child using Algorithm 2 ;

7 10 if 11 Replace 2

 10112 Replace worst individual by Child in New-Pop; 8 Mut-child = Mutation(Child); 9 Evaluate Mut-child using Algorithm 2; Mut-Child better than child then nd worst individual by Mut-child in New-Pop 12 Pop=New-Pop;

Algorithm 2 :

 2 EvaluateConf(c, O, S, ψ, σ, ω, A) ∈ R Input : Configuration c, Scenario O, Seeds S Fitness func. ψ, Aggregation func. σ, ω, Algorithm A Output: evaluate quality of one configuration 1 begin 2

HVFig. 4 :HVFig. 5 :

 45 Fig. 4: Scenario with linear geometry and 3, 5, 7 objectives

L

  γ = 4.0 03D 0.89376(1.8281e-05) ↑ 0.89419(1.3943e-05) ↑0.89419(1.2420e-05) ↓0.89328(1.2611e-05) ↑0.89376(1.8281e-05) 0.89377(1.8035e-05) 05D 0.97715(1.1496e-05) ↑ 0.97754(8.2635e-06) ↑0.97753(1.0646e-05) ↓0.97682(1.0701e-05) ↑0.97715(1.1496e-05) ↑0.97716(1.2953e-05) 07D 0.99482(5.2873e-05) ↑ 0.99500(3.4305e-05) ↑0.99494(6.3347e-05) ↓0.99479(4.2571e-05) ↑0.99482(5.2873e-05) 0.99483(4.2165e-05) DTLZ3 03D 0.92623(8.8163e-04) ↓0.91765(6.6721e-03) ↓0.92330(1.9209e-03) 0.92632(6.7448e-04) ↑0.92623(8.8163e-04) 0.92643(1.0158e-03) 05D 0.98925(3.7990e-04) ↓0.66321(4.5688e-01) ↓0.98887(4.2865e-04) 0.98921(3.2035e-04) ↑ 0.98925(3.7990e-04) 0.98907(4.3676e-04) 07D 0.99051(3.6548e-02) ↓0.67700(4.3885e-01) ↓0.96755(1.1470e-01) ↑ 0.99781(3.7914e-04) ↑0.99051(3.6548e-02) 0.93313(2.0597e-01)

TABLE I :

 I Weighted and unconstrained scalarizing functions.

	SF	Minimize g(f (x); λ) :=	Parameter values
	Chebyshev (CHE) [22]	maxi λi|yi|	-
	Achievement	max	y i λ i

  Pareto front shapes. TableIshows six SFs compatible with Pareto dominance which can handle convex, linear and concave Pareto front shapes. Both CHE and ASF generate at least weakly Pareto optimal solutions and do not require additional parameters. The AASF, PBI, EWC and VADS functions can avoid the generation of weakly Pareto optimal solutions by introducing additional model parameters. However, these parameters are very sensitive to the Pareto front shape and the MOP's dimensionality.

			y			t > 0
	Distance Scaling (VADS) [26]	λ λ	•	y y	t	suggested t = 100

different

TABLE II :

 II The parameters settings used by MOEA/D. The mark p in column H means that the original set of weight vectors generated by SLD is pruned in order to obtain the desirable population size.

	m	n	H	popsize	N F E	B
	3	7	15	136	60,000 27
	5	9	6p	180	80,000 36
	7	11	4	210	90,000 42

TABLE III :

 III EVOCA's recommendation for different scenarios using 3, 5 and 7 objectives

	Geo	mean-mean	median-median	best-worst	worst-worst	mean-worst	median-worst
	Convex	AASF (α = 5.3727) EWC (p = 7.2)	CHE	AASF (α = 1.5305) AASF (α = 0.6977) AASF (α = 1.4065)
	Concave	VADS (p = 11.9)	PBI (θ = 8.2)	PBI (θ = 2.6)	VADS (p = 6.3)	VADS (p = 11.9)	VADS (p = 12.1)
	Linear	PBI (θ = 15.9)	PBI (θ = 10.4)	AASF (α = 0.0469)	PBI (θ = 4.2)	PBI (θ = 8.4)	PBI (θ = 11.3)

TABLE IV :

 IV Results reported by EVOCA with different robust measures for different scenarios (2.8481e-07) 0.99994(4.0573e-07) ↑ 0.99994(3.3306e-16) ↓0.99994(9.7262e-07) ↓0.99994(4.9379e-07) ↑0.99994(2.5289e-07) DTLZ1 03D 0.99681(3.2364e-05) 0.99681(2.5263e-05) ↑ 0.99683(1.1418e-05) 0.99681(2.3043e-05) 0.99680(2.7473e-05) 0.99681(2.3141e-05) 05D 0.99995(3.9159e-06) 0.99995(3.7555e-06) ↑ 0.99996(5.0619e-07) ↓0.99995(3.9932e-06) 0.99995(2.3360e-06) ↑0.99995(3.7648e-06) 07D 0.99999(2.5451e-06) 0.99999(3.7277e-06) ↑ 0.99999(4.3969e-07) ↓0.99997(1.1753e-05) ↓0.99998(4.1686e-06) 0.99999(2.9046e-06) 03D 0.92821(2.0129e-06) 0.92821(1.2526e-06) ↑ 0.92821(1.3764e-06) ↑ 0.92821(1.0095e-06) ↑0.92821(2.0129e-06) 0.92821(2.0137e-06) 05D 0.99017(1.2045e-06) 0.99017(9.6148e-07) ↑ 0.99017(6.5248e-07) 0.99017(1.0589e-06) ↑0.99017(1.2045e-06) 0.99017(1.1618e-06) 07D 0.99863(1.5886e-06) ↑0.99863(1.0217e-06) ↑ 0.99863(7.1644e-07) 0.99863(2.2111e-06) ↑0.99863(1.5886e-06) 0.99863(1.9602e-06)

	Sc Geo	Dim	mean-mean	median-median	best-worst	worst-worst	mean-worst	med-worst
			03D 0.99949(7.6989e-06) ↑ 0.99994(2.4944e-08) ↑0.9968(7.6713e-07) ↓0.99937(1.6849e-05) ↓0.99916(2.2555e-05) ↓0.99935(1.4162e-05)
		L γ = 0.3	05D 0.99978(3.6113e-06) ↑ 1.0000(4.0000e-08) ↑0.99999(4.4408e-16) ↓0.99976(5.9568e-06) ↓0.99973(1.1476e-05) ↓0.99976(6.4822e-06)
			07D 0.99982(8.5306e-06) ↑ 1.0000(1.7950e-08) ↑ 1.0000(0.0000e+00) ↓0.99981(8.6050e-06) ↓0.99979(1.0696e-05) ↓0.99981(9.7633e-06)
	Convex	L γ = 0.5	03D 0.99778(3.6635e-06) ↓0.99758(4.4021e-07) ↓0.99557(3.0983e-06) ↓0.99768(9.3282e-06) ↓0.99748(1.5878e-05) ↓0.99767(8.3159e-06) 05D 0.99976(9.1288e-06) ↑ 0.99997(5.1207e-08) ↑0.99993(2.5377e-07) ↓0.99973(1.0707e-05) ↓0.99968(1.9492e-05) ↓0.99973(1.0752e-05) 07D 0.99981(2.0659e-05) ↑ 0.99999(3.3306e-16) ↑0.99999(4.4457e-06) ↓0.99978(2.2967e-05) ↓0.99975(2.4902e-05) ↓0.99977(2.8233e-05)
			03D 0.99335(2.0098e-04) ↓0.0000(0.0000e+00) ↓0.99010(3.0459e-04) 0.99334(1.6041e-04)	0.99340(1.3607e-04)	0.99342(7.1454e-05)
		DTLZ3 -1	05D 0.99991(1.4226e-05) ↓0.0000(0.0000e+00) ↓0.99978(1.0693e-05) ↓0.99989(2.1022e-05) ↓0.99983(4.3120e-05) ↓0.99989(2.0928e-05)
			07D 0.99996(1.4068e-05) ↓0.0000(0.0000e+00) ↓0.99983(2.1933e-04) ↓0.99995(2.5509e-05) ↓0.99994(2.3587e-05) ↓0.99995(2.0053e-05)
			03D 0.97481(9.1966e-07) 0.97481(3.1155e-07) ↓0.97481(1.1706e-06) 0.97481(2.3795e-07)	0.97481(2.6042e-07) 0.97481(3.3757e-07)
	L γ = 1.0 07D 0.99994Concave 05D 0.99886(1.4922e-07) 0.99886(1.4083e-07) ↑ 0.99886(4.1831e-07) 0.99886(2.8511e-07) Linear L γ = 2.0	0.99886(1.5205e-07)	0.99886(1.3743e-07)

ACKNOWLEDGMENT

The first author acknowledges support from IPN project no. SIP 20181450, CONACyT and CINVESTAV-IPN to pursue graduate studies in Computer Science. The last author gratefully acknowledges support from the CONACyT grant no. 2016011920 (Investigación en Fronteras de la Ciencia 2016).


Sreeja S Nair
email: sreeja.nair@lip6.fr

Gustavo Petri
email: gustavo.petri@arm.com

Marc Shapiro
email: marc.shapiro@acm.org

Invariant Safety for Distributed Applications

Keywords: Replicated data, Consistency, Automatic verification, Distributed application design, Tool support

We study a proof methodology for verifying the safety of data invariants of highlyavailable distributed applications that replicate state. The proof is (1) modular: one can reason about each individual operation separately, and (2) sequential: one can reason about a distributed application as if it were sequential. We automate the methodology and illustrate the use of the tool with a representative example.

INTRODUCTION

A distributed application often replicates its data to several locations, and accesses the closest available replica. Examples include social networks, multi-user games, co-operative engineering tools, collaborative editors, source control repositories, or distributed file systems. To ensure availability, an update must not synchronise across replicas; otherwise, when a network partition occurs, the system will block. Asynchronous updates may cause replicas to diverge or to violate the data invariants of the application.

To address the first problem, Conflict-free Replicated Data Types (CRDTs) [START_REF] Shapiro | Conflict-free replicated data types[END_REF] have mathematical properties to ensure that all replicas that have received the same set of updates converge to the same state [START_REF] Shapiro | Conflict-free replicated data types[END_REF]. To ensure availability, a CRDT replica executes both queries and updates locally and immediately, without remote synchronisation. It propagates its updates to the other replicas asynchronously.

There are two basic approaches to update propagation: to propagate operations, or to propagate states. In the former approach, an update is first applied to some origin replica, then sent as an operation to remote replicas, which in turn apply it to update their local state. Operation-based CRDTs require the the message delivery layer to deliver messages in causal order, exactly once; the set of replicas must be known.

In the latter approach, an update is applied to some origin replica. Occasionally, one replica sends its full state to some other replica, which merges the received state into its own. In turn, this replica will later send its own state to yet another replica. As long as every update eventually reaches every replica transitively, messages may be dropped, re-ordered or duplicated, and the set of replicas may be unknown. Replicas are guaranteed to converge if the set of states, as a result of updates and merge, forms a monotonic semi-lattice [START_REF] Shapiro | Conflict-free replicated data types[END_REF]. Due to these relaxed requirements, state-based CRDTs have better adoption [START_REF] Baquero | Composition in state-based replicated data types[END_REF]. They are the focus of this work.

As a running example, consider a simple auction system. The state of an auction consists of status, a set of bids, and a winner. This state is replicated at multiple servers; CRDTs ensures that all replicas eventually converge. Users at different locations can start an auction, place bids, close the auction, declare a winner, inspect the local replica, and observe if a winner is declared and who it is. All replicas will eventually agree on the same auction status, same set of bids and the same winner.

However, the application may also require to maintain a correctness property or invariant over the data. An invariant is an assertion on application data that must evaluate to true in every state of every replica. For instance, the auction's invariant is that: when the auction is closed, there is a winner; there is a single winner; and the winner's bid is the highest.

Such an invariant is easy to ensure in a sequential system, but concurrent updates might violate it. In this case, the application would need to synchronise some updates between replicas, in order to maintain the invariant. For instance, in the absence of sufficient synchronisation, a replica might close the auction and declare a winner, while concurrently a user at a different replica is placing a higher bid.

This problem has been addressed before, by stating correctness rules and proof obligations; however, previous work considers only the operation-based approach [START_REF] Gotsman | Cause I'm Strong Enough: Reasoning about consistency choices in distributed systems[END_REF][START_REF] Marcelino | Bringing hybrid consistency closer to programmers[END_REF][START_REF] Najafzadeh | The CISE tool: Proving weaklyconsistent applications correct[END_REF].

In this paper, we propose a proof methodology for applications that use statebased CRDTs. We exploit the properties of state-based CRDTs to reason about a concurrent system in a sequential manner. We have also developed a tool named Soteria, to automate our proof rule. Soteria detects concurrency bugs and provides counterexamples.

SYSTEM MODEL

An application consists of state, some operations, a merge function, and an invariant. The state is replicated at any number of replicas. A client chooses any arbitrary replica for its next operation, called the origin replica for that operation. A replica occasionally sends its state to some other replica, which the receiving replica merges into its own state. In summary, the state of any given replica changes, either by executing an update operation for which it is the origin, or by merging the state received from a remote replica. Each replica is sequential. A merge is the only point where a replica observes concurrent operations submitted to other replicas. Each replica executes An application invariant is an assertion over state. The invariant must evaluate true in every state of every replica. Despite being evaluated against local state, an invariant is in effect global, since it must be true at all replicas, and replicas eventually converge.

Figure 1 depicts the evolution of state in our auction application. Each line represents a replica, time progressing from left to right. A box represents local state. A curved arrow represents an update operation, labelled with the operation name. A diagonal arrow shows propagation (labelled with the propagated state), merged at the receiving replica.

We assume here that application state is a composition of CRDTs. This is not a limitation, since many basic CRDT types have been proposed, which extend familiar sequential data types with a concurrency semantics [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF].

The CRDT convergence rules, or lattice rules, are the following. The state of a replica progresses monotonically with time. The set of states forms a semi-lattice, i.e., is equipped with a partial order and a least-upper-bound function. A state transition represents the execution of either an update operation or merge. An update is an inflation, i.e., the resulting state is no less (in the partial order) than the previous one. Merge computes a state which is the least-upper-bound of the current local state and the received remote state.

We write σ for the current state of a replica, σ new for the new state after an operation or merge, and σ ′ for the incoming state from a remote replica. In the Boogie specification language (used by our tool, described in Section 4), we denote operation execution with the keyword call, an assertion with assert, and assumptions with assume. An operation op executes only if its precondition is true. Thus, we can write the above lattice conditions as follows:

• An update operation op is an inflation.

call σ new = op(σ ) assert σ new ≥ σ • Merge is a least upper bound call σ new = merge(σ , σ ′ ) assert σ new ≥ σ ∧ σ new ≥ σ ′ # upper bound assert ∀σ * , σ * ≥ σ ∧ σ * ≥ σ ′ =⇒ σ * ≥ σ new # least
Let us illustrate with the auction example (for simplicity we consider a single auction). Its state is as follows:

• Status: the status of an auction can move from its initial state, INVALID (under preparation), to ACTIVE (can receive bids) and CLOSED (no more bids accepted), such that INVALID < ACTIVE < CLOSED. • Winner: The winner of the auction. It is either ⊥, initially, or the bid with the highest amount. In case of multiple bids with same amount, the bid with the lowest id wins. It is ordered such that ∀b ∈ Bids, ⊥ ≤ b. • Bids: Set of bids placed (initially empty)

-BidId: A unique identifier for each bid placed -Placed: A boolean flag to indicate whether the bid has been placed or not, ordered TRUE > FALSE. It is enabled once when the bid is placed. Once placed, a bid cannot be withdrawn. -Amount: An integer representing the amount of the bid; this cannot be modified once the bid is created. Figure 1 illustrates how the auction state evolves over time; status, bids and winner are represented by a circle, rectangle and a star respectively. The application state is geo-replicated at data centres in Australia, Belgium and China.

We now specify the merge operation for an auction. We denote the receiving replica state σ = (status, winner, Bids); the received state is denoted σ ′ = (status ′ , winner ′ , Bids ′ ). In the absence of any extra synchronisation, it is possible to violate the invariant, as the following execution scenario illustrates. Alice from Australia starts an auction, by setting its status to ACTIVE (green in the figure) in the Australian replica. Henceforth, the auction can receive bids. The Australian replica sends its updated state to the other two replicas, which merge it into their own states. Now Bob in Belgium places a bid for $100 (blue). This update is sent to other replicas. Charlie, in China, sees the auction and Bob's bid. He updates the China state with a higher bid of $105 (dotted blue), which is sent to both other replicas. However, due to a network failure, the remote replicas do not receive this update (red dotted lines). Meanwhile Alice, unaware of Charlie's bid, closes the auction and declares Bob's bid as the winner. Later when the network heals, the updated states are sent and merged. The auction is new closed, and contains Bob's $100 bid and Charlie's $105 bid. Unfortunately, Bob's bid is the winner, violating the application invariant.

In the next section, we can discuss how to ensure invariants of applications build on top of the system model we described.

PROVING INVARIANTS

As explained earlier, each replica executes a sequence of state transitions, due either to a local update, or to a merge incorporating remote updates. Thus, concurrency can be observed only through merge.

Let us call safe state a replica state that satisfies the invariant. Assuming the current state is safe, any update (local or merge) must result in a safe state. To ensure this, every update is equipped with a precondition that disallows any unsafe execution. 1Thus, a local update executes only when, at the origin replica, the current state is safe and its precondition currently holds. Similarly, merge executes only with two safe states that together satisfy a merge precondition.

Formally, an update u (an operation or a merge), mutates the local state σ , to a new state σ new = u(σ ). To preserve the invariant, Inv, we require that

σ ∈ Pre u =⇒ u(σ ) ∈ Inv
To illustrate local preconditions, consider an operation close_auction(w: BidId), which sets auction status to CLOSED and the winner to w. The developer may have written a precondition such as status = ACTIVE, because closing an auction doesn't make sense otherwise. In order to ensure the invariant that the winner has the highest amount, one needs to strengthen it with the clause is_highest(Bids, w), defined as ∀b ∈ Bids : b.Amount ≤ w.Amount.

To illustrate merge precondition, consider a CRDT whose state is the pair of integers, σ = (n, m) ∈ N × N. It has two operations, incn and incm, that respectively increment n or m by 1, and a merge function:

merge(σ , σ ′ ) = (max(n, n ′ ), max(m, m ′ ))
We wish to maintain the invariant that their sum is no more than 10:

Inv ≜ (n + m) ≤ 10
The precondition of incn is Pre incn ≜ (n + m) ≤ 9; similarly for incm. Starting from a safe state (4, 5), two replicas may independently increment to states (5, 5) and (4, 6) respectively. Both are safe. However, merging them would violate the invariant. Therefore, merge(σ , σ ′ ) must have precondition

Pre merge ≜ max(n, n ′ ) + max(m, m ′ ) ≤ 10
Since merge can happen at any time, it must be the case that its precondition is always true, i.e., it constitutes an additional invariant. Now our global invariant consists of two parts: first, the application invariant, and second, the precondition of merge. We can now state our proof rule informally as follows:: 2State-Based Safety Rule. Define the precondition of merge to be the weakestprecondition of merge, for the application invariant. The initial state must satisfy, and each local update or merge operation must preserve, the conjunction of: (i) the application invariant, and (ii) the precondition of merge. □ In our Boogie notation, each operation can be verified as follows:

assume Inv ∧ Pre merge ∧ Pre op call σ new = op(σ ) assert Inv ∧ Pre merge
The case of the merge function can be verified with the following condition:

assume Inv ∧ Inv ′ ∧ Pre merge call σ new = merge(σ , σ ′ ) assert Inv ∧ Pre merge
Note that there are two copies of state, the unprimed local state of the replica applying the merge, and the primed state received from a remote replica. Inv ′ denotes that σ ′ preserves the invariant Inv.

Applying the proof rule

Let us apply the proof methodology to the auction application. Its invariant is the following conjunction:

(1) A bid is placed only when status is ACTIVE.

(2) And: Once a bid is placed, its amount does not change.

(3) And: There is no winner until status is CLOSED. (4) And: There is a single winner, the bid with the highest amount (breaking ties using the lowest identifier). Computing the weakest-precondition of each update operation, for this invariant, is obvious. For instance, as discussed earlier, close_auction(w: BidId) gets precondition is_highest(Bids, w), because of Invariant Term 4 above.

Despite local updates to each replica preserving the invariant, Figure 1 showed that it is susceptible of being violated by merging. This is the case if Bob's $100 bid in Belgium wins, even though Charlie concurrently placed a $105 bid in China; this occurred because status became CLOSED in Belgium while still ACTIVE in China. The weakest-precondition of merge for Term 4 expresses that, if status in either states is CLOSED, the winner should be the bid with the highest amount in both the states. Therefore, merge(σ , σ ′ ) must have the following additional precondition:

status = CLOSED =⇒ is_highest ( Bids , winner ) ∧ is_highest ( Bids ′ , winner ) ∧ status ′ = CLOSED =⇒ is_highest ( Bids , winner ′ ) ∧ is_highest ( Bids ′ , winner ′ )
Furthermore, the code for merge uses Term 2, for which its weakest-precondition is as follows: These two merge preconditions now strengthen the global invariant, in order to preserve safety in concurrent executions. Let us now consider how this strengthening impacts the local update operations. Since starting the auction doesn't modify any bids, this operation trivially preserves it. Placing a bid might violate it, if the auction is concurrently closed in some other replica; conversely, closing the auction could violate it, if a higher bid is concurrently placed in a remote replica. Thus, the auction application is safe when executed sequentially, but is unsafe when updates are concurrent. This indicates the specification has a bug, which we now proceed to fix.

Concurrency Control for Invariant Preservation

As we discussed earlier, the preconditions of operations and merge are strengthened in order to preserve the invariant. This provides a sequentially safe specification. An application must also preserve the precondition of merge in order to ensure concurrent safety. Violating this indicates the presence of a bug in the specification. In that case, the developer needs to strengthen the application by adding appropriate concurrency control mechanisms, ie., the operations that fail to preserve the precondition of merge might need to synchronise. The required concurrency control mechanisms are added as part of the state in our model. The modified application state is now composed of the CRDTs that represents the state and the concurrency control mechanism. Together, it behaves like a composition of state-based CRDTs. The whole state should now ensure the lattice conditions described in section 2.

Recall that in the auction example, placing bids and closing the auction were not preserving the precondition of merge. This requires strengthening the specification by adding a concurrency control mechanism to restrict these operations. We can enforce them to be strictly sequential, thereby avoiding concurrency at all. But this will affect the availability of the application.

A concurrency control can be better designed with the workload characteristics in mind. For this particular application, we know that placing bids are very frequent operations than closing an auction. Hence we try to formulate a concurrency control like a readers-writer lock. In order to realise this we distribute tokens to each replica. As long as a replica has the token, it can allow placing bids. Closing the auction requires recalling the tokens from all replicas. This ensures that there are no concurrent bids placed and thus a winner can be declared, respecting the application safety.

The entire specification of the auction application can be seen in Figure 2. The shaded lines in blue indicate the effect of adding concurrency control to the state.

An alternative approach to our treatment of concurrency control could be to consider the invariant as a resource in the style of Concurrent Separation Logic [START_REF] O'hearn | Resources, concurrency, and local reasoning[END_REF]. In this case, access to the application state, described through a separation logic invariant, is guarded by a concurrency control mechanism (typically some form of a lock). However, this approach is tied to separation logic reasoning, where assertions act as resources, and allows one to distinguish local from global resources. We consider that this was not essential for the kind of proofs that we conduct, but it might be more promising when verifying client programs of our data types.

AUTOMATING THE VERIFICATION

In this section we present a tool to automate the verification of invariants as discussed in the previous sections. Our tool, called Soteria is based on the Boogie [START_REF] Barnett | Boogie: A modular reusable verifier for object-oriented programs[END_REF] verification framework. The input to Soteria is a specification of the application written in Boogie, an intermediate verification language.

A specification in Soteria will consist of the following parts:

• State: a declaration of the state. It can be a single CRDT or a composition of CRDTs.

• Comparison function:

The programmer provides a comparison function (annotated with keyword @gteq) that determines the partial order on states. • Operations: The programmer provides the implementation of the operations and their respective preconditions, Pre op . Operations are encoded either imperatively as Boogie procedures or declaratively as postconditions. • Merge function: The special merge operation is distinguished from the other operations (with annotation @merge). The programmer must provide a precondition to merge that is strong enough to prove the invariant. • Application Invariant: The programmer provides the invariant (with keyword @invariant) to be verified by the tool as a Boogie assertion over the state.

In addition, Boogie often requires additional information such as: • User-defined data types, • Constants to declare special objects such as the origin replica me, or to bound the quantifiers, • Axioms for inductive functions over aggregate data structures, for instance, to compute the maximum of a set of values, • Loop invariants. The specification of the auction application can be seen in Figure 2.

Verification

The verification of a specification is performed in multiple stages; in order:

(1) Syntactic checks: validates the specification for syntactical errors and checks whether the pre/post conditions are sound. (2) Compliance check: checks whether the specification provides all the elements explained earlier.

(3) Convergence check: checks whether the specification respects the properties of a state-based CRDT, ie., each operation inflates the state and merge is the least upper bound. (4) Safety check: verifies the safety of the application invariant, as discussed in section 3. This stage is divided further into two sub-stages:

• Sequential safety: whether each individual operation (or merge) upholds the invariant. If not, the designer needs to strengthen the precondition of the corresponding operation (or merge) • Concurrent safety: whether every operation (and merge) upholds the precondition of merge. Note that, while this check relates to the concurrent behaviour of state-based CRDTs, the check itself is completely sequential, ie., it does not require reasoning about operations performed by other processes. This check ensures that the invariant remains safe during concurrent operation. If this check fails, the application needs stronger concurrency control. Each check in Soteria3 generates counterexamples when the verification fails. These counterexamples might guide the developer in debugging the specification according to the verification steps.

RELATED WORK

Several works have concentrated on the formalisation and specification of eventually consistent systems [START_REF] Burckhardt | Principles of Eventual Consistency, volume 1 of Foundations and Trends in Programming Languages[END_REF][START_REF] Burckhardt | Replicated data types: Specification, verification, optimality[END_REF][START_REF] Sivaramakrishnan | Declarative programming over eventually consistent data stores[END_REF]. A number of works concentrate on the specification and correct implementation of CRDTs [START_REF] Gomes | A framework for establishing strong eventual consistency for conflict-free replicated datatypes[END_REF][START_REF] Jagadeesan | Eventual consistency for CRDTs[END_REF]. Our work also verifies the CRDT (lattice) conditions, but additionally verifies an arbitrary application invariant over a replica state.

Gotsman et al. [START_REF] Gotsman | Cause I'm Strong Enough: Reasoning about consistency choices in distributed systems[END_REF] provides a proof methodology for proving invariants of CRDTs that propagate operations. The associated tools [START_REF] Marcelino | Bringing hybrid consistency closer to programmers[END_REF][START_REF] Najafzadeh | The CISE tool: Proving weaklyconsistent applications correct[END_REF] performs the check using an SMT solver as the backend and Nair and Shapiro [START_REF] Nair | Improving the "Correct Eventual Consistency" tool[END_REF] discusses some concurrency control suggestions by using the counterexamples generated by the failed proofs. Gotsman et al. [START_REF] Gotsman | Cause I'm Strong Enough: Reasoning about consistency choices in distributed systems[END_REF] assume that the underlying network ensures causal consistency, and their methodology requires reasoning about concurrent behaviours. This requires checks for each pair of operations in the application (reflected as stability verification conditions). Gotsman et al. [START_REF] Gotsman | Cause I'm Strong Enough: Reasoning about consistency choices in distributed systems[END_REF] uses an abstract notion of tokens as concurrency control mechanisms. The operations acquire tokens in order to preserve the application invariant.

In contrast, Soteria focuses on state-based CRDTs. We check convergence by verifying the lattice conditions of section 2 and that because of the rules shown in section 3, we can reduce the problem of verifying the invariant to sequential proof obligations. This is reflected by the fact that all of our proofs are standard pre/post conditions checks using the Boogie framework. Boogie framework. In contrast with Gotsman et al. [START_REF] Gotsman | Cause I'm Strong Enough: Reasoning about consistency choices in distributed systems[END_REF], Soteria includes concrete specification of concurrency control as part of the application state.

To the best of our knowledge, ours is the first attempt in automated verification of invariants of state-based CRDTs.

CONCLUSION

We have presented a proof methodology to verify invariants of state-based CRDT implementations guaranteeing: (1) that the implementation satisfies the lattice conditions of state-based CRDTs [START_REF] Baquero | Composition in state-based replicated data types[END_REF], and (2) that the implementation satisfies programmer provided invariant reducing the problem to checking that each operation of the data type satisfies a precondition of the merge function of the state.

We implemented Soteria, a tool sitting on top of the Boogie verification framework, to specify the implementation, its invariant and validate it.

In future work, we plan to automate concurrency control synthesis. The synthesised concurrency control can be analysed and adjusted dynamically to minimise the cost of synchronisation. Another direction for future work can be to decouple the update propagation mechanism of CRDT from the proof rule resulting in a generic proof rule to verify distributed systems.

Figure 1 :

 1 Figure 1: Evolution of state in an auction application

  merge (( status , winner , Bids ), ( status ′ , winner ′ , Bids ′ )): status new := max ( status , status ′ ) if winner ′ ⊥ then winner new := winner ′ else winner new := winner ∀b , if b ∈ Bids ∩ Bids ′ then Bids new .b. placed := Bids .b. placed ∨ Bids ′ .b. placed else if b ∈ Bids ′ then Bids new .b. placed := Bids ′ .b. placed else Bids new .b. placed := Bids .b. placed ∀b , if b ∈ Bids then Bids new .b. amount := Bids .b. amount else Bids new .b. amount := Bids ′ .b. amount

  ∀b ∈ Bids∩Bids ′ , Bids .b. amount = Bids ′ .b. amount

Technically, this is at least the weakest-precondition of the update for the invariant. It strengthens any a priori precondition that the developer may have set.

We omit the full formalisation and the proof of soundness for brevity.

The tool along with some sample specifications can be accessed at https://github.com/sreeja/soteria_tool.

ACKNOWLEDGMENTS

The authors would like to thank the anonymous reviewers for their comments which helped in improving this paper. This research is supported in part by the RainbowFS project (Agence Nationale de la Recherche, France, number ANR-16-CE25-0013-01) and by European H2020 project 732 505 LightKone (2017-2020).