Anne-Elisabeth Falq
email: anne-elisabeth.falq@lip6.fr

Pierre Fouilhoux
email: pierre.fouilhoux@lip6.fr

Safia Kedad-Sidhoum
email: safia.kedad_sidhoum@cnam.fr

Dominance inequalities for scheduling around an unrestrictive common due date

Keywords: scheduling, integer programming, common due date, dominance properties

The problem considered in this work consists in scheduling a set of tasks on a single machine, around an unrestrictive common due date to minimize the weighted sum of earliness and tardiness. This problem can be formulated as a compact mixed integer program (MIP). In this article, we focus on neighborhood-based dominance properties, where the neighborhood is associated to insert and swap operations. We derive from these properties a local search procedure providing a very good heuristic solution. The main contribution of this work stands in an exact solving context: we derive constraints eliminating the non locally optimal solutions with respect to the insert and swap operations. We propose linear inequalities translating these constraints to strengthen the MIP compact formulation. These inequalities, called dominance inequalities, are different from standard reinforcement inequalities. We provide a numerical analysis which shows that adding these inequalities significantly reduces the computation time required for solving the scheduling problem using a standard solver.

Introduction

The scheduling problem studied in this work falls within the just-in-time scheduling field. In this framework, each task has a due date, and any deviation from this due date is penalized. From one hand, the tardiness is penalized to model the customer dissatisfaction, on the other hand, the earliness is penalized to model the induced storage costs . The reader can refer to the seminal surveys of [START_REF] Baker | Sequencing with earliness and tardiness penalties: A review[END_REF] and [START_REF] Kramer | A unified heuristic and an annotated bibliography for a large class of earliness-tardiness scheduling problems[END_REF] for the early results in this field.

We consider a set J of n tasks with fixed processing times (p j ) j∈J ∈ R J + , to be non-preemptively processed on a single machine. These tasks share a common due date d for which tasks should be preferably completed. In this work, we assume that the due date is unrestrictive, i.e. d p j . In a just-in-time framework, tasks completing before or after d will therefore incur penalties according to unit earliness (resp. tardiness) penalties (α j ) j∈J ∈ R J + (resp. (β j ) j∈J ∈ R J + ). A schedule is defined by the task completion times denoted by (C j ) j∈J . Using [x] + to denote the positive part of x ∈ R, the earliness (resp. tardiness) of a task j ∈ J is given by [d -C j ] + (resp. [C j -d] + ). Given parameters d, p, α and β, the Unrestrictive Common Due Date Problem (UCDDP) aims at finding a schedule minimizing the total penalty:

j∈J α j [d-C j ] + + β j [C j -d] +
This criterion is non-regular, i.e. is not a non-increasing function of completion time C j for any j ∈ J. Moreover, the penalty function is not a linear function of the completion times.

If all task penalties are equal, i.e. α j = β j for any task j ∈ J, the UCDDP is solvable in polynomial time [START_REF] Kanet | Minimizing the average deviation of job completion times about a common due date[END_REF]. If task penalties are symmetric, i.e. α j = β j for all j ∈ J, the problem is NP-hard [START_REF] Hall | Earliness-tardiness scheduling problems, I: weighted deviation of completion times about a common due date[END_REF]. Therefore, the problem that we consider with arbitrary α and β coefficients is also NP-hard.

A wide range of variables can be used to formulate a single machine scheduling problem as a Mixed Integer Program (MIP): completion time variables, timeindexed variables, linear ordering variables, positional date and assignment variables [START_REF] Queyranne | Polyhedral approaches to machine scheduling[END_REF]. However, only few works based on linear programming are proposed for the UCDDP. Biskup & Feldmann (2001a) propose a compact MIP formulation based on disjunctive variables, which is only given to compare the proposed heuristic method to an exact method for small instances. [START_REF] Van Den Akker | Combining column generation and lagrangean relaxation to solve a single-machine common due date problem[END_REF] propose a formulation based on an exponential number of binary variables using column generation and Lagrangian relaxation. In [START_REF] Falq | Mixed integer formulations using natural variables for single machine scheduling around a common due date[END_REF], we propose a MIP based on natural variables, similar to completion time variables, together with a compact MIP based on partition variables. Furthermore, the UCDDP have been solved through several other dedicated exact methods like branch-and-bound algorithms (e.g. [START_REF] Sourd | New exact algorithms for one-machine earliness-tardiness scheduling[END_REF]) and dynamic programming methods (e.g. [START_REF] Hall | Earliness-tardiness scheduling problems, I: weighted deviation of completion times about a common due date[END_REF], [START_REF] Hoogeveen | Scheduling around a small common due date[END_REF], [START_REF] Tanaka | An exact algorithm for the single-machine total weighted tardiness problem with sequence-dependent setup times[END_REF]). Moreover, a benchmark is provided in Biskup & Feldmann (2001a), together with a heuristic method. An important remark is that the efficiency of the latter approaches comes from the exploitation of dominance properties. In particular, the dedicated Branch-and-Bound proposed in [START_REF] Sourd | New exact algorithms for one-machine earliness-tardiness scheduling[END_REF] for an exact solving of UCDDP, solves up to 1000-tasks instances of the benchmark provided in Biskup & Feldmann (2001a).

Among the dominance properties, we can distinguish the structural ones, which allow to restrict the solution set to the set of the schedules having a specific structure. These properties are already taken into account in the compact formulation that we propose in [START_REF] Falq | Mixed integer formulations using natural variables for single machine scheduling around a common due date[END_REF]. Indeed, thanks to their structure, the dominant schedules are encoded by only n binary variables resulting in a formulation whose size does not depend on the time horizon, in contrast with the time-indexed formulations. In the present work, we focus on another type of dominance property called neighborhood based dominance properties. We propose to model them by linear inequalities.

One important contribution of this article is to propose inequalities strengthening the linear formulation in a new way. They are different from standard reinforcement and symmetry-breaking inequalities. The standard reinforcement inequalities cut fractional points to improve the lower bound obtained from the linear relaxation, the symmetry-breaking inequalities cut integer points, -which may be optimal-to reduce the search space. In contrast, Dominance inequalities aims at cutting sub-optimal integer points. This article is organized as follows. In Section 2, we recall structural dominance properties leading to reformulate the UCDDP into a partition problem; then we adapt some neighborhood based dominance properties used on schedules for partitions. Section 3 translates in a linear way the constraints proposed in Section 2, which leads to a new compact linear formulation, enriched with dominance inequalities. Section 4 presents some experimental results to show the contribution of considering dominance in an exact solving process, as in a heuristic approach.

Dominance properties

A subset of solutions is said dominant if it contains at least one optimal solution, and strictly dominant if it contains all the optimal solutions. For brevity, a schedule will be said dominant (resp. strictly dominant) if it belongs to a dominant (resp. strictly dominant) set. In this article, we will only use dominance properties about solutions, even if there also exist dominance properties about instances or problems [START_REF] Jouglet | Dominance rules in combinatorial optimization problems[END_REF].

Structural dominance properties

In order to describe dominant schedules, we first provide some useful definitions. A task j ∈ J is said early (resp. tardy) if it completes at time C j d ( resp. C j > d). Among the early tasks, the one completing at d (if it exists), is called the on-time task. The early-tardy partition of a schedule is the pair (E, T ) where E (resp. T ) is the early (resp. tardy) task subset. Moreover, we will use α-ratio (resp. β-ratio) to designate α j /p j (resp. β j /p j ). We define a block as a feasible schedule without idle time between task execution, a d-schedule as a feasible schedule with an on-time task, and a d-block as a d-schedule which is also a block. A schedule is said V-shaped if early tasks are scheduled in non-decreasing order of their α-ratios and the tardy ones in non-increasing order of their β-ratios.

The following lemma gives dominance properties already known for the unrestrictive common due date problem with symmetric penalties [START_REF] Hall | Earliness-tardiness scheduling problems, I: weighted deviation of completion times about a common due date[END_REF]. These results have been extended to asymmetric penalties in [START_REF] Falq | Mixed integer formulations using natural variables for single machine scheduling around a common due date[END_REF], using the same task shifting and exchange proof arguments. Lemma 1 (Falq et al. ( 2021))

The set of V-shaped schedules is strictly dominant for the UCDDP. The set of d-blocks is dominant for the UCDDP. Moreover, if unit earliness and tardiness penalties are positive, i.e. (α, β) ∈ R * + J ×R * + J , then the set of blocks is strictly dominant.

Thanks to the d-block dominance, we can make some assumptions on parameters p, α, β without loss of generality. A task having a zero processing time can be inserted between two tasks of a schedule without impacting other tasks. In particular in a d-block, such a task can be inserted at the due date, which incurs no penalty variation. Hence, we assume that ∀ j ∈ J, p j ∈ R * + . Since the due date d is large enough, a task having a zero earliness (resp. tardiness) penalty can be added at the beginning (resp. at the end) of a d-block without incurring any cost variation. Hence, we assume that ∀ j ∈ J, (α j , β j ) ∈ (R * + ) 2 .

Since they form a dominant set, we only consider V-shaped d-blocks in the following. Note that in such schedules, tardy tasks are completely processed after d, since there is no straddling task, i.e. there is no task starting before d and completing after d. Therefore, the early (resp. tardy) task set can be referred to the left (resp. right) side of d.

In general, a V-shaped d-block cannot be encoded by its early-tardy partition. Indeed, if two early (resp. tardy) tasks have the same α-ratio (resp. β-ratio), one cannot determine which one is processed first, they can be sequenced arbitrarily. However, swapping two such tasks in a V-shaped d-block results in another Vshaped d-block having the same penalty. Furthermore, all the V-shaped d-block having the same early-tardy partition have the same penalty. Based on this remark, two V-shaped d-blocks having the same early-tardy partition will be said equivalent. We denote by ∼ this relation.

We define an ordered bi-partition of a set A as a couple (A 1 , A 2 ) where {A 1 , A 2 } is a partition of A. This is not only a partition into two subsets since the two subsets are not symmetric, i.e. (A 1 , A 2 ) (A 2 , A 1 ). Note that the early-tardy partition of a schedule is an ordered bi-partitions. More precisely, there is a one-to-one correspondence between equivalence classes for ∼ and ordered bi-partitions (E, T ) of J where E ∅. Indeed, the set of the early tasks of a d-block cannot be empty since it contains at least the on-time task. Let P * 2 (J) denote the set of such (E, T ), and P 2 (J) denote the set of any ordered bi-partition of J. In the sequel, we will only say partition to refer to an ordered bi-partition or to an early-tardy partition, when there is no ambiguity.

For any partition (E, T ) ∈ P * 2 (J), let f (E, T ) denote the penalty of the equivalence class (E, T ), that is the penalty of any V-shaped d-block of partition (E, T ). For sake of consistency, we extend the definition of f to P 2 (J) by setting f (∅, J) as the penalty of any V-shaped block starting at time d. From now on, our aim is to find a partition of J minimizing f , since the UCDDP can be formulated as follows.

(F 1 ) : min (E,T ) ∈ P * 2 (J)

f (E, T )

Neighborhood based dominance properties

Let us recall some definitions used in local search context (Aarts & Lenstra, 2003).

A neighborhood function N is a function which associates to any solution S a subset of solutions N(S), called the neighborhood of S. A solution of N(S) is called a neighbor of S. Moreover, a solution S is locally optimal with respect to minimizing function φ if φ(S) φ(S ), for any neighbor S ∈ N(S). If, on the contrary, there exists S ∈ N(S) such that φ(S ) < φ(S), S is dominated (by S ).

Given a neighborhood function, the set of locally optimal solutions always contains all optimal solutions, and is therefore a strictly dominant set. This statement can be seen as a generic dominance property denoted G. This kind of dominance property is commonly used in local search. Indeed, a step of a local search procedure consists in enumerating some of the neighbors of a given solution S, computing their values and, if a better solution is found, then moving to the best solution found at the current iteration, which is equivalent to discard S since it is dominated.

