Vincent Visseq

Pierre Alart
email: pierre.alart@univ-montp2.fr

David Dureisseix

High performance computing of discrete nonsmooth contact dynamics with domain decomposition

Keywords: discrete elements, multibody, parallel computing, distributed memory, message passing, FETI method, additive Schwarz method, LMGC90, granular material

We investigate two algorithms for solving large-scale granular problems using domain decomposition methods. These numerical schemes can be connected to classical domain decompositions, with and without overlapping, developed in continuum mechanics. The two algorithms are compared in terms of implementation and parallel performance. We focus the numerical study on communication procedures between processors, load balancing and scalability, topics particularly crucial in nonlinear evolutive problems. This is the accepted version of the following article:

Introduction

Numerical simulations of granular media, at the grain scale, leads to several issues: non smoothness of frictional contact behavior (when no regularization is used), multiplicity of solutions (for dense and a not so dynamical regime at macroscopic scale), and time consuming solvers (for large collections of grains and interactions). To address these issues, several answers had been proposed, for instance, several time stepping schemes [START_REF] Alart | Using nonsmooth analysis for numerical simulation of contact mechanics[END_REF][START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF], and parallelization of solvers [START_REF] Herrmann | Simulations of granular materials on different scales[END_REF][START_REF] Renouf | Conjugate gradient type algorithms for frictional multi-contact problems: Applications to granular materials[END_REF][START_REF] Renouf | A parallel version of non smooth contact dynamics algorithm applied to simulation of granular medium[END_REF][START_REF] Shojaaee | An adaptive hierarchical domain decomposition method for parallel contact dynamics simulations of granular materials[END_REF].

Here, we focus on non smooth contact dynamic (NSCD) method for the dynamical modeling and time stepping, and a domain decomposition method (DDM) as a support for distributed memory parallelization with message passing [START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF][START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF]. We discuss algorithms with two different proposals with a distribution of contact interactions among subdomains: a non overlapping FETI-like DDM, defining an interface gluing step thanks to the use of Lagrange multipliers, and a minimal overlapping Schwarz-like DDM with impulse coupling terms. In addition, two different implementations of the underlying communication scheme between processors are proposed. The numerical and parallel performances are compared and assessed on large scale 3D granular media simulations.

Up to our knowledge, the only solution methods coupling a domain decomposition method with non-regular implicit contact dynamics for granular models are available in [START_REF] Breitkopf | Modélisation parallèle des matériaux granulaires[END_REF][START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF][START_REF] Shojaaee | An adaptive hierarchical domain decomposition method for parallel contact dynamics simulations of granular materials[END_REF]. The approach described in [START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF] mainly uses an asynchronous distributed-memory Non-Linear Gauss-Seidel solver (close to the one used in [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF] for a synchronous solver and in [START_REF] Hoang | Domain decomposition with discrete element simulations using shared-memory parallel computing for railways applications[END_REF] for an asynchronous version, both for shared-memory architectures. They correspond to the algebraic partitioning of the reduced dynamics, i.e. a splitting of grains between processors, while herein, we rely on a splitting of interactions between processors). Good parallel efficiencies are obtained in [START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF] with the use of a tuned dynamic load balancing, data migration, and duplication of the whole grain database among processors to minimize information exchanges. It illustrates the fact that there is still a scope for efficiency gains using a load-balancing algorithm. This tendency is confirmed in the similar approach of [START_REF] Shojaaee | An adaptive hierarchical domain decomposition method for parallel contact dynamics simulations of granular materials[END_REF] for 2D problems addressing a large number of processors (up to 256). The distribution of interactions between processors has first been studied in [START_REF] Breitkopf | Modélisation parallèle des matériaux granulaires[END_REF] for a shared-memory architecture, smallsized 2D problems and the same kind of decomposition using a 'box method'; the common point with the Schwarz-NSCD method described in this article is the information exchange between processors, related to the grain global forces. Section 2 briefly recalls the background theory of the nonsmooth contact dynamics (NSCD) in the context of time-stepping schemes. Section 3 focuses on domain decomposition for such a problem, and proposes several alternatives, that are compared on a first test example. Section 4 is related to the parallel performances for a large scale problem on a supercomputing architecture. Finally, Section 5 concludes and discusses outlooks of this work.

Nonsmooth contact dynamics reference problem

With a time-stepping scheme on a dynamical analysis of a granular media, no event detection is performed. Once the solution is known at the beginning t -= t i of a time slab [t i , t i+1 ], whose known quantities are denoted with a superscript (-), the quantities at the end of the time slab (without a superscript) have to be determined.

Grain nonsmooth dynamics

Considering a rigid model, the dynamics of a grain over a time step is expressed with the vector equation [START_REF] Moreau | Numerical aspects of sweeping process[END_REF]:

M V -R = R d (1)
where the prescribed right-hand side is

R d = R D +M V -, V
is the velocity of the grain (it contains the translational degrees of freedom (dof), and the rotational ones), while R is the resultant impulse on the grain due to interactions with other grains and R D are the external prescribed impulses. The matrix M contains both the mass (for the translational dof) and the inertia (for the rotational dof).

A choice leading to a constant, and diagonal, matrix M consists in expressing the global coordinates of rotation vectors in the inertia eigenbasis of each grain, see Figure 1. The assembly of these equations (independent for each grain) for all the involved grains in a granular medium is formally written in the same way as equation ( 1), which then represents the assembly of the individual grain equations.

Contact interactions

For a unilateral contact, Moreau proved via a viability lemma [START_REF] Moreau | Numerical aspects of sweeping process[END_REF], that a velocity-impulse formulation may be equivalent to the classical complementary law between gap and contact force. It is then possible to admit velocity jumps and simultaneous contacts. The second interest is to handle, simultaneously, more complex interactions like friction, wear, adhesion, some of them often written with velocity. The interaction laws are not detailed in this paper. We simply summarize them in a formal relationship,

R(v, r) = 0, (2)
between the relative velocity at the contact point between the two interacting grains, v, and the impulse at the same contact point, r. The relationship R is usually composed of nonlinear and multivalued mappings between the previous two dual quantities. The assembly of the interaction-related quantities for all interactions is also written formally in the same way [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF]. Both v and r are expressed in the local coordinate basis to the contacts between the interacting grains, see Figure 1. Therefore, they are linked to the global kinematic and static quantities using a compatibility condition, H being the mapping matrix with a finite rotation:

v = H T V and R = Hr (3)

Reference problem and generic solver

Taking the dynamic equation ( 1) and the compatibility conditions (3) into account, the reduced dynamics involving 'material' quantities can be obtained:

P n t s S j S i e 1 e 2 e 3 y i x i z i y j x j z j G j G i Figure 1.
Coordinate basis; Rg = (e 1 , e 2 , e 3 ): global coordinate basis (for interface quantities), R i = (x i , y i , z i ): local (related to the grain S i ) coordinate basis (for grain dynamics), Rα = (n, t, s): local (related to an interaction α between two grains) coordinate basis (for interactions). Details of the local contact frame (n, t, s) at a contact α between two touching particles S i and S j at point P .

Therefore, they are linked to the global kinematic and static quantities using a compatibility condition, H being the mapping matrix with a finite rotation:

v = H T V and R = Hr (3)

Reference problem and generic solver

Taking the dynamic equation ( 1) and the compatibility conditions (3) into account, the reduced dynamics involving 'material' quantities can be obtained:

W r -v = -v d (4)
where W is the Delassus operator, W = H T M -1 H, and v d = H T M -1 R d . To close the problem, one adds the constitutive relation [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF], and the reference problem reads:

