Pascale Bendotti

Pierre Fouilhoux

Cécile Rottner

Sub-symmetry-breaking inequalities and application to the Unit Commitment Problem

come

Introduction

We consider an Integer Linear Program (ILP) of the form (ILP ) min c(x) | x ∈ X , with X ⊆ P(m, n) and c : P(m, n) → R

where P(m, n) is the set of m × n binary matrices. A symmetry is defined as a permutation π of the indices {(i, j) | 1 ≤ i ≤ m, 1 ≤ j ≤ n} such that for any solution matrix x ∈ X , matrix π(x) is also solution and has same cost, i.e., π(x) ∈ X and c(x) = c(π(x)). The symmetry group G of (ILP ) is the set of all such permutations. It partitions the solution set X into orbits, i.e., two matrices are in the same orbit if there exists a permutation in G sending one to the other. In [START_REF] Bendotti | Orbitopal fixing for the full (sub-)orbitope and application to the unit commitment problem[END_REF], symmetries arising in solution subsets of (ILP ) are called sub-symmetries.

In this article, we focus on (sub-)symmetry groups containing all sub-column permutations of a given solution submatrix. These symmetry groups are assumed to be known or previously detected [START_REF] Liberti | Reformulations in mathematical programming: automatic symmetry detection and exploitation[END_REF][START_REF] Berthold | Detecting orbitopal symmetries[END_REF].

Various techniques, so called symmetry-breaking techniques, are available to handle symmetries in (ILP ). The general idea is, in each orbit, to pick one solution, defined as the representative, and then restrict the solution set to the set of all representatives. The most common choice of representative is based on the lexicographical order. Column y ∈ {0, 1} m is said to be lexicographically greater than column z ∈ {0, 1} m if there exists i ∈ {1, ..., m -1} such that ∀i ′ ≤ i, y i ′ = z i ′ and y i+1 > z i+1 , i.e., y i+1 = 1 and z i+1 = 0. We write y ≻ z (resp. y z) if y is lexicographically greater than z (resp. greater than or equal to z). A technique is said to be full symmetry-breaking (resp. partial symmetry-breaking) if the solution set is exactly (resp. partially) restricted to the representative set.

Symmetry-breaking techniques can be based on specific branching and pruning rules during the B&B search [START_REF] Margot | Exploiting orbits in symmetric ILP[END_REF][START_REF] Ostrowski | Orbital branching[END_REF][START_REF] Bendotti | Orbitopal fixing for the full (sub-)orbitope and application to the unit commitment problem[END_REF]. Alternatively, variable aggregation can be a very efficient technique [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF]. In general, aggregating symmetric variables in (ILP ) yields a relaxation of (ILP ) [START_REF] Fischetti | Orbital shrinking: Theory and applications[END_REF]. When the integer decomposition property [START_REF] Baum | Integer rounding and polyhedral decomposition for totally unimodular systems[END_REF] holds, i.e., the aggregated solution can be disaggregated, aggregated (ILP ) has same optimal value as (ILP ).

Other symmetry-breaking techniques rely on full or partial symmetry-breaking inequalities. Techniques based on symmetry-breaking inequalities are flexible, i.e., at any node of the B&B tree, the branching rule can be derived from any linear inequality on the variables. Most techniques based on branching and pruning rules are either full symmetry-breaking or flexible. Variable fixing [START_REF] Kaibel | Orbitopal fixing[END_REF][START_REF] Bendotti | Orbitopal fixing for the full (sub-)orbitope and application to the unit commitment problem[END_REF] is both full symmetry-breaking and flexible. The size of the LP solved at each node of the branching tree is generally invariant with pruning and variable fixing techniques. This size is decreased (resp. increased) by the use of variable aggregation (resp. symmetry-breaking inequalities) techniques.

Symmetry-breaking inequalities can be derived from the linear description of the convex hull of an arbitrary representative set [START_REF] Friedman | Fundamental Domains for Integer Programs with Symmetries[END_REF]. In most works, each chosen representative x is lexicographically maximal in its orbit, i.e., x g(x), for each g ∈ G. The convex hull of the representative set is called the symmetry-breaking polytope [START_REF] Hojny | Polytopes associated with symmetry handling[END_REF]. When x is a matrix and symmetry group G acts on the columns of x, the symmetry-breaking polytope is called orbitope. Even if complete linear descriptions for symmetry-breaking polytopes may be hard to reach in general, integer programming formulations for these polytopes still yield full symmetry-breaking inequalities [START_REF] Hojny | Polytopes associated with symmetry handling[END_REF]. Instead of considering orbits of solutions, [START_REF] Liberti | Reformulations in mathematical programming: automatic symmetry detection and exploitation[END_REF][START_REF] Liberti | Stabilizer-based symmetry breaking constraints for mathematical programs[END_REF] introduce inequalities enforcing a lexicographical order within orbits of variables.

If symmetry group G is the symmetric group S n acting on the columns, the chosen representative x of an orbit may be such that its columns x(1), ..., x(n) are lexicographically non-increasing, i.e., for all j < n, x(j) x(j + 1). The convex hull of all m × n binary matrices with lexicographically non-increasing columns is called the full orbitope [START_REF] Kaibel | Packing and partitioning orbitopes[END_REF]. Full symmetry-breaking inequalities are introduced by Friedman [START_REF] Friedman | Fundamental Domains for Integer Programs with Symmetries[END_REF]:

m i=1 2 m-i x i,j ≥ m i=1 2 m-i x i,j+1 , ∀j ∈ {1, ..., n -1} (1)
As the 2 m-i term might cause numerical intractability, alternative inequalities featuring binary coefficients can be used, at the expense of losing the full symmetry-breaking property. An option is to use column inequalities introduced in [START_REF] Kaibel | Packing and partitioning orbitopes[END_REF]:

i k=1
x k,j ≥ x i,j+1 , ∀j ∈ {1, ..., n -1}

Another option is to use a partial-symmetry-breaking form of Friedman inequalities, as in [START_REF] Jans | Solving lot-sizing problems on parallel identical machines using symmetry-breaking constraints[END_REF][START_REF] Lima | Symmetry breaking in MILP formulations for Unit Commitment problems[END_REF]:

m i=1 x i,j ≥ m i=1 x i,j+1 , ∀j ∈ {1, ..., n -1} (3)
The latter inequalities enforce that the total number of ones in each column is non-increasing. Note that this does not guarantee lexicographically non increasing columns for the representatives.

For the particular case of packing (resp. partition) problems, i.e., problems whose solution matrix features at most (resp. exactly) one 1-entry in each row, a class of full symmetry-breaking inequalities is introduced in [START_REF] Kaibel | Packing and partitioning orbitopes[END_REF]. These inequalities lead to a complete linear description of two special cases of full orbitopes: the packing (resp. partitioning) orbitope, i.e., the convex hull of all m × n binary matrices with lexicographically non-increasing columns and with at most (resp. exactly) one 1-entry per row. For the full orbitope, an O(mn 3 ) extended formulation is given in [START_REF] Kaibel | Branched polyhedral systems[END_REF]. To the best of our knowledge, it has never been used in practice to handle symmetries. A complete linear description of the full orbitope in the x variable space seems hard to reach [START_REF] Loos | Describing Orbitopes by Linear Inequalities and Projection Based Tools[END_REF]. For the full orbitope restricted to 2-column matrices, a complete linear description in the x space is available [START_REF] Loos | Describing Orbitopes by Linear Inequalities and Projection Based Tools[END_REF].

The Unit Commitment Problem (UCP) has demonstrated to be a good candidate to apply symmetry-breaking techniques [START_REF] Ostrowski | Modified orbital branching for structured symmetry with an application to unit commitment[END_REF][START_REF] Lima | Symmetry breaking in MILP formulations for Unit Commitment problems[END_REF][START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF][START_REF] Bendotti | Orbitopal fixing for the full (sub-)orbitope and application to the unit commitment problem[END_REF]. The Min-up/min-down Unit Commitment Problem (MUCP) is to find a minimum-cost power production plan on a discrete time horizon for a set of production units. At each time period, the total production has to meet a forecast demand. Each unit must satisfy minimum up-time and down-time constraints besides featuring production and start-up costs. In real-world UCP, some more technical constraints have also to be taken into account, such as ramp constraints limiting the increase and decrease in generated power at each time period. Classical formulations of the MUCP feature two sets of binary variables: up variables x and start-up variables u. In practical instances, there are several sets of identical units, i.e., units with identical characteristics, which induce symmetries. Indeed, assuming a solution is expressed as a matrix where column j corresponds to the up/down trajectory of unit j over the time horizon, then any permutation of columns corresponding to identical units leads to another solution with same cost. Moreover, in some subproblems, there exist symmetries not contained in the symmetry group of the original problem. Indeed, units which are ready to start (resp. shut down) at a given time t can permute their plans from time t to T .

We propose a general framework to build full symmetry-breaking inequalities in order to handle sub-symmetries arising from solution subsets whose symmetry groups contain the symmetric group acting on some sub-columns. One additional variable per subset Q considered may be needed in these inequalities, depending whether variables x are sufficient to indicate that "x belongs to subset Q". The proposed framework is applied to derive such inequalities when the symmetry group is the symmetric group S n acting on the columns. It is also applied to derive inequalities breaking both symmetries and sub-symmetries in the MUCP. We present experimental results comparing these sub-symmetry-breaking inequalities to state-of-the-art symmetry-breaking formulations, such as the MUCP formulation featuring inequality (3) [START_REF] Lima | Symmetry breaking in MILP formulations for Unit Commitment problems[END_REF], aggregated (x, u) or aggregated interval MUCP formulations [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF]. When the MUCP is considered, the integer decomposition property holds for the (x, u) formulation and thus efficient aggregation techniques apply [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF]. When the ramp-constrained MUCP is considered, the integer decomposition property does not hold anymore for the (x, u) MUCP formulation, there the corresponding aggregated solutions can no longer be disaggregated. We show that our inequalities outperform all above mentioned formulations in the ramp-constrained case.

In Section 2, we describe the framework. In Section 3, an application to the symmetric group case is presented. In Section 4, the framework is applied to derive sub-symmetry-breaking inequalities dedicated to the MUCP. Experimental results on MUCP instances (with or without ramp constraints) are presented in Section 5.