• Operation-based neighborhoods

We call operation any (eventually partial) function, which maps a solution to another solution. In this work, we will consider neighborhood functions based on a set of operations. The neighborhood of a solution S is then the set of the solutions obtained by applying to S any operation defined on S. This kind of neighborhood functions allows to use the generic dominance property G in a different way. Instead of considering sequentially the neighborhood of each solution, we will consider sequentially each operation. For each one, any solution is compared to its neighbor obtained using this operation. This differentiates the solution-centered and the operation-centered point of view, as developed in the next section/point.

• The solution-centered and the operation-centered points of view Figures 1 and 2 illustrate these two points of view on the same set of solutions {A, B, . . . , N} represented by the blue points. We consider a given set of operations and the associated neighborhood function. An arrow is set from a solution X to a solution Y if X is compared to Y in order to determine whether X is dominated by Y.

In Figure 1, we focus on one solution, J, which is compared to all the solutions obtained by applying to J an operation defined on J. Since J is compared to all its neighborhood, i.e. {C, F, G, I, M, N}, one can determine if J is locally optimal. In Figure 2, we focus on one operation. All the solutions where this operation can be applied are compared with the obtained neighbors. Solutions A, B, C, G and N are not compared to others solutions, since the considered operation cannot be applied on these solutions. Conversely, since solutions D, E, F, H, I, J, K, L and M are compared to one of their neighbors, they might be discarded. For example if E is better than H, H can be discarded, no matter whether H is better than K. However, we cannot say that the non-discarded solutions are locally optimal, since only one neighbor is taken into account. To this end, all the operations must be considered. Let us now present the main ideas of a generic method to handle operationbased dominance property. First, find operations inducing an objective function variation that can be explicitly expressed. Then derive, for each operation, a constraint ensuring that any solution is not dominated by applying this operation. Using all the obtained constraints, only the dominant set of property G subsists as all the non-locally optimal solutions have been removed (where local is understood with respect to the neighborhood defined by the operations). In Section 2.3, we introduce two families of operations over the partitions, and follow this method, resulting in constraints (I u ), (I u ) and (S u,v ) translating Property 3.

• A • B • C • D • E • F • G • H • J • K • L • M • N • I Figure 1: Illustration of the solution-centered point of view • A • B • C • D • E • F • G • H • I • J • K • L • M • N

Operations for early-tardy partitions

In the case of problems where the solutions are partitions, two operations are commonly considered to define a neighborhood: the insertion, which consists in moving two elements, one in each subset and the swap which consists in swapping two elements of each subset. Let us define solutions which are locally optimal in the neighborhoods associated with these operations.

Definition 2 Let (E, T ) ∈ P * 2 (J). (i) (E, T ) is an insert local optimum if        ∀v ∈ T, f E, T f E ∪{v}, T \{v} , ∀u ∈ E, f E, T f E \{u}, T ∪{u} . (ii) (E, T ) is a swap local optimum if ∀u ∈ E, ∀v ∈ T, f E, T f E \{u}∪{v}, T ∪{u}\{v} .

Property 3

The set of insert locally optimal partitions, as well as the set of swap locally optimal partitions, is strictly dominant when minimizing f over P * 2 (J).

In terms of scheduling, the insert operation consists in removing a task j ∈ J from the early (resp. tardy) side and inserting j on the tardy (resp. early) side, as early (resp. as tardy) as possible according to its β-ratio (resp. α-ratio). The tasks scheduled before (resp. after) j are right-shifted (resp. left-shifted) by p j time units. The swap operation consists in sequentially inserting an early task on the tardy side and inserting another tardy task on the early side, or vice-versa, as described above. Let us denote by insert u (S) (resp. swap u,v (S)) the schedule obtained from a schedule S by the insert (resp. swap) operation on task u (resp. on tasks u, v).

Since the insert and swap operations are fundamentally defined on partitions (i.e. over equivalence classes of V-shaped d-blocks), these operations preserve the equivalence. Therefore we have ∀u ∈ J, S ∼ S ⇒ insert u (S) ∼ insert u (S ) and

∀(u, v) ∈ J 2 , S ∼ S ⇒ swap u,v (S) ∼ swap u,v (S )
, assuming that u and v are not on the same side in these schedules.

In order to provide an expression of the penalty variation induced by such an operation, it is convenient to choose a specific V-shaped d-block as representative of an equivalence class. More precisely, the chosen representative will depend on the inserted task (resp. on the swapped tasks). Let us introduce some notations, related to a given task u ∈ J. To describe the early side of a V-shaped d-block with regard to u, the set of remaining tasks J \{u} can be split according to their α-ratio into two subsets A(u) and Ā(u) defined as follows.

A(u) = i ∈ J α i p i > α u p u and Ā(u) = i ∈ J \{u} α i p i α u p u
Similarly, we introduce the two following subsets to describe the tardy side.

B(u) = i ∈ J β i p i > β u p u and B(u) = i ∈ J \{u} β i p i β u p u
Note that if u is early in a V-shaped d-block, early tasks belonging to A(u) are necessarily scheduled after u, because of their α-ratio. To illustrate, we can observe in Figure 3 that the tasks of A(u) ∩ E are scheduled after u in the schedule S.

Conversely, an early task of Ā(u) is not necessarily scheduled before u, when its α-ratio is the same as u. However, we will consider a representative where u is placed after all the early tasks of Ā(u), that is as tardy as possible according to its α-ratio. Similarly, in case where u is tardy, we will consider a representative where u is scheduled before all tardy tasks of B(u). Such a representative will be called a u-canonical V-shaped d-block. The schedule S (resp. S ) in Figure 3 represents the shape of a u-canonical V-shaped d-block when u is early (resp. tardy).

• Penalty variation induced by an insert operation

Given u ∈ J, let (E, T ) be a partition such that u ∈ E. We aim to express the variation of f induced by the insertion of u in T, i.e. between (E, T ) and (E \{u}, T ∪{u}). To this end, we consider a u-canonical representative S of (E, T ), and the V-shaped d-block S obtained from S by inserting u in T , as early as possible, that is just after tardy tasks of B(u). Note that S is thus a u-canonical representative of E \ {u}, T ∪ {u} . Let (e, t) (resp. (e , t )) denote the earliness and tardiness vector of tasks in S (resp. S ).

As we can observe in Figure 3, early tasks of A(u) and tardy tasks of B(u) are identically scheduled in S and S . The penalty variation induced by the insertion of u is then only due to the move of u and to the right-shifting of early tasks of Ā(u) and tardy tasks of B(u). Each task j ∈ Ā(u)∩E of S is postponed p u time units later in S , while staying early. We then have e j = e jp u . The earliness penalty of task j in S is therefore α j e j = α j e j -α j p u , which represents a reduction of α j p u compared to its earliness penalty in S. Summing up over Ā(u)∩E, we obtain a reduction of the earliness penalties of p u α Ā(u)∩ E , where x(I) = i∈I x i for any x ∈ R n and any I ⊆ [1..n], Similarly, the right-shifting of tasks in B(u)∩T induces an increase of the tardiness penalties of p u β B(u)∩T , since for each j ∈ B(u)∩T we have β j t j = β j (t j + p u ).

S Ā(u)∩E A(u)∩E B(u)∩T B(u)∩T d u S Ā(u)∩E A(u)∩E B(u)∩T B(u)∩T d u p u p u
Moreover, since e u = p A(u) ∩ E , removing u from the early side induces a reduction of its earliness penalty of α u p A(u)∩E . Similarly, introducing u on the tardy side induces an increase of its tardiness penalty of

β u t u = β u p B(u)∩T +p u since t u = 0.
Finally, the penalty variation between S and S , is given by the following expression.

∆ u (E, T ) = -α u p A(u)∩E + β u p B(u)∩T + p u + p u β B(u)∩T -α Ā(u)∩E
Since S and S are representative of (E, T ) and (E\{u}, T∪{u}) respectively, ∆ u (E, T ) is also the variation of f induced by the insertion of u in T . Property 4.(i) follows.

The insert operation of the tardy task u on the early side applied on schedule S provides schedule S. Note that this statement requires that S is u-canonical. This observation allows to establish that the penalty variation induced by inserting u on the early side in S is simply -∆ u (E, T ), and results in Property 4.(ii). Property 4

Let (E, T ) be a partition.

(i) For any u ∈ E, f (E \{u}, T ∪{u}) = f (E, T ) + ∆ u (E, T ). (ii) For any v ∈ T , f (E ∪{v}, T \{v}) = f (E, T ) -∆ v (E, T ).
Let us introduce, for a given task u ∈ J, the two following constraints.

u ∈ E ⇒ ∆ u (E, T ) 0 (I u ) u ∈ T ⇒ ∆ u (E, T ) 0 (I u )
Thanks to Property 4, (I u ) (resp. (I u )) discards every partition where u is early (resp. tardy), but which would have a lower penalty if u were tardy (resp. early). Note that the constraint ∆ u (E, T ) 0 (resp. ∆ u (E, T ) 0) might be not satisfied by an optimal solution where u is tardy (resp. early). Moreover, only considering constraints (I u ) and (I u ) for a given u is not sufficient to discard all insert-dominated partitions. Indeed, these constraints do not take into account the whole insert-neighborhood of the solutions as each one is only compared to its neighbor obtained by an insert operation on task u. It is thus needed to consider constraints (I u ) and (I u ) for every u ∈ J to translate the dominance of the insert local optimal solutions. In Section 3, we will explain how these constraints can be used in a linear formulation.

• Penalty variation induced by a swap operation

This section follows the same organization as the previous one. The objective is to obtain constraints translating the dominance of the swap locally optimal solutions.

Given (u, v) ∈ J 2 such that u v, let (E, T ) be a partition such that u ∈ E and v ∈ T , and S be a u and v-canonical representative of (E, T ). We denote by S the schedule obtained from S by swapping u and v so that S is also both u and v-canonical, that is by scheduling u after all the early tasks having the same α-ratio α u /p u and v before all the tardy tasks having the same β-ratio β v /p v . S is a representative of E \{u}∪{v}, T \{v}∪{u} . Let (e, t) (resp. (e , t )) denote the earliness and tardiness vector of tasks in S (resp. S ).

If α v p v < α u p u , all the early tasks of Ā(v) are scheduled before u in S, since the V-shaped property holds, but the early tasks of A(v)\{u} can be scheduled before or after u. This case is illustrated in Figure 4. Conversely, if α v p v α u p u , all the early tasks of A(v) are scheduled after u in S, but those of Ā(v) \ {u} can be scheduled before or after u. This case is illustrated in Figure 5.

A(v)∩ Ā(u)∩E S Ā(u)∩E A(u)∩E A(v)∩E Ā(v)∩E d | u S Ā(u)∩E ∪{v} A(u)∩E A(v)\{u}∩E Ā(v)∩E d | v
v p v < α u p u A(u)∩ Ā(v)∩E S Ā(u)∩E A(u)∩E A(v)∩E Ā(v)∩E d | u S Ā(u)∩E A(u)∩E ∪{v} A(v)∩E Ā(v)\{u}∩E d | v
v p v α u p u B(v)∩ B(u)∩T S B(v)∩T B(v)∩T B(u)∩T B(u)∩T B(v)∩T d | v S B(v)∩T ∪{u} B(v)∩T B(u)∩T B(u)∩T \{v} d | u
Figure 6: Tardy side variation induced by swapping u and v when

β v p v β u p u B(u)∩ B(v)∩T S B(v)∩T B(v)∩T B(u)∩T B(u)∩T d | v S B(v)∩T B(v)∩T ∪{u} B(u)∩T \{v} B(u)∩T d | u
Figure 7: Tardy side variation induced by swapping u and v when