W r -v = -v d R(v, r) = 0 (5)
The difficulty to solve this problem is at least two-folds: on one hand, the number of unknowns (number of interaction quantities r and v) may be large (for instance, an average of 2 10 6 unknowns for the 3D problem illustrating this paper), and the Delassus operator W is not well conditioned. On the other hand, the constitutive relation is nonsmooth (e.g. it is non linear, and not differentiable). To address the nonsmoothness issue, the NSCD (nonsmooth contact dynamics) method with a nonlinear Gauss-Seidel (NLGS) solver [START_REF] Moreau | Numerical aspects of sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF] is used. The convergence of such an algorithm is a difficult topic with few theoretical results [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF] because of the multiplicity of solutions and the lack of regularity of the equations. Based on numerical experiments, the NLGS algorithm reveals to be the most robust solver that supports the more general interactions.

To address the large size of the problem, two substructuring approaches are proposed.

DOMAIN DECOMPOSITION APPROACHES

Several domain decomposition approaches have been proposed for contact problems, see [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Avery | The FETI family of domain decomposition methods for inequality-constrained quadratic programming: Application to contact problems with conforming and nonconforming interfaces[END_REF][START_REF] Dostál | A theoretically supported scalable TFETI algorithm for the solution of multibody 3d contact problems with friction[END_REF][START_REF] Dostál | Scalable TFETI with optional preconditioning by conjugate projector for transient frictionless contact problems of elasticity[END_REF][START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF][START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF][START_REF] Avery | The FETI family of domain decomposition methods for inequality-constrained quadratic programming: Application to contact problems with conforming and nonconforming interfaces[END_REF] amongst others, but a few are concerned with multiple contacts in granular media (for an algebraic-like partition, see [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF], and for a geometric partitioning of the discrete granular domain, see [START_REF] Breitkopf | Modélisation parallèle des matériaux granulaires[END_REF][START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF][START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF][START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF]). Here, we focus on the comparison of two formulations, one similar to FETI approaches, the other closer to additive Schwarz approaches, and of two implementations of the underlying communication scheme.

W r -v = -v d ( 4
)
where W is the Delassus operator, W = H T M -1 H, and

v d = H T M -1 R d .
To close the problem, one adds the constitutive relation [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF], and the reference problem reads:

W r -v = -v d R(v, r) = 0 (5)
The difficulty to solve this problem is at least two-folds: on one hand, the number of unknowns (number of interaction quantities r and v) may be large (for instance, an average of 2 × 10 6 unknowns for the 3D problem illustrating this paper), and the Delassus operator W is not well conditioned. On the other hand, the constitutive relation is nonsmooth (e.g. it is non linear, and not differentiable). To address the nonsmoothness issue, the NSCD (nonsmooth contact dynamics) method with a non-linear Gauss-Seidel (NLGS) solver [START_REF] Moreau | Numerical aspects of sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF] is used. The convergence of such an algorithm is a difficult topic with few theoretical results [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF] because of the multiplicity of solutions and the lack of regularity of the equations. Based on numerical experiments, the NLGS algorithm reveals to be the most robust solver that supports the more general interactions.

To address the large size of the problem, two substructuring approaches are proposed, as described below.

Domain decomposition approaches