For a given solution subset Q, the symmetry group G Q of the corresponding subproblem is different from G and may contain symmetries undetected in G. In practice it is too expansive to compute the symmetry group for every subset Q ⊂ X . However for many problems, symmetries of G can be deduced from the problem's structure, and so can symmetries of G Q , for some particular solution subsets Q. In this case, symmetries of G Q are a priori known, and thus do not need to be computed. Such symmetries may be handled together with symmetries of G. In this section, we introduce subsymmetry-breaking inequalities designed to simultaneously handle symmetries and sub-symmetries in symmetric groups.

First, we provide an overview of sub-symmetries in integer programming.

Background on sub-symmetries

The sub-symmetry concepts introduced in [START_REF] Bendotti | Orbitopal fixing for the full (sub-)orbitope and application to the unit commitment problem[END_REF] are briefly recalled in this section. Consider a subset Q ⊂ X of solutions of (ILP ). The sub-symmetry group G Q relative to subset Q is defined as the symmetry group of subproblem min{cx | x ∈ Q}. Permutations in sub-symmetry group G Q are referred to as sub-symmetries.

Let {Q s ⊂ X , s ∈ {1, ..., q}} be a set of solution subsets. To each Q s , s ∈ {1, ..., q}, corresponds a sub-symmetry group G Qs . Let O s k , k ∈ {1, ..., o s }, be the orbits defined by G Qs on subset Q s , s ∈ {1, ..., q}, and O = {O s k , k ∈ {1, ..., o s }, s ∈ {1, ..., q}}. For given x ∈ P(m, n), let us define G(x) = Qs∋x G Qs , the set of all permutations π in q s=1 G Qs such that π can be applied to X.

Matrix x ′ is said to be in relation with x ∈ P(m, n) if there exist r ∈ N and permutations π 1 , ..., π r such that π k ∈ G(π k-1 ...π 1 (x)), ∀k ∈ {1, ..., r}, and x ′ = π 1 π 2 ...π r (x). The generalized orbit O of x with respect to {Q s , s ∈ {1, ..., q}} is thus the set of all x ′ in relation with x. By definition, for any generalized orbit O, there exist orbits σ 1 , ..., σ p ∈ O such that O = ∪ p s=1 σ s . The set of representatives {ρ(σ), σ ∈ O} is said to be orbit-compatible if for any generalized orbit O = ∪ p s=1 σ s , σ 1 , ..., σ p ∈ O, there exists j such that ρ(σ j ) = ρ(σ s ) for all i such that ρ(σ j ) ∈ σ s . Such a solution ρ(σ j ) is said to be a generalized representative of O.

Given x ∈ X and sets R ⊂ {1, ..., m} and C ⊂ {1, ..., n}, we consider submatrix (R, C) of x, denoted by x(R, C), obtained by considering columns C of x on rows R only. Symmetry group G Q is the sub-symmetric group with respect to (R, C) if it is the set of all permutations of the columns of x(R, C). If G Q is the sub-symmetric group with respect to (R, C) then subset Q is said to be sub-symmetric with respect to (R, C).

Consider a set S of solution subsets Q s , s ∈ {1, ..., q}, such that each subset Q s , s ∈ {1, ..., q}, is sub-symmetric with respect to (R s , C s ). For each orbit O s k , k ∈ {1, ..., o s } of G Qs , s ∈ {1, ..., q}, its representative x s k ∈ O s k is chosen to be such that submatrix x s k (R s , C s ) is lexicographically maximal, i.e., its columns are lexicographically non-increasing. Such x s k is said to be the lex-max of orbit O s k with respect to (R s , C s ).

Property 1 ([4]

). The set of representatives {x s k , k ∈ {1, ..., o s }, s ∈ {1, ..., q}} is orbit-compatible. The full sub-orbitope P sub (S) associated to S is the convex hull of binary matrices x such that for each s ∈ {1, ..., q}, if x ∈ Q s then the columns of x(R s , C s ) are lexicographically non-increasing.

Definition and validity

Consider a set S of solution subsets Q s , s ∈ {1, ..., q}, such that each subset Q s , s ∈ {1, ..., q}, is sub-symmetric with respect to (R s , C s ). Consider integer variable z s , s ∈ {1, ..., q}, such that z s = 0 if variable x ∈ Q s , and such that z s ≥ 1 if x ∈ Q s . For any x ∈ X , function Z associates x to a vector Z(x) such that z s , s ∈ {1, ..., q}, is the s th component of Z(x) denoted by Z s (x)

Note that in many cases, function Z is linear, i.e., each integer variable z s is a linear expression of variables x. In such cases, no additional variable z s is needed. In some cases where function Z is not linear, variable z s can be linearly expressed from variables x using only a few additional inequalities or integer variables.

Given c, c ′ ∈ C s such that c < c ′ , the sub-symmetry-breaking inequality, denoted by (Q s (c, c ′ )), is defined as follows.

x r1,c ′ ≤ z s + x r1,c

where

r 1 = min(R s ) (4)
For each orbit O s k , k ∈ {1, ..., o s }, of G Qs , s ∈ {1, ..., q}, the chosen representative is the lexmax of orbit O s k with respect to (R s , C s ). Then by Property 1, this set of representatives is orbit-compatible. In particular, solution set X can be restricted to the set of representatives by considering its intersection with the full sub-orbitope P sub (S). If x ∈ Q s , inequality (Q s (c, c ′ )) enforces that the first row of submatrix x(R s , C s ) is lexicographically non-increasing, hence the following result.

Lemma 1 (Validity). If x ∈ P sub (S), then (x, Z(x)) satisfies inequality (Q s (c, c ′ )) for each s ∈ {1, ..., q} and c, c ′ ∈ C s such that c < c ′ .
Note that an inequality similar to (4) applied to a row of R s distinct from r 1 may not be valid when used alongside with (4), as shown in Example 1.

Example 1. Let S = {Q 1 }, q = 1, where subset Q 1 is as follows.

Q 1 = x ∈ P(4, 3) ∩ X | 3 c=1 x 2,c = 3
Let us suppose the symmetry group of Q 1 is the sub-symmetric group with respect to submatrix ({3, 4}, {1, 2, 3}). Variable z 1 can be defined using equality

z 1 = 3 - 3 c=1 x 2,c . Note that z 1 = Z 1 (x) = 0 when 3 c=1 x 2,c = 3, i.e., x ∈ Q 1 , and is positive otherwise. Here the first row in R 1 is r 1 = min(R 1 ) = 3, thus given c, c ′ ∈ {1, 2, 3}, c < c ′ , inequality (Q 1 (c, c ′ )) is as follows x 3,c ′ ≤   3 - 3 j=1 x 2,j   + x 3,c (5)
Inequality (Q 1 (c, c ′ )) enforces that row 3 of a solution matrix x is lexicographically ordered, i.e., x 3,1 ≥ x 3,2 ≥ x 3,3 , whenever 3 c=1 x 2,c = 3. Now consider solutions x 1 , x 2 ∈ Q 1 : x 1 =     1 0 0 1 1 1 1 0 0 0 1 1     and x 2 =     1 0 0 1 1 1 0 0 1 1 1 0    
Inequality (5) cuts off solution x 2 from the feasible set. Inequality (6) corresponds to inequality (4) applied to row 4:

x 4,c ′ ≤   3 - 3 j=1 x 2,j   + x 4,c (6)
Inequality (6) would cut off x 1 . This shows that inequalities (5) and ( 6) cannot be used simultaneously.

Note that in the general case, inequalities (4) may only be partial-symmetry-breaking. Indeed, for given s ∈ {1, ..., q} and c, c 4) is not sufficient to select the lexmax representatives.