β v p v > β u p u
Similarly, we can distinguish two cases depending on the relative order of the β-ratios of u and v, as illustrated in Figures 6 and7.

As shown in Figures 4 and5, the earliness of u in S is e u = p A(u)∩E , while the tardiness of u in S is e u = p B(u)\{v}∩T +p u (Cf. Figures 6 and7). Note that v is removed from B(u) since v is not tardy in S , and therefore cannot contribute to the tardiness of u. Similarly, we have t v = p B(v)∩T + p v and e v = p A(v)\{u}∩E since u is not early in S .

Moreover, the tasks of A(u) ∩ E are identically scheduled in S and S only if

α v p v < α u p u .
In this case, tasks of Ā(u)∩ E are not consecutive in S since v separates them into two blocks: tasks of Ā(v)∩E which have been left-shifted by p v -p u time units, and tasks of A(v)∩ Ā(u)∩ E which have been right-shifted by p u time units (Cf. Figure 4). In the opposite case, i.e. if α v p v α u p u , tasks of A(u) ∩ E are not consecutive in S . Indeed, v separates them into two blocks: tasks of A(v) ∩ E which are identically scheduled in S and S ; and tasks of A(u)∩ Ā(v)∩E which are left-shifted by p v time units. Moreover, in that case tasks of A(v)∩E are left-shifted by p vp u time units (Cf. Figure 5). Note that, in the two previous paragraphs, as in Figures 4 to 7, we assumed that p v -p u 0. In the contrary case, i.e. if p v -p u < 0, tasks are not left-shifted by p v -p u time units but rather right-shifted by p up v time units.

From these observations, we can express the earliness penalty variation of all the tasks in E \ {u} by expressing, in each case, the penalty variation induced by a block shifting as described for the insert operation. The same method can be applied for the tasks of T \ {v}. The reader can refer to Figures 6 and7 for illustration. Finally, the penalty variation between S and S is given by the following expression.

∆ u,v (E, T ) = -α u p A(u)∩E + β u p B(u)\{v}∩T + p u -β v p B(v)∩T + p v + α v p A(v)\{u}∩E +          (-p u + p v ) α Ā(v)∩E -p u α A(v)∩ Ā(u)∩E if α v p v < α u p u (-p u + p v ) α Ā(u)∩E + p v α A(u)∩ Ā(v)∩E otherwise +          (p u -p v ) β B(v)∩T + p u β B(v)∩ B(u)∩T if β v p v β u p u (p u -p v ) β B(u)∩T -p v β B(u)∩ B(v)∩T otherwise Property 5 Let (E, T ) be a partition. For any (u, v) ∈ E ×T , f E \{u}∪{v}, T \{v}∪{u} = f (E, T ) + ∆ u,v (E, T ).
For a given pair of tasks (u, v) ∈ J 2 such that u v, let us introduce the following constraint.

(u, v) ∈ E ×T ⇒ ∆ u,v (E, T ) 0 (S u,v )
Thanks to Property 5, (S u,v ) discards every partition where u is early, v is tardy and which would have a lower penalty in the contrary case. As for the insert operation, constraint ∆ u,v (E, T ) 0 might be not satisfied by optimal partitions where u is not early or v is not tardy. Moreover, it is needed to consider constraints (S u,v ) for every (u, v) ∈ J 2 such that u v to translate the dominance of swap locally optimal solutions. In Section 3, we will explain how these constraints can be used in a linear formulation.

Neighborhood based dominance in linear programming

The dominance properties described in Section 2 can be used in a linear programming framework. In this section, we provide linear inequalities translating constraints (I u ), (I u ) and (S u,v ) for all tasks u and v.

A compact MIP translating F 1

We first recall the linear compact MIP formulation for the UCDDP given in [START_REF] Falq | Mixed integer formulations using natural variables for single machine scheduling around a common due date[END_REF]. In this formulation denoted by F 2 , a partition is encoded by a vector δ of n binary variables. Given δ ∈ {0, 1} J , the partition encoded by δ is

{ j ∈ J | δ j = 1}, { j ∈ J | δ j = 0} .
In other words, for each j ∈ J, δ j indicates whether j is early or not.

Moreover, F 2 also uses X, a vector of binary variables indexed by the set J < = {(i, j) ∈ J 2 | i < j}, in order to linearize products of type δ i δ j as proposed by [START_REF] Fortet | L'algèbre de Boole et ses applications en recherche opérationelle (Boole's algebra and its applications in operations research)[END_REF]. Lemma 6 [START_REF] Fortet | L'algèbre de Boole et ses applications en recherche opérationelle (Boole's algebra and its applications in operations research)[END_REF])

Let (δ, X) ∈ R J ×R J < .
If δ ∈ {0, 1} J and (δ, X) satisfies the following inequalities:

∀(i, j) ∈ J < , X i, j δ i -δ j (1)

∀(i, j) ∈ J < , X i, j δ j -δ i (2) ∀(i, j) ∈ J < , X i, j δ i +δ j (3) ∀(i, j) ∈ J < , X i, j 2-(δ i +δ j ) (4) then ∀(i, j) ∈ J < , X i, j = 1 if δ i δ j 0 otherwise , (1 -δ i )(1 -δ j ) = 2-(δ i +δ j )-X i, j 2 and
δ i δ j = δ i +δ j -X i, j 2 .
Let us consider the polyhedron P = (δ, X) ∈ R J ×R J < (1 -4) and the set of its integer points int δ (P) = P∩{0, 1} J ×{0, 1} J < . From Lemma 6, if a partition (E, T ) is encoded by δ, there exists a unique X such that (δ, X) ∈ int δ (P). We will say that (δ, X) encodes (E, T ).

To obtain an expression of the penalty of a partition from its encoding (δ, X), two orders on J are introduced. Let ρ and σ be two functions from [1..n] to J, such that α ρ(k)

p ρ(k) k∈[1..n] and β σ(k) p σ(k) k∈[1.
.n] are non-increasing A d-block is said ρ-σ-shaped if the early (resp. tardy) tasks are processed in decreasing order of ρ -1 (resp. increasing order of σ -1 ). Each equivalence class of V-shaped d-blocks admits a unique ρ-σ-shaped representative. This representative is used to provide the following expression of the penalty of the partition encoded by (δ, X) ∈ {0, 1} J ×{0, 1} J < .

g(δ, X) = j∈J α j          ρ -1 ( j)-1 k=1 p ρ(k) δ j +δ ρ(k) -X j,ρ(k) 2          + β j          σ -1 ( j)-1 k=1 p σ(k) 2-(δ j +δ σ(k) ) -X j,σ(k) 2 + p j (1-δ j )         
Finally, the compact formulation for the UCDDP provided in Falq et al. ( 2021) is the following. (F 2 ) : min

(δ,X) ∈ int δ (P) g(δ, X)
Formulation F 2 is a direct linear translation of F 1 . Indeed, there is a one to one correspondence between their solution sets, i.e. between int δ (P) and P 2 (J), and g(δ, X) = f (E, T ) for any (δ, X) encoding (E, T ).

3.2. Linear inequalities translating constraints (I u ), (I u ) and (S u,v )

• Linear inequalities translating constraints (I u ) and (I u ) for any u ∈ J Let u ∈ J. If δ ∈ {0, 1} J encodes a partition (E, T ), the penalty variation ∆ u (E, T ) can be expressed linearly from δ as follows.

∆ u (δ) = -α u i∈A(u) p i δ i + β u i∈B(u) p i (1-δ i )+ p u + p u i∈ B(u) β i (1-δ i ) - i∈ Ā(u) α i δ i Moreover, if δ ∈ {0, 1} J encodes (E, T ), we can translate constraint (I u ) as follows. (E, T ) satisfies (I u ) ⇔ u ∈ E and ∆ u (E, T ) 0 or u ∈ T ⇔ δ u = 1 and ∆ u (E, T ) 0 or δ u = 0
To unify these two cases into one inequality, we introduce the following constant which is an upper bound of -∆ u (δ) for any δ ∈ {0, 1} J .

M u = α u p A(u) -β u p u + p u α Ā(u) δ ∈ {0, 1} J such that δ v + (1-δ u ) = 2, i.e.
such that (δ u , δ v ) = (0, 1), this inequality might not be satisfied if M u,v 0, since -M u,v -2 M u,v in this case. To provide an inequality satisfied by every δ ∈ {0, 1} J such that (δ u , δ v ) (1, 0) we introduce the following constant.

M u,v =        M u,v if M u,v 0 M u,v /2 otherwise We then have -M u,v -2 M u,v and -M u,v -M u,v
. Therefore, the following inequality is satisfied by every δ ∈ {0, 1} J such that

δ v + (1-δ u ) ∈ {1, 2}, i.e. such that (δ u , δ v ) (1, 0). ∆ u,v (δ) -M u,v δ v + (1-δ u ) (7)
Conversely, for every

δ v + (1-δ u ) = 0, δ ∈ {0, 1} J such that i.e. such that (δ u , δ v ) =
(1, 0), inequality ( 7) is satisfied if and only if ∆ u,v (δ) 0. Finally, δ ∈ {0, 1} J satisfies (7) if and only if the partition encoded by δ satisfies (S u,v ). Property 8 follows.

Property 8

Let δ ∈ {0, 1} J and let (E, T ) be the partition encoded by δ. For any (u, v) ∈ J 2 such that u v, δ satisfies (7) if and only if (E, T ) satisfies constraint (S u,v ).

In the sequel, insert inequalities (5) and ( 6) and swap inequalities (7) will be called dominance inequalities.

Link between dominance properties and operations

In this section, we show how to benefit from the fact that each dominance inequality is based on an operation. When a vector δ ∈ {0, 1} J encoding a partition (E, T ) does not satisfy a dominance inequality, applying the corresponding operation to (E, T ) provides a partition with a strictly lower penalty. The following property, resulting from properties 4, 5, 7 and 8, formally states this result.

Property 9

Let δ ∈ {0, 1} J and let (E, T ) be the partition encoded by δ. For any u ∈ J, (i) δ does not satisfy (5) for u if and only if u ∈ E and f E\{u}, T∪{u} < f (E, T ), (ii) δ does not satisfy (6) for u if and only if u ∈ T and f E∪{u}, T\{u} < f (E, T ).

(iii) Moreover, for any (u, v) ∈ J 2 such that u v, δ does not satisfy (7) for (u, v) if and only if (u, v) ∈ E ×T and f E \{u}∪{v}, T \{v}∪{u} < f (E, T ).

Proof : Let us fix u ∈ J. From Property 7, δ does not satisfy inequality (5) if and only if (E, T ) does not satisfy constraints (I u ), which is equivalent to u ∈ E and ∆ u (E, T ) < 0.

Using Property 4, it is equivalent to u ∈ E and f E \{u}∪{v}, T \{v}∪{u}f (E, T ) < 0, which proves (i). The proofs of (ii) and (iii) follow the same scheme.

Property 9 will be used in Section 4 to propose a local search procedure. The following corollary, which directly derives from the negation of statements (i), (ii) and (iii), ensures that the solution provided by this local search procedure is an insert and swap local optimum.

Corollary 10

Let δ ∈ {0, 1} J and let (E, T ) be the partition encoded by δ.

(i) (E, T ) is an insert local optimum if and only if δ satisfies inequalities (5) and ( 6) for all u ∈ J. (ii) (E, T ) is a swap local optimum if and only if δ satisfies inequality (7) for all (u, v) ∈ J 2 such that u v.

The following section presents experimental results to assess the practical relevance of the dominance inequalities.

Numerical results

All experiments are carried out using a single thread with Intel(R) Xeon(R) X5677, @ 3.47GHz, and 144Gb RAM. Linear programs (LP) and MIP are solved with Cplex 12.6.3.0.