Several domain decomposition approaches have been proposed for contact problems, see [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Avery | The FETI family of domain decomposition methods for inequality-constrained quadratic programming: Application to contact problems with conforming and nonconforming interfaces[END_REF][START_REF] Dostál | A theoretically supported scalable TFETI algorithm for the solution of multibody 3d contact problems with friction[END_REF][START_REF] Dostál | Scalable TFETI with optional preconditioning by conjugate projector for transient frictionless contact problems of elasticity[END_REF][START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF][START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF][START_REF] Avery | The FETI family of domain decomposition methods for inequality-constrained quadratic programming: Application to contact problems with conforming and nonconforming interfaces[END_REF] amongst others, but a few are concerned with multiple contacts in granular media (for an algebraic-like partition, see [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF], and for a geometric partitioning of the discrete granular domain, see [START_REF] Breitkopf | Modélisation parallèle des matériaux granulaires[END_REF][START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF][START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF][START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF]). Here, we focus on the comparison of two formulations, one similar to FETI approaches, the other closer to additive Schwarz approaches, and of two implementations of the underlying communication scheme.

Splitting choice and grain multiplicity

The domain has to be split into subdomains in order to use parallel computing with data locality. This decomposition is performed as frequently as needed to take into account the migration of grains from one subdomain to another. Since the nonsmoothness may occur in interactions between grains, we choose to distribute interactions among subdomains as in [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF] (we proceed by distributing the middle points between the centers of mass of interacting grains, according to their coordinates, using an arbitrary regular underlying grid, see Figure 2). If a grain indexed with S i is connected with m i subdomains, m i is called its multiplicity number. The interface between two subdomains is defined as the set of these grains, that joins the subdomains. The nonsmoothness is therefore localized only within the subdomains. This modeling choice is identical to [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] and somehow the dual of that proposed in [START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF], where nonlinearities (contact on crack lips) are isolated at the interfaces. Note that a direct algebraic partitioning of the reference problem can also be chosen, leading to a dual partitioning and a different algorithm [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF]. Some advantages and disadvantages have been discussed in previous works [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF][START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF].

A grain is defined as belonging to (at least) one subdomain E if it is connected to an interaction that belongs to this subdomain. Therefore, a boolean matrix B E selecting kinematic degrees of freedom of grains belonging to subdomain E allows to define the grain velocities in this subdomain as:

V E = B E V (6)
With this definition of the mapping matrix, one can check that the diagonal matrix of the grain multiplicities is

n sd E B T E B E .

FETI-like domain decomposition (NSCDD algorithm)

partitioning and not a geometrical partitioning. For the elastic deformable model of the grains, this splitting can be performed with a classical mesh decomposer. This leads to a partition of unity over the inertia parameters, as:

ME = B E DM B T E , (7)
with:

D kl = 0 if k = l 1/m i if k = l (8)
for entries k related to the grain S i . The partition of unity property reads:

M = E B T E ME B E .
This topic is investigated in details in [START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF]. In each subdomain E, the problem is identical to the global one (with the subscript E), provided that a term arising from the inter-grain interface is added. It can be built from the interconnecting condition (on the velocity jumps of boundary grains) that has been added to 'glue' neighboring subdomains, where A ΓE is a signed boolean matrix with a finite rotation, to map the grain velocities V E to the global coordinate basis into which the null velocity jump on the grain interface is expressed:

n sd E=1 A ΓE V E = 0 (9)
Γ denotes the global interface of all the interface grains. Formally the previous summation is performed on all the subdomains (whose number is equal to n sd ); in a practical way, for a given grain interface, only the neighboring subdomains have to be considered. We then obtain a FETI-like formulation [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF] for the reference problem using a multiplier field F Γ and the notation ÂT

ΓE = H T E M -1 E A T ΓE , WE = H T E M -1 E H E : WE r E -v E -ÂT ΓE F Γ = -v d E R(v E , r E ) = 0 E = 1, . . . , n sd n sd E=1 A ΓE V E = 0 (10)
The reduced problem on (r

E ,v E ,F Γ ), with the notations f = E A ΓE M -1 E R d E , v d E = H T E M -1 E R d E and X = E A ΓE M -1 E A T ΓE and a partial condensation of the problem, reads: WE r E -v E -ÂT ΓE F Γ = -v d E R(v E , r E ) = 0 E = 1, . . . , n sd XF Γ - n sd E=1 ÂΓE r E = f (11)
One easily reformulates the interface equation as an incremental problem [START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF]: if F Γ is associated to a velocity field V with the dynamical equations, and if this last field is not continuous at the interface, the correction of the impulse field F Γ is ∆F Γ such that

X∆F Γ = n sd E=1 A ΓE V E = [V ] |Γ (12)
the last term being the residual on the interface, i.e. the velocity jump [V ] |Γ . As for many domain decomposition approaches, the goal is to be able to localize the same typical problem that is under consideration on each subdomain independently, while designing a suited coupling recovery algorithm between subdomains, i.e. on the interface.

Here, the algorithmic formulation described in Algorithm 1 has been implemented into the LMGC90 platform [START_REF] Dubois | LMGC90[END_REF] for time-evolution problems (N is the number of time steps). At each new time step of the incremental solving procedure, the mapping H and the contact graph have to be updated within a contact detection phase. Eventually, the domain could also be repartitioned according to the new contact graph.

WE rE -vE = -v d E R(v E , rE ) = 0 Update (r E , v E ) ← (r E ,

NSCDD interface solving procedures

Centralized communication scheme (CCS). At first glance, the interface gluing step ( 12) is defined as a global linear equation linking all the subdomains. It requires communications between the subdomains such that one process gathers all the velocity contributions to the vector of velocity jumps [V ] |Γ . The value of the Lagrange multipliers F Γ computed sequentially is then distributed such that subdomain E receives its minimal required data amount A T Γ E F Γ .

Decentralized communication scheme (DCS). Due to the structure of the interface operator X, extensively studied in [START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF], each distributed database (per process related to subdomain E * ) is sufficient to get the elementary contribution to the interface operator:

X Γ E * = n sd E=1 A Γ E * E M -1 E A T Γ E * E (13)
A Γ E * E is a signed boolean matrix, mapping grains of subdomain E to velocity jumps of the elementary interface Γ E * (restriction of the global interface to the boundary of subdomain E * ). Then, an elementary interface problem can be defined as:

X Γ E * ∆F Γ E * = n sd E=1 A Γ E * E V E (14)
Finally, the data gathering of

n sd E=1 A Γ E * E V E
on each process corresponds to data exchanges over an unstructured topology. That means the exchanges do not concern only the subdomains which are geometrically neighbors, since some large/elongated rigid bodies may go through several subdomains. A common example of such a situation is a granular medium confined by rigid walls as illustrated in Figure 3. A wall performs a connection and exchanges between several non-neighboring subdomains. With the computation of the assembling of local contributions, it is easy to show that this is the expected iterated vector:

∆F Γ = n sd E=1 B Γ E D E B T Γ E E ∆F Γ E (15)
B Γ E is a boolean matrix selecting interface grains among subdomains, B Γ E E is a boolean matrix selecting elementary interface grains among subdomains and D E is a diagonal matrix with value 1/m i for entries corresponding to grain i.

Schwarz-like domain decomposition

We define a Schwarz-like domain decomposition method based on fundamental properties of rigid body dynamics dealing with diffuse contact network. Considering the selecting operator B E , directly derived from partitioning strategy defined in section 3.1 in (6), dynamic equations (1) per subdomain, i.e. for all the grains within a given subdomain, read:

B E (M V ) = B E R + R d (16)
Considering an overlapping between subdomains, the restricted inertia matrix to a subdomain E can be defined as,

M E = B E M B T E ( 17
)
The overlapping is rendered by the lack of partition of the unity in [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF] in comparison with [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. Here the domain partitioning defines a minimal overlapping between subdomains, since the dynamical operators M E , for E = 1, ..., n sd duplicate the contribution of the interface grains. Then the left hand side reads:

B E (M V ) = M E V E (18)
In a similar way, the given impulse R d , which comes from velocity at the beginning of the time step and the external forces acting on grains, depends only on inner quantities of subdomain E:

B E R d = R d E (19)
As an interface grain supports contacts in more than one subdomain the coupling term between subdomains arises from contact impulses as:

R = n sd E * =1 B T E * H E * r E * (20)
When extracting resultant impulse quantities according to B E , B E R splits into internal and coupling terms such that:

B E R = R E + R ddm E with R E = H E r E (21)
and

R ddm E = E * =E B E B T E * R E * (22)
where B E B T E * defines a boolean matrix whose non-zero entries correspond to grains from subdomain E * belonging also to subdomain E.

The dynamics of the granular medium per subdomain is therefore expressed with the vector equation:

M E V E -H E r E = R d E + R ddm E ( 23
)
and the reduced problem, with

v ddm E = H T E M -1 E R ddm E , reads: W E r E -v E = -v d E -v ddm E R(v E , r E ) = 0 E = 1, . . . , n sd (24)
The algorithmic formulation of this method (see Algorithm 2) does not define the interface behavior but only data exchanges between subdomains over an unstructured topology (as for the NSCDD decentralized communication scheme). This differs also from the dual partitioning developed in [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF], for which a non smooth problem similar to [START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF] is defined on the global interface composed of contacts, as a standard subdomain. Impulses acting on the boundary grains of a subdomain related to contacts belonging to interconnected subdomains are then taken as boundary conditions without defining a global interface problem as in the NSCDD method. This structure of DDM algorithm has the flavor of Schwarz algorithms [START_REF] Lions | On the Schwarz alternating method. I[END_REF][START_REF] Cai | A minimum overlap restricted additive schwarz preconditioner and applications in 3d flow simulations[END_REF][START_REF] Gander | Overlapping schwarz waveform relaxation for the heat equation in n dimensions[END_REF]. This method will be called Schwarz-NSCD in the following.

= H T E M -1 E E * =E B E B E * R E * Compute (r E ,v E ) with n GS non-linear Gauss-Seidel iterations on W E rE -vE = -v d E -vddm E R(v E , rE ) = 0 Update (r E , v E ) ← (r E , vE ) end while

Update grain positions in parallel end for

Algorithm 2 may be interpreted as a Gauss-Seidel-like method with a special numbering. Indeed, choosing α, β as indexes of contacts belonging to subdomain E, the Non Linear Gauss-Seidel k + 1 iteration for contact α reads:

W αα r k+1 α -v k+1 α = -v d α -vddm α -β<α W αβ r k+1 β -β>α W αβ r k β R(r k+1 α , v k+1 α ) = 0 ( 25
)
If i is the index (in the global numbering of the sample) of the interface grain associated to a contact α, and γ E * the indexes of contacts of subdomain E * interconnected with E via the grain i, and choosing n GS = 1, vddm α reads:

vddm α = H T α M -1 i E * =E R i E * = H T α M -1 i E * =E γ E * H γ E * r k γ E * = E * =E γ E * W αγ E * r k γ E * (26)
The choice n GS = 1 means that the interface exchanges are performed after each Gauss-Seidel iteration; this is justified in previous works [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF] for iterative schemes with robust but slow convergence. Rewriting equations (25) using the

W αα r k+1 α -v k+1 α = -v d α -β<α W αβ r k+1 β -β>α W αβ r k β -E * =E γ E * W αγ E * r k γ E * R(r k+1 α , v k+1 α ) = 0 ( 27
) These last equations reveal an hybrid Non linear Gauss-Seidel/Jacobi algorithm for the calculus of (r E , vE ) since the last summation relates only to values of r γ E * at the previous iteration.

Comparison of the different approaches

The influence of the proposed communication schemes is studied regarding the CPU time percentage consumed during MPI (Message Passing Interface [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF]) exchanges (Table 1) with respect to the whole CPU time of a simulation. The test consists of a sample with 16000 spheres submitted to an isotropic compaction, over 500 time steps (Figure 4). Increasing the number of subdomains with a fixed problem size leads to increase the interface size and to reduce the local problem size on each subdomain.

Results (obtained on the supercomputer HPC@LR) presented in Table 1 show clearly the gain we may obtain considering DCS compared to CCS. Decentralized communication scheme (similar for NSCDD-DCS or Schwarz-NSCD) decreases drastically the time consumed in MPI exchanges in comparison to centralized one, such that time percentage consumed in exchanges between processors stays quite small for a reasonable domain partitioning. An heuristic rule for 3D dense granular material may be to chose a domain partitioning such that each domain contains at least 10 3 grains.

PERFORMANCES AND SCALABILITY RESULTS

Numerical example

Sample preparation. A dense packing composed of 2.10 5 spheres is prepared in two phases. First a layer-by-layer deposition model based on simple geometrical rules is performed [START_REF] Voivret | Space-filling properties of polydisperse granular media[END_REF]. The particles are deposited sequentially on a substrate. Each new particle is placed at the lowest possible position on the free surface as a function of its diameter. This procedure leads to a random packing in which each particle supports at least three contacts. To avoid long range ordering a small polydispersity in size is used.

Once performed this geometrical process, the second phase is a compaction consisting of an isotropic compression inside a prismatic frame of dimensions L 0 × l 0 × h 0 in which three walls are fixed, and three walls are subjected to a compressive stress σ 0 . The gravity g and friction coefficients µ between the particles and the walls are set to zero during the compression in order to avoid force gradients and obtain an isotropic dense packing, see Figure 5(a).

Triaxial test. The isotropic sample is then subjected to a vertical compression by downward displacement of the top wall with a constant velocity v z for a constant confining stress σ 0 acting on the lateral walls, see Figure 5(b). The friction coefficient µ between particles is set to 0.35 and to zero with the walls. Since we are interested in a quasistatic macroscopic behavior, the shear rate should be such that the kinetic energy supplied by shearing is negligible compared to the static pressure. The cumulative strain components ε λ are defined as:

ε λ = λ λ0 dλ λ = ln 1 + ∆λ λ 0 , λ ∈ {h, L, l} (28)
where h 0 , L 0 and l 0 are the initial height, width and length of the simulation box, respectively, and ∆h = h -h 0 , ∆L = L -L 0 and ∆l = l -l 0 are the corresponding cumulative displacements. The volumetric strain reads:

ε V = V V0 dV V = ln 1 + ∆V V 0 , (29)
where V 0 is the initial volume and ∆V = V -V 0 is the total volume change. The cumulative shear strain is defined as:

ε V ≡ ε h -ε l . (30)
The solid fraction ρ is defined by ρ = V p /V , with V p the volume occupied by particles.

Mechanical validation. The simulation is performed using the domain partitioning depicted in Figure 6 up to ε V = 0.65. The domain partitioning and the NSCDD algorithm does disturb neither the macroscopic behavior nor the micromechanical statistical analysis of the granular medium, as shown in [START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF], whereas increasing the number of particles allows to avoid boundary effects and improves the statistical analysis one can perform. Due to the indeterminacy of the problem, we monitored both the quality of the solution obtained considering classical convergence criteria, into each subdomain [START_REF] Cambou | Micromechanics of Granular Materials[END_REF] and over the interface [START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF], and the mean/max interpenetration quantities: mean interpenetration remained below 0.02% of the min radius of the sample and maximum In the following we choose to study the numerical performances of the only NSCDD algorithm, with decentralized communication scheme, for a large granular sample.

Performances and scalability results

Numerical example

Sample preparation. A dense packing composed of 2 × 10 5 spheres is prepared in two phases. First a layer-by-layer deposition model based on geometrical rules is performed [START_REF] Voivret | Space-filling properties of polydisperse granular media[END_REF]. The particles are deposited sequentially on a substrate. Each new particle is placed at the lowest possible position on the free surface as a function of its diameter. This procedure leads to a random packing in which each particle supports at least three contacts. To avoid long range ordering a small polydispersity in size is used.

Once performed this geometrical process, the second phase is a compaction consisting of an isotropic compression inside a prismatic frame of dimensions L 0 × l 0 × h 0 in which three walls are fixed, and three walls are subjected to a compressive stress σ 0 . The gravity g and friction coefficients µ between the particles and the walls are set to zero during the compression in order to avoid force gradients and obtain an isotropic dense packing, see Figure 5(a).

Triaxial test. The isotropic sample is then subjected to a vertical compression by downward displacement of the top wall with a constant velocity v z for a constant confining stress σ 0 acting on the lateral walls, see Figure 5(b). The friction coefficient µ between particles is set to 0.35 and to zero with the walls. Since we are interested in a quasistatic macroscopic behavior, the shear rate should be such that the kinetic energy supplied by shearing is negligible compared to the static pressure. The cumulative strain components ε λ are defined as:

ε λ = λ λ0 dλ λ = ln 1 + ∆λ λ 0 , λ ∈ {h, L, l} (28)
where h 0 , L 0 and l 0 are the initial height, width and length of the simulation box, respectively, and ∆h = h -h 0 , ∆L = L -L 0 and ∆l = l -l 0 are the corresponding cumulative displacements. The volumetric strain reads: where V 0 is the initial volume and ∆V = V -V 0 is the total volume change.

ε V = V V0 dV V = ln 1 + ∆V V 0 , (29)
The cumulative shear strain is defined as:

ε V ≡ ε h -ε l . (30)
The solid fraction ρ is defined by ρ = V p /V , with V p the volume occupied by particles.

Mechanical validation. The simulation is performed using the domain partitioning depicted in Figure 6 up to ε V = 0.65. The domain partitioning and the NSCDD algorithm does disturb neither the macroscopic behavior nor the micromechanical statistical analysis of the granular medium, as shown in [START_REF] Visseq | Dense granular dynamics analysis by a domain decomposition approach[END_REF], whereas increasing the number of particles allows to avoid boundary effects and improves the statistical analysis one can perform. Due to the indeterminacy of the problem, we monitored both the quality of the solution obtained considering classical convergence criteria, into each subdomain [START_REF] Cambou | Micromechanics of Granular Materials[END_REF] and over the interface [START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF], and the mean/max interpenetration quantities: mean interpenetration remained below 0.02% of the min radius of the sample and maximum interpenetration remained below 10% of the minimum radius of the sample. Figure 7 depicts one map of the particles velocities at ε V over a slice of the sample. We see clearly the topology of the shear bands associated to localization phenomena. Figure 8 shows the contact network we deal with, composed of more than 2 × 10 6 force components.

Scalability in nonsmooth evolutive problems

We choose to quantify the parallel performances of the NSCDD algorithm by the usual quantities: the speedup S p , which is the ratio of the run time of the non-parallel version on a single processor to the run time of the parallel version on N p processors and the efficiency E p = S p /N p . Results are obtained setting the parameter n GS equals to 1 and for 3 typical states of the sample during its evolution. Complementary informations on those states are given in Table 2.

• "Loose state" (ρ = 0.567) occurs in the preparation phase which is an essential stage for providing an isotropic sample. More precisely it is obtained by deposit routines, at the beginning of the isotropic compaction phase. This state is characterized by a relatively small number of active contacts in comparison with the number of particles. Moreover a frictionless contact law leads to a quick convergence; the reduced dynamic problem is therefore quite cheap to solve. The goal is to provide a random close packing (RCP) sample.

• "RCP state" (ρ = 0.648 in our case) is used to start triaxial test, where the number of contacts is maximal. Due to confinement, the number of NSCDD iterations is large because of the indeterminacy of the problem amplified by the friction.

• "Residual state" (from ρ = 0.617 to ρ = 0.612 in our case) corresponds to a plateau of the shear stress evolution for which the shear bands appearance in the sample implies the decrease of the number of contacts in comparison with the RCP state.

Figure 9 depicts evolution of the speedup as a function of the number of processors (equal to the number of subdomains). Results obtained strongly depend on the state of the granular material. For the "loose state" the speedup diverges rapidly from the linear graph with slope equal to 1. Instead that, for the "RCP state" and the "Residual state", the speedup behavior is almost linear but not with the same slope (apart from one value related to the "RCP state" for 36 processors which illustrates variability phenomenon when using computer center). Speedup culminates to a value of 88 for 96 subdomains in "RCP state". This shows that the NSCDD method is particularly relevant for dense granular material. Table 3 shows that the three analyzed states differ as regards to the repartition, among the main stages of the algorithm, of CPU elapsed time percentage. Strong differences appear for "loose state" compared to the two others: partitioning stage becomes non negligible, as regards to solving procedure of reduced dynamics, whereas the number of contacts N C and the NSCDD iterations It NSCDD are smaller and the number of migrating grains is larger (see Table 2). In all cases the partitioning stage is re-done every 10 time steps to detect migrating grains. More precisely a grain is declared migrant when it supports contacts in a subdomain where it did not stand previously.

Table 4 shows that the geometric partitioning proposed, based on a regular underlying grid, does not assure a perfect load balancing. Indeed the number of bodies belonging to a subdomain varies up to 23%. For interface quantities this phenomena is even more present, as the number of interface grains per subdomain N m>2 E and the number of interconnected subdomains for a given subdomain N Isd E vary up to 73.8% and 91.4% respectively. This implies strong fluctuations on the CPU time percentage consumed in MPI exchanges during the iterative process for solving non smooth dynamics. Nevertheless, the contacts, which are the main unknowns of dense granular problems, are well distributed among subdomains. Obtaining a good load balancing for both grains and contacts is quite a hard topic when considering a simple geometric partitioning but may be attempted using automatic graph partitioner such as ParMETIS [START_REF] Karypis | ParMETIS: Parallel graph partitioning and sparse matrix ordering[END_REF], or other approaches such as orthogonal recursive bisection that proved their efficiency in [START_REF] Koziara | A distributed memory parallel multibody contact dynamics code[END_REF][START_REF] Shojaaee | An adaptive hierarchical domain decomposition method for parallel contact dynamics simulations of granular materials[END_REF].

Conclusions

This article focused on numerical and parallel performances of a domain decomposition method suited to granular dynamics analysis. Two versions, namely a FETI-like approach and a Schwarz-like approach, have been implemented into the LMGC90 software platform and tested on typical 3D simulations of granular dynamics, and on a modern supercomputing hardware platform. Even without a dedicated load balancing, the speedup of the simulations exhibits good parallel behavior, for large scale granular systems. To do so, both contact detection and solver have been parallelized, up to 100 processors used for a problem with an average of 2 × 10 6 frictional contact interaction unknowns. Next step will consist in a dynamic load balancing strategy, with a distributed database dynamical update to improve efficiency. On the other hand, additional strategies may decrease the computational cost by model adaptation and homogenization, to tackle even larger simulations, while maintaining a good parallel efficiency. These topics are currently under development.

Figure

  Figure1: Coordinate basis; R g = (e 1 , e 2 , e 3 ): global coordinate basis (for interface quantities), R i = (x i , y i , z i ): local (related to the grain S i ) coordinate basis (for grain dynamics), R α = (n, t, s): local (related to an interaction α between two grains) coordinate basis (for interactions). Details of the local contact frame (n, t, s) at a contact α between two touching particles S i and S j at point P .

Figure 2 :

 2 Figure 2: Geometrical partitioning of the discrete domain: inner grains are plotted empty and interface grains striped.

Figure 3 :

 3 Figure 3: A rigid wall supporting contacts in a long range (a whole side of the 16000 spheres sample partitioned in twelve subdomains), is involved in the local interface of six subdomains.

Figure 4 :

 4 Figure 4: Sample of 16000 spheres submitted to isotropic compaction; tested decompositions to study the time spent in communications depending on the communication scheme. For convenience walls are not plotted on these figures.

Figure 5 .

 5 Figure 5. Isotropic compaction (a) and triaxial loading (b)

Copyright c 2013 Figure 5 :

 20135 Figure 5: Isotropic compaction (a) and triaxial loading (b)

Figure 6 :

 6 Figure 6: Domain partitioning of the sample at ε V = 0 (a) and ε V = 0.7 (b) and underlying structure of corner grains at ε V = 0.7 (c). For convenience (a) and (b) gray scale is restricted to [1, 4] and (c) gray scale is restricted to [3, 7]; m = 1: inner grain; m = 2: face grain; m > 2: corner grain. The configuration (b) is split into 100 subdomains in order to adapt the partitioning to the aspect ratio.

Figure 7 :

 7 Figure 7: Velocity component along y axis at ε V = 0.5

Figure 8 :

 8 Figure 8: Force contact network at ε V = 0.5; full sample (a) and zoom (b). Only normal components are depicted.

Figure 9 :

 9 Figure 9: Speedup and efficiency as functions of the number of processors.

  1: Coordinate basis; R g = (e 1 , e 2 , e 3 ): global coordinate basis (for interface quantities), R i = (x i , y i , z i ): local (related to the grain S i ) coordinate basis (for grain dynamics), R α = (n, t, s): local (related to an interaction α between two grains) coordinate basis (for interactions). Details of the local contact frame (n, t, s) at a contact α between two touching particles S i and S j at point P .

  vE ) Compute RE and correct the velocity on interface grains to get A ΓE VE In sequential, but may be possibly parallelized (DCS version) Compute ∆F Γ as: X∆F Γ =

	n sd E=1 A ΓE VE and update interface
	impulses F Γ
	end while
	Update grain positions in parallel
	end for

Table 1 :

 1

	CCS	NSCDD-DCS/Schwarz-NSCD

Comparison of elapsed CPU time percentage consumed during MPI exchanges for centralized (NSCDD-CCS) and decentralized (NSCDD-DCS and Schwarz-NSCD) communication schemes; isotropic compaction of a 16000 spheres sample. n sd Partitioning grid (n x × n y × n z ) NSCDD-

Table 2 :

 2 Main characteristic values of the 3 states used to study the speedup behavior of the simulations; N C is the mean number of contacts per time step, It NSCDD is the mean number of NSCDD iterations per time step and N MG is the mean number of migrating grains per time step (i.e. whose subdomain is changing).

	State	N C	It NSCDD N MG
	"Loose"	339957	701	38
	"RCP"	669665	11978	24
	"Residual" 567957	9993	29

Table 3 :

 3 Time consuming percentage of the main steps of the NSCDD algorithm; Partitioning stands for the data gathering, sequential routines of partitioning of the sample, scattering of local interface informations and treatment of migrating grains; Contact detection relates to the parallel construction of H E operators; Inloop 1 stands for parallel calculus routines for iterative solving of the reduced dynamics (Gauss-Seidel iterations and interface updating); Inloop 2 stands for MPI exchanges during the iterative solving of the reduced dynamics; Update stands for parallel calculus routines for updating velocities and generalized configuration of grains.

	State	Partitioning Contact detection Inloop 1 Inloop 2 Update
	"Loose"	11.0%	2.6%	66.7%	10.7%	9%
	"RCP"	4.4%	0.7%	76.7%	15.6%	3.9%
	"Residual"	2.0%	0.2%	82.1%	12.5%	3.2%

Table 4 :

 4 Characteristic values per subdomains ("RCP state"); N S E stands for the number of grains belonging to subdomain E; N C E stands for the number of contacts; N Isd E stands for the number of subdomains interconnected with a considered subdomain E; N m>2 E stands for the number of corners grains; T MPI Exch stands for the percentage of MPI exchange time consumed when solving reduced dynamics; minimum, maximum and variation percentage.

		Min	Max Variation percentage
	N S E N C E N Isd E N m>2 E	2086 6514 20 445	2633 6795 52 965	23.2% 4.2% 91.4% 73.8%
	T MPI Exch 3.84% 26.7%	149%

This first DDM considers a non overlapping partition of the sample. For consistency with the rigid model of the grains, the masses and moments of inertia are distributed among the neighboring subdomains according to their multiplicity number. More precisely the distribution of masses and inertia is an algebraic

Acknowledgements. This work was realized with the support of HPC@LR, a Center of Competence in High-Performance Computing from the Languedoc-Roussillon region, funded by the Languedoc-Roussillon region, Europe and University Montpellier 2 Sciences et Techniques. The HPC@LR Center is equipped with an IBM hybrid Supercomputer.


Asma Atamna

Anne Auger

Nikolaus Hansen

Augmented Lagrangian Constraint Handling for CMA-ES-Case of a Single Linear Constraint

We consider the problem of minimizing a function f subject to a single inequality constraint g(x)  0, in a black-box scenario. We present a covariance matrix adaptation evolution strategy using an adaptive augmented Lagrangian method to handle the constraint. We show that our algorithm is an instance of a general framework that allows to build an adaptive constraint handling algorithm from a general randomized adaptive algorithm for unconstrained optimization. We assess the performance of our algorithm on a set of linearly constrained functions, including convex quadratic and ill-conditioned functions, and observe linear convergence to the optimum.

Introduction

Evolution strategies (ESs) are derivative-free continuous optimization algorithms that are now well-established to solve unconstrained optimization problems of the form min x f (x), f : R n ! R, where n is the dimension of the search space. The state-of-theart ES, the covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF], is especially powerful at solving a wide range of problems and particularly ill-conditioned problems [START_REF] Hansen | Impacts of Invariance in Search: When CMA-ES and PSO Face Ill-Conditioned and Non-Separable Problems[END_REF][START_REF] Auger | Experimental Comparisons of Derivative Free Optimization Algorithms[END_REF]. It typically exhibits linear convergence. The default CMA-ES algorithm implements comma selection where the best solution is not preserved from one iteration to the next one (contrary to plus selection). Comma selection is an important feature of CMA-ES that entails robustness of the algorithm to various types of ruggedness including noise.

Linear convergence being a central aspect of an ES in the unconstrained case, a

(1 + 1)-ES using an adaptive augmented Lagrangian constraint handing-to deal with a single inequality constraint-has been introduced in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF] with the motivation to obtain a linearly converging algorithm. Empirical results show the linear convergence of the algorithm on the sphere and moderately ill-conditioned ellipsoid functions, subject to one linear constraint. In [START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF], the authors present a variant of the previous (1 + 1)-ES with augmented Lagrangian constraint handling and study theoretically its linear convergence using a Markov chain approach. In both mentioned works, the step-size is adapted using the 1/5th success rule [START_REF] Kern | Learning Probability Distributions in Continuous Evolutionary Algorithms -A Comparative Review[END_REF] while the covariance matrix is fixed to the identity. On ill-conditioned problems, however, adapting the covariance matrix is crucial. It is hence natural to wonder whether it is possible to design a CMA-ES variant with augmented Lagrangian constraint handling. The algorithms presented in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF], however, use plus selection and can thus a priori not be used directly to design such a variant.