′ ∈ C s such that c < c ′ , inequality (Q s (c, c ′ )) only enforces that the first row of submatrix x(R s , C s ) is lexicographically non-increasing when x ∈ Q s . In the case when x r1,c ′ < x r1,c , then sub-columns x(R s , {c ′ }) ≺ x(R s , {c}). Otherwise, when x r1,c ′ = x r1,c , inequality (
To enforce a lexicographical order, subsequent rows of submatrix x(R s , C s ) should be considered until a tie-break row is found. It is shown in the next section that inequalities (Q s (c, c ′ )) for all s ∈ {1, ..., q} and c < c ′ ∈ C s enforce that x ∈ P sub (S) provided a tie-break condition on set S is fulfilled.

Full symmetry-breaking sufficient condition

In this section, we introduce a condition for inequalities (4) to be full symmetry-breaking.

For each s ∈ {1, ..., q}, consider

R s = {r s 1 , ..., r s |Rs| } and C s = {c s 1 , ..., c s |Cs| }, where r s 1 < ... < r s |Rs| and c s 1 < ... < c s |Cs| .
For given s ∈ {1, ..., q} and any two columns

c s l-1 , c s l ∈ C s , if there is a solution x ∈ Q s such that columns c s l-1 and c s l are equal from row r s 1 to row r s k-1 , it must be ensured that row r s k is lexicographically non increasing, i.e., x r s k ,c s l-1 ≤ x r s k ,c s l .
The key idea is to exhibit another set Q p ∈ S for quartet (Q s , k, l, x), such that Q p contains x and is sub-symmetric with respect to (R p , C p ), where the first row of R p is r s k and C p contains columns c s l-1 and c s l . Then inequality

(Q p (c s l-1 , c s l )) will ensure that x r s k ,c s l-1 ≥ x r s k ,c s l .
For each quartet (Q s , k, l, x), the existence of such a subset Q p in S will be ensured by tie-break condition (C), defined as follows:

(C)    ∀s ∈ {1, ..., q}, ∀k ∈ {2, ..., |R s |}, ∀l ∈ {2, ..., |C s |} If x ∈ Q s such that x r s k ′ ,c s l-1 = x r s k ′ ,c s l , ∀k ′ ∈ {1, ..., k -1}, then there exists p ∈ {1, ..., q} such that x ∈ Q p , C p ⊇ {c s l-1 , c s l } and r s k = min(R p ) If condition (C) holds, inequalities (Q s (c s l-1 , c s l ))
, ∀s ∈ {1, ..., q}, ∀l ∈ {2, ..., |C s |} exactly restrict the solution set to the representative set X ∩ P sub (S). They are therefore full symmetry breaking, with respect to the sub-symmetries defined by S. This gives the idea of the proof for the following theorem.

Theorem 1. If condition (C) holds, then:

(i) (x, Z(x)) satisfies (Q s (c s l-1 , c s l )), ∀s ∈ {1, ..., q}, ∀l ∈ {2, ..., |C s |} (ii) x ∈ P sub (S) are equivalent.
For general set S, condition (C) may not hold. Fortunately, it will be shown that we can construct from S another set S satisfying (C) and such that P sub ( S) = P sub (S).

The idea is to divide each Q s , s ∈ {1, ..., q} in smaller subsets such that for each row r s k ∈ R s and each column c s l ∈ C s , there is a subset Q, which is sub-symmetric with respect to (R, C) = ({r s k , ..., r s |Rs| }, {c s l-1 , c s l }). Set S is defined as follows.

S = Q s (k, l) | s ∈ {1, ...,
Q s (k, l) = x ∈ Q s | x r,c s l-1 = x r,c s l , ∀r ∈ {r s 1 , ..., r s k-1 }
Note that for solution x ∈ Q s such that columns c s l-1 and c s l are equal from row r s

1 to row r s k-1 , the set exhibited for quartet (Q s , k, l, x) is Q s (k, l). Note also that Q s (1, l) = Q s , l ∈ {2, ..., |C s |}.
We thus have the following result:

Lemma 2. Set S satisfies (C) and is such that P sub ( S) = P sub (S).

Proof. The symmetry group of

Q s (k, l) is the sub-symmetric group with respect to (R, C) = ({r s k , ..., r s |Rs| }, {c s l-1 , c s l }). Thus if some solution x ∈ Q s is such that columns c s l-1 and c s l are equal from row r s 1 to row r s k-1 , then subset Q s (k, l) contains x and is such that C ⊇ {c s l-1 , c s l } and min(R) = r s k . Condition (C) is therefore satisfied by S.
It can be readily checked that the full sub-orbitopes defined by S and S are the same.

It follows, from Theorem 1, that inequalities (Q(c, c ′ )), c < c ′ ∈ C, Q ∈ S are full symmetrybreaking with respect to the sub-symmetries defined by S.

Corollary 1. If for each Q ∈ S, (x, Z(x)) satisfies inequality (Q(c, c ′ )), ∀c < c ′ ∈ C, then x ∈ P sub (S).
We can then consider S instead of S. This implies to add one inequality at least (resp. one variable at most), per subset Q ∈ S, i.e., O(qmn) inequalities at least (resp. variables at most).

Example. Referring to Example 1, S = Q 1 (1, l), Q 1 (2, l), l ∈ {2, 3} . For each l ∈ {2, 3}, Q 1 (1, l) = Q 1 as for any s, Q s (k, l) = Q s whenever k = 1. We also have Q 1 (2, l) = x ∈ Q 1 | x 3,l-1 = x 3,l . For each l ∈ {2, 3}, z l associated to subset Q 1 (2,
l) can be expressed as follows:

z l = 2z 1 + (x 3,l-1 -x 3,l ). Indeed, when z 1 = 0, inequality (5) becomes x 3,l-1 ≤ x 3,l . Thus, z l = 0 if x 3,l-1 = x 3,l and z l ≥ 1 otherwise. When z 1 = 1, z l ≥ 1.
Hence the following inequalities are full symmetry-breaking:

x 3,l-1 ≤ 3 -3 j=1 x 2,j + x 3,l ∀l ∈ {2, 3} x 4,l-1 ≤ 6 + x 3,l-1 -x 3,l -2 3 j=1 x 2,j + x 4,l ∀l ∈ {2, 3}
In Sections 3 and 4, inequalities (4) are built in a more straightforward way, in the sense that set S already satisfies condition (C) in the two applications studied.

Application to the symmetric group case

In this section, we apply the framework of Section 2 to any problem whose symmetry group G is the symmetric group S n acting on the columns. The collection S S of subsets considered will lead to inequalities restricting any solution x ∈ X to be in the full orbitope. These inequalities feature variables z which can be explicitly expressed from x with O(mn) linear inequalities. Here, the sub-symmetries considered are restrictions of symmetries' actions to solution subsets.

A complete linear description of the 2-column full orbitope, featuring additional integer variables, is proposed in [START_REF] Loos | Describing Orbitopes by Linear Inequalities and Projection Based Tools[END_REF]. In the general n-column case, we show that these inequalities can also be derived using the framework described in Section 2, and can be used as full symmetry-breaking inequalities.

We consider S S = Q i,j , i ∈ {0} ∪ {1, ..., m -1}, j ∈ {2, ..., n} , where

Q i,j = x ∈ X | x i ′ ,j-1 = x i ′ ,j ∀i ′ ∈ {1, ..., i} .
Subset Q i,j is the set of feasible solutions such that columns j -1 and j are equal from row 1 to row i. Note that Q 0,j = X . The symmetry group of Q i,j is then the sub-symmetric group with respect to (R i , {j -1, j}) where R i = {i + 1, ..., m}. It can be readily checked that in this case, S already satisfies condition (C).

Let variable z i,j be such that z i,j = 0 if x ∈ Q i,j and 1 otherwise. Note that for all j ∈ {2, ..., n}, Q 0,j = X , thus z 0,j = 0, ∀x ∈ X . Note also that X ∩ P sub (S S ) is a subset of the full orbitope.

Thus, given that the columns of any x ∈ X ∩ P sub (S S ) are in a non-increasing lexicographical order, function Z is such that Z(x) = z, where z satisfies the following linear inequalities.

             z 1,j-1 = x 1,j-1 -x 1,j ∀j ∈ {2, ..., n} (7a)
z i,j-1 ≤ z i-1,j-1 + x i,j-1 ∀i ∈ {2, ..., m}, j ∈ {2, ..., n} (7b)
z i,j-1 + x i,j ≤ 1 + z i-1,j-1 ∀i ∈ {2, ..., m}, j ∈ {2, ..., n} (7c)
x i,j-1 ≤ z i,j-1 + x i,j ∀i ∈ {2, ..., m}, j ∈ {2, ..., n} (7d)

z i-1,j-1 ≤ z i,j-1 ∀i ∈ {2, ..., m}, j ∈ {2, ..., n} (7e)
Constraint (7a) sets variable z 1,j-1 to 1 whenever columns j -1 and j are different and in a non-increasing lexicographical order on row 1, and to 0 when they are equal. Constraints (7b) and (7c) set variable z i,j-1 to 0 when z i-1,j-1 = 0 and columns j -1 and j are equal on row i. Constraint (7d) sets variable z i,j-1 to 1 if columns j -1 and j are different and in a non-increasing lexicographical order on row i. Constraint (7e) sets z i,j-1 to 1 when variable z i-1,j-1 = 1, i.e., when columns j -1 and j are different before row i.

For each i ∈ {0, ..., m-1} and j ∈ {2, ..., n} inequality (4) is inequality (Q i,j (j -1, j)) as follows:

x i+1,j ≤ z i,j-1 + x i+1,j-1 ∀i ∈ {1, ..., m}, ∀j ∈ {2, ..., n}
It ensures that if columns j -1 and j of x are equal from row 1 to i, then row i + 1 is in a non-increasing lexicographical order. Note that if z i-1,j-1z i,j-1 = -1 then necessarily x i,j = 0. Thus inequality ((Q i,j (j -1, j))) can be lifted to

x i,j ≤ (2z i-1,j-1 -z i,j-1 ) + x i,j-1 (8)
In the special case when n = 2, by replacing variable z i,j by y i,j where z i,j = 1 -i i ′ =1 y i,j , for each i ∈ {1, ..., m}, j ∈ {1, 2}, inequalities (7a)-( 8) yield the complete linear description of the 2-column full orbitope proposed in [START_REF] Loos | Describing Orbitopes by Linear Inequalities and Projection Based Tools[END_REF].

In the general n-column case, inequalities (7a)-( 8) are still full symmetry-breaking (by Theorem 1), and then can be used in practice to restrict the feasible set to any full orbitope. In this case, O(mn) additional variables and constraints are needed.

Application to the Unit Commitment Problem

The framework of Section 2 is now applied to the Unit Commitment Problem, which features many sub-symmetries non detected by the symmetry group G.

Given a discrete time horizon T = {1, ..., T }, a demand for electric power D t is to be met at each time period t ∈ T . Power is provided by a set N of n production units. At each time period, unit j ∈ N is either down or up, and in the latter case, its production is within [P j min , P j max ]. Each unit must satisfy minimum up-time (resp. down-time) constraints, i.e., it must remain up (resp. down) during at least L j (resp. ℓ j ) periods after start up (resp. shut down). Each unit j also features three different costs: a fixed cost c j f , incurred each time period the unit is up; a start-up cost c j 0 , incurred each time the unit starts up; and a cost c j p proportional to its production. The Min-up/min-down Unit Commitment Problem (MUCP) is to find a production plan minimizing the total cost while satisfying the demand and the minimum up and down time constraints. The MUCP is strongly NP-hard [START_REF] Bendotti | On the complexity of the unit commitment problem[END_REF].

In the real-world Unit Commitment Problem (UCP), some more technical constraints have also to be taken into account, such as ramp constraints or reserve requirement constraints, and the start-up costs are an exponential function of the unit downtime. From a combinatorial point of view, the MUCP is the core structure of the UCP. In this section, we study the MUCP with and without ramp constraints.

For each unit j ∈ N and time period t ∈ T , let us consider three variables: x t,j ∈ {0, 1} indicates if unit j is up at time t; u t,j ∈ {0, 1} whether unit j starts up at time t; and p t,j ∈ R is the quantity of power produced by unit j at time t. Without loss of generality we consider that L j , ℓ j ≤ T . Formulation F (x, u) for the MUCP is as follows [START_REF] Rajan | Minimum up/down polytopes of the unit commitment problem with start-up costs[END_REF][START_REF] Ostrowski | Modified orbital branching for structured symmetry with an application to unit commitment[END_REF][START_REF] Bendotti | The min-up/min-down unit commitment polytope[END_REF].

(F (x, u)) min x,u,p n j=1 T t=1 c j f x t,j + c j p p t,j + c j 0 u t,j s. t. t t ′ =t-L j +1 u t ′ ,j ≤ x t,j ∀j ∈ N , ∀t ∈ {L j , ..., T } (9)
t t ′ =t-ℓ j +1 u t ′ ,j ≤ 1 -x t-ℓ j ,j ∀j ∈ N , ∀t ∈ {ℓ j , ..., T } (10)
u t,j ≥ x t,j -x t-1,j ∀j ∈ N , ∀t ∈ {2, ..., T } (11)
P j min x t,j ≤ p t,j ≤ P j max x t,j ∀j ∈ N , ∀t ∈ T (12)

n j=1 p t,j ≥ D t ∀t ∈ T (13)
x t,j , u t,j ∈ {0, 1} ∀j ∈ N , ∀t ∈ T

For convenience, we will also use variable w t,j = x t-1,jx t,j + u t,j , indicating whether unit j shuts down at time t.

Symmetries and sub-symmetries in the UCP

Symmetries in the MUCP (and in the UCP) arise from the existence of groups of identical units, i.e., units with identical characteristics (P min , P max , L, ℓ, c f , c 0 , c p ). The instance is partitioned into types h ∈ {1, ..., H} of n h identical units. The unit set of type h is denoted by N h = {j h 1 , ..., j h n h }. The solutions of the MUCP can be expressed as a series of binary matrices. For a given type h, we introduce matrix x h ∈ P(T, n h ) such that entry x h t,k corresponds to variable x t,j h k , where j h k is the index of the k th unit of type h, k ∈ {1, ..., n h }. Column j of matrix x h corresponds to the up/down plan relative to the j th unit of type h. Similarly, we introduce matrices u h and p h . The set of all feasible x = (x t,j ) t∈T ,j∈N is denoted by X MUCP . Note that any solution matrix x (resp. u, p) can be partitioned in H matrices x h (resp. u h , p h ). Since all units of type h are identical, their production plans can be permuted, provided that the same permutation is applied to matrices x h , u h and p h . Thus, the symmetry group G contains the symmetric group S n h acting on the columns of x h , for each unit type h. Consequently, for each type h, feasible solutions x h can be restricted to be in the T × n h full orbitope. As binary variables u are uniquely determined by variables x, breaking the symmetry on x variables will break the symmetry on u variables. Note that this restriction to the T × n h full orbitope for each type h can possibly be done using inequalities from Section 3 featuring z variables.

There are also other sources of symmetry, arising from the possibility of permuting some subcolumns of matrices x h . For example, consider two identical units. Suppose at some time period t, these two units are down and ready to start up. Then their plans after t can be permuted, even if they do not have the same up/down plan before t.

More precisely, a unit j ∈ N is ready to start up at time t ∈ {1, ..., T } if and only if ∀t ′ ∈ {tℓ j , ..., t -1}, x t ′ ,j = 0. Similarly, a unit j ∈ N k is ready to shut down at time t ∈ {1, ..., T } if and only if ∀t ′ ∈ {t -L j , ..., t -1}, x t ′ ,j = 1.

Sub-symmetry-breaking inequalities for the MUCP

For each time period t ∈ {1, ..., T } and any two consecutive units j h k , j h k+1 of type h, k ∈ {1, ..., n h -1}, consider the following subsets of X MUCP :

Q t k,h = x ∈ X MUCP | x t ′ ,j = 0, ∀t ′ ∈ {t -ℓ h , ..., t -1}, ∀j ∈ {j h k , j h k+1 } Q t k,h = x ∈ X MUCP | x t ′ ,j = 1, ∀t ′ ∈ {t -L h , ..., t -1}, ∀j ∈ {j h k , j h k+1 }
where ℓ h (resp. L h ) is the minimum down (resp. up) time of units of type h.

Note that

Q t k,h and Q t k,h are different from subsets Q i,j defined in Section 3. Actually, Q t,j h k+1 ⊂ Q t k,h and Q t,j h k+1 ⊂ Q t k,h . Let G Q t k,h and G Q t k,h
be the sub-symmetry groups associated to

Q t k,h and Q t k,h , t ∈ {1, ..., T }, h ∈ {1, ..., H}, k ∈ {1, ..., n h -1}. The sub-symmetries in G Q t k,h (resp. G Q t k,h
) are called start-up sub-symmetries (resp. shut-down sub-symmetries). Most of these sub-symmetries are not detected in the symmetry group of the MUCP.

Groups

G Q t k,h
and

G Q t k,h
contain the sub-symmetric group associated to the submatrix defined by rows and columns ({t, ..., T }, {j h k , j h k+1 }). Applying results from Section 2, variables z t k,h and z t k,h , indicating whether x ∈ Q t k,h and x ∈ Q t k,h respectively, can be directly derived from variables x and u:

z t k,h = x t-ℓ h ,j ′ + t-1 t ′ =t-ℓ h +1 u t ′ ,j ′ + x t-ℓ h ,j + t-1 t ′ =t-ℓ h +1 u t ′ ,j z t k,h = 1 -x t-L h ,j ′ + t-1 t ′ =t-L h +1 w t ′ ,j ′ + 1 -x t-L h ,j + t-1 t ′ =t-L h +1 w t ′ ,j
where j = j h k and j ′ = j h k+1 for sake of clarity.

Consider S MUCP = Q t k,h , Q t k,h , t ∈ {1, .
.., T }, h ∈ {1, ..., H}, k ∈ {1, ..., n h -1} . In this case, set S directly satisfies condition C.

For each h ∈ {1, ..., H}, k ∈ {1, ..., n h -1} and t ∈ {1, ..., T }, inequalities (Q t k,h (j, j ′ )) and (Q t k,h (j, j ′ )), where j = j h k and j ′ = j h k+1 , are as follows.

x t,j ′ ≤ x t-ℓ h ,j ′ + t-1 t ′ =t-ℓ h +1 u t ′ ,j ′ + x t-ℓ h ,j + t-1 t ′ =t-ℓ h +1 u t ′ ,j + x t,j x t,j ′ ≤ 1 -x t-L h ,j ′ + t-1 t ′ =t-L h +1 w t ′ ,j ′ + 1 -x t-L h ,j + t-1 t ′ =t-L h +1 w t ′ ,j + x t,j
Strengthening symmetry-breaking inequalities Inequalities (Q t k,h (j, j ′ )) and (Q t k,h (j, j ′ )) can be further strengthened, using the relationship between variables x and u.

First note that by definition of variables w:

x t,j ′ -x t-ℓ h ,j ′ + t-1 t ′ =t-ℓ h +1 u t ′ ,j ′ = u t,j ′ - t t ′ =t-ℓ h +1 w t ′ ,j ′ x t,j + 1 -x t-L h ,j + t-1 t ′ =t-L h +1 w t ′ ,j = -w t,j + 1 + t t ′ =t-L h +1 u t ′ ,j
As if u t,j ′ = 1 (resp. w t,j = 1), then t t ′ =t-ℓ h +1 w t ′ ,j ′ = 0 (resp.

t t ′ =t-L h +1 u t ′ ,j = 0)
, the following Start-Up-Ready and Shut-Down-Ready inequalities are valid and stronger than inequalities

(Q t k,h (j, j ′ )) and (Q t k,h (j, j ′ )). u t,j ′ ≤ x t-ℓ h ,j + t-1 t ′ =t-ℓ h +1 u t ′ ,j + x t,j (15)
w t,j ≤ 1 -x t-L h ,j ′ + t-1 t ′ =t-L h +1 w t ′ ,j ′ + 1 -x t,j ′ (16)
Note that for any h ∈ {1, ..., H} and k ∈ {1, ..., n h -1}, Q

1 k,h = Q 1 k,h = X MUCP .
As condition (C) is satisfied by S MUCP , any x = (x 1 , ..., x H ) satisfying inequalities ( 15) and ( 16) is such that x h is in the T × n h full orbitope, h ∈ {1, ..., H}. Hence inequalities ( 15) and ( 16) ensure in particular that any solution x h is in the full orbitope.

Sub-symmetry-breaking inequalities for the ramp-constrained MUCP

In the real-world UCP, each unit j must also feature ramp-up (resp. ramp-down) constraints, i.e., the maximum increase (resp. decrease) in generated power from time period t to time period t + 1 is RU j (resp. RD j ). Moreover, if unit i starts up at time t (resp. shuts down at time t + 1), its production at time t cannot be higher than SU j (resp. SD j ).

For each unit j ∈ N and time period t ∈ {2, ...T }, ramp constraints can be formulated as follows:

p t,j -p t-1,j ≤ RU j x t-1,j + SU j u t,j (17)
p t-1,j -p t,j ≤ RD j x t,j + SD j w t,j (18)
The MUCP formulation including ramp constraints can be further strengthened with valid inequalities as proposed in [START_REF] Ostrowski | Tight mixed integer linear programming formulations for the unit commitment problem[END_REF][START_REF] Pan | A polyhedral study of the integrated minimumup/-down time and ramping polytope[END_REF]. As the aim of this article is to compare symmetry-breaking techniques, we will only consider the classical MUCP formulation ( 9) -( 14) with ramp-constraints ( 17) -( 18), as done in [START_REF] Ostrowski | Modified orbital branching for structured symmetry with an application to unit commitment[END_REF][START_REF] Knueven | Generating cuts from the ramping polytope for the unit commitment problem[END_REF]. When ramp-constraints are considered, the symmetry group of set Q t k,h still contains the subsymmetric group associated to the submatrix defined by rows and columns ({t, ..., T }, {j h k , j h k+1 }). Therefore, inequalities [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF] can still be used.

However the symmetry group of set Q t k,h no longer contains the sub-symmetric group associated to the submatrix defined by rows and columns ({t, ..., T }, {j h k , j h k+1 }). Indeed, if two identical units have been up for at least L h time periods at time t -1, they may produce distinct power values at time t -1 and thus, because of ramp constraints, their up/down trajectories from time t to T cannot be permuted. Therefore, inequalities ( 16) can no longer be used.

Note that when two identical ramp-constrained units are ready to shut down, there still exist some sub-symmetries that could be exploited. These sub-symmetries are more intricate because they depend, for example, on the quantity of power produced by both units, or on the time of their last start-up.

Experimental results

In this section, we compare various formulations for the MUCP with or without ramp constraints. Some symmetry-breaking techniques need to interfere with the branching process. These are typically implemented using a callback instruction which deeply affects the performance of commercial solvers like Cplex. Consequently in our computational comparison, we only consider symmetrybreaking techniques that do not require the use of a callback.

Aggregated formulations for the UCP

In [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF], the authors propose to break symmetries of the UCP by aggregating variables corresponding to identical units. This method is shown to outperform existing symmetry-breaking inequalities.

• Aggregated (x, u) formulation In the case of the MUCP, variables x, u of formulation [START_REF] Hojny | Polytopes associated with symmetry handling[END_REF][START_REF] Jans | Solving lot-sizing problems on parallel identical machines using symmetry-breaking constraints[END_REF][START_REF] Kaibel | Branched polyhedral systems[END_REF][START_REF] Kaibel | Orbitopal fixing[END_REF][START_REF] Kaibel | Packing and partitioning orbitopes[END_REF][START_REF] Knueven | Generating cuts from the ramping polytope for the unit commitment problem[END_REF] are aggregated into variables x t,h = j∈N h x t,j ∈ {0, ..., n h } (resp. u t,h = j∈N h u t,j ∈ {0, ..., n h }) indicating how many units of type h are up (resp. start up) at time t. Variables p t,h = j∈N h p t,j ∈ R is the total amount of power produced at time t by units of type h. Formulation (9-14) becomes

A-( x, u) min x, u, p H h=1 T t=1 c h f x t,h + c j p p t,h + c j 0 u t,h s. t. t t ′ =t-L h +1 u t ′ ,h ≤ x t,h ∀h ∈ N h , ∀t ∈ {L h , ..., T } (19) t t ′ =t-ℓ h +1 u t ′ ,h ≤ n h -x t-ℓ h ,h ∀h ∈ N h , ∀t ∈ {ℓ h , ..., T } (20)
u t,h ≥ x t,h -x t-1,h ∀h ∈ N h , ∀t ∈ {2, ..., T } (21)
P h min x t,h ≤ p t,h ≤ P h max x t,h ∀h ∈ N h , ∀t ∈ T (22) H h=1 p t,h ≥ D t ∀t ∈ T (23) x t,h , u t,h ∈ {0, ..., n h } ∀h ∈ N h , ∀t ∈ T (24)
When aggregating variables corresponding to h identical units, one must ensure that the aggregated production plan, satisfing ( 19)- [START_REF] Ostrowski | Modified orbital branching for structured symmetry with an application to unit commitment[END_REF], can be disaggregated into h feasible production plans satisfying ( 9)- [START_REF] Kaibel | Orbitopal fixing[END_REF]. Inequalities ( 9)-( 12) have the integer decomposition property [START_REF] Baum | Integer rounding and polyhedral decomposition for totally unimodular systems[END_REF], i.e., any integer solution ( x, u, p) of formulation ( 19)-( 24) can be disaggregated into an integer solution (x, u, p) of formulation ( 9)- [START_REF] Knueven | Generating cuts from the ramping polytope for the unit commitment problem[END_REF]. A disaggregation algorithm for the MUCP is proposed in [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF].

When ramp constraints are considered in formulation ( 9)-( 14), the integer decomposition property is lost. Examples of aggregated solutions which cannot be disaggregated are given in [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF].

• Aggregated interval formulation As the integer decomposition property depends on the formulation considered, an interval-based formulation is introduced in [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF] for the ramp-constrained MUCP. For each unit j ∈ N , for each interval {t 0 , ..., t 1 -1} of size t 1t 0 ≥ L j , variable y t0,t1 j = 1 if and only if unit j starts up at time t 0 , remains up on interval {t 0 , ..., t 1 -1} and shuts down at time t 1 . For each time period t ∈ T , variable p t0,t1 t,j represents the quantity of power produced by unit j at time t if y t0,t1 j = 1, and p t0,t1 t,j = 0 otherwise. The formulation is as follows. min y,p n j=1 {t0,...,t1-1}∈Yj

c j t0,t1 y t0,t1 j + c j p t1-1 t=t0 p t0,t1 t,j s. t. A t0,t1 j p t0,t1 j ≤ b t0,t1 j y t0,t1 j ∀j ∈ N , ∀{t 0 , ..., t 1 -1} ∈ Y j (25)
{t 0 , ..., t 1 -1} ∈ Y j s.t. t ∈ {t 0 , ..., t 1 + ℓ j } y t0,t1 j ≤ 1 ∀j ∈ N , ∀t ∈ T (26)
j∈N {t0,...,t1-1}∈Yj p t0,t1 t,j

≥ D t ∀t ∈ T (27)
y t0,t1 j ∈ {0, 1} ∀j ∈ N , ∀{t 0 , ..., t 1 -1} ∈ Y j (28)
where Y j = {{t 0 , ..., t 1 -1} ∈ T × T | t 1t 0 ≥ L j } and where P t0,t1 j = {p t0,t1 j ∈ R T + | A t0,t1 j p t0,t1 j ≤ b t0,t1 j } is the feasible production polytope of unit j, if unit j starts up at time t 0 , remains up on interval {t 0 , ..., t 1 -1} and shuts down at time t 1 . In the case of the ramp-constrained MUCP, P t0,t1 j is defined by power limits (29) and ramp constraints (30)-(33):

P j min ≤ p t0,t1 t,j ≤ P j max (29) p t0,t1 t,j -p t0,t1 t-1,j ≤ RU j ∀t ∈ {t 0 + 1, ..., t 1 -1} (30)
p t0,t1 t-1,j -p t0,t1 t,j ≤ RD j ∀t ∈ {t 0 + 1, ..., t 1 -1} (31)
p t0,t1 t0,j ≤ SU j (32)

p t0,t1 t1-1,j ≤ SD j (33)
Inequalities ( 25)-( 26) have the integer decomposition property, thus variables y t0,t1 j (resp. p t0,t1 t,j ) can be aggregated into variables y t0,t1 h = j∈N h y t0,t1 j and p t0,t1 t,h = j ∈ N h p t0,t1 t,j , leading to aggregated formulation Int( y):

Int( y) min y, p H h=1 {t0,...,t1-1}∈Y h c h t0,t1 y t0,t1 h + c h p t1-1 t=t0 p t0,t1 t,h s. t. A t0,t1 h p t0,t1 h ≤ b t0,t1 h y t0,t1 h ∀h ∈ {1, ..., H}, ∀{t 0 , ..., t 1 -1} ∈ Y h (34) {t 0 , ..., t 1 -1} ∈ Y j s.t. t ∈ {t 0 , ..., t 1 + ℓ j } y t0,t1 j ≤ n h ∀h ∈ {1, ..., H}, ∀t ∈ T (35) H h=1 {t0,...,t1-1}∈Y h p t0,t1 t,h ≥ D t ∀t ∈ T (36)
y t0,t1 h ∈ {0, ..., n h } ∀h ∈ {1, ..., H}, ∀{t 0 , ..., t 1 -1} ∈ Y h (37)

Experimental settings

In this section, we compare various symmetry-breaking formulations for the MUCP with or without ramp-constraints. As shown in [START_REF] Ostrowski | Modified orbital branching for structured symmetry with an application to unit commitment[END_REF], neither Friedman inequalities (1) nor column inequalities (2) are competitive with respect to the classical UCP formulation when solved by Cplex.

On the opposite, the weaker form of Friedman inequality (3) has been shown in [START_REF] Lima | Symmetry breaking in MILP formulations for Unit Commitment problems[END_REF] to outperform Default Cplex.

Hence the following formulations for the MUCP are compared: -F (x, u): (x, u)-formulation (9)-( 14) -A-( x, u):

Aggregated ( x, u)-formulation (19)-( 24) (only when disaggregation applies) -Int( y):

Aggregated interval formulation (34)-(37) -W (x, u):

(x, u)-formulation (9)-( 14) with weaker Friedman inequalities (3) -F (x, u, z): (x, u)-formulation (9)-( 14) with variables z, inequalities (7a)-(7e) and subsymmetry-breaking inequalities (8) -LF (x, u):

(x, u)-formulation (9)-( 14) with sub-symmetry-breaking inequalities ( 15)-( 16). Formulation F (x, u, z) is obtained from the classical MUCP formulation F (x, u) by a direct use of the inequalities given in Section 3. As seen in Section 4, taking into account sub-symmetries in the MUCP leads to formulation LF (x, u) featuring lifted symmetry breaking-inequalities ( 15) and ( 16), namely Start-up-ready and Shut-down-ready inequalities, in place of inequalities (7a)- [START_REF] Friedman | Fundamental Domains for Integer Programs with Symmetries[END_REF]. Note that the start-up and shut-down sub-symmetries of the MUCP are not handled by formulations F (x, u), W (x, u) and F (x, u, z).

Formulations F (x, u), W (x, u), F (x, u, z) and LF (x, u) feature O(nT ) variables while formulation A-( x, u) (resp. Int( y)) features O(HT ) (resp. O(T 2 H)) variables, where H is the number of groups of identical units.

For the ramp-constrained MUCP, inequalities ( 17)-( 18) enforcing ramp-constraints are added to formulations F (x, u), W (x, u), F (x, u, y) and LF (x, u). Aggregated formulation A-( x, u) can no longer be used, as its solutions cannot be disaggregated [START_REF] Knueven | Exploiting identical generators in unit commitment[END_REF]. Note also that in this context, Start-up-ready inequalities are adjoined to LF (x, u), but Shut-down-ready inequalities cannot.

In formulation Int( y), the production limit constraint (29) is always included in inequalities (34) defining feasible productions. In the ramp-constrained case, the ramp constraints (30)-(33) are also included.

All experiments are performed using Cplex 12.8 C++ API on 8 threads of a PC with a 64 bit Intel Core i7-6700 processor running at 3.4GHz, and 32 GB of RAM memory. The UCP instances are solved until optimality (defined within 10 -7 of relative optimality tolerance) or until the time limit of 3600 seconds is reached.

Instances

We generate MUCP instances as follows.

For each instance, we generate a "2-peak per day" type demand with a large variation between peak and off-peak values: during one day, the typical demand in energy has two peak periods, one in the morning and one in the evening. The amplitudes between peak and off-peak periods have similar characteristics to those in the dataset from [START_REF] Carrion | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF].

We consider the parameters (P min , P max , L, ℓ, c f , c 0 , c p ) of each unit from the dataset presented in [START_REF] Carrion | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF]. We draw a correlation matrix between these characteristics and define a possible range for each characteristic. In order to introduce symmetries in our instances, some units are randomly generated based on the parameters correlations and ranges. Each unit generated is duplicated d times, where d is randomly selected in [1, n F ] in order to obtain a total of n units. The parameter F is called symmetry factor, and can vary from 2 to 4 depending on the value of n. Note that these instances are generated along the same lines as literature instances considered in [START_REF] Bendotti | The min-up/min-down unit commitment polytope[END_REF], but with different F factors.