The numerical experiments are performed on the instance benchmark proposed by Biskup & Feldmann (2001a), available online on OR-Library [START_REF] Biskup | Common due date scheduling[END_REF]. For each number of tasks n ∈ {10, 20, 50, 100, 200}, ten triples (p, α, β) of N * n 3 are given. For each one, we assume that d = p(J), so that the due date is unrestrictive. For the sake of comparison, we additionally construct instances with n ∈ {60, 80} (resp. n ∈ {120, 150, 180}) by only considering the first n tasks of the previous 100-task (resp. 200-task) instances. Unless otherwise specified, the gap, time and number of nodes presented in the following tables are average values over the ten instances for a given n and the time limit is set to 3600 seconds.

To measure the improvement induced by the insert or swap inequalities, we compare the four following formulations.

F 2 : the formulation defined in Section 3, only with inequalities (1-4) F i : the formulation obtained from F 2 by adding ( 5) and ( 6) for all u ∈ J F s : the formulation obtained from F 2 by adding (7) for all (u, v) ∈ J 2 s.t. u v F i+s : the formulation obtained from F s by adding ( 5) and ( 6) for all u ∈ J For a given formulation F, we distinguish two settings: a setting with all available features, that is using Cplex default, denoted by F d , and a setting with less Cplex features, denoted by F l . Two types of features are disabled in this setting: the cut generation, which produces reinforcement inequalities and adds them to the formulation, and the primal heuristic procedures. The cut generation is disabled in order to measure the impact of the dominance inequalities on the linear relaxation value of the formulation F 2 , rather than their impact on the linear relaxation value of a strengthened formulation. The primal heuristic procedures have been disabled to focus on the lower bound since we have other methods to quickly obtain good feasible solutions (Cf. Section 4.3).

This results in eight formulation settings:

F 2 l , F i l , F s l , F i+s l , F 2 d , F i d , F s d and F i+s d .
For each one, inequalities (1-4), as well as inequalities (5-7) when included, are added initially.

Let us recall that only the δ variables need to be integer in F 2 . Indeed, from Lemma 6, if δ ∈ {0, 1} J , inequalities (1-4) ensures that X ∈ {0, 1} J < . It is also the case for F i , F s and F i+s . Therefore, unless otherwise specified, δ variables are set as binary variables, while X variables are set as continuous variables. Consequently, the branching decisions only involve δ.

Solving MIP formulations to optimality

Table 1 provides the results obtained by solving MIP to optimality, using the eight formulation settings. Each line corresponds to the ten instances of a given size n. More precisely, Table 1 entries are the following.

#opt : the number of instances solved to optimality within the time limit time : the average running time in seconds over the instances solved to optimality #nd : the average number of nodes, except the root node, in the search tree, over the instances solved to optimality For a given formulation setting, we choose to stop the run at a line of the table if less than 5 over ten instances are solved to optimality. For the subsequent lines, we report a "-" in the table.

Using formulation setting F 2 l , the ten n-task instances are solved to optimality within the time limit for n up to 50. In contrast, using F i l , it is the case for n up to 60, using F s l for n up to 120 and using F i+s l for n up to 150. Within approximately 5 minutes, F 2 l solves 50-task instances, F i l 60-task instances, F s l 100-task instances and F i+s l 120-task instances. This computation time decrease is due to a drastic reduction in the number of nodes. For example, for n = 50 the number of nodes goes from more than 53 000 for F 2 l to only 31 for F i+s l . With this latter formulation setting, the number of nodes is low, it is at most 200, even for large size instances. However, the time limit is reached for some 180-and 200-task instances, since the size of the linear program solved at each node is large.

In light of the four first columns of Table 1, we can conclude that, with less Cplex features, adding insert and swap inequalities significantly reduces the number of nodes and hence the computation time. More precisely, adding only swap inequalities is better than adding only insert inequalities, but adding both of them provides the best performance.

The four last columns show the same improvement in terms of computation time and number of nodes when Cplex default features are used.

Let us now focus on the 4th and 5th columns to compare the impact of the dominance inequalities and the impact of Cplex default features. For small instances, i.e. n ∈ {10, 20}, Cplex default features allow to solve the problem at the root node (Cf. F 2 d columns). However, from n = 50, the number of nodes grows fast, so that no 60-task instance can be solved within 3600 seconds. Conversely, we already noticed that adding swap and insert inequalities limits the number of nodes (Cf. F i+s l columns), so that the ten 150-task instances are solved within 3600 seconds. Finally, adding the swap and insert inequalities provides better results than adding Cplex default features.

Up to size 60, F i+s d solves all instances at the root node, and is faster than F i+s l . For larger instances, except 200-task instances, F i+s l and F i+s d solves the problem in similar computation times, even if F i+s l explores a smaller number of nodes: for example, it is two times smaller for n = 150. For n = 200, F i+s d solves 4 over 10 instances, while F i+s l only solves 1 over 10 instances. To conclude, F i+s l and F i+s d offer comparable performances, so that, for both settings, the formulations providing the best results are the ones with insertion and swap inequalities.

Lower bound obtained at the root node

To further investigate the impact of dominance inequalities, we focus in this section on the root node of the search tree for different formulation settings. More precisely, we compare the different lower bounds obtained at the root node.

In the Cplex framework, setting the node limit to 0 allows to only solve the root node of a MIP: the branch-and-bound algorithm is stopped before the first branching. If Cplex default features are activated, the preprocessing is applied and the cuts are added before the algorithm stops. For a given formulation setting F, the corresponding run with the node limit set to 0 is denoted by F-rn. This results in eight runs:

F 2 l -rn, F i l -rn, F s l -rn, F i+s l -rn, F 2 d -rn, F i d -rn, F s d -rn, F i+s d -rn.
Note that, in the Cplex framework, solving F-rn is different from solving the linear relaxation of F, denoted by F-lp. Indeed, F-lp is obtained by setting δ variables as continuous variables, which desactivates most of the Cplex features. In particular, the reinforcement cuts cannot be added since they are not valid for the relaxed formulation. Similarly, the inference procedure on the binary variables cannot be applied. We run the four linear relaxations F 2 -lp, F i -lp, F s -lp and F i+s -lp. Surprisingly, the obtained values are the same for these four relaxations. In other words, adding insert and swap inequalities does not improve the linear relaxation value. Therefore, we only present in Table 2 the results for F 2 -lp.

To measure the quality of the nine different lower bounds obtained, we compute, when it is possible, the optimality gap, i.e. (OPT -LB)/OPT where OPT denotes the optimal value and LB the lower bound. When the optimal value is not known, we compute a gap using the best upper bound that we get U B, i.e. (U B-LB)/U B. Such gaps are indicated with a "*" in Table 2. For each of the nine runs, the entries of Table 2 are the following. L-gap : the average optimality gap of the lower bound obtained at the root node time : the average running time in seconds over the ten instances

The obtained lower bound is exactly the same using F 2 -lp, F 2 l -rn or F s l -rn. We deduce that with less Cplex features and without insert inequalities, setting the δ variables as binary or continuous variables, provides the same lower bound. Moreover, this lower bound is quite weak, since the average optimality gap is larger than 40% even for the 10-task instances. The computation times using F 2 -lp and F 2 l -rn are similar: 2 seconds for the 100-task instances and about 12 minutes for the 500-task instances. The computation time required for F s l -rn is larger: almost 20 seconds for the 100-task instances and 47 minutes for the 500-task instances.

The lower bound obtained when only considering the insert inequalities is slightly better when the δ variables are set as binary variables for n ∈ {10, 20}. Indeed, the average optimality gap is 33% instead of 41% when n = 10, and 66% instead of 68% when n = 20 (Cf. F 2 -lp and F i l -rn columns). The computation time using F i l -rn is comparable to the computation time using F 2 -lp and F 2 l -rn.

The lower bound obtained when considering both insert and swap inequalities, is significantly better when the δ variables are set as binary variables. Indeed, the average optimality gap is smaller than 39% for any value of n, and it is equal to 0 for n = 10 (Cf. F i+s l -rn column). The computation time for F i+s l -rn is between those for F i l -rn and F s l -rn: 14 seconds for the 100-task instances and about 30 minutes for the 500-task instances.

The lower bound provided using F 2 d -rn, is better than the one obtained using F 2 l -rn, that is with less Cplex features. Indeed, the average optimality gap is 7% instead of 41% when n = 10, and 46% instead of 94% when n = 120. However, the lower bound is weaker that the one obtained for F i+s l , whose optimality gap is 0 for n = 10 and 38% for n = 120. Moreover, the computation times using F 2 d -rn increases fast with the increase of n so that the root node cannot be solved within one hour for sizes larger than 120.

Combining Cplex features with insert inequalities gives almost the same results (Cf. F i d -rn column). Conversely, combining Cplex features with swap inequalities gives better results (Cf. F s d -rn column). In particular, the average computation time is reduced so that instances up to size 200 can be solved. Moreover, the optimality gap is less than 22% for all solved instances. Finally, using F i+s d -rn gives even better results, the average optimality gap does not exceed 15%, even for 200task instances, which are solved in 418 seconds, instead of 1200 using F s d -rn.

In a nutshell, combining insert and swap inequalities is the best to obtain a lower bound at the root node. Not using Cplex features allows its fast computation (Cf. F i+s l -rn column). Conversely, using them allows to obtain a better lower bound at the expense of the computation time (Cf. F i+s d -rn column).

Using swap and insert inequalities to obtain an upper bound

In this section, we propose two upper bounds on the optimal value. The first one is derived from the fractional solution obtained at the root node by a simple rounding procedure. The second one is obtained by applying in addition a local search procedure.

We derive an integer solution (δ, X) by rounding a fractional solution ( δ, X), as follows.

Note that this algorithm can be seen as a local search procedure for the neighborhood associated to the insert and swap operations. Moreover, this procedure can be applied to any integer solution. Particularly, by sake of comparison we apply it to the solutions obtained by the heuristic "Heur II" provided by Biskup & Feldmann (2001a).

We finally compare the upper bounds given by the six following heuristic solutions.

BF : the solution obtained by the Biskup and Feldmann heuristic BF+ : the solution obtained by applying Insert_swap_improvement to BF R1 : the solution obtained by rounding the fractional solution of F 2 -lp R1+ : the solution obtained by applying Insert_swap_improvement to R1 R2 : the solution obtained by rounding the fractional solution of F i+s d -rn R2+ : the solution obtained by applying Insert_swap_improvement to R2

In the sequel, we will use the same notation for both a heuristic solution and its value, which provides an upper bound on the optimal value. To measure the quality of these upper bounds, Table 3 presents their optimality gap denoted by U-gap, i.e. (U B-OPT )/OPT where OPT denotes the optimal value and U B the upper bound. The Biskup and Feldmann heuristic provides a solution in less than 1 second. Applying rounding and Insert_swap_improvement to a fractional solution provides a solution in less than 1 second for instances up to size 200. Therefore, the time needed to obtain R1 and R1+ (resp. R2 and R2+) is essentially the computation time required to solve F 2 -lp (resp. F i+s d -rn) given in Table 2.

As shown in Table 3, BF is a good upper bound. Indeed, its optimality gap is smaller than 0.35% for instance sizes larger than 50. However, this bound is improved by Insert_swap_improvement: the optimality gap of BF+ is smaller than 0.02% for all the instances. With an optimality gap larger than 170%, R1 is a very weak upper bound, while R1+, with an optimality gap smaller than 0.01%, is very good, and even slightly better than BF+. With an optimality gap smaller than 17%, R2 is a better upper bound than R1, and R2+ is exactly the same as R1+.