In this context, we consider the constrained problem of minimizing f : R n ! R subject to a single inequality constraint g(x)  0, g : R n ! R. More formally, we write min x f (x) subject to g(x)  0 .

(1)

We bring to light that the algorithms previously presented in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF] derive from a more general framework that seamlessly allows to build an adaptive constraint handling algorithm from a general adaptive stochastic search method. We then naturally apply this finding to build a (µ/µ w , )-CMA-ES variant with adaptive augmented Lagrangian constraint handling. We opted for using the median success rule step-size adaptation (MSR) [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF] because it is an extension of the 1/5th success rule algorithm used in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF]. We then test the resulting algorithm-the (µ/µ w , )-MSR-CMA-ES with augmented Lagrangian constraint handling-on a set of functions, including convex quadratic as well as ill-conditioned functions, subject to one linear inequality constraint.

The rest of this paper is organized as follows: we introduce some basics about augmented Lagrangian in Section 2. Then, we define the general framework and apply it to the (µ/µ w , )-MSR-CMA-ES in Section 3. We present our empirical results in Section 4 and conclude with a discussion in Section 5.

Notations

We introduce here the notations that are not explicitly defined in the rest of the paper. We denote R + the set of positive real numbers and R + > the set of strictly positive real numbers. N > is the set of natural numbers without 0. x 2 R n is a column vector, x | is its transpose, and 0 2 R n is the zero vector. ||x|| denotes the Euclidean norm of x and ⇠ equality in distribution. (µ/µ w , ) denotes comma selection with weighted recombination and (1 + 1) denotes plus selection with one parent and one offspring. I n⇥n 2 R n⇥n is the identity matrix. x i is the ith component of vector x. The derivative with respect to x is denoted r x . Finally, 1 {A} returns 1 if A is true and 0 otherwise.