In order to determine which symmetry-breaking technique performs best with respect to the number of rows and columns of matrices in feasible set X , we consider various instance sizes n ∈ {20, 30, 60} and T ∈ {48, 60}, and various symmetry factors F ∈ {2, 3, 4}. For each size (n, T ) and symmetry factor F , we generate a set of 20 instances. Symmetry factor F = 4 is not considered for instances with a small number n of units (n = 20 or 30), as it leads to very small sets of identical units.

Table 1 provides some statistics on the instances characteristics. For each instance, a group is a set of two or more units with the same characteristics. Each unit which has not been duplicated is a singleton. The first and second entries column-wise are the number of singletons and groups. The third entry is the average group size and the fourth entry is the maximum group size. Each entry row-wise corresponds to the average value obtained over 20 instances with same size (n, T ) and same symmetry factor F .

The ramp-constrained MUCP instances considered are the same as in the non-ramp-constrained case, with additional ramp characteristics RU j = P j max -P j min 3

, RD j = P j max -P j min 2 and SU j = SD j = P j min . Note that a sign "-" in the column entry corresponding to the CPU time means that no instance could be solved within the time limit.

Results for the MUCP

It is clear from Table 2 that aggregated (x, u) formulation A-( x, u) outperforms by far all the other formulations. This could be explained by the reduced size of aggregated formulation A-( x, u), but also by the good performance of Cplex on ILP featuring integer variables (with bounds greater than 1). This efficiency will certainly be preserved any time the integer decomposition property holds for an (x, u) formulation of the UCP. Aggregated interval formulation Int( y) is in average one or even two order of magnitude slower than F (x, u), F (x, u, z) and LF (x, u). Formulations F (x, u, z) and W (x, u) are always outperformed by F (x, u) and LF (x, u). Formulations F (x, u) and LF (x, u) are quite comparable on (n, T ) = (20, 48) instances. Interestingly, on (n, T ) = (20, 96) instances, LF (x, u) is better than F (x, u). Otherwise, when n is larger (i.e., n ≥ 30), and when T = 96 or when F = 2, F (x, u) outperforms LF (x, u), probably because more sub-symmetry-breaking inequalities are added in this case. When the horizon size is smaller (i.e., T=48) and when F ∈ {3, 4}, fewer sub-symmetry-breaking inequalities are needed, which allows formulation LF (x, u) to outperform F (x, u) for n = 30 and n = 60.