Finally, BF+, R1+ and R2+ are very good upper bounds. However it is worth noticing that even if the computation time to obtain BF+ is about 1 second, the bound is obtained without any guarantee, since no lower bound is provided. Conversely, the computation time to obtain R2+ is larger: 25 seconds for n = 100 and about 7 minutes for n = 200, but a lower bound is provided. R2+ is then guaranteed to be at 14% of the optimal value for n = 100, and at 15% for n = 200 (Cf. L-gap of F i+s d -rn in Table 2). R1+ is a compromise between BF+ and R2+. Indeed, for instances up to size 200, R1+ is provided in less than 20 seconds together with a lower bound, but the guarantee obtained from this lower bound is quite weak (97% for n = 200, Cf. F 2 l -rn in Table 2).

Insert and swap operations use cases

Insert and swap operations can be used in different ways. Table 4 presents the best way to use them depending on the expected solution quality.

-To obtain an upper bound: apply rounding and Insert_swap_improvement to the fractional solution given by F 2 -lp. (Cf. R1+ column in Table 4).

-To obtain an upper bound with a better guarantee than the one obtained with R1+: apply rounding and Insert_swap_improvement to the fractional solution given by F i+s d -rn. (Cf. R2+ column in Table 4). -To obtain a 5%-approached solution: use F i+s d , setting the gap limit to 5%. (Cf. F i+s d -5% column in Table 4). -To obtain an exact solution: use F i+s d . (Cf. F i+s d column in Table 4). Table 4 sums up the performance of the four above mentioned use cases. To measure the performances on the 200-task instances, no time limit is fixed. The entries of Table 4 4 to offer an overview.

Table 4 shows that the number of nodes is lowered by 37.0% while the computation time is only lowered by 10.8% on average for n = 200. In addition, for the six 200-task instances where F i+s d reaches the time limit, only less than 100 nodes are explored. The limit for solving F i+s l is thus the size or the difficulty of the LPs solved at each node, rather than the number of nodes.

Trying to address this issue, we implemented a separation algorithm for the insert and swap inequalities using a callback function. The time needed to solve 50-task instances using this separation algorithm and Cplex features was 1513 seconds with 925 nodes in average. We observe that 98% of the computation time is used by the UserCut Callback to add 71 inequalities in average. This is not surprising since the separation algorithm consists in simply evaluating the terms of inequality ( 5) and ( 6) for the n possible tasks u, and the terms of inequality (7) for the n 2 possible couples (u, v), which results in an O(n 3 ) procedure.

Providing a faster separation algorithm could reduce the computation time, but the branching scheme, and then the number of nodes, would be the same. Since this number of nodes is quite large compared to the performance of F i+s d (which solves all 50-task instances at the root node), we conclude that adding dominance inequalities through a separation procedure reduces their impact. Indeed, when initially added, the dominance inequalities allow to the Cplex presolve phase to fix some variables to 0 or 1. The number of LPs variables is then reduced and the value obtained at each node is improved. When the δ variables are set as continuous variables, this presolve is not executed. It is then consistent with the observation that adding dominance inequalities in this latter case has no impact (Cf. Section 4.2).

Conclusion

In this work, we propose a new way to use neighborhood-based dominance properties, which results in a new kind of reinforcement inequalities. In contrast with the commonly used reinforcement inequalities, which cut fractional points, these inequalities cut non locally optimal solutions. In particular, for the compact formulation F 2 , we provide linear inequalities cutting all the solutions which are not insert and swap locally optimal.

From a practical point of view, we show that adding insert and swap inequalities greatly improves performances of F 2 . Indeed, instead of 50-task instances, we can now solve up to 150-task instances to optimality within one hour.

Insert and swap inequalities can also be used to provide a heuristic solution which is slightly better than the one proposed by Biskup & Feldmann (2001a). For instances up to size 200, this heuristic solution is obtained in less than 20 seconds. A lower bound providing a 15% gap can also be obtained in less than 420 seconds.

We observe that insert and swap inequalities do not improve the linear relaxation value of the compact formulation F 2 . However, used in conjunction with Cplex features, they allow to improve the lower bound obtained at the root node. Two issues follow. Firstly, for a version of F 2 reinforced by cuts or by branching decisions, do dominance properties improve the linear relaxation value? Secondly, which procedure implemented in the Cplex features take advantage of the insert and swap inequalities? Addressing these issues requires an appropriate experimental framework.

Moreover, this work could be extended to other problems where the solutions can be encoded by partitions (any kind of partition, not necessarily ordered bipartitions). For instance, inequalities similar to insert and swap inequalities could be used in a generalization of UCDDP to a parallel machine framework. Indeed, if the tasks share a common due date,the dominant schedules can be encoded by ordered 2m-partitions, where m is the number of machines. This is true even if the common due date, the processing times and the unit earliness and tardiness penalties depend on the machine. Beyond the scheduling field, such inequalities could also be used in the maximum cut problem [START_REF] Karp | Reducibility among combinatorial problems[END_REF] or in a maximum k-cut problem [START_REF] Frieze | Improved approximation algorithms for max kcut and max bisection[END_REF].

For other combinatorial problems where solutions do not have a partition structure; some neighborhood-based dominance inequalities could also be designed using appropriate operations.

Figure 2 :

 2 Figure 2: Illustration of the operation-centered point of view

Figure 3 :

 3 Figure 3: Illustration of the insert operation of an early task u on the tardy side

Figure 4 :

 4 Figure 4: Early side variation induced by swapping u and v when α v p v < α u p u

Figure 5 :

 5 Figure 5: Early side variation induced by swapping u and v when α v p v

F

  are the following.

	L-gap : the average optimality gap of the provided lower bound
	U-gap : the average optimality gap of the provided upper bound
	time : the average running time in seconds
	#nd : the average number of nodes except the root node
	New experiments are conducted for the results reported in F i+s d -5% and F i+s

d columns when n = 200. These results are gathered with the previously obtained results in Table

Table 2 :

 2 Effect of adding insert and swap inequalities on lower bounds

		2 -lp	F 2 l -rn	F i l -rn	F s l -rn	F i+s l -rn	F 2 d -rn	F i d -rn	F s d -rn	F i+s d -rn
	n	L-gap time	L-gap time	L-gap time	L-gap time	L-gap time	L-gap time	L-gap time	L-gap time	L-gap time
	10	41% 0	41% 0	33% 0	41% 0	0%	0	7%	1	5%	1	0%	1	0%	0
	20	68% 0	68% 0	66% 0	68% 0	12% 1	28% 2	28% 2	6%	1	2%	0
	50	86% 0	86% 1	86% 1	86% 6	28% 6	42% 27	41% 31	17% 5	11% 3
	60	89% 0	89% 1	89% 1	89% 7	36% 7	41% 91	41% 95	22% 9	16% 5
	80	92% 1	92% 1	92% 1	92% 11	34% 8	43% 345	43% 359	21% 28	15% 10
	100	93% 2	93% 2	93% 2	93% 19	35% 14	45% 1091	44% 1152	21% 62	14% 25
	120	94% 3	94% 4	94% 11	94% 31	38% 15	46% 3189	46% 3192	22% 133	16% 52
	150	96% 6	96% 13	96% 15	96% 60	34% 29	-	-	-	-	22% 352	15% 130
	180	96% 12	96% 19	96% 23	96% 98	34% 49	-	-	-	-	22% 766	15% 274
	200	97% 19	97% 25	97% 31	97% 126	39% 72	-	-	-	-	22% 1204	15% 418
	500	99%* 722	99%* 698	99%* 742	99%* 2820	36%* 1870	-	-	-	-	-	-	-	-

Table 3 :

 3 Comparison of different heuristics providing an upper bound

	to obtain:	an upper bound		a lower bound		a 5%-approximation	an exact solution
	use:		R1+			R2+		F i+s d -5%	F i+s d
	n	L-gap U-gap time	L-gap U-gap time	time	#nd	time	#nd
	50	86% 0.00% <1	11% 0.00% 3	8	0	4	34
	100	93% 0.00%	2	14% 0.00% 25	160	114	165	141
	200	97% 0.01% 20	15% 0.01% 418	7420	928	8317 1474
	500	-	-(99%) 778	-	-	-	-	-	-	-

Table 4 :

 4 Different ways of using insert and swap inequalities

Since we have ∆ u (δ) -M u for any δ ∈ {0, 1} J , the following inequality is satisfied by every δ ∈ {0, 1} J such that (1-δ u ) = 1.

Conversely, for every δ ∈ {0, 1} J such that (1-δ u ) = 0, inequality ( 5) is satisfied if and only if ∆ u (δ) 0. Considering these two cases, δ ∈ {0, 1} J satisfies (5) if and only if the partition encoded by δ satisfies (I u ). Property 7.(i) follows.

Similarly, to translate constraint (I u ), we introduce the following constant which is an upper bound of -∆ u (δ) for any δ ∈ {0, 1} J .

Conversely, if δ u = 0, δ satisfies (6) if and only if -∆ u (δ) 0. Considering these two cases, δ ∈ {0, 1} J satisfies (6) if and only if if and only if the partition encoded by δ satisfies (I u ). Property 7.(ii) follows.

Property 7

Let δ ∈ {0, 1} J and let (E, T ) be the partition encoded by δ. For any u ∈ J, (i) δ satisfies inequality (5) if and only if (E, T ) satisfies constraint (I u ).

(ii) δ satisfies inequality (6) if and only if (E, T ) satisfies constraint (I u ).

• Linear inequalities translating constraint (S u,v ) for any (u, v) ∈ J 2

Let (u, v) ∈ J 2 such that u v. If δ ∈ {0, 1} J encodes a partition (E, T ), the penalty variation ∆ u,v (E, T ) can be expressed linearly from δ as follows.

To unify these cases into one inequality, we introduce the following constant M u,v which is an upper bound on -∆ u,v (δ) for any δ ∈ {0, 1} J .

In particular, this inequality is satisfied for every δ ∈ {0, 1} J such that (δ u , δ v ) = (0, 0) or (1, 1). For

By construction, (δ, X) satisfies inequalities (1-4) (Cf. Lemma 6). It is thus a solution of F 2 , and g(δ, X) is an upper bound of the optimal value. However, it is not necessarily a solution for F i , F s and F i+s formulations, since (δ, X) does not necessarily satisfy the insert and swap inequalities.

In order to transform such a solution into a solution satisfying the dominance inequalities, we can iteratively apply the operation associated to each violated dominance inequality, until all of them are satisfied. Algorithm 1 presents a way to implement this procedure that we call Insert_swap_improvement. From Property 9, if an insert (resp. a swap) inequality is not satisfied, applying the appropriate insert (resp. swap) operation provides a strictly better solution. Therefore, each solution is considered at most once in this procedure. Since the number of solutions is finte, the Insert_swap_improvement procedure finishes.

The returned solution is an insert and swap local optimum, since it satisfies all dominance inequalities (Cf. Corrolary 10).

Insert_swap_improvement

input: δ ∈ {0, 1} J output: δ encoding an insert and swap local optimum


Fernando Casas
email: fernando.casas@uji.es

Philippe Chartier
email: philippe.chartier@inria.fr

Alejandro Escorihuela-Tomàs

Yong Zhang

May

Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations

Keywords: Composition methods, projection on the real-axis, pseudo-symmetry, pseudo-symplecticity

In this paper, we are concerned with the construction and analysis of a new class of methods obtained as double jump compositions with complex coefficients and projection on the real axis. It is shown in particular that the new integrators are symmetric and symplectic up to high orders if one uses a symmetric and symplectic basic method. In terms of efficiency, the aforementioned technique requires fewer stages than standard compositions of the same orders and is thus expected to lead to faster methods.

Introduction

Given a differential equation

̇ ≡ = ( ), (0) = 0 , (1)
composition methods constitute a powerful technique to raise the order of a given integrator applied to (1) with time-step , as high as might be required, by considering expressions of the form