Augmented Lagrangian Methods

Augmented Lagrangian methods are constraint handling approaches that transform the constrained optimization problem into an unconstrained one where an augmented Lagrangian is optimized [START_REF] Hestenes | Multiplier and Gradient Methods[END_REF][START_REF] Powell | A Method for Nonlinear Constraints in Minimization Problems[END_REF].

The augmented Lagrangian consists of a Lagrangian L and a penalty function, with L : R n+1 ! R defined as

L(x, ) = f (x) + g(x) (2)
for the objective function f subject to one constraint g(x)  0, where 2 R is the Lagrange factor. The Lagrangian encodes the KKT stationarity condition which states that, given some regularity conditions are satisfied (constraint qualifications), if x ⇤ 2 R n is a local minimum of the constrained problem, then there exists a constant ⇤ 2 R + , called the Lagrange multiplier, such that

r x f (x ⇤ ) + ⇤ r x g(x ⇤ ) | {z } rxL(x ⇤ , ⇤ ) = 0 ,
where we assume here that f and g are differentiable at x ⇤ . A penalty function is combined with the Lagrangian L to create the augmented Lagrangian h. There exist different ways to construct the augmented Lagrangian and we refer to [START_REF] Nocedal | Numerical Optimization[END_REF] for a deeper discussion about this topic. In this work, we use the following augmented Lagrangian