Results for the ramp-constrained MUCP

Recall that aggregated formulation A-( x, u) can no longer be used in this context.

Tables 4 and5 provide, for each formulation and each group of 20 instances, the exact same column entries as those in Tables 2 and3.

First note that the ramp constraints make the MUCP instances much harder to solve by Cplex in general, as the CPU times in Tables 4 and5 are much larger than those in Tables 2 and3. For example, the integrality gap is in average 10 times larger for ramp-constrained problems on (n, T ) = (60, 48) and F = 2 instances. instance with n > 20 but one. Moreover, on n = 30 instances, and on (n, T ) = (60, 96) instances, the root node cannot be processed at all within the time limit for formulation Int( y); the number of nodes explored is 0 and the optimality gap is 100%. Formulation LF (x, u) is more efficient than all considered formulations. Except on (n, T ) = (20, 96) and F = 2 instances where F (x, u, z) manages to solve to optimality two instances more than LF (x, u), formulation LF (x, u) is able to solve a larger number of instances to optimality than all considered formulations. In particular, it manages to solve to optimality two of the large-size instances (i.e., (n, T ) = (60, 96)), while other formulations do not reach optimality on any of these instances. Moreover, formulation LF (x, u) solves 52 instances to optimality among the 80 instances with n = 30, while F (x, u, z) or F (x, u) (resp. W (x, u)) only manages to solve to optimality 18 (resp. 24) of them. Among the 80 instances with n = 20, formulation LF (x, u) solves 57 instances to optimality, while F (x, u, z) (resp. W (x, u), F (x, u)) only manage to solve to optimality 49 (resp. 42, 34) of them. Formulation Int( y) solves to optimality only 20 of them. Formulation LF (x, u) also globally improves the solving time. For example, on instances of size (n, T ) = (60, 48) and F = 3, the average CPU time of formulation LF (x, u) is 1450 seconds, while this number increases to 3422 (resp. 2527, 2689) for F (x, u) (resp. F (x, u, z), W (x, u)).