= 1 • 2 • ⋯ • , ( 2
)
where the coefficients 1 , 2 , … , are appropriately chosen so as to satisfy some universal algebraic conditions [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Murua | Order conditions for numerical integrators obtained by composing simpler integrators[END_REF][START_REF] Chartier | An algebraic theory of order[END_REF]. It is known in particular that if is of order , i.e. satisfies

( 0 ) -( 0 ) = ( +1 ),
where denotes the exact flow of (1), then will be at least of order + 1 (i.e., local error + 2) if the following two conditions are satisfied

( ) ∑ =1 = 1 and ( ) ∑ =1 +1 = 0. (3)
Given that these two equations have no real solution for odd and arbitrary , a series of authors (e.g. [START_REF] Suzuki | Fractal decomposition of exponential operators with applications to manybody theories and Monte Carlo simulations[END_REF][START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]) suggested to start from a second-order method and to consider symmetric compositions only, i.e., schemes with coefficients satisfying the additional condition

+1-= , = 1, … , .
This has led to so-called triple-jump compositions ( = 3, 3 = 1 ) obtained by iterating the process described above to construct a sequence of symmetric methods with even orders (see, e.g., [START_REF] Hairer | Geometric Numerical Integration[END_REF] pp. 44).

In spite of its simplicity, the triple-jump rationale leads to inefficiencies for high orders as compared to methods obtained by solving directly the order conditions [START_REF] Hairer | Geometric Numerical Integration[END_REF]. On top of this, it also suffers from the occurrence of negative time-steps, although this fact is not specific to triplejump methods and concerns all composition or splitting methods of orders higher than two. This, of course, is a severe limiting factor for equations where the vector field (usually an operator) is not reversible, the prototypical example of which being the heat equation. To circumvent this difficulty, several authors have suggested to use complex time-steps (or complex coefficients) in the context of parabolic equations [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF][START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF]. One indeed easily sees that, already for = 2, solutions of equations ( ) -( ) exist in ℂ.

Generally speaking, suppose that is an integrator of order , denoted  [ ] in the sequel for clarity, and consider the composition (2) with = 2,

 [ +1] =  [ ] 1 • [ ] 2 .
(4)

Then, if the coefficients verify conditions ( ) -( ), that is to say if

1 = ̄ 2 ≡ = 1 2 + 2 sin( 2 +1 +1 ) 1 + cos( 2 +1 +1 ) for -2 ≤ ≤ 2 -1 if is even -+1 2 ≤ ≤ -1 2 if is odd , (5)
then (4) results in a method of order + 1, which can subsequently be used to generate recursively higher order composition schemes by applying the same procedure. The choice = 0,

= [ ] ∶= 1 2 + 2 sin +1 1 + cos +1 = 1 2 + 2 tan 2( + 1) = 1 2 cos 2( +1) 2( +1) (6)
gives the solutions with the smallest phase. If we start with a symmetric method of order 2,  [2] , and apply composition (4) with corresponding coefficients (6), we can construct the following sequence of methods:

 [2] ⟶  [3] ⟶  [4] ⟶  [5] ⟶  [6] ,
all of which have coefficients with positive real part [START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF]. The final method of order 6 involves 16 evaluations of the basic scheme  [2] . By contrast, there are composition methods of order 6 (both with real and complex coefficients) involving just 7 evaluations of [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]. It is thus apparent that this direct approach does not lead to cost-efficient high-order schemes. One should remark that the composition (4) does not provide a time-symmetric method, i.e.,  [ +1] -• [ +1] is not the identity map, even if  [ ] happens to be symmetric. As we have mentioned before, symmetry allows to raise the order by two at each iteration by considering the triple-jump composition

 [2] [BCCM13,
 [2 +2] =  [2 ] 1 • [2 ] 2 • [2 ] 1 (7)
starting from a symmetric method. Apart from the real solution, the complex one with the smallest phase is

1 = e ∕( +1) 2 1∕( +1) -2 e ∕( +1) , 2 = 1 -2 1 , ( 8
)
and symmetric methods up to order 8 with coefficients having positive real part are possible if one starts with a symmetric second-order scheme 1 . These order barriers has been rigorously proved in [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF].

The simple third-order scheme (4) corresponding to = 2 has been in fact rediscovered several times in the literature [START_REF] Bandrauk | Improved exponential split operator method for solving the time-dependent Schrödinger equation[END_REF][START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF][START_REF] Chambers | Symplectic integrators with complex time steps[END_REF][START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF][START_REF] Suzuki | General theory of higher-order decomposition of exponential operators and symplectic integrators[END_REF]. In particular, it was shown in [START_REF] Chambers | Symplectic integrators with complex time steps[END_REF] that the method, when applied to the two-body Kepler problem, behaves indeed as a fourthorder integrator, the reason being attributed to the fact that the main error term in the asymptotic expansion is purely imaginary. In this note we elaborate further the analysis and provide a comprehensive study of the general composition (4), paying special attention to the qualitative properties the method shares with the continuous system (1). In addition, we show how it is possible to combine compositions and a trivial linear combination to raise the order, while still preserving the qualitative properties of the basic integrator up to an order higher than of the method itself.

Composition and pseudo-symmetry or pseudo-symplecticity

In what follows, we will assume for convenience that all values of in (1) lie in a compact set where the function is smooth. Before starting the analysis, it is worth recalling the notions of adjoint method and symplectic flow.

The adjoint method * of a given method is the inverse map of the original integrator with reversed time step -, i.e., * ∶= -1 -. A symmetric method satisfies * = [Cha15, HLW06]. The vector field in (1) is Hamiltonian if there exists a function ( ) such that = ∇ ( ), where = ( , ) and is the basic canonical matrix. Then, the exact flow of (1) is a symplectic transformation, ′ ( ) ′ ( ) = for ≥ 0 [BC16, SSC94]. It then makes sense introducing the following definitions, taken from [START_REF] Chartier | Reversible B-series[END_REF] and [START_REF] Aubry | Pseudo-symplectic Runge-Kutta methods[END_REF]: Definition 1 Let be a smooth and consistent integrator:

1. it is pseudo-symmetric of pseudo-symmetry order if for all sufficiently small , the following relation holds true:

* = + ( +1 ), (9)
where the constant in the -term depends on bounds of derivatives of on .

2. it is pseudo-symplectic of pseudo-symplecticity order if for all sufficiently small , the following relation holds true when is applied to a Hamiltonian system:

( ′ ) ′ = + ( +1 ), ( 10
)
where the constant in the -term depends on bounds of derivatives of on .

Remark 1 A symmetric method is pseudo-symmetric of any order ∈ ℕ, whereas a method of order is pseudo-symmetric of order ≥ . A similar statement holds for symplectic methods.

As a first illustration of Definition 1, let us consider again a symmetric 2nd-order method  [2] and form the composition

[3] =  [2] • [2] ̄ with = 1 2 + √ 3
6 . Then, if the vector field under consideration is real-valued, its real part

ℜ( [3] ) = 1 2 [3] + [3] = 1 2  [2] • [2] ̄ +  [2] ̄ • [2] .
1 It is actually possible to reach order 14 if, in the construction, one uses formula (7) alternatively with coefficients

1 , 2 and coefficients ̄ 1 , ̄ 2 [CCDV09].
is a method of order 4 and pseudo-symmetric of pseudo-symmetry order 7. This result is a consequence of the fact that ( [3] ) * = [3] and the following general statement, which lies at the core of the construction procedure described in this paper.

Proposition 1 Let be any consistent smooth method for equation (1) and consider the new method = 1 2 + * .

Assume also that is pseudo-symmetric of order . Then is of pseudo-symmetry order 2 + 1. If is furthermore of pseudo-symplecticity order , then is of pseudo-symplecticity order min(2 + 1, ).

Proof: By assumption, there exists a smooth function ( , ) ↦ ( ), defined for all in a compact set and for all sufficiently small real , such that * = + +1 or

-1 -= + +1 or -1 = -+ (-) +1 -, ( 11
) so that = + 1 2 +1 .
Composing the third relation of (11) from the left by , we obtain

id = • -+ (-) +1 ′ • -⋅ -+ ( 2( +1) ), (12)
where the -term depends on bounds of the derivatives of and on . Similarly, composing the second relation of ( 11) from the right by -, we get

id = • -+ +1 • -. ( 13
)
As a consequence, we have

+1 • -= (-) +1 ′ • -⋅ -+ ( 2( +1) ).
We are then in position to write

• -= + 1 2 +1 • -+ 1 2 (-) +1 - = • -+ 1 2 (-) +1 ′ • -⋅ -+ 1 2 +1 • -+ ( 2( +1) ) = id + ( 2( +1) ),
which proves the first statement. Now, if is in addition of pseudo-symplecticity order , then its adjoint * is also of pseudo-symplecticity order , so that relation (11) leads to

+ ( +1 ) = ( * ) * = ( ′ + +1 ′ ) ( ′ + +1 ′ ) = + ( +1 ) + +1 ( ′ ) ′ + ( ′ ) ′ + ( 2( +1) ),
which implies that

+1 ( ′ ) ′ + ( ′ ) ′ = ( min(2( +1), +1) ).
As an immediate consequence, we have that

( ′ ) ′ = + ( min(2( +1), +1) )
which proves the second statement. □

This result can be rendered more specific as follows:

Proposition 2 Let  [2 ] be a smooth method of order 2 ≥ 2 and pseudo-symmetry order ≥ 2 +1.

Let us consider the composition method

[2 +1] =  [2 ] 1 •  [2 ] 2 , ( 14
)
where the coefficients 1 and 2 satisfy both relations 1 + 2 = 1 and 2 +1 1

+ 2 +1 2 = 0. Then the method ̂ = 1 2 [2 +1] + [2 +1] ( 15
)
is of order

2 + 1 if = 2 + 1, 2 + 2 if ≥ 2 + 2 (16)
when the vector field in ( 1) is real, and of pseudo-symmetry order

2 + 1 if = 2 + 1, min( , 4 + 3) if ≥ 2 + 2. ( 17
)
If in addition, is a (real) Hamiltonian vector field and  [2 ] is of pseudo-symplecticity order , then ̂ is of pseudo-symplecticity order

min( , 2 + 1) if = 2 + 1, min( , , 4 + 3) if ≥ 2 + 2. ( 18
)
Remark 2 Note that in Proposition 2, one has necessarily ≥ 2 + 1. This can be seen straightforwardly by a direct computation of

 [2 ] -• [2 ] ( ) with  [2 ] ( ) = ( ) + 2 +1 ( ) + ( 2 +2
).

Proof: Noticing that 1 and 2 are complex conjugate and (1) is real, and taking into account that  [2 ] is of pseudo-symmetry order , we have

[2 +1] =  [2 ] 2 •  [2 ] 1 = ( [2 ] 2 ) * + ( +1 ) • ( [2 ] 1 ) * + ( +1 ) = ( [2 ] 2 ) * • ( [2 ] 1 ) * + ( +1 ) = ( [2 +1] ) * + ( +1 ). (19)
Moreover, by construction, [2 +1] is at least of order 2 + 1, so that

[2 +1] + ( 2 +2 ) = [2 +1] = ( [2 +1] ) * + ( 2 +2 ), (20)
and altogether

̂ = + ( max(2 +2, +1) ). (21)
Now, since the pseudo-symmetry order of [2 +1] is at least 2 + 1, the method

= 1 2 [2 +1] + ( [2 +1] ) *
is, according to Proposition 1, of pseudo-symmetry order 4 + 3 and of pseudo-symplecticity order min(4 + 3, ). The first (16), second (17) and third (18) statements on orders then follow from (21). □

In the Appendix we provide an alternative proof of Proposition 2 based on the Lie formalism, which allow us, in addition, to generalise the previous result on pseudo-symplecticity to other geometric properties the continuous system may possess (such as in volume preserving flows, isospectral flows, differential equations evolving on Lie groups, etc.).

Notice that, according with Proposition 2, if we start from = 1, that is to say from a basic symmetric ( = +∞) and symplectic ( = +∞) method of order 2, we get a method of order 4 that is pseudo-symmetric and pseudo-symplectic of order 7 just by considering the simple composition (14) and taking the real part of the output at each time step. If this technique is applied to a symmetric and symplectic method of order 4, i.e. with = 2, then ̂ is of order 6 and pseudo-symmetric and pseudo-symplectic of order 11.

Let us consider, in particular, the 4 ℎ -order symmetric scheme (7) with = 2 as basic scheme. Then, the resulting 6 ℎ -order integrator ̂ only requires the evaluation of 6 second-order methods  [2] , whereas the corresponding 6 ℎ -order scheme obtained by the triple-jump technique involves 9 evaluations. This number is reduced to 7 by considering general compositions of  [2] [BCCM13]. If we take this 6 ℎ -order composition of 7 schemes as the basic method  [6] , the resulting integrator of order 8, ̂ , involves the evaluation of 14  [2] , whereas 15 evaluations are required by pure composition methods. Notice that ̂ is pseudo-symmetric and pseudo-symplectic of order 15, so that for values of sufficiently small, it preserves effectively the symmetry up to round-off error while the drift in energy for Hamiltonian systems is hardly noticeable.

Families of pseudo-symplectic methods

There is another possibility to increase the order, though, and it consists in applying the technique of Proposition 2 recursively. Thus, if denote by ̂ (1) ≡ ̂ the method of eq. ( 15), we propose to apply the following recurrence:

For = 2, 3, … Φ ( ) = ̂ ( -1) [2 ] • ̂ ( -1) ̄ [2 ] ̂ ( ) = 1 2 Φ ( ) + Φ ( ) (22)
where [2 ] is given by (6). Then, according with Proposition 2, it is possible to raise the order up to the pseudo-symmetry order of the underlying basic method  [2 ] . Thus, in particular, the maximum order one can achieve by applying this technique to the basic symmetric method  [2] is 7, whereas if we start with a basic symmetric method of order 4,  [4] , the maximum order is 11. It is 15 from a symmetric method  [6] of order 6 and so on and so forth.

To give an assessment of the computational cost of the methods obtained by applying this type of composition, we notice that the computation of Φ ( ) and Φ ( ) required to form ̂ ( ) by ( 22) at the intermediate stages can be done in parallel, whereas at the final stage it only requires taking the real part. Thus, the method of order 6 constructed recursively from  [2] only requires the effective computation of 4 basic methods  [2] .

Starting from a symmetric second-order method  [2] , say Strang splitting for instance, it is important to monitor the sign of the real part of all coefficients involved in the previous iteration. It is immediate to see that in the recursive construction

 [2] → ̂ (1) → ̂ (2) → ̂ (3)
envisaged in the recurrence (22), the basic method  [2] is used with the following coefficients

= 1 ∶ [2] , ̄ [2] = 2 ∶ [4] [2] , ̄ [4] [2] , [4] ̄ [2] , ̄ [4] ̄ [2] = 3 ∶ [6] [4] [2] , [6] ̄ [4] [2] , [6] [4] ̄ [2] , [6] ̄ [4] ̄ [2] , ̄ [6] [4] [2] , ̄ [6] ̄ [4] [2] , ̄ [6] [4] ̄ [2] , ̄ [6] ̄ [4] ̄ [2]
Given the expression of [ ] (see (6)), these coefficients have arguments of the form

2 ∑ =1 ± 1 2 + 1 = 2 ± 1 3 ± 1 5 ± ⋯ ± 1 2 + 1 , = 1, 2, 3,
so that their maximum argument is

2 3 ∑ =1 1 2 + 1 .
For all the coefficients to have positive real parts, a necessary and sufficient condition is thus that

3 ∑ =1 1 2 + 1 ≤ 1.
It clearly holds for methods ̂ (1) , ̂ (2) and ̂ (3) , of respective orders 4, 6 and 7, since 1∕3+1∕5+1∕7 = 71∕105. Similarly, starting form a symmetric method of order 4 having real or complex coefficients with maximum argument 4 , the condition becomes

2 4 + 4 ∑ =1 1 2 + 3 = 2 4 + 1888 3465 ≤ 1.
For instance, suppose that in (1) can be split as ( ) = ( ) + ( ), so that the exact -flows [ ] and [ ] corresponding to and , respectively, can be computed exactly. Then, the following composition

 [4] = [ ] 1 • [ ] 1 • [ ] 2 • [ ] 2 • [ ] 3 • [ ] 2 • [ ] 2 • [ ] 1 • [ ] 1 (23)
with

1 = 1 10 - 1 30 , 2 = 4 15 + 2 15 , 3 = 4 15 - 1 5 and 1 = 2 = 3 = 4 = 1 4
provides a 4th-order symmetric scheme (see [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF]). Taking (23) as basic method we get max =1,2,3 Arg( ) = arccos (4∕5) so that 2 4 + 1888 3465 < 0.409666 + 1888 3465 < 0.96 < 1 and thus all methods

 [4] → ̂ (1) → ̂ (2) → ̂ (3) → ̂ (4)
of respective orders 4, 6, 8, 10 and 11 obtained by the procedure ( 22) have all their coefficients with positive real parts. As far as the part is concerned, the maximum argument is less than 0.55 2 .

Numerical experiments

In this section we illustrate the previous results on several numerical examples, comprising Hamiltonian systems and partial differential equations of evolution previously discretised in space.

Harmonic oscillator

We first consider the simple harmonic oscillator, with Hamiltonian

= ( ) + ( ) = 1 2 2 + 1 2 2 .
If we denote by ( ) the exact matrix evolution associated with the Hamiltonians = , and , i.e., ( ( ), ( )) = ( )( (0), (0)) , then Table 1: Main term in the truncation error, degree of symmetry and symplecticity for schemes ̂ ( ) obtained from the basic leapfrog integrator for the simple harmonic oscillator. 2 stands for the 2 × 2 identity matrix.

( ) = cos( ) sin( ) -sin( ) cos( ) , ( ) = 1 0 1 , ( ) = 1 0 -1 , ( ) - • --2 det ̂ (1) 0 -1
respectively. We take as basic symmetric (and symplectic) scheme the leapfrog/Strang splitting:

 [2] = ( ∕2) ( ) ( ∕2) (24)
and compute the first three iterations in (22). In Table 1 we collect the main term in the truncation error for the resulting integrators ̂ ( ) , = 1, 2, 3. We also check their time-symmetry and the preservation of the symplectic character of the approximate solution matrix by computing its determinant (a 2 × 2 matrix is symplectic iff det( ) = 1). One can observe that these results are in agreement with the previous estimates.

Next we take initial conditions ( , ) = (2.5, 0), integrate until the final time = 10 4 with  [2] , ̂ (1) , and ̂ (2) and compute the relative error in energy along the evolution. The result is depicted in Figure 1. We see that for ̂ (1) and ̂ (2) the error in energy is almost constant for a certain period of time, and then there is a secular growth proportional to ( 7 ).

Kepler problem

Next, we consider the two-dimensional Kepler problem with Hamiltonian

( , ) = ( ) + ( ) = 1 2 - 1 . ( 25
)
Here = ( 1 , 2 ), = ( 1 , 2 ), = , is the gravitational constant and is the sum of the masses of the two bodies. Taking = 1 and initial conditions

1 (0) = 1 -, 2 (0) = 0, 1 (0) = 0, 2 (0) = √ 1 + 1 - , ( 26
)
if 0 ≤ < 1, then the solution is periodic with period 2 , and the trajectory is an ellipse of eccentricity . Note that the gradient function must here be implemented carefully so as to be analytic for complex values of = 2 1 + 2 2 . Here, we define it using the following determination of the complex logarithm (analytic on the complex plane outside the negative real axis):

∀ ( , ) ∈ ℝ 2 s.t. + ∉ ℝ -, ( + ) = log | + | + 2 arctan + | + | .
As a consequence, the analytic continuation of the function

1∕ 3 = 1∕( 2 1 + 2 2 ) 3∕2 writes exp - 3 2 ( + ) , where = ℜ( 2 1 + 2 2 ) and = ℑ( 2 1 + 2 2 ).
Here, as with the harmonic oscillator, we take as basic method the 2nd-order Strang splitting

 [2] = [ ] ∕2 • [ ] • [ ] ∕2 , ( 27
)
where [ ] (respectively, [ ] ) corresponds to the exact solution obtained by integrating the kinetic energy ( ) (resp., potential energy ( )) in (25). We take = 0.6, integrate until the final time = 20 with Strang and the schemes obtained by the recursion (22) with = 1, 2, 3 for several time steps and compute the relative error in energy at the final time. Figure 2 show this error as a function of the inverse of the step size 1∕ to illustrate the order of convergence: order 2 for Strang, order 4 for ̂ (1) and order 6 for ̂ (2) . For ̂ (3) , and contrary to what happens to the harmonic oscillator, the observed numerical order is higher than expected, varying between 7 and 8. We do not have at present a theoretical explanation for this phenomenon. Figure 2 (right) depicts the time evolution of this error when the final time is = 10 4 .

The semi-linear reaction-diffusion equation of Fisher

Our third test-problem is the scalar equation in one-dimension

( , ) = Δ ( , ) + ( ( , )), (28)
with periodic boundary conditions on the interval [0, 1]. Here ( ) is a nonlinear reaction term.

For the purpose of testing our methods, we take Fisher's potential [START_REF] Sari | Encyclopedia of Applied and Computational Mathematics, chapter Fisher's equation[END_REF] ( ) = (1 -)

as considered for example in [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF]. The splitting corresponds here to solving, on the one hand, the linear equation with the Laplacian (as ), and on the other hand, the non-linear ordinary differential equation

( , ) = ( , )(1 -( , )),
with initial condition ( , 0) = 0 ( ), whose analytical solution is given by the well-defined (for small enough complex time ) formula ( , ) = 0 ( ) + 0 ( )(1 -0 ( )) (e -1) 1 + 0 ( ) (e -1) .

Here we aim to solve Eq. ( 28) with periodic boundary conditions on the interval [0, 1], and initial condition 0 ( ) = sin(2 ). Numerically, the interval is discretised on a uniform grid, i.e., = ∕ , = 0, … , -1, ∈ ℕ, and ( , ) is approximated by Fourier pseudo-spectral methods. In this way we construct a vector with components ( ) ≈ ( -1 , ), = 1, 2, … , . If we denote by the numerical solution computed by a certain integrator with step size , and by ∕2 the corresponding numerical solution computed by the same integrator with step size ∕2, then the quantity ∶= ‖ -∕2 ‖ ∞ is a good indicator of the convergence order. Numerical simulations were carried out in quadruple precision (with Intel Fortran) such that roundoff errors are suppressed. Figure 3 shows the successive errors , at final time = 10, of the methods obtained with the sequence (22) with the Strang splitting as the basic method  [2] (left) and the fourth order scheme  [4] given by (23) (right) with different time steps = 0.1∕2 , = 1, … , 7. One can clearly observe that the convergence order matches the previous analysis with a slightly better performance for the highest order, analogously to the Kepler problem. Figure 4 shows the successive errors versus the number of basic integrators in each case.

The semi-linear complex Ginzburg-Landau equation

Our final test problem is the complex Ginzburg-Landau equation on the domain ( , ) ∈ [-100, 100]× [0, 100],

( , ) = Δ ( , ) + ( , ) -| ( , )| 2 ( , ), (29)
with = 1 + 1 , = 1 -3 and initial condition ( , 0) = 0 ( ). Here, , 1 and 3 denote real coefficients. In physics, the Ginzburg-Landau appears in the mathematical theory used to model superconductivity. For a broad introduction to the rich dynamics of this equation, we refer to [vS95]. Here, we will use the values 1 = 1, 3 = -2 and = 1, for which plane wave solutions establish themselves quickly after a transient phase (see [START_REF] Winterbottom | Oscillatory pattern formation with a conserved quantity[END_REF]). In addition, we set 0 ( ) = 0.8 cosh( -10) 2 + 0.8 cosh( + 10) 2 , so that the solution can be represented in Figure 5. To apply the composition methods presented in previous sections, it seems natural to split equation (29) as

( , ) = (1 + 1 )Δ ( , ) + ( , ), (30)
whose solution is ( , ) = e e (1+ 1 )Δ 0 ( ) for ≥ 0, and

( , ) = -(1 -3 )| ( , )| 2 ( , ) (31)
with solution is for ≥ 0 ( , ) = e -(1-3 ) ∫ 0 ( , ) 0 ( ) = e -2 log(1+2 0 ( )) 0 ( ).

Here we have first solved the equation for ( , ) ∶= | ( , )| 2 , given by

( , ) = -2 2 ( , ),
with solution ( , ) = 0 ( ) 1 + 2 0 ( ) .

Considering now as a complex variable with positive real part does not raise any difficulty for the first part, since e e (1+ 1 )Δ is well-defined. More care has to be taken for the second part, since ↦ | | 2 is not a holomorphic function, and this prevents us from solving (29) in its current form.

As a consequence, we first rewrite (29) as a system for ( ( , ), ( , )) where ( , ) = ℜ( ( , )) and ( , ) = ℑ( ( , )):

⎧ ⎪ ⎨ ⎪ ⎩ ( , ) = Δ ( , ) -1 Δ ( , ) + ( , ) -( 2 ( , ) + 2 ( , ))( ( , ) + 3 ( , )) ( , ) = 1 Δ ( , ) + Δ ( , ) + ( , ) -( 2 ( , ) + 2 ( , ))(-3 ( , ) + ( , )) (32)
and now solve it for complex time ∈ ℂ with ℜ( ) ≥ 0. Observing that

-1 -3 3 -1 = 3 -1 and 1 -1 1 1 = 1 -1 , with 1 = 0 0 ̄ , 3 = - 0 0 -̄ , = 1 1 and -1 = -2 1 2 1 2 -2
, system (32) can be rewritten as

⎧ ⎪ ⎨ ⎪ ⎩ ̃ ( , ) = Δ ̃ ( , ) + ̃ ( , ) - ̃ ( , ) ̃ ( , ) ̃ ( , ) = ̄ Δ ̃ ( , ) + ̃ ( , ) -̄ ̃ ( , ) ̃ ( , ) (33)
where ̃ ( , ) = 4 ̃ ( , ) ̃ ( , ) and where

̃ ̃ = 1 2 -1 1 - .
It is not difficult to see that the exact solution of the second part of (33) is given by

̃ ( , ) = ̃ 0 ( )e -2 log(1+2 ̃ 0 ( )) ̃ ( , ) = ̃ 0 ( )e -̄ 2 log(1+2 ̃ 0 ( )) (34)
where ̃ 0 ( ) is now defined as ̃ 0 ( ) ∶= 4 ̃ 0 ( ) ̃ 0 ( ). Note that here, by convention, the logarithm refers to the principal value of log( ) for complex numbers: if = ( + ) = e with -< ≤ , then

log ∶= ln + = ln | | + arg = ln(| + |) + 2 arctan + √ 2 + 2 .
Since log( ) is not defined for ∈ ℝ -, this means that the solution ( ̃ ( , ), ̃ ( , )) is defined only as long as 1 + 2 ̃ 0 ( ) ∉ ℝ -. Finally, the solution ( ( , ), ( , )) is of the form

( , ) = 0 ( ) (e -( , ) +e -̄ ( , ) ) 2 -0 ( ) (e -( , ) -e -̄ ( , ) ) 2 ( , ) = 0 ( ) (e -( , ) -e -̄ ( , ) ) 2 + 0 ( ) (e -( , ) +e -̄ ( , ) ) 2
where ( , ) ∶= log(1 + 2 ̃ 0 ( )) = log(1 + 2 0 ( )) with 0 ( ) = 2 0 ( ) + 2 0 ( ).

Denoting

= ( 1 , … , ) ∈ ℝ and = ( 1 , … , ) ∈ ℝ , we eventually have to numerically solve the following system: where is the diagonal matrix with , = 2 + 2 . Equation ( 29) is solved with periodic boundary conditions on the interval [-100, 100]. Now, in the previous example, the interval is discretised on a uniform grid, i.e., = ∕ , = 0, … , -1, ∈ ℕ with = 512, and ( , ) is approximated by Fourier pseudo-spectral methods. The successive errors ∶= ‖ -∕2 ‖ ∞ are shown also here to confirm the convergence order. Figure 6 shows the successive errors, at final time = 10, of the schemes obtained by applying the sequence (22) from the basic Strang splitting and the fourth-order scheme (23) with = 0.1∕2 , = 1, … , 7. The observed order of convergence matches the previous analysis with a slightly better performance for the highest order. Figure 7 shows the successive errors versus the number of basic integrators. where is the Lie derivative associated with ,

̇ = -1 + -( + 3 ) ̇ = 1 + + -(-3 + )
= ∑ ≥1 ( ) . (35)
Analogously, for a given integrator one can associate a series of linear operators so that For the adjoint integrator * = -1 , one clearly has

( * ( )) = exp -(-) [ ]( ).
This shows that is symmetric if and only if ( ) = 1 + 3 3 + ⋯, and in particular, that symmetric methods are of even order.

An integrator  [2 ] of order 2 ≥ 2 can be associated with the series

Φ( ) = exp + 2 +1 2 +1 + 2 +2 2 +2 + ⋯ (36)
for certain operators . Then, the adjoint method ( [2 ] ) * has the associated series

Φ * ( ) = exp + 2 +1 2 +1 -2 +2 2 +2 + ⋯ .
In consequence,  [2 ] is pseudo-symmetric of order ≥ 2 + 1.

(i) Let us analyse first the case > 2 + 1. Then, 2 +2 = ⋯ = = 0 in (36) and the series of operators associated with the composition [2 +1] =  [2 ] 1 •  [2 ] 2 is Ψ( ) = Φ( 2 ) Φ( 1 ) ≡ exp( ( )), with 0 = 1 , 2 = 3 + 1 24 [ , [ , 1 ]], etc. In general, 2 is a linear combination of the operators { , 2 +1 , 2 +2 , …} and their nested Lie brackets. In addition, 2 ( ) + 2 (-) = ( 4 +4 ), so that we can write  = + 2 +3 + ( 4 +4 ) = exp 2 +3 + ( 4 +4 ) and ( ) = exp 2 exp 2 +3 exp 2 + ( 4 +4 ), whence the following statements follow at once:

• Method (38) is of order 2 + 2, since ( ) = exp( ) + ( 2 +3 ).

• Since ( ) only contains even powers of (up to +1 +1 ), then e 2 e 2 +3 e 2 is a symmetric composition and is pseudo-symmetric of order min( , 4 + 3).

• Let us suppose that scheme (38) is applied to a Hamiltonian system and that  [2 ] is of pseudo-symplecticity order . Since is an operator in the free Lie algebra generated by { , 2 +1 , 2 +2 , …}, clearly the composition e 2 e 2 +3 e 2 is symplectic (at least up to terms ( )). As a matter of fact, this can be extended to any geometric property the differential equation (1) has: volume-preserving, unitary, etc., as long as the basic scheme  [2 ] preserves this property up to order .

Finally, in view of ( 19)-( 21) and recalling that ≥ 2 + 2, the same considerations apply if we take the complex conjugate instead of the adjoint, i.e., to the scheme

̂ = ℜ( [2 +1] ) = 1 2 [2 +1] + [2 +1] . (39)
(ii) We analyse next the case = 2 + 1. Then 2 +2 ≠ 0 in (36) and, if 1 and 2 verify equations (37), then ( ) read ( ) = + 2 +2 0 + ( 2 +3 ) with

0 = ( 2 +2 1 + 2 +2 2 ) 2 +2 + 1 2 ( 2 2 +1 1 -1 2 +1 2 )[ , 2 +1 ].
Notice that, whereas has non-vanishing real and imaginary parts. In any event, the same procedure as in the previous case can be carried out, leading to the conclusion that method (38) is still of order 2 + 2.

The situation is different, however, for method (39), since relations (19)-( 21) do not provide further information. We have to analyse instead In consequence, ̂ is a method of order 2 + 1, pseudo-symmetric of order 2 + 1 and pseudosymplectic of order min( , 2 + 1).

Figure 1 :

 1 Figure 1: Error in energy along the integration for the harmonic oscillator taking (24) as the basic integrator in the sequence (22).

Figure 2 :

 2 Figure 2: Left figure: Relative error in energy vs. the inverse of the step size after approximately 3.183 periods ( = 20) for the Kepler problem for the schemes obtained by the recurrence (22). Right figure: Evolution of this error along the integration.

Figure 3 :

 3 Figure 3: Successive errors versus time step for Eq. (28) of the composition methods starting from the Strang scheme (left) and the fourth order scheme  [4] (right).

Figure 4 :

 4 Figure 4: Successive errors versus number of basic integrators for Eq. (28) of the composition methods starting from the Strang scheme (left) and the fourth order scheme [4] (right).

Figure 5 :

 5 Figure 5: Colormaps of the amplitude | ( , )| 2 (left) and real part ℜ( ( , )) (right) of the solution of (29).

Figure 6 :

 6 Figure 6: Successive errors versus time step for Eq. (29) of the composition methods starting from the Strang scheme (left) and the fourth order scheme  [4] (right).

Figure 7 :

 7 Figure 7: Successive errors versus number of basic integrators for Eq. (29) of the composition methods starting from the Strang scheme (left) and the fourth order scheme [4] (right).

(

  ( )) = exp( ( ))[ ]( ), with ( ) = ∑ ≥1 for all functions [BCM08]. The integrator is of order if 1 = and = 0 for 2 ≤ ≤ .

  ( 2 +3 ). Noting thatℜ exp ( ) = + ℜ ( ) + ( 4 +4 ) = + 2 +2 ℜ( 0 ) + ( 2 +3 )then we can writeℜ(Ψ( )) = exp 2 exp 2 +2 ℜ( 0 ) + ( 2 +3 ) exp 2 + ( 4 +4).

Acknowledgements

The work of the first three authors has been supported by Ministerio de Economía y Competitividad (Spain) through project MTM2016-77660-P (AEI/FEDER, UE). PC acknowledges funding by INRIA through its Sabbatical program and thanks the University of the Basque Country for its hospitality. This work was initiated during YZ's visit at the University of Rennes, IRMAR.

Appendix

In this Appendix we provide an alternative proof of Proposition 2 via Lie formalism. This allows us not only to gain some additional insight into the structure of the methods, but also to generalize the result on pseudo-symplecticity to other properties of geometric character, very often related to Lie groups, the differential equation may possess.

To begin with, if is the exact flow of the equation (1), then for each infinitely differentiable map , the function ( ( )) admits an expansion of the form [Arn89, SSC94]

where ( ) can be formally determined by applying the Baker-Campbell-Hausdorff formula [START_REF] Blanes | A Concise Introduction to Geometric Numerical Integration[END_REF] as

Here [⋅, ⋅] denotes the usual Lie bracket. Clearly, the order of [2 +1] is 2 + 1 if

so that 2 = ̄ 1 ≡ is given by eq. ( 5) (with = 2 ). In that case we can write

whereas for the adjoint method one has

In particular

2 ) 2 +3 . The series Ψ( ) can also be written as

where ( ) is determined by applying the symmetric BCH formula [START_REF] Blanes | A Concise Introduction to Geometric Numerical Integration[END_REF] as

By the same token,

Consider now the method

Clearly, its associated series of operators,

can be expressed as

where  = 1 2 exp ( ( )) + 1 2 exp (-(-)) .

By expanding, we have