h(x, , !) = f (x) + ( g(x) + ! 2 g 2 (x) if + !g(x) 0 2 2! otherwise , (3)
where ! > 0 is a penalty factor. The same augmented Lagrangian was used for the first time within an ES in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF]. The function h is minimized successively with respect to x, and and ! are updated so that approaches the Lagrange multiplier ⇤ and ! favors feasible solutions. By adapting , the penalty factor ! does not have to grow to infinity to achieve convergence, unlike with quadratic penalty function methods [START_REF] Nocedal | Numerical Optimization[END_REF]. Let x opt be the optimum of the constrained problem in (1) and let opt be the corresponding Lagrange multiplier. If f and g are differentiable at x opt , then for all ! > 0,

r x h(x opt , opt , !) = r x f (x opt ) + max(0, opt + !g(x opt ))r x g(x opt ) = 0 .

A General Framework for Adaptive Augmented Lagrangian Constraint Handling

In [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF] and [START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF], the authors present two (1 + 1)-ESs with an augmented Lagrangian constraint handling approach for the optimization problem in (1). The algorithms derive from a general framework for building a constraint handling adaptive algorithm. This framework starts with a randomized adaptive algorithm for minimizing an unconstrained function f : R n ! R: the randomized adaptive algorithm can be identified by the sequence of its states s t at iteration t that are iteratively computed from an update function F such that