As there is an important variability in the computation time for instances with same size (n, T ) and same F , we introduce the improvement score. For given formulations F 1 and F 2 , the improvement score I of F 1 with respect to F 2 is as follows.

I = 2 CP U (F 2 ) -CP U (F 1 ) CP U (F 2 ) + CP U (F 1 )
The improvement score I is a performance ratio comparing formulation CPU times pairwise. Table 6 provides, for each formulation F ∈ { F (x, u, z), LF (x, u) }, the average improvement score of F with respect to F (x, u) on each group of 20 instances. Formulation Int( y) is not included as on most instance groups, it solves no instance to optimality.

Formulation LF (x, u) outperforms all other formulations. In particular, even if on (20, 96) and F = 2 instances the average CPU time of LF (x, u) is slightly higher than F (x, u, z), the average improvement score of LF (x, u) is more important. This shows that LF (x, u) have higher CPU time than F (x, u, z) on instances on which the difference in CPU time is not very significant with respect to Cplex's CPU time. On the opposite, LF (x, u) have low CPU time on instances on which this difference in CPU time represents an important improvement. Note that formulation W (x, u) appears to perform better than F (x, u, z) on T = 48 instances. Recall that W (x, u) is only partialsymmetry-breaking. Thus, when T is smaller, the number of feasible columns featuring a given number of 1-entries is also smaller. On the opposite, when T = 96, the number of one-entries is not a very discriminating indicator among symmetric columns. Therefore W (x, u) is not able to break as much symmetries, and F (x, u, z) globally performs better.

For example, on (n, T ) = (60, 48), F = 4 instances, the improvement score of LF (x, u) is 109%, while it is 45.8% for F (x, u, z) and 27.9% for W (x, u). On (n, T ) = (30, 48), F = 3 instances, the improvement score of LF (x, u) is 114%, while it is -19.3% for F (x, u, z) and 25.3% for W (x, u). On (n, T ) = (20, 96), F = 3 instances, this number increases to 47.6% for LF (x, u) (resp. 22.9% for F (x, u, z), -11.9% for W (x, u)).

Conclusion

We propose a framework to build sub-symmetry-breaking inequalities, in order to handle the symmetries arising from a collection of sub-symmetric solution subsets. These inequalities may require to introduce one additional variable z per solution subset considered. Depending on the subset structure, variable z may only be a linear expression of variables x, and therefore does not need to be added to the model as an additional variable. The derived sub-symmetry-breaking inequalities are full symmetry-breaking under a mild condition. If this condition is not satisfied anyway, a new collection of sub-symmetric subsets can be constructed such that the derived inequalities are full symmetry-breaking.

Our experimental results for the MUCP show that aggregation of the classical formulation is a very efficient technique to handle symmetries and sub-symmetries arising in the MUCP. When ramp constraints are taken into account in the MUCP, disaggregation is no longer possible. Our sub-symmetry-breaking inequalities can still be used and outperform all other formulations. Sub-symmetry-breaking inequalities are always applicable as the solution subsets considered can capture the specific conditions under which the symmetries hold. On the opposite, aggregated formulations require specific conditions to be applicable.

One perspective is to use the framework proposed in this article to derive new sub-symmetrybreaking inequalities for "ready to shut down" sub-symmetries in the ramp-constrained case. Another perspective is to apply the proposed framework to other problems featuring sub-symmetric solution subsets such as covering problems, or bin packing variants where one item can be placed in multiple bins. It would also be useful to study how the presented framework could be automated, so that sub-symmetric subsets are automatically detected and variables z automatically constructed.

  q}, k ∈ {1, ..., |R s |}, l ∈ {2, ..., |C s |} where for each s ∈ {1, ..., q}, for each l ∈ {2, ..., |C s |}, for each k ∈ {1, ..., |R s |},

Table 1 :

 1 Instance characteristics

	Tables 2 and 3 provide, for each formulation and each group of 20 instances:
	#opt:	Number of instances solved to optimality,
	#nodes:	Average number of nodes,
	gap:	Average optimality gap,
	CPU time:	Average CPU time in seconds.

Table 2 :

 2 Formulation Int( y) is still the less efficient formulation. It does not solve to optimality any

	Instances	Formulation #opt	#nodes gap (%)	CPU time
	(20,48) F = 2	F (x, u)	20	1271 0	2.6
		A-( x, u)	20	0 0	0.2
		Int( y)	16	205 667 0.005 02	781.6
		W (x, u)	20	4809 0	13.7
		F (x, u, z)	20	3838 0	23.2
		LF (x, u)	20	1915 0	6.6
	F = 3	F (x, u)	20	806 0	2.6
		A-( x, u)	20	0 0	0.3
		Int( y)	18	152 948 0.001 57	572.1
		W (x, u)	20	1600 0	4.4
		F (x, u, z)	20	683 0	6.7
		LF (x, u)	20	271 0	3.5
	(20,96) F = 2	F (x, u)	20	148 942 0	267.3
		A-( x, u)	20	0 0	0.7
		Int( y)	6	13 180 9.368 57	2977.5
		W (x, u)	18	110 644 0.000 15	459.6
		F (x, u, z)	18	118 877 0.000 18	497.8
		LF (x, u)	19	52 881 0.000 13	215.1
	F = 3	F (x, u)	18	29 418 0.032 71	376.7
		A-( x, u)	20	2360 0	8.2
		Int( y)	7	32 859 0.168 55	2574.8
		W (x, u)	19	79 864 0.000 61	357.1
		F (x, u, z)	18	39 694 0.000 92	458.3
		LF (x, u)	19	19 831 0.000 85	229.1
	(30,48) F = 2	F (x, u)	20	3207 0	5.5
		A-( x, u)	20	0 0	0.2
		Int( y)	13	615 979 0.002 11	1483.4
		W (x, u)	18	67 553 0.000 18	370.7
		F (x, u, z)	18	41 669 0.000 21	391.5
		LF (x, u)	20	3894 0	19.4
	F = 3	F (x, u)	19	157 903 0.000 01	189.1
		A-( x, u)	20	6 0	0.4
		Int( y)	13	267 844 0.007 72	1548.7
		W (x, u)	19	108 689 0.000 02	283.9
		F (x, u, z)	19	30 450 0	249.5
		LF (x, u)	20	6044 0	18.5
	(30,96) F = 2	F (x, u)	20	39 633 0	131.2
		A-( x, u)	20	21 0	0.5
		Int( y)	0	21 783 2.002 12	-
		W (x, u)	18	54 860 0.000 29	679.8
		F (x, u, z)	18	20 062 0.000 83	758.5
		LF (x, u)	19	17 040 0	299.5
	F = 3	F (x, u)	20	15 653 0	53.5
		A-( x, u)	20	215 0	0.6
		Int( y)	0	30 267 4.389 68	-
		W (x, u)	13	324 913 0.000 09	1400.7
		F (x, u, z)	17	123 943 0.000 03	907.1
		LF (x, u)	19	146 639 0.000 03	299.6

Performance indicators relative to the comparison of six formulations for MUCP instances with symmetries and n = 20, 30

Table 6 :

 6 Improvement score w.r.t. F (x, u) W (x, u) F (x, u, z) LF (x, u) Improvement scores of formulations W (x, u), F (x, u, z) and LF (x, u) w.r.t formulation F (x, u) for ramp-constrained MUCP instances with symmetries

	(20,48) F = 2	-6.61%	-14.2%	83.7%
	F = 3	-11.3%	-12.2%	111%
	(20,96) F = 2	-13%	9.05%	23.9%
	F = 3	-11.3%	22.9%	47.6%
	(30,48) F = 2	20.8%	-17.8%	89.4%
	F = 3	25.3%	-19.3%	114%
	(30,96) F = 2	-15.8%	-19.4%	40.4%
	F = 3	-11.8%	-7.63%	76.4%
	(60,48) F = 2	26.5%	26.2%	47.8%
	F = 3	35.5%	30.9%	104%
	F = 4	27.9%	45.8%	109%
	(60,96) F = 2	0%	0%	0%
	F = 3	0%	0%	0%
	F = 4	0%	0%	4.72%


Omar A Mures

Emilio J Padrón

Bruno Raffin

Leveraging the Power of Big Data Tools for Large Scale Molecular Dynamics Analysis

Keywords: Molecular Dynamics, Big Data, MapReduce, Flink, High Performance Data Analysis, Parallelization

Parallel Molecular Dynamics simulations are generating atom trajectories of growing sizes and complexity. Analyzing these trajectories is expensive computationally and time consuming. One reason is the lack of tools that enable the computational biologist to easily implement the analysis while ensuring reduced processing times exploiting the benefits of parallel architectures. In this paper, we present a comparison between two parallel analytics frameworks based on the Map/Reduce paradigm: HiMach, a dedicated framework for trajectory analysis based on MPI, and Flink, a Big Data analytics framework. Both frameworks enable to hide the complexity of parallel code creation to the programmer, providing significant performance gains compared to a sequential execution.

I. Introduction

R ECENTLY in the HPC vendor and science com- munity a new tendency is on the rise, since both fields face similar issues, these two worlds need to come together in order to solve bigger and more complex science problems. As a consequence, base technologies themselves are becoming more closely related. Nowadays scientists are taking advantage more than ever of the evolving computing resources available to them. From their local workstations to high performance clusters in large data centers, they rely on these tools to discover new phenomena in their respective fields.

To accomplish the aforementioned tasks, there is an inherent need to simplify the creation of analysis code that can exploit these heterogeneous resources. Big Data analytics tools based on the MapReduce paradigm [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF], such as Hadoop [START_REF]Hadoop[END_REF], enable to take advantage of commodity compute clusters for analysing very large data sets generated from the web or business applications. But when trying to analyse complex scientific datasets they can fall short from a programmers perspective and performance wise.

Molecular Dynamics (MD) is a computer simulation technique for studying physical interactions of atoms and molecules [START_REF] Berni | Studies in molecular dynamics. i. general method[END_REF]. It is a type of N-body simulation [START_REF] John | The complexity of n-body simulation[END_REF], which simulates a dynamic system of particles, typically under the influence of physical forces, such as gravity. Atoms and molecules interact during a preset period of time, allowing the scientist to observe the dynamical evolution of the system.

The most common approach determines the trajectories of atoms and molecules by solving Newton's equations of motion for a system of interacting particles. The forces and energies of said particles are computed using inter-atomic potentials or molecular mechanics force fields. This method is applied today in chemical physics, material science and the modeling of biomolecules.

Although plenty of progress has been made in the creation of efficient parallel code for Molecular Dynamics simulations, code for analysis of the results that they produce in parallel is scarce. As a consequence, scientist have to rely on serial tools to analyze the aforementioned data. In this paper we focus on the analysis of atom trajectories generated from Molecular Dynamics simulations. These datasets can be very large and time consuming to analyze, sometimes reaching thousands of gigabytes in size. The size of trajectories is expected to keep on growing further as exascale computers [START_REF] Germann | Exascale computing and what it means for shock physics[END_REF] become available, aggravating this need for parallel analysis code. Current serial approaches force the scientist to filter trajectories sometimes randomly in order to find the phenomena of interest, this is far from a perfect approach, and if not resolved, could hurt the scientists' opportunities to perform new discoveries and rendering the purpose of creating more precise simulations useless.

Existing analysis libraries such as MDAnalysis [START_REF] Michaud-Agrawal | Mdanalysis: A toolkit for the analysis of molecular dynamics simulations[END_REF] or VMD [START_REF] Humphrey | Vmd: visual molecular dynamics[END_REF] rely on traditional parallelization approaches based on MPI and/or OpenMP that are often complex to deploy at large scale and use efficiently for computational biologists. In this paper, we investigate how the Map/Reduce paradigm could be leveraged to enable the analysis of large trajectories. We have performed a comparison between two parallel analytics frameworks based on the Map/Reduce paradigm: HiMach, a framework for trajectory analysis based on MPI and VMD, and Flink, a new generation Big Data analytics engine. Both frameworks hide the complexity of parallel code creation to the user, providing significant performance gains compared to sequential executions.

A brief background overview of the subject is given in section II, the proposed solution is explained in section III and the results obtained are presented in section IV.

II. Background

MapReduce is a programming model for processing and generating large data sets in parallel, on commodity clusters. It relies on a simple programming model that uses two main functions, Map and Reduce. They are inspired by the map and reduce functions that were commonly used in functional programming, although their purpose in the MapReduce framework is not the same as in their original forms.

These two functions (see Figure 1) have to be defined by the programmer in order to obtain parallelism for data-intensive tasks. The Map function consumes key-value pairs and produces one or more intermediate pairs:

M ap(k1, v1) → list(k2, v2)
The Reduce function, processes a certain key and a list of values associated to it, usually outputting a list with less values than the input list:

Reduce(k2, list(v2)) → list(v3)
For example, a typical MapReduce workload that counts words in a text, will Map (extract) the words in all lines of the text emitting (k2, v2), where k2 is the word and v2 is 1. This will allow to Reduce (summarize) the results, counting the number of occurrences of each word. In order to do this, each reducer uses the input (k2, list(v2)); where k2 is the word and list(v2) is the list of values associated to k2 (all ones), to sum up the values, obtaining the aforementioned count.

The MapReduce system performs the processing by marshaling the distributed servers, running the required tasks in parallel, managing all communications and data movement between the various parts of the system, while providing redundancy and fault tolerance. As such, a in a single-threaded environment, MapReduce will usually not be faster than a traditional implementation; any performance gains are usually only seen in multi-threaded scenarios.

From a performance standpoint, MapReduce has several limitations. First, as we have seen in Figure 1, after each Map phase, results are written to disk as Map Output Files (MOFs) even if they can fit in main memory, and will be later read by the Reduce phase in order to compute the final results. This is good from a fault tolerance point of view, but when trying to achieve high performance it presents a clear problem, as avoiding disk writes and reads could increase performance substantially. Furthermore, until every map task has finished, the reduce phase cannot start; even if the required data has already been obtained. Second, since the paradigm is quite simple, it presents almost no opportunities for the framework to perform automatic optimizations. While a experimented MapReduce programmer can use these two abstractions to model complicated algorithms, it is clear that this paradigm can fall short in some complex scenarios where more flexibility is required, apart from the fact that there is not enough semantic information in this two constructs in order to infer properties about the user defined functions.

Several efforts have been made to mitigate these problems, for instance the parallel trajectory analysis framework HiMach [START_REF] Tu | A scalable parallel framework for analyzing terascale molecular dynamics simulation trajectories[END_REF], is a MapReduce like system where the user writes sequential trajectory analysis code and the system carries out the parallel analysis of the trajectory automatically. This solution avoids the use of MOFs, as well as, using high performance MPI code that can take advantage of In-finiBand interconnects, improving performance substantially. Although this framework solves some of the issues that have been presented earlier, the parallelization paradigm offered still suffers from the same issues as MapReduce.

Since for complex analysis the MapReduce method is not optimal, a more advanced paradigm that relies on Directed Acyclic Graphs (DAGs) [START_REF] Thulasiraman | Graphs: Theory and Algorithms[END_REF] has been chosen in order to allow the user to express more complex analytics tasks. This new approach uses Parallelization Contracts (PACTs) [START_REF] Alexandrov | Mapreduce and pact-comparing data parallel programming models[END_REF], which can be thought of as a generalization of MapReduce in order to ease the expression of many operations, which in turn leads to better expressiveness and enables optimizations to be inferred automatically. The chosen framework is the result of the Stratosphere project [START_REF] Alexandrov | The stratosphere platform for big data analytics[END_REF], an open source project from the Apache Foundation called Flink [12]. To our Knowledge Stratosphere/Flink has never been used for analyzing trajectories.

The usage of such framework will allow computational biologists to more efficiently analyze massive trajectories, that until now had to be simplified, which is not ideal. With the presented framework, a scientist will be able to leverage in a simple way the power of new generation MapReduce tools, letting them worry about what analysis they want to perform instead of how to parallelize or deal with simplifying the data so it is usable.

III. MDReduce Framework

Since our main objective is easing the creation of parallel code by the scientist, initially we have chosen the Python [START_REF]Python[END_REF] language for its ease of use and plethora of external libraries. Python has been amply utilized by the scientific community, with libraries for efficient numerical analysis like Numpy [START_REF] Van Der Walt | The numpy array: a structure for efficient numerical computation[END_REF]. This library is a key component of MDAnalysis, one of the tools that allows scientists to efficiently create parallel Molecular Dynamics analysis code using our system. This are the key technologies that were used in order to implement the MDReduce framework in conjunction with Flink and HiMach.

A. Architecture