s t+1 = F f (s t , U t+1 ) , (4)
where the superscript indicates the function being minimized and where (U t ) t2N> is a sequence of independent identically distributed (i.i.d.) random vectors. For instance, in the case of a (1 + 1)-ES in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF], the state is a vector of the search space (current estimate of the optimum) and a step-size. We assume that the state s t of the algorithm includes a vector X t 2 R n which typically encodes the current estimate of the optimum at iteration t. Note that the transition function F above includes a step where candidate solutions are sampled from the current state s t and the random vector U t+1 , and evaluated on the objective function f .

From the adaptive algorithm above, we construct an algorithm with adaptive constraint handling to take into account a single constraint in the following way: we add to the state of the algorithm two scalars t and ! t that correspond respectively to the Lagrange factor and the penalty factor of the augmented Lagrangian h at iteration t. Therefore, the state at iteration t is s t 0 = [s t , t , ! t ]. The objective function used at each iteration to evaluate a candidate solution X i t+1 is now

h ( t ,!t) (X i t+1 ) := h(X i t+1 , t , ! t ) , ( 5
)
where h is the augmented Lagrangian defined in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF]. Finally, the update of the state s t 0 of the adaptive algorithm with augmented Lagrangian constraint handling takes place in two steps: first, s t is updated via

s t+1 = F h ( t ,! t ) (s t , U t+1 ) , (6)
where candidate solutions are now evaluated on h ( t ,!t) instead of f . Then, the parameters t and ! t of h are updated. In [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF], t is updated according to

t+1 = max(0, t + ! t g(X t+1 )) , (7)
while in [START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF], the authors use the following update

t+1 = t + ! t g(X t+1 ) . (8)
For ! t , the following update is used in both [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF] and [START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF]

! t+1 = 8 > < > : ! t 1/4 if ! t g 2 (X t+1 ) < k 1 |h(Xt+1, t ,!t) h(Xt, t ,!t)| n or k 2 |g(X t+1 ) g(X t )| < |g(X t )| ! t 1 otherwise , (9)
for some constants > 1, k 1 , k 2 2 R + . Based on these examples, we introduce some general update functions G and G ! for the updates of t and ! t defined implicitly via

t+1 = G g (( t , ! t ), X t+1 ) (10)
! t+1 = G (f,g) ! ((X t , t , ! t ), X t+1 ) . (11)
The superscript in G and G ! indicates that the function value is used in the update.

3.1 The (µ/µ w , )-MSR-CMA-ES with Adaptive Augmented Lagrangian

We now apply the general framework sketched above to the covariance matrix adaptation evolution strategy (CMA-ES) with median success rule step-size adaptation (MSR). We start by presenting the algorithm for the unconstrained case then we give the updates of the augmented Lagrangian parameters t and ! t .

The (unconstrained) CMA-ES with MSR The original CMA-ES with MSR is given in Algorithm 1, without the highlighted parts. The algorithm proceeds iteratively: at each iteration t, candidate solutions (offspring) X i t+1 , i = 1, • • • , , are sampled according to Line 5, where X t 2 R n is the current estimate of the optimum (mean vector), t 2 R + is the step-size, and

U i t+1 2 R n , i = 1, • • • , , are i.i.d
. random vectors sampled from the normal distribution N (0, C t ), with mean 0 2 R n and covariance matrix C t 2 R n⇥n . The offspring are ordered according to their fitness (f -value in the unconstrained case) in Line 6, where i : is the index of the ith best offspring. The µ best offspring (parents) are then recombined (Line 7) to create the new mean vector X t+1 , where the weights

w i > 0, i = 1, • • • , µ, satisfy w 1 > • • • > w µ and P µ i=1 w i = 1.
The step-sized t is adapted in Lines 8 to 11 using the MSR step-size adaptation [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF]. MSR is a success-based step-size adaptation method which extends the well-known 1/5th success rule step-size adaptation [START_REF] Kern | Learning Probability Distributions in Continuous Evolutionary Algorithms -A Comparative Review[END_REF], used with plus selection, to comma selection. The step-size is adapted depending on "success", where the success is defined as the median offspring X m( ) t+1 (fitness-wise) of the current population being better than the jth best offspring X j: t of the previous population. In practice, we choose j to be the 30th percentile-the value for which the median success probability is roughly 1/2 on the sphere function with optimal step-size [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF]. The number K succ of offspring better than

X j: t is computed in Line 8. Note that K succ /2 is equivalent to h(X m( ) t+1 , t , ! t )  h(X j:
t , t , ! t ). Therefore, we define the success measure z t in Line 9 such that z t 0 if and only if X m( ) t+1 is successful. z t is cumulated in q t+1 (Line 10) and, finally, t is updated in Line 11: it increases in the presence of success (q t+1 > 0) and decreases otherwise in order to increase the probability of success.

The covariance matrix C t is adapted with CMA [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] in Lines 12 and 13. The update is a combination of the so-called rank-one-update and rank-µ-update. A detailed discussion on CMA can be found in [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF].

Finally, the jth best offspring is updated in Line 17. Therefore, the state of the algorithm in the unconstrained case is

s t = (X t , t , q t , p t , C t , X j: t ) .
The constrained (µ/µ w , )-MSR-CMA-ES with adaptive augmented Lagrangian As explained in the general framework, the fitness f is replaced with the augmented Lagrangian h in the constrained case. The parameters t and ! t are adapted in Lines 15 and 16 in Algorithm 1, where changes in comparison to the unconstrained case are highlighted in gray.

The Lagrange factor t is adapted in Line 15. It is increased when the new solution X t+1 is unfeasible and decreased otherwise, unless it is zero. The derivation of this update is discussed in details in [START_REF] Nocedal | Numerical Optimization[END_REF].

For the penalty parameter ! t , we use the original update proposed in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF] for the

(1 + 1)-ES with augmented Lagrangian. The update rule is given in Line 16. ! t is increased either when (i) the augmented Lagrangian h does not change "enough" after t and ! t are updated to avoid stagnation. This is translated by the first inequality where

! t g 2 (X t+1 ) ⇡ |h(X t+1 , t + t , ! t + ! t ) h(X t+1 , t , ! t )| is compared to the change in h due to the change in X t , |h(X t+1 , t , ! t ) h(X t , t , ! t )|.
! t is also increased when (ii) the change in the value of the constraint function is not large enough (second inequality in Line 16). To prevent an unnecessary ill-conditioning of the problem, ! t is decreased whenever conditions (i) and (ii) are not satisfied.

Empirical Results

We evaluate Algorithm 1 on the sphere function (f sphere ), two ellipsoid functions (f elli ) with condition numbers ↵ = 10 2 , 10 6 , f 2 sphere , f 0.5 sphere , the different powers function

Algorithm 1 (µ/µ w , )-MSR-CMA-ES with Augmented Lagrangian Constraint Handling