Once the reader has been acquainted with the general concepts that have been used in the implementation of MDReduce, it is interesting to take a look at Figure 2, that presents a more accurate picture of the system. In the current MDReduce implementation, the Core module contains all the classes related to the system internals. First, the Configuration class takes care of parsing either a configuration file (it uses the YAML syntax to store several settings) or the arguments that the user passes to the class constructor. It will take care of initializing the requested back-end (Flink or HiMach), and other trivial setup tasks. Directly related to this class are the available Backends. These classes are in charge of dealing with Flink and HiMach internals, these are helper classes for initialing these frameworks and running the requested jobs on them.

API

Fig. 2 MDReduce system architecture.

The main supported operations that the user can employ to analyze trajectories in parallel are the following:

In addition, the Analysis module contains all the classes related to the Molecular Dynamics analysis methods included in the MDReduce framework. First, an Algorithm base class has been created, so that developing new analysis scripts is more intuitive for programmers. Second, several analysis classes have been included to perform the analyses that will be detailed in subsection III-B. These two classes, WaterCount and DihedralAngles offer an insight for programming analysis codes using the MDReduce platform.

Finally, through the MDReduce API, external programs that use the existing algorithms or perform new types of analyses that can be created by scientists and run on the platform. In order to ease the integration of the framework with external programs, output can be consumed by and external client, or stored in data files. These files can then be analyzed in local workstations, since the amount data should have been substantially reduced, and should now be usable on lesser machines. As can be seen in Figure 2, Zeppelin [START_REF]Zeppelin[END_REF] a new framework that enables interactive data analytics has also been integrated with the system. This allows scientists to submit jobs interactively, and dynamically explore the resulting data. This is a convenient tool to ease the creation and representation of different analyses, since it allows the user to interactively modify code and visualize analysis results.

B. Case Study

To test the feasibility of the proposed system for Molecular Dynamics analyses, we choose two typical tasks for a computational biologist. First, how to obtain a water count in a GLIC channel [START_REF] Aryal | Hydrophobic gating in ion channels[END_REF], [START_REF] Corringer | Atomic structure and dynamics of pentameric ligand-gated ion channels: new insight from bacterial homologues[END_REF] using MDReduce. Second, another common analysis performed in trajectories consists in computing the dihedral angles [START_REF] Ricarda | A prokaryotic perspective on pentameric ligand-gated ion channel structure[END_REF] in the protein backbone. As seen in Figure 3, the objective of this analysis is obtaining the amount of water ions that are comprised between certain residues in the ion channel. The final objective is obtaining a water ion count for every frame in the trajectory, since the user will vary the simulation parameters in order to isolate the causes of the closing and opening of the GLIC ion channel. This structure has both open and closed states, the transition mechanism between different states is unknown for now. Understanding this process is crucial, as empowering molecules such as anesthetics and alcohols are believed to affect transition barriers or the relative free energy of states.

To investigate these phenomena, extensive simulations of the GLIC channel with different starting conditions and parameters have been performed. This allows computational biologists to draw conclusions about the pore hydration, organization and the pore lining helix. To characterize the structure of this helix, the first part is performing this analysis. The number of water molecules in the pore zone between the two hydrophobic residues I9 and I16 are counted. After this the scientist can infer which region has been more affected by pore dewetting.

This type of analysis can be easily characterized in MDReduce using either the Flink or HiMach backends. Since this is an operation that needs to be computed for every frame, this translates easily to the MapReduce paradigm. One just needs to write a mapper that will perform the water count operation for every frame. The Figure 4 represents how this analysis can be performed using the Flink back-end.

The user has to override Flink's Map function, with the operations that will be applied to the trajectory. First, the frame number corresponding to the instance of mapper is selected, then it is processed with MDAnalysis in order to align it to a reference frame and obtain the water count in the ion channel using atom selections. As can be seen, this code is completely straightforward requiring slight modifications to be run in parallel. Water count code example in Flink.

B.2 Dihedral Angles

α ઠ β α α-carbon Φ Ψ I9' I16' α-carbon Φ Ψ N N N N Fig. 5
Graphical representation of the dihedral angle calculation operation in the protein.

As seen in Figure 5, the objective of this analysis is obtaining a way to visualize energetically allowed regions for the backbone dihedral angles of amino acid residues in the protein structure. As mentioned earlier, understanding the reasons why the transition between the open and closed states of the ion channel is crucial for several topics. This second analysis code is the missing link to study this phenomena that depends on the pore organization and lining of the M2 helix. TO help in the description of the state of the pore, the tilt angles are computed and split in two angles: the radial and lateral tilt angles of the M2 helix around the center of mass in respect to the channel axis.

The second analysis, can be performed in MDReduce using either both available back-ends. Since this is an operation that again, needs to be computed for every frame, it is easily adapted to the MapReduce paradigm. The user just needs to write a mapper that will perform the required operations for every frame. The Figure 6 represents how this analysis can be performed using the HiMach backend.

IV. Results

To test the proposed system we have used the Pluton cluster, a supercomputing facility situated in the Universidade da Coruña. In the Table II, the hardware that the cluster possesses can be seen. • Broadcast: This strategy utilizes a single file that contains the trajectory, that is then distributed to every map task and accessed in par-allel.

• Split: The trajectory in this strategy is split in one file per frame, which is then accessed independently by every map task in parallel.

All strategies were then executed using both backends available in MDReduce to see which of them fares better for each type of analysis.

Figure 7 shows the job execution times of the water count analysis. It is kind of difficult to see which approach fares better, since with every approach the execution time as we use more than one processor descends greatly. The fastest time for this workload, was achieved by the HiMach back-end with the broadcast strategy. This is due to the fact that the HiMach runtime uses MPI and the Infiniband capabilities of our cluster in an optimal way. Bandwidth with each back-end of the water count job.

To further see what is happening in this analysis code in Figure 8, we can see the bandwidth achieved by all the strategies. Here the reader can see better which back-end is the fastest. This results mimic the ones obtained before, but paint a clearer picture about the performance that each strategy achieves. As was mentioned earlier, the best performance is achieved by HiMach followed closely by Flink. The first thing that one can deduce examining this graph, is that clearly the split access pattern is less efficient than the broadcast pattern. This was predicted, since reading through small files generally causes lots of seeks and stresses the underlying file system, which results in an inefficient data access pattern.

Fig. 9 Speedup with each back-end of the water count job.

In addition, as to ascertain the efficiency of the tested frameworks, in Figure 9 we can see the speedups of all the different code versions. Here we can observe other interesting behavior. Both frameworks scale well but as the number of processors is increased a couple of discrepancies can be observed, and the scaling weakens. In addition, we can see that with the split strategy both back-ends have close speedups. The performance difference is more significant for the broadcast access pattern, probably due to a higher impact of the network performance that HiMachbetter use than Flink.

V. Future Work

Since several inefficiencies in the current implementation are believed to be related to the usage of the Python interpreter from the Java Flink processes, it could be interesting to test the ZipPy [START_REF] Zhang | Accelerating iterators in optimizing ast interpreters[END_REF] python implementation in our system. This Python 3 implementation uses Truffle [START_REF] Würthinger | One vm to rule them all[END_REF] to generate optimized JVM bytecode from Python scripts. This could certainly improve the efficiency of the analysis performed using the MDReduce framework.

In addition, it could also be interesting to extend our study to include a Spark [START_REF] Zaharia | Fast and interactive analytics over hadoop data with spark[END_REF] back-end, today probably the most popular Dta analytics framework.

Moreover, to solve one of the main issues with the split data access pattern, other storage solutions like HBase [START_REF] Mehul | Hadoop-hbase for large-scale data[END_REF] could be explored. This column oriented data store has the potential to avoid the bottleneck that the plethora of small files in a long trajectory could present to the MDReduce framework.

VI. Conclusions

Using MDReduce has proven to improve not only processing times, since it allows users to write almost sequential code that can then be exploited for parallel analysis, but also in terms of programming efficiency. It is clear that using the introduced system, is minimally more complex than writing sequential analysis code for straightforward tasks.

The usage of Flink, can yield performance advantages due to incorporated job optimizer, but also makes programming more complex analysis algorithms possible. It also allow the user to take advantage of the Hadoop ecosystem, for instance allowing the users to visualize data dynamically using Zeppelin.

Our initial performance tests have shown that Hi-Mach, due to its simplicity, is the best option for running straightforward analysis tasks. However the analysis we have tested are simple map/reduce schemes. They did not allow us to probe the capabilities of Flink designed to support more complexe parallelization schemes.

Fig. 1 Typical

 1 Fig.1Typical MapReduce job workflow.

B. 1 Fig. 3

 13 Fig. 3Graphical representation of the water counting operation in an ion channel.

3 u

 3 Fig. 4

Fig. 6

 6 Fig. 6 Dihedral angle code example in HiMach.

Fig. 7

 7 Fig. 7Run time with each back-end of the water count job. Legend should be move elsewhere (inside the graph or above it) to have bigger graphs

Fig. 8

 8 Fig. 8

TABLE I

 I Supported operations. Note: The HiMach backend only supports the Map and Reduce operations.

	Operation Description
		Takes one element and
	Map	produces zero, one, or more
		elements.
	Reduce	Combines a group of elements into a single element.
		Evaluates a function for each
	Filter	element and retains those
		that fulfill it.
		Aggregates a group of
	Aggregate	values into a single value.
		Minimum, maximum, and sums.
		Joins two data sets by
	Join	creating all pairs of elements that are equal on
		their keys.
	Union	Produces the union of two data sets.
	Sort	Sorts dataset depending on a certain field.

TABLE II

 II Hardware used for testing. The analysis algorithms explained in subsection III-B, have been used to run the performance tests carried out in this section.

	Hardware	compute-0-16
	CPU Model	2 x Intel Xeon E5-2660
	Cores/Threads	16/32
	CPU Speed/Turbo 2.20 GHz/3 GHz
	Memory	64 GB DDR3 1600 Mhz
	Disk	1 x HDD 1 TB SATA3 7.2K
	Networks	InfiniBand FDR
	Add performance number for sequential execution (pure
	MDAnalyis and/or VMD code)

These tests have tested several different strategies, that access data with multiple I/O patterns, in order to discern the best I/O access pattern. The tested patterns are:

Computer Architecture Group, Universidade da Coruña, e-mail: omar.alvarez@udc.es

Computer Architecture Group, Universidade da Coruña, e-mail: emilioj@udc.es

University GrenobleAlpes, INRIA, e-mail: bruno.raffin@inria.fr