0 given n 2 N>, = 2 1/n , k1 = 3, k2 = 5, µ, 2 N>, j = 0.3 , 0  wi < 1, µ X i=1 wi = 1, µeff = 1/ µ X i=1 w 2 i , c = 0.3, d = 2 2/n, cc = 4 + µeff/n n + 4 + 2µeff/n c1 = 2 (n + 1.3) 2 + µeff , cµ = min ✓ 1 c1, 2 µeff 2 + 1/µeff (n + 2) 2 + µeff ◆ 1 initialize X0 2 R n , 0 2 R + > , C0 = In⇥n, t = 0, q0 = 0, p0 = 0, constrained problem // true if the problem is constrained, false otherwise 2 if constrained problem 3 initialize 0 2 R, !0 2 R + > 4 while stopping criteria not met 5 X i t+1 = Xt + tU i t+1 , U i t+1 ⇠ N (0, Ct), i = 1, • • • , // sample candidate solutions 6 Extract indices {1 : , • • • , : } of ordered candidate solutions such that ( h(X 1: t+1 , t, !t)  . . .  h(X : t+1 , t, !t) if constrained problem f (X 1: t+1 )  . . .  f (X : t+1 ) otherwise 7 Xt+1 = µ X i=1 wiX i: t+1 = Xt + t µ X i=1 wiU i: t+1 // recombine µ best candidate solutions 8 Ksucc = 8 > < > : P i=1 1 {h(X i t+1 , t ,!t)h(X j: t , t ,!t)} if constrained problem P i=1 1 {f (X i t+1 )f (X j: t )} otherwise 9 zt = 2 ✓ Ksucc 2 ◆ // compute success measure 10 qt+1 = (1 c )qt + c zt 11 t+1 = t exp ✓ qt+1 d ◆ // update step-size 12 pt+1 = (1 cc)pt + p cc(2 cc)µeff ✓ Xt+1 Xt t ◆ // cumulation path for CMA 13 Ct+1 = (1 c1 cµ)Ct + c1pt+1p | t+1 + cµ µ X i=1 wi ✓ X i t+1 Xt t ◆ ✓ X i t+1 Xt t ◆ | // update covariance matrix 14 if constrained problem 15 t+1 = max(0, t + !tg(Xt+1)) // update Lagrange factor 16 !t+1 = 8 > < > : !t 1/4 if !tg 2 (Xt+1) < k1 |h(X t+1 , t ,!t) h(Xt, t ,!t)| n or k2|g(Xt+1) g(Xt)| < |g(Xt)| !t 1 otherwise // update penalty factor 17 X j: t+1 = Xt + tU j: t+1
// update jth best solution 18 t = t + 1 (f diff pow ), and the Rosenbrock function (f rosen ), with one linear inequality constraint. The functions are defined in Table 1. We consider the case where the constraint is active at the optimum x opt , i.e. g(x opt ) = 0. We choose the optimum to be at x opt = (10, • • • , 10) | and construct the constraint function, g(x) = b | x + c, so that the KKT stationarity condition is satisfied at

x opt with opt = 1. Therefore, b = r x f • (x opt ) | and c = r x f • (x opt )x opt ,
for each function. Note that all considered functions are differentiable at

x opt = (10, • • • , 10) | . Name Definition f ↵ sphere (x) 1 2 P n i=1 x 2 i ↵ felli(x) 1 2 P n i=1 ↵ i 1 n 1 x 2 i Name Definition fdiff pow(x) pP n i=1 |xi| 2+4 i 1 n 1 frosen(x) P n 1 i=1 10 2 (x 2 i xi+1) 2 + (xi 1) 2
Table 1: Definitions of the tested functions, where f sphere := f 1 sphere .

For the step-size and the covariance matrix adaptation, we use the Python implementation of CMA-ES whose source code can be found at [1], with the default parameter setting detailed in [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF]. We run the algorithm 11 times in n = 10, with X 0 sampled uniformly in [START_REF] Auger | Experimental Comparisons of Derivative Free Optimization Algorithms[END_REF][START_REF] Auger | Experimental Comparisons of Derivative Free Optimization Algorithms[END_REF] n , 0 = 1, 0 = 5, and ! 0 = 1. The results are presented for one run in Figures 1 (f sphere , f 2 sphere , and f 0.5 sphere ) and 2 (f elli with ↵ = 10 2 , 10 6 , f diff pow , and f rosen ). On the left column of each figure are graphs of the evolution of the distance to the optimum ||X t x opt ||, the step-size t , the distance to the Lagrange multiplier | t opt |, and the penalty factor ! t in log-scale. On the right column of the figures are graphs representing the evolution of the coordinates of the mean vector X t .

Graphs on the right column of Figures 1 and2 show the overall convergence of the algorithm to x opt . We also observe linear convergence of X t to x opt , as well as linear convergence of t to opt and t to 0 (left column of Figures 1 and2). Moreover, ||X t x opt ||, | t opt |, and t decrease at the same rate. On the other hand, the penalty factor ! t is observed to converge to a stationary value after a certain number of iterations. We sometimes observe a stagnation in graphs of kX t x opt k due to numerical precision.

The largest convergence rate (when excluding the initial adaptation phase) is observed on f sphere and the smallest one on f 0.5 sphere , where there is a factor of approximately 1.5 between the two convergence rates. However, there is some variance in the empirical convergence rate. In particular, on 11 performed runs we observe the highest variance in the empirical convergence rate for f elli with ↵ = 10 6 , f diff pow , and f rosen .

On f elli with ↵ = 10 6 , f diff pow , and f rosen , we observe a stagnation of X t in the early stages of the algorithm (left column in Figure 2). The reason is that the adaptation of the covariance matrix takes longer on ill-conditioned problems. This explains the slow convergence of some coordinates of X t to 10 (right column in Figure 2). Once the covariance matrix is adapted, the convergence occurs.

When comparing 11 single runs of Algorithm 1 to the (1 + 1)-ESs with augmented Lagrangian in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF] (not shown for space reasons) on constrained f sphere , f elli (in n = 10), it appears that on f sphere , Algorithm 1 needs approximately up to 1.5 times more function evaluations than algorithms in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF] to reach a distance to the optimum of 10 4 . On f elli with ↵ = 10 2 , however, Algorithm 1 is faster and needs approximately 1.3 times less function evaluations to reach the same distance, with ↵ = 10 6 , Algorithm 1 is around 167 times faster to reach a target of 15 (this large difference is due to the adaptation of the covariance matrix). Fig. 1: Single runs of (µ/µ w , )-MSR-CMA-ES with augmented Lagrangian on f sphere (top row), f 2 sphere (middle row), and f 0.5 sphere (bottom row) in n = 10. The optimum x opt = (10, • • • , 10) | . Left: evolution of the distance to the optimum, the distance to the Lagrange multiplier, the penalty factor, and the step-size in log-scale. Right: evolution of the coordinates of X t .

Discussion

Linear convergence is a key aspect of ESs in both unconstrained and constrained optimization scenarios. As stated in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF], the minimum requirement for a constraint handling Fig. 2: Single runs of (µ/µ w , )-MSR-CMA-ES with augmented Lagrangian on f elli with ↵ = 10 2 (first row), f elli with ↵ = 10 6 (second row), f diff pow (third row), and f rosen (fourth row) in n = 10. The optimum x opt = (10, • • • , 10) | . Left: evolution of the distance to the optimum, the distance to the Lagrange multiplier, the penalty factor, and the step-size in log-scale. Right: evolution of the coordinates of X t .

ES is to converge linearly on convex quadratic functions with a single linear constraint. On the other hand, an algorithm for constrained optimization should be able to tackle ill-conditioned problems. Having that in mind, we proposed a (µ/µ w , )-CMA-ES with an augmented Lagrangian approach for handling one inequality constraint, where the choice of the augmented Lagrangian constraint handling was motivated by the promising results of its implementation for the (1 + 1)-ESs with 1/5th success rule in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF].

Moreover, we showed that our algorithm-as well as (1 + 1)-ESs with augmented Lagrangian constraint handling in [START_REF] Arnold | Towards au Augmented Lagrangian Constraint Handling Approach for the (1 + 1)-ES[END_REF][START_REF] Atamna | Analysis of Linear Convergence of a (1 + 1)-ES with Augmented Lagrangian Constraint Handling[END_REF]-is an instance of a more general framework for building an adaptive constraint handling algorithm from a general adaptive algorithm for unconstrained optimization.

Experiments on linearly constrained convex quadratic functions, as well as ill-conditioned functions (including the ellipsoid and Rosenbrock functions), showed linear convergence of our algorithm to the unique optimum of the constrained problem.

?? Research centre Saclay-Île-de-France, TAO team, lastname@lri.fr

Acknowledgments This work was supported by the grant ANR-2012-MONU-0009 (NumBBO) of the French National Research Agency.