Ralph O Suarez

Olivier Commowick

Xavier Tomas-Fernandez

Sanjay Prabhu

Simon K Warfield

AUTOMATED DETECTION OF WHITE MATTER FIBER TRACTS

Keywords: automation, diffusion, tractography, DTI, fiber segmentation, STAPLE

White matter (WM) fiber bundles of the brain can be delineated by diffusion tractography utilizing anatomical regions-of-interest (ROI). These ROIs can specify seed regions in which tract generation algorithms are initiated. Interactive identification of such anatomical ROIs enables the detection of the major WM fiber tracts, but suffers from inter-rater and intra-rater variability, and is time consuming. We developed and compared three techniques for automated delineation of ROIs for the detection of two major WM fiber tracts in 12 healthy subjects. Tracts identified automatically were compared quantitatively to reference standard tracts derived from carefully hand-drawn ROIs. Based on comparative performance of the experimental techniques, a multi-subject label fusion algorithm was found to generate tracts most consistent with the reference standard.

INTRODUCTION

Diffusion imaging is a magnetic resonance imaging (MRI) technique that measures 3D diffusion properties of water molecules in the brain. Because the presence of microstructural barriers to otherwise uniform diffusion gives shape to the orientation diffusion function, inferences can then be made regarding the underlying white matter (WM) microstructure [START_REF] Conturo | Tracking neuronal fiber pathways in the living human brain[END_REF][START_REF] Westin | Processing and visualization for diffusion tensor MRI[END_REF][START_REF] Mori | Imaging cortical association tracts in the human brain using diffusion-tensorbased axonal tracking[END_REF]. In addition, scalar quantities derived from the diffusion orientation, such as fractional anisotropy and mean diffusivity, have important applications in medicine as these measures reflect regional axon myelination and fiber connectivity characteristics [START_REF] Goldberg-Zimring | Diffusion tensor magnetic resonance imaging in multiple sclerosis[END_REF].

Many important clinical applications for WM tractography involve the segmentation of specific fiber bundles of interests (as opposed to whole-brain evaluations). For example, presurgical applications typically involve a focused assessment of the major fiber bundles distinctly located within or in the vicinity of the surgical site [START_REF] Yogarajah | Defining Meyer's loop-temporal lobe resections, visual field deficits and diffusion tensor tractography[END_REF]; evaluation of disease progression in neurodegenerative disorders similarly involves a regional assessment of WM connectivity and axonal integrity, often within the fiber bundles specific to clinical symptoms [START_REF] Goldberg-Zimring | Diffusion tensor magnetic resonance imaging in multiple sclerosis[END_REF].

A common approach used for the segmentation of important WM pathways initially involves determination of specific sets of regions-of-interest (ROIs) in the subject anatomy. These ROIs are then used for tractography seeding, or as inclusion and/or exclusion regions. To achieve the segmentation of specific WM fiber bundles, a technique currently under wide use involves the manual drawing of tractography seeding ROIs by a trained expert [START_REF] Wakana | Fiber tract-based atlas of human white matter anatomy[END_REF][START_REF] Van Hecke | A tracking-based diffusion tensor imaging segmentation method for the detection of diffusion-related changes of the cervical spinal cord with aging[END_REF]. While generally considered to be accurate, this procedure is timeconsuming, requires extensive training, and presents opportunity for human fatigue, error, or inter-rater variability to impact the segmentation outcome. As such, the atomization of WM fiber bundle segmentation is a topic that is increasingly addressed in the literature [START_REF] Van Hecke | A tracking-based diffusion tensor imaging segmentation method for the detection of diffusion-related changes of the cervical spinal cord with aging[END_REF][START_REF] Hua | Tract probability maps in stereotaxic spaces: analyses of white matter anatomy and tract-specific quantification[END_REF][START_REF] Zhang | Automated fiber tracking of human brain white matter using diffusion tensor imaging[END_REF].

In this article, we present a comparative study of three experimental methodologies which we designed for automated segmentation of the corpus callosum and corticospinal WM tracts in individual subject anatomies. Tractography segmentations obtained by the experimental methods are compared against reference standard segmentations derived from manual placement of seeding ROIs by a trained expert, and are judged based on the rootmean-square difference between the resulting probabilistic tract density images. Based on our findings in this study, we recommend an automated algorithm derived from the Simultaneous Truth and Performance Level Estimation framework (STAPLE) [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF].

One should note, however, that our study does not address other important questions still under research in the field, such as, determination of the most robust algorithms for either WM fiber tractography or for nonrigid subject-tosubject anatomical registration. Accordingly, we recognize that more advanced tractography and/or nonrigid alignment algorithms may be implemented with our automated fiber segmentation procedure in the future.

MATERIALS AND METHODS

Subjects

We enrolled an age-matched study group made up of 12 young volunteers (mean age: 16.8 yrs, SD: 1.1 yrs), with no history of neurobiological disorders; 5 were male and 7 female. In compliance with Institutional Review Board (IRB) guidelines, institutional study approval, and signed individual participant consent forms were obtained.

MR Imaging

All of the study participants underwent MRI scanning at three Tesla (3T) using a Magnetom Trio scanner (Siemens, USA) to obtain high-quality T1-weighted anatomical and diffusion tensor images. Anatomical images: Highresolution 3D Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) T1-weighted images were obtained using a 32-channel head coil; we used a 24cm FOV, 1.0mm contiguous slice thickness, sagittal slices covering the entire head, TR/TE=1410ms/2.27ms, matrix 256x256, TI=800 ms, and a flip angle=9 o . Diffusion images: High spatial resolution diffusion weighted tensor images (DTI) were acquired using 30 directions and 5 baselines, B value=1000; we used isotropic 2x2x2 mm 3 voxels (24cm FOV, matrix 128x128, 2mm thick contiguous slices covering the entire brain). In order to minimize eddy current distortion, a double balanced gradient scheme was applied. To achieve the smallest possible TE in order to minimize geometric distortion and signal loss, a parallel acquisition with an acceleration factor of 2 was used with a 32-channel head coil.

Reference standard ROI placement

We evaluated the capability of automatically reproducing corpus callosum (CC) and corticospinal tract (CST) fiber bundle segmentations in all of the study volunteers. Initially, CC and CST seeding ROIs were manually draw directly on the particular subject's red, green, and blue (RGB) colorcoded images representing local tensor orientation [START_REF] Douek | MR color mapping of myelin fiber orientation[END_REF]. For delineation of the CC fiber bundle, one ROI was drawn encompassing approximately the 5 central sagittal slices outlining the RGB red-colored lateral radiations of the corpus callosum [START_REF] Catani | A diffusion tensor imaging tractorgrphy atals for vitual in vivo dissections[END_REF]. The CST bundle was based on two seeding ROIs: 1) approximately four axial slices of the bluepurple colored inferior-to-superior radiations of the posterior limb of the internal capsule, and 2) approximately four axial slices the blue-purple radiations of the corticospinal tracts traversing the pons (see Figure 1).

Tractography

We applied a previously unreported single tensor tractography scheme which we configured in order to merge log-Euclidean tensor interpolation [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF], combined tensor deflection [START_REF] Lazar | White matter tractography using diffusion tensor deflection[END_REF], and primary eigenvector trajectory orientation to stochastically identify highly oriented trajectories in diffusion tensor volumes. All of our fiber bundle segmentations were determined based on an ROI seeding strategy which made no requirement for exclusion regions. Standard criteria of low fractional anisotropy and rapid trajectory angle changes were used to terminate unrealistic trajectories [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF]. Prior to tractography, each subject's diffusion images were rigidly registered and resampled in order to match their T1-weighted anatomical volumes.

Automated methods

We designed three experimental methods for automated projection of manually drawn seeding ROIs from a reference template anatomy onto individual subject anatomies of interest. These transformations are based on nonrigid alignment of skull-stripped, T1-weighted grayscale images. The nonrigid alignment algorithm we used applies a block matching approach together with a local correlation coefficient as the similarity measure, and features an innovative outlier rejection scheme which effectively rejects local pairings that disagree completely with the rest of the transformation [START_REF] Commowick | Title: Création et Utilisation d'Atlas Anatomiques Numériques pour la Radiothérapie (Design and Use of Anatomical Atlases for Radiotherapy)[END_REF][START_REF] Garcia | A robust and efficient block-matching framework for non-linear registration of thoracic images[END_REF]. The resulting deformation fields from nonrigid alignment of T1-weighted images are subsequently applied to manually drawn ROI volumes, in effect projecting them onto that subject's diffusion image for subsequent tractography generation. The automated methods we designed are defined as follows:

1) Single-subject projection: we tested the procedure of projecting hand-drawn seeding ROIs from a single template anatomy onto the subject anatomy of interest. To test the reliability of this approach, we performed 11 similar experiments for each of the 12 study participant (leave-oneout). We refer to the single-subject projection method as SSproj from this point forward.

2) Atlas projection: this procedure is based on projecting atlas-drawn template seeding ROIs onto the subject anatomy of interest. These projections were done based on nonrigid alignment of the atlas and the T1-weighted subject anatomy. Atlases were generated by first performing nonrigid alignment of 11 individual skull-stripped, T1-weighted subject anatomies (leave-one-out), followed by computation of an average T1-weighted volume. Deformation fields determined from the alignment of T1-weighted images were then applied to the corresponding diffusion images, similarly aligning all diffusion images for computation of an average diffusion volume. In this way, we created geometrically unbiased RGB and T1-weighted atlases, for manual placement of seeding ROIs and subsequent atlas-tosubject alignments, respectively. We refer to this atlas projection method as Atlsproj from this point forward.

3) ROI label map consensus: this automated procedure is based on the determination of a local consensus ROI region in the subject anatomy. Initially, we project manually drawn seeding ROIs from 11 template anatomies onto the subject of interest, where we then perform a label map consensus analysis based on the STAPLE framework [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF]. STAPLE statistical evaluation of the underlying consensus across individually projected template ROIs consequently results in a subject-specific region which we designate as that subject's tractography seeding ROI. We refer to the ROI map consensus method as MPcns from this point forward.

Comparative evaluation of experimental results

Fiber bundle segmentations derived from each of the experimental methods were judged based on how similar the resulting probabilistic tract density images were to those generated by manually drawn ROIs in the same subject. Voxel values in probabilistic tract density images were determined based on the number of fiber trajectories traversing each voxel, divided by the total number of trajectories in the entire volume. As such, these images simultaneously quantify both the spatial extent, and also the density of the 3D fiber bundles produced by the tractography (see Figure 2). We use the root-mean-squared (RMS) difference between these images as quantitative measure of tract segmentation error. It follows that tract density images demonstrating perfect correspondence to the reference standard will yield RMS error equal to zero. This metric is comparable to Jaccard's Similarity Index, but offers the advantage of accounting for continuous intensity mismatch, in addition to the binary overlap.

RESULTS

Group-level results

CC segmentation outcomes: the lowest mean RMS error in the study group was observed for MPcns (Mean: 3.25E-5, SD: 0.61E-5); the second lowest error was observed for Atlsproj (Mean: 4.14E-5, SD: 0.71E-5); the highest RMS error was observed for SSproj (4.24E-5, 0.87E-5). ANOVA analysis revealed no statistically significant difference in the mean RMS error across all of the experimental methods we tested.

CST segmentation outcomes: CST segmentations generally yielded higher mean RMS error than was observed for CC (approximately one order of magnitude greater); this was true across all of the methods we tested. The lowest error we observed was for MPcns (10.01E-5, 1.49E-5); the second lowest for SSproj (12.44E-5, 2.16E-5); and the highest RMS error was observed for Atlsproj (14.51E-5, 2.36E-5). ANOVA analysis revealed no statistically significant difference in the mean RMS error across all of the experimental methods we tested.

Subject-level results

All three of the automated segmentation methods we designed were able to delineate CC and CST fiber bundles in 12 out of the 12 study participants, Jaccard's Similarly Indices of 0.49 or greater in every case. The resulting fiber bundles generally differed however, as was reflected in the calculated values of RMS error, and 3D visualizations (see Figure 3). CC segmentation outcomes: MPcns resulted in the lowest RMS error, in 11 out of 12 subjects; the second lowest error was observed for Atlsproj, in 5 out of 12 subjects; the highest error was observed for SSproj, in 6 of 12 subjects (see Figure 4).

CST segmentation outcomes: MPcns yielded the lowest error in 12 out of 12 subjects; the second lowest error was observed for SSproj, in 12 out of 12 subjects; the highest RMS error was observed for Atlsproj, in 12 out of 12 subjects (see Figure 4).

DISCUSSION

We designed and tested three algorithms for automated segmentation of CC and CST WM fiber bundles. Our goal was to evaluate the feasibility of reproducing reference standard segmentations derived from the manual drawing of seeding ROIs by a trained expert. Performance of the automatic methods was quantified by computing the RMS difference in the resulting probabilistic tract density images. All three of the proposed methods were able to reproduce CC and CST, in 12 out of 12 subjects, and all generally yielded comparable segmentation outcomes in group-level assessment. However, we demonstrated lower error from the MPcns procedure in every case we tested. Based on these findings, we recommend our Simultaneous Truth and Performance Level Estimation (STAPLE) procedure for the reliable automation of WM fiber bundle segmentation.

The results we presented should be viewed relative to important limitations of our study. The main limitation derives from the lack of a true validation method. Our findings were based exclusively on the ability to reproduce fiber segmentations derived from manually drawn ROIs; however, it is not clear how accurate this reference standard should be considered. Furthermore, our study was done with healthy volunteers; therefore, it remains to be seen how the recommended automated methods will perform in patient populations suffering from structural abnormalities.

In this study, we performed nonrigid alignments based on the T1-weighted modality; we note that other modalities should be evaluated in future work. In particular, increased reliability may be achieved by nonrigid aliments based on fractional anisotropy, mean diffusivity, or full tensor volumes. Additionally, a more exhaustive list of prominent WM fiber bundles should be explored in future study, including: the cingulum bundle, the optical radiation tracts, and the arcuate fasciculus.

Fig. 1 .

 1 Fig. 1. Hand-draw seeding ROIs (white) overlaid on the RGB image of a representative study participant. Left panel: a single seeding ROI drawn for tract generation of the corpus callosum (CC). Right panel: depiction of two ROIs used for tract generation of corticospinal tract (CST): the posterior limb of the internal capsule and the inferior-superior radiations of the pons.

Fig. 2 .

 2 Fig. 2. Probabilistic tract density images of the corpus callosum (CC) and the corticospinal tract (CST) in a representative study participant.

Fig. 3 .

 3 Fig. 3. Comparison of 3D subject tractography results for CC (top row) and CST (bottom row). Left panel: reference standard fiber bundles from manually drawn seeding ROIs. Right panel: fiber bundles generated from automated SSproj, Atlsproj, and MPcns methodologies on the same representative subject.

Fig. 4 .

 4 Fig. 4. RMS error in corpus callosum (CC) (top panel) and corticospinal tract segmentations (CST) (lower panel) as a function of subject of interest, error bars represent the standard deviation in the mean RMS error from SSproj. Notice that the lowest RMS error was in every case obtained by the Mpcns procedure.


H Le Dret

A Raoult

Homogenization of hexagonal lattices

Keywords: Graphene sheet, homogenization, Γ-convergence, Cauchy-Born rule. AMS Subject Classification: 74Q05, 74Q15, 74K35, 49J45

We characterize the macroscopic effective mechanical behavior of a graphene sheet modeled by a hexagonal lattice of elastic bars, using Γ-convergence.

Introduction

We consider a graphene sheet modeled by a hexagonal network of elastic bars, see [START_REF] Caillerie | Discrete homogenization in graphene sheet modeling[END_REF]. We are interested in deriving an equivalent continuum mechanics model for the deformations of the sheet by means of a homogenization procedure when the rest lengths of the bars go to 0, using Γ-convergence techniques in order to obtain rigorous convergence results. We are not concerned here with electronic properties of graphene, nor quantum or relativistic effects that occur in graphene. Our analysis also applies to other hexagonal networks of elastic springs or trusses of elastic bars. There is a comprehensive body of work on the homogenization of discrete networks, see for instance [START_REF] Alicandro | A general integral representation result for continuum limits of discrete energies with superlinear growth[END_REF][START_REF] Blanc | From molecular models to continuum mechanics[END_REF][START_REF] Braides | From Discrete to Continuum: a Variational Approach[END_REF][START_REF] Braides | Continuum limits of discrete systems without convexity hypotheses[END_REF][START_REF] Ming | Cauchy-Born rule and the stability of crystalline solids: Static problems[END_REF]. Let us mention that in a recent independent work [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF] on stochastic lattices leading to nonlinearly elastic models for polymers in the homogenization limit, a convergence theorem is proved which also applies to hexagonal lattices. This theorem gives rise to a formula for the homogenized energy density that is equivalent to the one we prove here. The methods used and the context are however quite different from ours.

We start with a careful presentation of the discrete problem: how to select the nodes in a hexagonal lattice that belong to the graphene sheet under consideration, how to impose a condition of place on part of the sheet, how to apply forces. We are thus able to treat realistic cases of graphene sheets by tackling all these aspects, which are frequently set aside in the literature, where often only infinite crystals without realistic boundary conditions or applied forces are considered, see 1 UPMC Univ Paris 06, UMR 7598 LJLL, Paris, F-75005 France 2 Laboratoire MAP5, UMR CNRS 8145, Université Paris Descartes, Paris, France 1 for example the critical remarks made by Ericksen [START_REF] Ericksen | On the Cauchy-Born rule[END_REF] on this subject. We do this however at the expense of a nonnegligible amount of notation and special cases that have to be checked separately in the sequel.

We next rewrite the problem as a sequence of problems in the calculus of variations, in the same spirit as [START_REF] Braides | From Discrete to Continuum: a Variational Approach[END_REF] and other works in the literature, by replacing the discrete displacements of the atoms in the sheet by functions defined on a domain. A hexagonal network is a complex lattice, and consequently, we have to introduce two independent functions to that effect, one of which is continuous piecewise affine, and the other is piecewise constant.

The original discrete problem is thus recast as a sequence of problems in the calculus of variations in which a functional, depending on a small parameter ε that represents the interatomic distance, is minimized over a set of admissible functions. The core of the article in Section 4 is then to characterize the limit of this sequence of minimization problems when ε tends to 0. We first prove uniform estimates, which is not immediate since the energy densities are not uniformly coercive, but vanish on roughly speaking half of the domain. Then, in a long series of technical lemmas, we prove the Γ-convergence of the sequence of functionals toward a functional of the calculus of variations, and we give a formula for the limit homogenized energy density. The general organization of the Γ-convergence argument is inspired by, but not a consequence of [START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF], with several simplifications on the one hand, and on the other hand several arguments that are entirely specific to the discrete to continuum limiting process that we consider here.

In Section 5, we establish a few properties of the homogenized energy density: material frame indifference, material symmetry, non convexity properties.

Finally, in Section 6 we show some numerical results which lead to several interesting observations in relation with the Cauchy-Born rule.

Part of the results of this article were announced in [START_REF] Dret | Homogenization of hexagonal lattices[END_REF].

Setting of the problem

In this article, we consider sheets, i.e., two-dimensional structures, that deform in three-dimensional Euclidean space. It is better for ulterior purposes to keep the two spaces separate, in the spirit of the Lagrangian description of continuum mechanics, rather than identifying the two-dimensional space with a particular plane of three-dimensional space. We thus consider R 2 equipped with a Euclidean structure and an orthonormal basis (e 1 , e 2 ), and R 3 also equipped with a Euclidean structure and an orthonormal basis (e 1 , e 2 , e 3 ). We can thus measure lengths both in R 2 and in R 3 .

Let ω be a bounded open connex subset of R 2 of class C 1 for simplicity. We assume that ω contains the reference configuration of a hexagonal lattice of nodes linked together by elastic bars of length ε 0 > 0, where ε 0 > 0 is the interatomic distance between nearest neighbors in graphene at equilibrium under zero loading. The sheet is subjected to boundary conditions of place, to be made precise later on.

It is also subjected to dead loading forces applied to the nodes and consequently deforms in R 3 .

Let us first describe the global, scale 1, hexagonal lattice in R 2 . We introduce the three vectors

s = √ 3e 1 , t = √ 3 2 e 1 +
3 2 e 2 and p = 1 3 (s + t).

In the description we use, the lattice is comprised of two types of nodes: The type 1 nodes that occupy points is + jt with (i, j) ∈ Z 2 , and the type 2 nodes that occupy points is + jt + p, again with (i, j) ∈ Z 2 , see Figure 1 below. The hexagonal lattice is thus a complex lattice, a superposition of two simple Bravais lattices which are translates of each other, shown with different dashed lines below. We are following here the standard description of such complex lattices. The hexagonal nature of the sheet is not yet apparent. We now assume that the internal energy of the sheet only derives from chemical bonds that join nearest neighboring type 1 and type 2 nodes. We model these bonds by bars. There are thus three types of bars: Type 1 bars parallel to sp, type 2 bars parallel to tp, and type 3 bars parallel to p, see Figure 2 below. This classification of bars is only for labeling reasons, all bars are physically equivalent. We scale the global, scale 1, complex lattice by the factor ε 0 . The sheet in its reference configuration is the largest subset of scaled closed hexagons, the union of which is contained in ω, see Figure 3 below. It is easily seen that the distance of ∂ ω to the set of nodes is at most 2ε 0 for ε 0 small enough. Moreover, it should be noted that most nodes are connected to three neighbors by bars, except for some nodes next to the boundary of ω that are connected to only two neighbors by bars.

The actual relative scales of the characteristic lengths of the domain ω and the interatomic distance are quite different from what is shown on Figure 3 in the situations we are interested in. Indeed the length of carbon-carbon bonds in graphene is 0.142 nm, whereas it is now possible to produce graphene at sizes up to 76 cm. Even though the interatomic distance is a fixed number ε 0 , the basis of our analysis will consist in embedding the problem into a similarly defined family of problems indexed by a sequence ε that tends to 0, which is quite reasonable given the above orders of magnitude. We will adopt this point of view from now on, all objects and quantities defined as above, with ε 0 replaced by ε.

Let us now turn to the mechanical side of the model. We first describe the deformations of the sheet. Let N ε 1 denote the set of integer pairs (i, j) ∈ Z 2 such that the type 1 node located at point ε(is + jt) belongs to the sheet and likewise N ε 2 for type 2 nodes located at ε(is + jt + p). The set N ε α thus indexes type α nodes in the sheet.

Under the action of applied loads and boundary conditions, each node of type α occupies a new equilibrium position in space χ ε α (i, j), given by a mapping

χ ε α : N ε α → R 3 , α = 1, 2.
This is a Lagrangian description of the deformations. We need to impose a condition of place on part of the sheet. In continuum mechanics, this is usually done on part of the boundary. Now of course, there is no particular reason here why any node would fall exactly on ∂ ω. It is however conceivable that part of the sheet could be bonded to a rigid substrate. We thus pick another regular open set ω 0 ⊂ ω of R 2 and impose the conditions

χ ε 1 (i, j) = (ε(is + jt); 0) if ε(is + jt) ∈ ω0 , (1)
on type 1 nodes, where the notation (x; 0) stands for the point in R 3 the first two coordinates of which are those of x in R 2 and the third is 0, which defines an arbitrary embedding of R 2 into R 3 . For ε small enough, condition (1) is non empty. Every type 2 node in the sheet is connected to either two or three type 1 nodes. We impose

χ ε 2 (i, j) = (ε(is + jt + p); 0) (2)
whenever the corresponding two or three type 1 nodes belong to ω0 . For instance, in Figure 3, there is no condition of place imposed on any type 2 node and a condition a place imposed on the two upper left type 1 nodes.

Roughly speaking, all nodes located in ω0 at a distance larger than ε of ω \ ω 0 are submitted to a condition of place. The set ω 0 , or more accurately the subset (ω 0 ; 0) of R 3 , thus models a rigid substrate to which the sheet is attached in a natural state. We could also impose a given deformation on ω 0 , for instance a prestressed compressive state. We now describe the energy of the sheet. As said earlier, we assume that each chemical bond is modeled by a bar. For simplicity, we assume that each bar acts as an elastic spring of stiffness κ > 0 and natural length ε. Thus, if a particular bar B k has deformed length k , then the elastic energy stored in this bar is given by

E ε k = κ( k -ε) 2 . (3)
The convergence analysis can be carried out for more general elastic energies. We retain expression (3) because it is the simplest frame indifferent energy such that the reference configuration is a natural state. As such, it is used in the chemical literature where experimental numerical values for κ can be found (for instance κ = 326 N.m -1 , see [START_REF] Odegard | Equivalent continuum modeling of nanostructured materials[END_REF]). When undergoing a deformation (χ ε 1 , χ ε 2 ), the sheet stores an elastic energy

E ε (χ ε 1 , χ ε 2 ) = n ε b ∑ k=1 E ε k ,
where n ε b is the total number of bars in the sheet. Note that bars both extremities of which are bonded to ω0 do not contribute to the total elastic energy. Only bars with at least one extremity not subjected to the condition of place are susceptible to length change under sheet deformation.

Deformed bar lengths are expressed using the relative displacements of nodes in space. A typical deformed bar length k thus assumes the form

k = |χ ε 1 (i 1 , j 1 ) -χ ε 2 (i 2 , j 2 )|,
where |•| denotes the Euclidean norm in R 3 and the integer pairs (i 1 , j 1 ) and (i 2 , j 2 ) correspond to each end of the bar.

To complete the description of the mechanical setting, we impose external dead loading forces on all nodes in the sheet. This includes the nodes that are bonded to the rigid substrate, even though they do not contribute to the minimization of the energy in the force term. We are thus given a function f : ω → R 3 , which we assume to be continuous and independent of ε, such that the external force acting on a node is ε 2 f (x), where x is the location of the node in the reference configuration. The ε 2 factor is irrelevant for the actual sheet, for which it is a constant. Since our plan is to let ε tend to 0, it will turn out to be the right scaling factor to yield a finite nonzero limit contribution.

The corresponding energy term reads

L ε (χ ε 1 , χ ε 2 ) = ε 2 ∑ (i, j)∈N ε 1 f (ε(is+ jt))• χ ε 1 (i, j)+ ∑ (i, j)∈N ε 2 f (ε(is+ jt + p))• χ ε 2 (i, j) ,
where • denotes the scalar product in R 3 . We consequently end up with a total energy for the sheet of the form

E ε (χ ε 1 , χ ε 2 ) = E ε (χ ε 1 , χ ε 2 ) -L ε (χ ε 1 , χ ε 2 ).
The equilibrium deformed configuration of the sheet minimizes the total energy among all possible deformations (χ ε 1 , χ ε 2 ) satisfying condition [START_REF] Alicandro | A general integral representation result for continuum limits of discrete energies with superlinear growth[END_REF]. The existence of such minimizers is obvious. Note that there is no uniqueness of minimizers.

Continuous formulation

In order to derive a limit continuous model when ε goes to 0, we replace the discrete unknowns χ ε α by unknown functions defined on ω, while keeping exactly the same values of the energy, see [START_REF] Alicandro | A general integral representation result for continuum limits of discrete energies with superlinear growth[END_REF][START_REF] Braides | From Discrete to Continuum: a Variational Approach[END_REF][START_REF] Meunier | Elastic limit of square lattices with three point interactions[END_REF]. We consider all the type 1 nodes in the sheet and the union of the equilateral triangles of edge length √ 3ε that they define. We call ωε the union of these triangles, see Figure 4. We assume that ω ε is included in ω. Since ω is regular, we see that meas (ω \ ω ε ) ≤ Cε for some constant C independent of ε.

We denote by L the Z-lattice generated by s and t. Let T ε be the triangulation of R 2 defined by the lattice εL. The type 1 nodes of the sheet are the vertices of this triangulation that belong to ωε , see Figure 4. We thus define a piecewise affine, R 3 -valued function ϕ ε on ωε by declaring that ϕ ε (ε(is + jt)) = χ ε 1 (i, j) for all (i, j) ∈ N ε 1 . Due to condition (1), we have

ϕ ε (x) = (x; 0) (4)
at all type 1 nodes of the sheet belonging to ω0 . So far ϕ ε is only defined on ωε .

To reformulate the discrete minimisation problem as a problem in the calculus of variations, we need to extend this function to ω in a controlled way.

For this, we add triangles that cover ∂ ω in a single or double layer as depicted in Figure 4 below (this is possible for ε small enough since we have chosen the largest union of hexagons included in ω and ω is regular). We call boundary triangle any triangle in ωε that touches ∂ ωε either on an edge or at a vertex, and exterior triangle any triangle in the added layer.

∂ ω ωε

A reflex angle "Small" triangles of ω ε ∂ in gray, see below. Proposition 3.1 There exists an extension operator A ε from the set of continuous piecewise affine functions on T ε ∩ ωε to the set of continuous piecewise affine functions on

T ε ∩ ω such that ∇(A ε ψ) L 2 (ω) ≤ C ∇ψ L 2 (ω ε )
where C is independent of ε, and if ψ is globally affine on the union of boundary triangles, then A ε ψ is the same affine function on ω \ ωε .

Proof. The proof is by inspection of all possible cases, which we skip here. It works mostly by extension from one boundary triangle with an edge in ∂ ωε to the exterior triangle sharing the same edge by the same affine function as in the boundary triangle, and then interpolation to fill out the remaining exterior triangles. There is a difficulty for reflex angles as indicated in Figure 4, where such an extension is inconsistent. In this case, we take the next type 1 node in ωε , located on the angle bisector and use these two values to define the value at the missing exterior node. We see that the gradients in the two exterior triangles are controlled by the gradients in the four boundary triangles.

Of course, when two reflex angle nodes are located next to each other, we pick one such extension for all the corresponding triangles, for example the leftmost one, and if a reflex angle node is located next to an acute or obtuse angle node, the reflex node extension takes precedence.

Finally, if a double layer of triangles is needed, as in the lower left corner of Figure 4, we just repeat the extension procedure.

In order to alleviate the notation, we will identify ϕ ε with its extension A ε ϕ ε in the sequel.

The deformations of type 2 nodes are taken into account via a piecewise constant deviation vector γ ε defined on ω. In any triangle that contains a type 2 node, which we call a full triangle (including exterior triangles), we let

γ ε (x) = χ ε 2 (i, j) -χ ε 1 (i, j),
with the exception of triangles with a type 2 node but no leftmost type 1 node, such as the two exterior triangles on the left, bottom of Figure 4, where we define

γ ε (x) = χ ε 2 (i, j) -χ ε 1 (i + 1, j)
instead. Finally, in triangles that do not contain a type 2 node, which we call empty triangles, we let γ ε (x) = 0. Due to condition (2), we have

γ ε (x) = ε(p; 0), (5)
for all x belonging to any full triangle all of which two or three type 1 nodes belong to ω0 , except for those with no leftmost type one node where

γ ε (x) = ε(p -s; 0). ( 6
)
If we know γ ε and ϕ ε , then we recover χ ε 1 and χ ε 2 . It is then a simple matter to express the deformed lengths 1 , 2 and 3 of bars of type 1, 2 and 3 respectively issuing from a type 2 node with these new unknowns. For most type 2 nodes, they read

1 = |ε∂ s ϕ ε (x) -γ ε (x)|, 2 = |ε∂ t ϕ ε (x) -γ ε (x)| and 3 = |γ ε (x)|, (7)
where ∂ u ϕ = Dϕ(u) is the directional derivative of ϕ in direction u, and x is any point in a full triangle containing these particular three bars (or two bars as at the top or on the right of Figure 4). The bars issuing from a type 2 node with no leftmost type 1 node require a special treatment. The type 1 and type 2 bars that attach them to their neighboring type 1 nodes have deformed lengths

1 = |γ ε (x)|, 2 = |ε∂ t-s ϕ ε (x) -γ ε (x)|, (8)
x being any point in the corresponding triangle. Formulas ( 7)-( 8) are the key to rewriting the discrete energy as a functional of the calculus of variations. Note that all the quantities they involve are piecewise constant.

From now on, all computations in the reference configuration will be performed in the oblique coordinate system based on s and t. It should be noticed that this entails a change of unit area as compared to the original orthogonal system in R 2 .

Let Y = T f ∪ T e be the unit area parallelogram obtained from the reference full triangle T f with vertices (0, 0), (1, 0) and (0, 1) in the oblique coordinate system, and the reference empty triangle T e with vertices (1, 0), (1, 1) and (0, 1). We use Y as the unit cell of our homogenization procedure. We recall that ωε denotes the union of all full and empty triangles the vertices of which are the type 1 nodes of the sheet. This constitutes the bulk of the sheet. We also introduce the triangles one of the vertices of which is a type 2 node linked by bars to only two type 1 nodes, and the two other vertices are these type 1 nodes, which we call small triangles (depicted in gray in Figure 4). We denote their union by ωε ∂ since it is close to part of the boundary. It has measure of the order of ε at most. We also assume that ω ε ∂ ⊂ ω. We denote by A(ε) the space of functions on ω that are of the form A ε ψ, with ψ piecewise affine continuous on T ε ∩ ωε and that satisfy condition [START_REF] Blanc | From molecular models to continuum mechanics[END_REF]. Likewise, we denote by C(ε) the space of functions on ω that are piecewise constant and zero in empty triangles on T ε ∩ ω and that satisfy conditions ( 5) and ( 6). We let

V (ε) = A(ε) ×C(ε).
Taking all of the above considerations into account, we can thus rewrite the elastic energy as an integral

E ε (χ ε 1 , χ ε 2 ) = I ε (ϕ ε , γ ε )
where for all (ψ, δ ) ∈ V (ε),

I ε (ψ, δ ) = ω ε W ε ε -1 x, Dψ(x), δ (x) dx + ω ε ∂ Z ε x, Dψ(x), δ (x) dx.
The first stored energy density W ε : R 2 × L (R 2 ; R 3 ) × R 3 → R, where L (R 2 ; R 3 ) denotes the space of linear mappings from R 2 into R 3 , is defined by

W ε (y, g, τ) = 2κ g(s) - τ ε -1 2 + g(t) - τ ε -1 2 + |τ| ε -1 2 (9) if y ∈ T f + L and W ε (y, g, τ) = 0, (10)
if y ∈ T e + L. The stored energy density is thus Y -periodic in the variable y. The factor 2 in front of equation ( 9) is because the reference full triangle is of area 1 2 . The second energy density

Z ε : ω ε ∂ × L (R 2 ; R 3 ) × R 3 → R is given by Z ε (x, g, τ) = 6κ g(u) - τ ε -1 2 + |τ| ε -1 2 , (11)
where u = s,t or ts depending on the orientation of the small triangle to which x belongs. The factor 6 in front of equation ( 11) is because reference small triangles are of area 1 6 . The energy terms corresponding to Z ε will have negligible contribution in the limit ε → 0.

We keep the force term as a discrete sum, noting that we can rewrite

L ε (χ ε 1 , χ ε 2 ) = F ε (ϕ ε , γ ε ) with F ε (ψ, δ ) = ε 2 ∑ (i, j)∈N ε 1 f (ε(is + jt)) • ψ(ε(is + jt)) + ∑ (i, j)∈N ε 2 f (ε(is + jt + p)) • (ψ(ε(is + jt)) + δ (ε(is + jt + p))) . (12)
In order not to complicate even more an already cumbersome notation in terms of indices, we wrote the above relation as though there were no small triangle with no leftmost type 1 node. It will become clear that these triangles play a negligible role in the sequel.

We now have a total energy functional defined for all (ψ, δ ) ∈ V (ε) by

J ε (ψ, δ ) = I ε (ψ, δ ) -F ε (ψ, δ ).
Finally, we extend the sheet energy functional to the space

H = L 2 (ω; R 3 ) × L 2 (ω; R 3 ) by letting J ε (ψ, δ ) = +∞, whenever (ψ, δ ) / ∈ V (ε).
It is clear that we have rephrased the equilibrium of the sheet as a problem in the calculus of variations:

Find (ϕ ε , γ ε ) ∈ H such that J ε (ϕ ε , γ ε ) = inf (ψ,δ )∈H J ε (ψ, δ ). ( 13
)
Our objective now is to let ε → 0 and find a limit problem that describes the asymptotic behavior of the continuous sheet deformation ϕ ε and deviation vector γ ε . This is a periodic nonlinear variational homogenization problem, see [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF][START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF], with several differences compared with the classical case: the energy functionals are +∞ outside of subspaces of H that depend on ε and the densities are not coercive on the unit cell Y .

Gamma-limit of the energies

We use Γ-convergence theory to study the asymptotic behavior of the minimization problem (13) when ε → 0. Let us briefly recall what Γ-convergence is about. Let (X, d) be a metric space and (J n ) n∈N a sequence of functionals X → R = R∪{+∞}. The sequence J n is said to Γ-converge to a functional J for the topology of X when n → +∞, if the following two conditions are satisfied: i) For all x ∈ X and all sequences x n → x in X, we have lim inf n→+∞ J n (x n ) ≥ J(x).

ii) For all x ∈ X, there exists a sequence y n → x such that J n (y n ) → J(x).

The two main virtues of Γ-convergence are first a compactness result, in that every sequence has a Γ-convergent subsequence, and second a result concerning minimizers or almost minimizers that states that if the minimizers of J n belong to a compact subset of X independent of n, then their limit points minimize J, see [START_REF] Maso | An Introduction to Γ-Convergence[END_REF]. The second result shows that Γ-convergence is the right tool to deal with sequences of problems in the calculus of variations.

For our purposes here, we are interested in computing Γ-lim ε→0 J ε in the strong topology of H, where ε denotes a sequence that tends to 0. It should be noted that the specific form of the energy introduced earlier plays next to no role in the ensuing analysis, and the Γ-convergence result holds true for more general energies, defined for instance on W 1,p (ω; R 3 ) with 1 < p < +∞.

A priori bounds

We start with an a priori bound. Let H 1 ω 0 (ω; R 3 ) denote the space of H 1 deformations ψ that satisfy the boundary condition of place ψ(x) = (x; 0) in ω 0 . Proposition 4.1 Let (ψ ε , δ ε ) be a sequence in H such that J ε (ψ ε , δ ε ) ≤ M for some M independent of ε. Then there exists C independent of ε such that

ψ ε H 1 (ω;R 3 ) ≤ C and δ ε L 2 (ω;R 3 ) ≤ Cε. In particular, δ ε L 2 (ω;R 3 ) → 0.
Proof. We first observe that the finiteness of the energy implies that (ψ ε , δ ε ) ∈ V (ε). In particular, ψ ε is piecewise affine and we have

Dψ ε L 2 (ω ε ) ≥ c Dψ ε L 2 (ω) ,
for some constant c > 0 independent of ε, by Proposition 3.1.

It is easy to see that for all S, T, z ∈ R 3 , we have

(|S -z| -1) 2 + (|T -z| -1) 2 + (|z| -1) 2 ≥ 1 6 (|S| 2 + |T | 2 + |z| 2 ) -6, (14)
without paying particular attention to optimizing the constants on the right. We let µ ε = δ ε ε . In view of equations ( 9), [START_REF] Ericksen | On the Cauchy-Born rule[END_REF] and estimate [START_REF] Dret | Homogenization of hexagonal lattices[END_REF], the elastic energy part can be estimated from below as

I ε (ψ ε , δ ε ) ≥ C 1 ω ε f |Dψ ε | 2 dx +C 2 ω ε ∂ |Dψ ε (u(x))| 2 dx +C 3 ω |µ ε | 2 dx -C 4 , (15)
with C 1 , C 2 and C 3 strictly positive and where ωε f denotes the union of all full triangles in ωε and u(x) = s,t or ts depending on to which small triangle of ω ε ∂ the point x belongs. Note that there is no problem extending µ ε by 0 to ω.

Let us consider each empty triangle in ωε . An empty triangle is either flanked by at least two full triangles, which is the general case, or by at least one full triangle and one small triangle. Let us start with the first case. Assuming the triangles are arranged as in Figure 6, we have

T ε,1 f T ε,2 f T ε e
∂ t-s ψ ε |T ε e = ∂ t-s ψ ε |T ε,1 f and ∂ t ψ ε |T e = ∂ t ψ ε |T ε,2 f
. Therefore,

T ε e |Dψ ε | 2 dx ≤ C T ε e |∂ t-s ψ ε | 2 + |∂ t ψ ε | 2 dx = C T ε,1 f |∂ t-s ψ ε | 2 dx +C T ε,2 f |∂ t ψ ε | 2 dx ≤ C T ε,1 f |Dψ ε | 2 dx +C T ε,2 f |Dψ ε | 2 dx,
and the same is true for the other two possible configurations. We perform a similar estimate for empty triangles flanked by one full triangle and one small triangle, we sum over all empty triangles and we obtain

ω ε \ω ε f |Dψ ε | 2 dx ≤ C ω ε f |Dψ ε | 2 dx + ω ε ∂ |Dψ ε (u(x))| 2 dx . ( 16
)
Putting estimates [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF] and ( 16) together, we obtain

I ε (ψ ε , δ ε ) ≥ C ω ε |Dψ ε | 2 dx + ω |µ ε | 2 dx -C ≥ C( Dψ ε 2 L 2 (ω) + µ ε 2 L 2 (ω) ) -C , (17)
with C > 0.

Let us now consider the force term F ε (ψ ε , δ ε ) given by formula [START_REF] Friesecke | Validity and failure of the Cauchy-Born hypothesis in a two-dimensional mass-spring lattice[END_REF]. For any full triangle T ε f , let

x ε i j = ε(is + jt) be the lower left vertex of T ε f and y ε i j = ε(is + jt + p) ∈ T ε
f . We use the fact that δ ε ∈ C(ε), hence piecewise constant and zero on empty triangles. Therefore we have

T ε f f (x) • δ ε (x) dx = ε 2 2 f (y ε i j ) • δ ε (y ε i j ) + T ε f ( f (x) -f (y ε i j )) • δ ε (x) dx, with | f (x) -f (y ε i j )| ≤ η( f , ε)
where η( f , •) denotes the modulus of continuity of the continuous function f . Therefore, it follows from the Cauchy-Schwarz inequality that

ε 2 ∑ (i, j)∈N ε 2 f (y ε i j ) • δ ε (y ε i j ) ≤ C δ ε L 2 (ω) , (18)
with

C = 2 f L 2 (ω) + η( f , 1)(meas ω) 1 2 .
We proceed similarly for the other terms. For instance, for the type 1 node terms, we write

T ε f f (x) • ψ ε (x) dx = ε 2 2 f (x ε i j ) • ψ ε (x ε i j ) + T ε f ( f (x) -f (x ε i j )) • ψ ε (x) dx + T ε f f (x ε i j ) • (ψ ε (x) -ψ ε (x ε i j )) dx. ( 19
)
Using the fact that

ψ ε (x) -ψ ε (x ε i j ) = Dψ ε (x)(x -x ε i j
) and using the Cauchy-Schwarz inequality again, we obtain

ε 2 ∑ (i, j)∈N ε 1 f (x ε i j ) • ψ(x ε i j ) + ∑ (i, j)∈N ε 2 f (y ε i j ) • ψ(x ε i j ) ≤ C ψ ε H 1 (ω) .
We finally obtain an estimate of the form

|F ε (ψ ε , δ ε )| ≤ C( ψ ε H 1 (ω) + δ ε L 2 (ω) ) = C( ψ ε H 1 (ω) + ε µ ε L 2 (ω) ). ( 20
)
We now note that ψ ε (x) = (x; 0) on a fixed subset of ω 0 by condition (4). Thus we apply Poincaré's inequality, with a constant that does not depend on ε and deduce from estimates [START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF] and ( 20) that there exists a constant C such that

ψ ε H 1 (ω) ≤ C and µ ε L 2 (ω) ≤ C. ( 21
)
The proof is concluded by noticing that δ ε = ε µ ε .

See [START_REF] Meunier | Elastic limit of square lattices with three point interactions[END_REF][START_REF] Raoult | Quasiconvex envelopes in nonlinear elasticity[END_REF] for a slightly different treatment of the force terms in a different context.

Note that, even though W ε vanishes on empty triangles, the fact that (ψ ε , δ ε ) ∈ V (ε) enabled us to recover uniform coercivity. Corollary 4.2 For any Γ-convergent subsequence, we have

Γ-lim ε→0 J ε (ψ, δ ) = +∞, if ψ / ∈ H 1 ω 0 (ω; R 3 ) or δ = 0.
Proof. Indeed, if Γ-lim ε→0 J ε (ψ, δ ) < +∞, it follows that there exists a sequence (ψ ε , δ ε ) → (ψ, δ ) in H that satisfies the hypothesis of Proposition 4.1. Therefore, we have δ = 0 and ψ ∈ H 1 (ω; R 3 ). Moreover, as ψ ε |ω 0 → (id; 0) strongly in L 2 loc (ω 0 ; R 3 ) by condition (4) and the fact that ψ ε is piecewise affine, we have that ψ ∈ H 1 ω 0 (ω; R 3 ).

Corollary 4.3

The minimizers of problem (13) for ε > 0 remain in a compact subset of H independent of ε.

Proof. Let (ϕ ε , γ ε ) be a minimizer. We let ψ ε (x) = (x; 0) for all x ∈ ω, and δ ε (x) = 0 in empty triangles, δ ε (x) = (ε p; 0), or exceptionally δ ε (x) = (ε(ps); 0), in full triangles. By the minimization property, we have

J ε (ϕ ε , γ ε ) ≤ J ε (ψ ε , δ ε ). By construction I ε (ψ ε , δ ε ) = 0. Moreover F ε (ψ ε , δ ε ) → 2 ω f (x)
• (x; 0) dx, thus is bounded. Therefore, (ϕ ε , γ ε ) satisfies the hypothesis of Proposition 4.1, and we apply Rellich's theorem to get the result. Corollary 4.3 shows that the strong topology of H is a topology in which the computation of the Γ-limit, if at all possible, will provide information on the asymptotic behavior of the minimizers.

The limit problem

Let us now introduce the function W 0 : R 2 × L (R 2 ; R 3 ) → R defined by

W 0 (y, g) = inf τ∈R 3 W ε (y, g, τ). ( 22
)
Note that this function no longer depends on ε. It is still Y -periodic and vanishes for y ∈ T e + L. We similarly define

Z 0 (x, g) = inf τ∈R 3 Z ε (x, g, τ). (23)
Even though we will not use it in the sequel, let us mention that

Z 0 (x, g) = 12κ |g(u)| 2 -1 + 2
, where u = s,t or ts depending on the orientation of the boundary triangle to which x belongs. Note that the infimum in both formulas ( 22) and ( 23) is attained.

We then define a reduced elastic energy I ε 0 : L 2 (ω; R 3 ) → R by

I ε 0 (ψ) = ω ε W 0 ε -1 x, Dψ(x) dx + ω ε ∂ Z 0 x, Dψ(x) dx, if ψ ∈ A(ε), I ε 0 (ψ) = +∞ if ψ ∈ L 2 (ω; R 3 ) \ A(ε).
Similarly, we define a reduced energy functional J ε 0 : L 2 (ω; R 3 ) → R by

J ε 0 (ψ) = I ε 0 (ψ) -F ε (ψ, 0) if ψ ∈ A(ε), J ε 0 (ψ) = +∞ if ψ ∈ L 2 (ω; R 3 ) \ A(ε).
The following proposition gives the connection between the functionals J ε and J ε 0 .

Proposition 4.4 For any subsequence such that both J ε and J ε 0 are Γ-convergent, we have,

(Γ-lim ε→0 J ε )(ψ, 0) = (Γ-lim ε→0 J ε 0 )(ψ), (24)
for all ψ ∈ H 1 ω 0 (ω; R 3 ), where the second Γ-limit is meant with respect to the strong topology of L 2 (ω; R 3 ).

Proof. We first show that

(Γ-lim ε→0 J ε )(ψ, 0) ≥ (Γ-lim ε→0 J ε 0 )(ψ). (25)
Let (ψ ε , δ ε ) → (ψ, 0) in H be a sequence that achieves the Γ-limit on the left of inequality (25). If the latter is +∞, there is nothing to prove. Let us assume that J ε (ψ ε , δ ε ) ≤ M for some M. Thus for all ε, (ψ ε , δ ε ) belongs to V (ε) and

J ε (ψ ε , δ ε ) ≥ J ε 0 (ψ ε ) -F ε (0, δ ε ), since F ε (ψ ε , δ ε ) = F ε (ψ ε , 0) + F ε (0, δ ε ).
Then, looking at estimate [START_REF] Raoult | Quasiconvex envelopes in nonlinear elasticity[END_REF] in the proof of Proposition 4.1, we see that

F ε (0, δ ε ) → 0. Thus, (Γ-lim ε→0 J ε )(ψ, 0) ≥ lim inf ε→0 (J ε 0 (ψ ε ) -F ε (0, δ ε )) ≥ (Γ-lim ε→0 J ε 0 )(ψ).
We now prove the reverse inequality. We can assume that ψ is such that (Γlim ε→0 J ε 0 )(ψ) < +∞. Let ψ ε ∈ A(ε) be sequence that achieves this Γ-limit. Since the infimum in formulas ( 22) and ( 23) is attained, we can construct δ ε ∈ C(ε) such that I ε (ψ ε , δ ε ) = I ε 0 (ψ ε ). Indeed, we simply define δ ε to take a minimizing value for W ε or Z ε that we select in any full triangle and δ ε = 0 in empty triangles. What remains to be seen to make sure that we are in C(ε) is that this construction gives rise to a function δ ε that satisfies conditions ( 5)- [START_REF] Braides | Continuum limits of discrete systems without convexity hypotheses[END_REF].

To do this, let us first consider the case of a type 2 node bonded to ω0 with three attached type 1 nodes in ω0 . In this case, Dψ ε = (id; 0) in the full triangle to which the type 2 node belongs and δ ε (x) = ε p achieves the minimum of W 0 , which is 0. Next, when there are only two attached type 1 nodes that both are in ω0 , we are dealing with a small triangle, and the energy Z ε only involves the directional derivative in direction u, where u is s, t or ts. This directional derivative is equal to u because both type 1 nodes are bonded. Therefore, the choices δ ε (x) = ε(ps) when u = ts and δ ε (x) = ε p when u = t or u = s achieve the minimum of Z 0 which is also equal to 0. Now we have

I ε (ψ ε , δ ε ) -F ε (ψ ε , 0) = I ε 0 (ψ ε ) -F ε (ψ ε , 0) = J ε 0 (ψ ε )
which is bounded by assumption. Arguing as in the proof of Proposition 4.1, we thus obtain that δ ε → 0 in L 2 . Consequently, F ε (0, δ ε ) → 0, see estimate [START_REF] Raoult | Quasiconvex envelopes in nonlinear elasticity[END_REF], and therefore

J ε (ψ ε , δ ε ) → (Γ-lim ε→0 J ε 0 )(ψ),
with (ψ ε , δ ε ) → (ψ, 0) in H, which completes the proof.

Let us now turn to the bulk of the proof, which is the computation of the Γlimit in the right-hand side of equation ( 24). For this, we introduce a homogenized energy density on L (R 2 ; R 3 ) defined by

W hom (g) = inf k∈N 1 k 2 inf θ ∈A(kY ) kY W 0 (y, g + Dθ (y)) dy , (26)
where A(kY ) denotes the set of continuous piecewise affine functions on the mesh defined on kY by Y +L and that vanish on ∂ (kY ). Finally, we define the functionals

J 0 (ψ) = ω W hom (Dψ(x)) dx -2 ω f • ψ dx if ψ ∈ H 1 ω 0 (ω; R 3 ) and J 0 (ψ) = +∞ if ψ ∈ L 2 (ω; R 3 ) \ H 1 ω 0 (ω; R 3 ), and
I 0 (ψ) = ω W hom (Dψ(x)) dx for ψ ∈ H 1 ω 0 (ω; R 3 ). Theorem 4.5 We have (Γ-lim ε→0 J ε 0 ) = J 0 .

The main convergence proof

The proof of Theorem 4.5 is long and technical and we break it into a series of lemmas. We borrow the global architecture of this proof from [START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF]. The detail is however quite different in places.

The force terms

Let us first deal with the force terms once and for all.

Proposition 4.6 Let ψ ε ∈ A(ε) be such that ψ ε → ψ strongly in L 2 (ω; R 3 ) and that is uniformly bounded in H 1 (ω; R 3 ). Then we have

F ε (ψ ε , 0) → 2 ω f • ψ dx.
Proof. We use the notation of the proof of Proposition 4.1. We have

F ε (ψ ε , 0) = ε 2 ∑ (i, j)∈N ε 1 f (x ε i j ) • ψ ε (x ε i j ) + ∑ (i, j)∈N ε 2 f (y ε i j ) • ψ ε (x ε i j ) .
Let us consider the type 1 nodes that are the lower left vertex of a full triangle T ε f and denote by M ε 1 ⊂ N ε 1 the set of indices corresponding to such nodes. For the corresponding terms, we have (viz. equation ( 19) reproduced here for the reader's convenience)

T ε f f (x) • ψ ε (x) dx = ε 2 2 f (x ε i j ) • ψ ε (x ε i j ) + T ε f ( f (x) -f (x ε i j )) • ψ ε (x) dx + T ε f f (x ε i j ) • (ψ ε (x) -ψ ε (x ε i j )) dx. Now T ε f ( f (x) -f (x ε i j )) • ψ ε (x) dx ≤ η( f , ε) T ε f |ψ ε (x)| dx, so that ∑ (i, j)∈M ε 1 T ε f ( f (x) -f (x ε i j )) • ψ ε (x) dx ≤ η( f , ε) ω |ψ ε (x)| dx ≤ Cη( f , ε) ψ ε L 2 (ω) → 0 when ε → 0. Similarly T ε f f (x ε i j ) • (ψ ε (x) -ψ ε (x ε i j )) dx ≤ Cε T ε f |Dψ ε (x)| dx so that ∑ (i, j)∈M ε 1 T ε f f (x ε i j ) • (ψ ε (x) -ψ ε (x ε i j )) dx ≤ Cε ω |Dψ ε (x)| dx ≤ Cε Dψ ε L 2 (ω) → 0 when ε → 0. Finally ∑ (i, j)∈M ε 1 T ε f f (x) • ψ ε (x) dx → 1 2 ω f (x) • ψ(x) dx since ψ ε → ψ strongly in L 2 and f 1 ∪T ε f f 2 weakly in L 2 .
It is clear that the type 1 nodes that are not of the previous kind have a negligible contribution in the limit ε → 0.

We proceed in the same fashion for the node 2 terms.

We can thus from now on concentrate on the elastic energy terms only.

Bound from below

We obtain the bound from below, i.e., condition i) of the definition of Γ-convergence in a series of Lemmas. The bound is obtained in increasing generality, starting with ψ affine in Lemmas 4.8 and 4.10, then ψ piecewise affine in Lemma 4.11 and finally general ψ in Proposition 4.13. The first step of the proof is actually a bound from above in a special case.

Lemma 4.7 Let U be a bounded open regular subset of R 2 and η > 0. Let g ∈ L (R 2 ; R 3 ), a ∈ R 3 and χ(x) = gx + a. For any sequence ε → 0, there exists a sequence χ ε such that χ ε is piecewise affine on

T ε ∩ U, χ ε → χ strongly in L 2 (U; R 3 ), χ ε -χ ∈ H 1 0 (U; R 3 ) and lim sup U W 0 ε -1 x, Dχ ε (x) dx ≤ U W hom (Dχ(x)) dx + η.
Proof. Let us choose k ∈ N and θ ∈ A(kY ) such that

W hom (g) ≤ 1 k 2 kY W 0 (y, g + Dθ (y)) dy ≤ W hom (g) + η meas (U)
.

Let U ε k be the union of all εkL-translates of the cell εkY included in U. Clearly, meas (U \U ε k ) → 0 when ε → 0. We extend θ to R 2 by kY -periodicity and set

χ ε (x) = χ(x) + εθ ε -1 x if x ∈ U ε k , χ ε (x) = χ(x)
otherwise. By construction, χ ε is piecewise affine on the triangulation restricted to U. It is also obvious that χ ε → χ strongly in L 2 (U; R 3 ) and that χ εχ ∈ H 1 0 (U; R 3 ). We have

U W 0 ε -1 x, Dχ ε (x) dx = U W 0 ε -1 x, g + Dθ ε -1 x dx = U ε k W 0 ε -1 x, g + Dθ ε -1 x dx + U\U ε k W 0 ε -1 x, g dx. Now U\U ε k W 0 ε -1 x, g dx ≤ C meas (U \U ε k ) → 0 when ε → 0. Moreover U ε k W 0 ε -1 x, g + Dθ ε -1 x dx = ∑ εkY cells ε 2 kY W 0 y, g + Dθ (y) dy ≤ meas (U ε k ) W hom (g) + η meas (U) ,
hence the Lemma by letting ε → 0.

We now start on the bound from below strictly speaking.

Lemma 4.8 Let O be a bounded open subset of R 2 , let g ∈ L (R 2 ; R 3 ), a ∈ R 3 and ψ(x) = gx + a. For any sequences ε → 0, ψ ε such that ψ ε is piecewise affine on

T ε ∩ O, ψ ε → ψ strongly in L 2 (O; R 3 ) and ψ ε -ψ ∈ H 1 0 (O; R 3 ), we have lim inf O W 0 ε -1 x, Dψ ε (x) dx ≥ O W hom (Dψ(x)) dx. (27)
Proof. Let M > 0 be such that Ō ⊂ M Y + x 0 for some x 0 ∈ R 2 . Let us assume without loss of generality that x 0 = 0. For all ε, we choose the largest natural integer k(ε) such that k(ε)εY ⊂ MY . For ε small enough, we have Ō ⊂ k(ε)εY , moreover meas (MY \ k(ε)εY ) → 0 when ε → 0. Let η > 0 be given and χ ε be the sequence given by Lemma 4.7 applied on the open set U = M Y \ Ō. We define

ψ ε = ψ ε in O and ψ ε = χ ε in MY \ Ō.
We have

O W 0 ε -1 x, Dψ ε (x) dx = I ε 1 -I ε 2 + I ε 3 , (28)
where

I ε 1 = k(ε)εY W 0 ε -1 x, D ψ ε (x) dx, I ε 2 = MY \O W 0 ε -1 x, D ψ ε (x) dx, I ε 3 = MY \k(ε)εY W 0 ε -1 x, D ψ ε (x) dx.
Due to the boundary conditions and the construction of χ ε , we have that ψ εψ ∈ H 1 0 (k(ε)εY ; R 3 ). By rescaling, the function

θ ε (y) = ε -1 ψ ε (εy) -ψ(y) belongs to H 1 0 (k(ε)Y ; R 3
) and is piecewise affine on the mesh defined by L. It is thus a competing test-function in the definition (26) of the homogenized density W hom . Therefore, we have

I ε 1 = ε 2 k(ε)Y W 0 y, g + Dθ ε (y) dy ≥ ε 2 k(ε) 2 W hom (g) = k(ε)εY W hom (Dψ(x)) dx.
(29) Moreover, by construction

lim sup I ε 2 ≤ MY \O W hom (Dψ(x)) dx + η. ( 30
)
Finally, since meas (MY \ k(ε)εY ) → 0 and the differential of ψ ε is controlled in this set, we have

I ε 3 → 0 (31)
when ε → 0. Putting equation ( 28) and estimates (29), ( 30) and ( 31) together, we obtain the Lemma.

We now need to obtain the bound from below (27) without the homogeneous boundary condition. This is classically done by using De Giorgi's slicing method. The slicing method does not work here because we need piecewise affine functions, and multiplying a piecewise affine function by a cut-off function destroys its piecewise affine character. We thus introduce a discrete version of the slicing argument, see also [START_REF] Alicandro | A general integral representation result for continuum limits of discrete energies with superlinear growth[END_REF]. First a preparatory lemma. Lemma 4.9 Let Π ε denote the piecewise affine interpolation operator associated with the triangulation T ε . There exists a constant C independent of ε such that for all finite unions U ε of triangles of T ε , all R 3 -valued functions ψ that are piecewise affine on T ε and all [0, 1]-valued functions θ , we have

Π ε (θ ψ) L 2 (U ε ;R 3 ) ≤ C ψ L 2 (U ε ;R 3 ) . (32) Moreover, if θ is W 1,∞ , we have U ε D(Π ε (θ ψ)) 2 dx ≤ C U ε |Dψ| 2 + Dθ 2 L ∞ |ψ| 2 dx . ( 33
)
Proof. Let T ε be a generic triangle in T ε with vertices S ε i . We denote by λ ε i the associated barycentric coordinates. In T ε , we thus have

Π ε (θ ψ)(x) = 3 ∑ i=1 λ ε i (x)θ (S ε i )ψ(S ε i ).
We introduce two functions Q ε 1 :

T ε ×[0, 1] 3 ×(R 3 ) 3 → R + and Q ε 2 : T ε ×(R 3 ) 3 → R + by Q ε 1 (T ε , θ i , v i ) = 3 ∑ i=1 λ ε i (x)θ i v i 2 L 2 (T ε ;R 3 ) , and
Q ε 2 (T ε , v i ) = 3 ∑ i=1 λ ε i (x)v i 2 L 2 (T ε ;R 3 )
.

Clearly, we have

Π ε (θ ψ) 2 L 2 (U ε ;R 3 ) = ∑ T ε ∈U ε Q ε 1 (T ε , θ (S ε i ), ψ(S ε i )), ψ 2 L 2 (U ε ;R 3 ) = ∑ T ε ∈U ε Q ε 2 (T ε , ψ(S ε i )).
Now every generic triangle is of the form x ε + εT f or x ε + εT e for some x ε ∈ L (for brevity, we assume the former), so that

Q ε 1 (T ε , θ i , v i ) = ε 2 3 ∑ i=1 λ i (x)θ i v i 2 L 2 (T f ;R 3 ) = ε 2 Q 1 1 (T f , θ i , v i ),
and

Q ε 2 (T ε , v i ) = ε 2 3 ∑ i=1 λ ε i (x)v i 2 L 2 (T f ;R 3 ) = ε 2 Q 1 2 (T f , v i ). Both functions Q 1 1 (T f , •, •) and Q 1 2 (T f , •) are continuous and quadratic with respect to (v i ) ∈ (R 3 ) 3 . Moreover, Q 1 2 (T f , v i ) = 0 if and only if v i = 0. Therefore, if we set C 2 = max θ i ∈[0,1] 3 ∑ |v i | 2 =1 Q 1 1 (T f , θ i , v i ) Q 1 2 (T f , v i )
we obtain estimate (32).

For estimate (33), we work on the reference scaled triangle εT f . For any affine function χ, we have

|Dχ| ≤ C ε |Dχ(εs)| + |Dχ(εt)|
for some constant C independent of ε. We apply this to χ = Π ε (θ ψ). We have

Dχ(εs) = Π ε (θ ψ)(εs) -Π ε (θ ψ)(0) = (θ ψ)(εs) -(θ ψ)(0) = θ (εs)(ψ(εs) -ψ(0)) + (θ (εs) -θ (0))ψ(0) = εθ (εs)Dψ(s) + (θ (εs) -θ (0))ψ(0). Therefore |D(Π ε (θ ψ))(εs)| ≤ ε(|Dψ(s)| + Dθ L ∞ |ψ(0)|) and thus |D(Π ε (θ ψ))| ≤ C(|Dψ| + Dθ L ∞ |ψ(0)|).
To conclude, it suffices to show that for any affine function ψ, we have

|ψ(0)| 2 ≤ C T f |ψ(x)| 2 dx
as the conclusion will follow by a simple scaling. But this is obvious, for if it were not true, we would have a sequence ψ n of affine functions tending to 0 in L 2 (T f ), but such that |ψ n (0)| = 1. This is impossible since we are in a finite dimensional space, in which L 2 convergence implies pointwise convergence.

We now perform the slicing step to establish estimate (27) without boundary conditions. Lemma 4.10 Let O be a bounded open subset of R 2 , let g ∈ L (R 2 ; R 3 ), a ∈ R 3 and ψ(x) = gx + a. For any sequences ε → 0, ψ ε such that ψ ε is piecewise affine on

T ε ∩ O, ψ ε is bounded in H 1 (O; R 3 ) and ψ ε → ψ strongly in L 2 (O; R 3 ), we have lim inf O W 0 ε -1 x, Dψ ε (x) dx ≥ O W hom (Dψ(x)) dx. Proof. Let N be an integer. Let O 0 O and r = d(O 0 , O) > 0. For k = 1, . . . , 2N + 1, we let O k = {x ∈ O, d(x, O 0 ) < kr 2N+1 }. For each i = 0, . . . , N -1, we pick a Lipschitz function θ i such that 0 ≤ θ i ≤ 1, θ i = 1 on O 2i+1 , θ i = 0 on O \ O 2i+3 and |Dθ i | ≤ C(2N+1) r .
We let

ψ ε i = ψ + Π ε θ i (ψ ε -ψ) . By construction, ψ ε i is piecewise affine on T ε ∩ O and such that ψ ε i -ψ ∈ H 1 0 (O; R 3 )
. By Lemma 4.9, ψ ε i → ψ strongly in L 2 (O; R 3 ). We can thus apply Lemma 4.8 to conclude that lim inf

O W 0 ε -1 x, Dψ ε i (x) dx ≥ O W hom (Dψ(x)) dx.
For ε small enough, we have that

Π ε θ i (ψ ε -ψ) = ψ ε -ψ in O 2i and that Π ε θ i (ψ ε -ψ) = 0 in O \ O 2i+3 . Therefore O W 0 ε -1 x, Dψ ε i (x) dx = O 2i W 0 ε -1 x, Dψ ε (x) dx + O 2i+3 \O 2i W 0 ε -1 x, Dψ ε i (x) dx + O\O 2i+3 W 0 ε -1 x, g dx ≤ O W 0 ε -1 x, Dψ ε (x) dx + O 2i+3 \O 2i W 0 ε -1 x, Dψ ε i (x) dx + O\O 2i+3 W 0 ε -1 x, g dx
For the last term, we have that

O\O 2i+3 W 0 ε -1 x, g dx ≤ C meas (O \ O 0 ).
For the second term, we see that all the triangles in which ψ ε i = ψ ε and ψ ε i = ψ are included in O 2i+3 \ O 2i for ε small enough, hence we use Lemma 4.9 and obtain

O 2i+3 \O 2i W 0 ε -1 x, Dψ ε i (x) dx ≤ C O 2i+3 \O 2i |g| 2 +|Dψ ε | 2 + (2N + 1) 2 r 2 |ψ ε -ψ| 2 dx.
Summing on the N slices, we thus obtain

1 N N ∑ i=0 O W 0 ε -1 x, Dψ ε i (x) dx ≤ O W 0 ε -1 x, Dψ ε (x) dx + C N O\O 0 (|g| 2 + |Dψ ε | 2 ) dx + CN r 2 O\O 0 |ψ ε -ψ| 2 dx +C meas (O \ O 0 ).
Now when ε → 0, we have O\O 0 |ψ ε -ψ| 2 dx → 0 by hypothesis. The corresponding term thus disappears and we are left with

O W hom (Dψ(x)) dx ≤ lim inf O W 0 ε -1 x, Dψ ε (x) dx + C N +C meas (O \ O 0 ).
We conclude by letting first N → +∞ and then O 0 → O.

The next Lemmas are the last steps in the proof of the bound from below.

Lemma 4.11 Let ψ ∈ H 1 ω 0 (ω; R 3 ) be piecewise affine. Then, for all sequences ψ ε in L 2 (ω; R 3 ) such that ψ ε → ψ strongly in L 2 (ω; R 3 ), we have that

lim inf I ε 0 (ψ ε ) ≥ I 0 (ψ).
Proof. If the left-hand side is +∞, there is nothing to prove. Hence we can assume that I ε 0 (ψ ε ) is bounded, which entails that ψ ε is bounded in H 1 and piecewise affine on T ε and we have

I ε 0 (ψ ε ) ≥ ω ε W 0 ε -1 x, Dψ ε (x) dx = p ∑ i=1 ω i ∩ω ε W 0 ε -1 x, Dψ ε (x) dx,
where ω i denotes a partition of ω in each part of which ψ is affine. By Lemma 4.10, we have lim inf

ω i ∩ω ε W 0 ε -1 x, Dψ ε (x) dx ≥ ω i W hom (Dψ(x)) dx,
for all ω i ω i , hence the Lemma.

To treat the case of a general limit function ψ, we need an intermediate technical result. Lemma 4.12 The function W hom is locally Lipschitz.

Proof. We first show that W 0 is locally Lipschitz. Let l 1 , l 2 ∈ L (R 2 ; R 3 ) and τ ∈ R 3 be such that W 0 (y, l 2 ) = W 1 (y, l 2 , τ). We have

W 0 (y, l 1 ) -W 0 (y, l 2 ) ≤ W 1 (y, l 1 , τ) -W 1 (y, l 2 , τ) ≤ C(1 + |l 1 | + |l 2 |)|l 1 -l 2 |.
Let us now take g, h ∈ L (R 2 ; R 3 ) and define ψ(x) = gx and χ(x) = hx. By Lemma 4.7, for all η > 0, there exists a subsequence ψ ε ∈ H 1 (Y ; R 3 ) piecewise affine on T ε and such that Dψ ε is bounded in L 2 , ψ ε → ψ strongly in L 2 and Y W 0 ε -1 x, Dψ ε (x) dx → W hom (g) + δ for some 0 ≤ δ ≤ η. We let

χ ε = h + ψ ε -g. Clearly, χ ε is bounded in H 1 , χ ε → h strongly in L 2 ,
and χ ε is piecewise affine on T ε . We thus have lim inf

Y W 0 ε -1 x, Dχ ε (x) dx ≥ W hom (h),
by Lemma 4.10. Therefore

W hom (h) -W hom (g) ≤ lim inf Y W 0 ε -1 x, Dχ ε (x) dx -lim Y W 0 ε -1 x, Dψ ε (x) dx + δ . Now since Y W 0 ε -1 x, Dχ ε (x) -W 0 ε -1 x, Dψ ε (x) dx ≤ C Y (1 + |Dψ ε | + |Dχ ε |)|Dψ ε -Dχ ε | dx ≤ C|g -h|
by Cauchy-Schwarz and the H 1 bound, we obtain the Lemma by letting η → 0.

We now are in a position to conclude the bound from below. Proposition 4.13 Let ψ ∈ H 1 ω 0 (ω; R 3 ). Then, for all sequences ψ ε in L 2 (ω; R 3 ) such that ψ ε → ψ strongly in L 2 (ω; R 3 ), we have that

lim inf I ε 0 (ψ ε ) ≥ I 0 (ψ).
Proof. Let P ε denote the orthogonal H 1 (ω; R 3 )-projection on A(ε), which is a finite dimensional, hence closed affine subspace. By standard finite element theory arguments, we have that P ε ψ → ψ strongly in H 1 when ε → 0. We pick a sequence ψ k of piecewise affine functions in H 1 ω 0 (ω; R 3 ) that converges strongly to ψ. For k fixed, we thus also have

P ε ψ k → ψ k strongly in H 1 when ε → 0. It follows that ψ ε -P ε ψ + P ε ψ k → ψ k strongly in L 2 (ω; R 3 ), so that by Lemma 4.11 lim inf I ε 0 (ψ ε -P ε ψ + P ε ψ k ) ≥ I 0 (ψ k ).
Again, we can assume that ψ ε ∈ A(ε) and ψ ε is bounded in H 1 , otherwise there is nothing to prove and thus both ψ ε and ψ ε -P ε ψ + P ε ψ k belong to A(ε). We can thus write

I ε 0 (ψ ε ) = I ε 0 (ψ ε -P ε ψ + P ε ψ k ) + I ε 0 (ψ ε ) -I ε 0 (ψ ε -P ε ψ + P ε ψ k ),
with

I ε 0 (ψ ε ) -I ε 0 (ψ ε -P ε ψ + P ε ψ k ) ≤ ω ε W 0 ε -1 x, Dψ ε (x) -W 0 ε -1 x, D(ψ ε -P ε (ψ -ψ k )) dx + ω ε ∂ Z 0 x, Dψ ε (x) -Z 0 x, D(ψ ε -P ε (ψ -ψ k )) dx ≤ C 1 + Dψ ε L 2 + DP ε (ψ -ψ k ) L 2 DP ε (ψ -ψ k ) L 2 ≤ C ψ -ψ k H 1 ,
by Cauchy-Schwarz and since an orthogonal projection has norm 1. It follows from the previous estimates that

lim inf I ε 0 (ψ ε ) ≥ I 0 (ψ k ) -C ψ -ψ k H 1 .
We now let k → +∞. Since W hom is locally Lipschitz, it is continuous and obviously with quadratic growth at most, therefore I 0 (ψ k ) → I 0 (ψ), while ψψ k H 1 → 0, and the proof is complete.

Recovery sequence

In order to complete the proof of Theorem 4.5, we need to establish condition ii) or the bound from above or recovery sequence, in the definition of Γ-convergence.

For this we first need a refinement of Lemma 4.7.

Lemma 4.14 Let ψ ∈ H 1 ω 0 (ω; R 3 ) be piecewise affine on ω and let η > 0. There exists a sequence ψ ε such that ψ ε → ψ strongly in L 2 (ω; R 3 ) and

lim sup I ε 0 (ψ ε ) ≤ I 0 (ψ) + η.
Proof. Because of the finite energy bound, the sequence ψ ε to be constructed must be in A(ε) that is to say piecewise affine functions on T ε ∩ ω such that ψ ε = A ε (ψ ε |ω ε ) plus the boundary conditions. Let ω i , i = 1, . . . , p, be a partition of ω such that ψ is affine on ω i with differential g i . Since ψ is affine on ω 0 , we may assume that ω0 ⊂ ω1 without loss of generality. As in Lemma 4.7, for each i, we choose k i ∈ N and θ

i ∈ A(k i Y ) such that W hom (g i ) ≤ 1 k 2 i k i Y W 0 (y, g i + Dθ i (y)) dy ≤ W hom (g i ) + η meas (ω) .
Let ω ε k i be the union of all εk i L-translates of the cell εk i Y included in ω ε ∩ ω i from which we remove all the boundary triangles. Clearly, meas (ω i \ ω ε k i ) → 0 when ε → 0. Again we set

ψ ε (x) = ψ(x) + εθ i x ε if x ∈ ω ε k i , ψ ε (x) = Π ε ψ(x), if x ∈ ω ε \ p i=1 ω ε k i , and
ψ ε = A ε (Π ε ψ |ω ε )
on the rest of ω.

Let us check that the condition of place is satisfied. Given that ω0 ⊂ ω1 , for all x in ω ε k 1 , we see that the choice θ 1 = 0 achieves the minimum since 0 ≤ W hom (g 1 ) ≤

1 k 2 1
k 1 Y W 0 (y, g 1 ) dy = 0. Thus ψ ε (x) = (x; 0) on ω ε k 1 . On ( ωε ∩ ω0 ) \ ω ε k 1 , we also have ψ ε (x) = (x; 0) by the second part of the definition of ψ ε . In particular, at all type 1 nodes belonging to ω0 , we have ψ ε (x) = (x; 0) and ψ ε satisfies the condition of place.

By construction, we thus have ψ ε is in A(ε) and ψ ε → ψ strongly in L 2 (ω; R 3 ).

We have

I ε 0 (ψ ε ) = p ∑ i=1 ω ε k i W 0 x ε , g i + Dθ i x ε dx + ω ε \ p i=1 ω ε k i W 0 x ε , D Π ε ψ(x) dx + ω ε ∂ Z 0 (x, DA ε (Π ε ψ |ω ε )) dx.
Now D Π ε ψ and DA ε (Π ε ψ |ω ε ) are both uniformly controlled and thus

ω ε \ p i=1 ω ε k i W 0 x ε , D Π ε ψ(x) dx + ω ε ∂ Z 0 (x, DA ε (Π ε ψ |ω ε )) dx → 0 when ε → 0. Finally ω ε k i W 0 x ε , g i + Dθ i x ε dx ≤ meas (ω ε k i ) W hom (g i ) + η meas (ω) ,
hence the Lemma by letting ε → 0.

Lemma 4.15 Let ψ ∈ H 1 ω 0 (ω; R 3 ) and let η > 0. There exists a sequence ψ ε such that ψ ε → ψ strongly in L 2 (ω; R 3 ) and

lim sup I ε 0 (ψ ε ) ≤ I 0 (ψ) + η.
Proof. We proceed by piecewise affine approximation as in the proof of Proposition 4.13 using the bound given by Lemma 4.14.

by letting τ = -τ. Thus W 0 (y, gS) = W 0 (y, g).

The two rotations R and S generate the group C 6 , which is thus a subgroup of the material symmetry group of W 0 .

For k ∈ N * , we let

W kY (g) = 1 k 2 inf θ ∈A(kY ) kY W 0 (y, g + Dθ (y)) dy. ( 34
)
The lack of symmetry of the integer multiples of the unit cell as well as numerical evidence show that we cannot expect W kY to be right-R invariant for k fixed. We will expand on this in the next section. Right-R invariance is however recovered in the k → +∞ limit as we now proceed to show.

Proposition 5.2 The material symmetry group of W hom contains the circular group C 6 .

Proof. For all q ∈ N * , we consider the hexagons H q of maximal size inscribed in kY with k = 2q, see Figure 7 below. We let A(H q ) be the set of piecewise affine functions on the lattice that vanish on ∂ H q and we likewise define

W H q (g) = 1 3q 2 inf
θ ∈A(H q ) H q W 0 (y, g + Dθ (y)) dy.

H q y q 2qY

Figure 7: The hexagon H q . Also figured is the microstructure around the center y q to make R-invariance clear.

Let H be the hexagon of unit area such that H q = √ 3q H. We have

W H q (g) = inf ψ∈A(( √ 3q) -1 ) H W 0 (( √ 3q) -1 x, Dψ(x)) dx,
where A(( √ 3q) -1 ) denotes the set of piecewise affine continuous functions on the scaled lattice on H that satisfy the boundary condition of place ψ(x) = g(x) on ∂ H.

We use the general Γ-convergence result of Proposition 4.16, which also applies for this boundary condition and deduce that

W H q (g) → inf θ ∈H 1 0 (H;R 3 )
H W hom (g + Dθ (x)) dx when q → +∞.

By Corollary 4.17, we know that W hom is quasiconvex, which means that inf

θ ∈H 1 0 (H;R 3 ) H W hom (g + Dθ (x)) dx = W hom (g).
We thus have an alternate representation formula for the homogenized energy density,

W hom (g) = lim q→+∞ 1 3q 2 inf θ ∈A(H q ) H q W 0 (y, g + Dθ (y)) dy .
Now it is apparent from Figure 7 that H q is invariant by rotation of R around its center node y q , and moreover that we have W 0 (y, g) = W 0 (R(yy q ) + y q , g) = W 0 (R(yy q ) + y q , gR) by Proposition 5.1 for the second equality. Therefore, for all θ ∈ A(H q ), we have H q W 0 (y, gR + Dθ (y)) dy = H q W 0 (R(yy q ) + y q , (g + Dθ (y)R -1 )R) dy = H q W 0 (z, (g + D θ (z))R) dz with the change of variables z = R(yy q ) + y q , θ (z) = θ (y) = H q W 0 (z, g + D θ (z)) dz by Proposition 5.1 again. Therefore W H q (gR) = W H q (g) hence W hom (gR) = W hom (g).

Right-S invariance can be seen directly on the integer multiples of Y . Indeed, for all θ ∈ A(kY ), we have Again, R and S generate C 6 , which completes the proof.

It is very likely that the material symmetry group of W hom is exactly C 6 . We will give some numerical evidence in the next section. We now turn to other properties of interest of W 0 and W hom .

Proposition 5.3 Let g ∈ L (R 2 ; R 3 ). We denote by A 1 , A 2 and A 3 the images of the three vertices of the reference triangle T f by the linear mapping g. We have i) W 0 (0, g) = 0 if and only if the radius r g of the circumcircle of the triangle A 1 A 2 A 3 is such that r g ≤ 1. If r g < 1, then the infimum in the definition of W 0 is attained at a vector τ that is out of the range of g.

ii) W 0 (0, g) > 0 if and only if the radius r g of the circumcircle of the triangle A 1 A 2 A 3 is such that r g > 1. The infimum in the definition of W 0 is attained at a vector τ that belongs to the range of g.

In the degenerate cases, we consider that the radius r g is equal to 0

if A 1 = A 2 = A 3 , |A 1 -A 3 | if A 1 = A 2 = A 3
and the permutations thereof, and +∞ if the three points are aligned and distinct.

Proof. We only deal with the nondegenerate cases, i.e., when the three points are affinely independent. The remaining cases are left to the reader.

It is easier to go back to the initial formulations in terms of nodes, so that A i denote the position of the three reference type 1 nodes in space, and we let M denote the position in space of the reference type 2 node attached to the three type 1 nodes. We thus seek to minimize the function

F : M → 3 ∑ i=1 (|M -A i | -1) 2
over R 3 and recover τ = M -A 1 afterwards.

It is clear that the minimum is attained. Assume first that r g ≤ 1. We can thus pick α such that α 2 + r 2 g = 1. Placing x on the normal to the plane spanned by the triangle and passing through the circumcenter, at a distance ±α of the plane, we see that F(x) = 0. Thus in this case, W 0 (0, g) = 0. Moreover, if r g < 1, then α = 0 and the two minimizing type 2 node positions are out of plane.

Assume now that r g > 1 and that the minimum is attained out of plane. The vectors (M -A i ) i=1,2,3 are thus linearly independent. Moreover, since F is differentiable out of the plane, it follows that ∇F(M) = 0. This means that 0

= 3 ∑ i=1 |M -A i | -1 |M -A i | (M -A i ).
But at least one of the terms |M -A i | -1 is nonzero, since r g > 1, contradiction. Therefore, the minimum is attained in the plane, and since there is no point in the plane at distance 1 of all three vertices, it follows that W 0 (0, g) > 0 in this case.

Corollary 5.4 Let g ∈ L (R 2 ; R 3 ) be such that r g ≤ 1. Then W hom (g) = 0.

The meaning of Corollary 5.4 is that the graphene sheet shares with nonlinear membranes a property of degeneracy under compression. However, compression for a nonlinear membrane is expressed in terms of the singular values of the deformation gradient being smaller than one, see [START_REF] Dret | The nonlinear membrane model as variational limit of nonlinear three-dimensional elasticity[END_REF], which is not the case for the condition on the circumcircle radius, which is a lot more subtle. In fact, we have that W hom vanishes on the rank-1-convex hull of the set of g such that r g ≤ 1.

Corollary 5.5 The function W 0 is not convex with respect to g. Proof. Indeed, the set of linear mappings g such that r g ≤ 1 is not convex. Consider for instance the segment g λ = (λ id + (1λ ) j; 0) for λ ∈ [0, 1], where id(s) = s, id(t) = t and j(s) = -t and j(t) = -s. It is easy to check that r g λ > 1 for all λ ∈ ] 1 3 , 2 3 [ and r g λ ≤ 1 elsewhere on the segment.

6 Numerical study and remarks on the Cauchy-Born rule

The representation formula for W hom is not very explicit, and neither is that for W 0 . Therefore, it is not easy to obtain quantitative information on W hom . We thus resort to numerical simulation and numerical observations lead us to a few remarks concerning the Cauchy-Born rule.

Let us briefly recall that the Cauchy-Born rule a priori concerns atomic crystals that are constrained by their own dimension, for instance a 2d crystal undergoing 2d deformations, which is not our case here. There is however a connection, as we will see shortly, both in the 2d-3d and 2d-2d cases. In the case of a complex lattice, such as here, the Cauchy-Born rule dictates that when a crystal is submitted to a homogeneous boundary condition of place ψ(x) = g(x) with g linear, the type 1 nodes should globally deform according to g and the type 2 nodes should relax the local energy inside each cell, see [START_REF] Ming | Cauchy-Born rule and the stability of crystalline solids: Static problems[END_REF][START_REF] Friesecke | Validity and failure of the Cauchy-Born hypothesis in a two-dimensional mass-spring lattice[END_REF]. The second part of the Cauchy-Born rule is exactly what W 0 does.

Interestingly enough, we observe numerically that for some g, the Cauchy-Born rule seems to apply in our case, whereas for other values of g, it appears to fail. We should emphasize that the hexagonal lattice is not especially engineered toward exhibiting this behavior, which emerges naturally.

We thus compute W kY (g) for various values of k and g. Here also, it is more convenient to work directly on the initial discrete node formulation rather than on the continuous formulation, especially since we do not have a closed form expression for W 0 . Energy minimization is performed by using the Polak-Ribière nonlinear conjugate gradient algorithm, which works quite well. Since the problem is non convex, we use random initial conditions to start the algorithm in order to avoid getting trapped at a local minimum as much as possible. All computations are performed with Scilab (http://www.scilab.org/). We use orthogonal Cartesian coordinates instead of the oblique system used in the theoretical discussion and 3 × 2 matrices F as deformation gradients instead of differentials g ∈ L (R 2 ; R 3 ) for numerical convenience. The value of κ is always 1. In this case, it is observed that W kY (F) does not depend on k and that the equilibrium configuration for each k is the periodic repetition of the Cauchy-Born configuration in one cell, where the four type 1 nodes follow F and the type 2 node assumes the position that minimizes W 0 in the plane. The numerical results for such values of the gradient F strongly suggest that W hom (F) = W 0 (y, F) where y is any point in the reference full triangle T f . Note that the initial condition for the gradient algorithm is always random and out of plane.

Next we show a non Cauchy-Born configuration that is achieved for Several interesting things are observed in this case. First of all, most type 2 nodes appear to be inplane with their attached type 1 nodes, indicating a strictly positive value of W 0 in most cells. The type 1 nodes however do not follow F in the interior but arrange themselves in 3d folding patterns as seen on Figure 9. Finally, W kY (F) is strictly decreasing with respect to k as shown in Figure 10. This indicates that taking the infimum with respect to k in formula ( 26) is necessary and that we have W hom (F) < W 0 (y, F) for this value of F, as opposed to the case considered in the previous computation. It should be noted that the necessity of minimizing over integer multiples of the unit cell in nonconvex homogenization was first proved in [START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF] in the continuous case. Recent examples in the continuous and discrete cases were given in [START_REF] Barchiesi | New counterexamples to the cell formula in nonconvex homogenization[END_REF]. As a side note regarding the Cauchy-Born rule in the context of atomic crystals, we observe that the same computation with the same F restricted to planar deformations yields a planar non Cauchy-Born configuration. Hence, the Cauchy-Born rules also appears to fail for a 2d-hexagonal crystal, see [START_REF] Friesecke | Validity and failure of the Cauchy-Born hypothesis in a two-dimensional mass-spring lattice[END_REF].

F =   1 -2 √ 3 3
We now test material symmetry. Figure 11 below shows W kY (FQ(θ )) with the same F as before, Q(θ ) ∈ SO(2) of angle θ ∈ [0, π], plotted against θ (this is for k = 50). Again, an interesting Cauchy-Born vs. non Cauchy-Born phenomenon is observed. For energies below approximately 0.13, the minimizing node configuration is the planar Cauchy-Born configuration. Hence, the corresponding parts of the graph do not depend on k (and this is a strong indication that the material symmetry group is exactly C 6 ). This behavior is presumably related to the low-energy Cauchy-Born regime evidenced in [START_REF] Conti | Sufficient conditions for the validity of the Cauchy-Born rule close to SO(n)[END_REF] for crystals confined in their own dimension. See also the remarks on material symmetry made in [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF], section 4.3.

In the Cauchy-Born regime of the previous plot, the convergence of the algorithm is fairly fast, with between 440 and 500 iterations of gradient. Energies above 0.13 exhibit non Cauchy-Born, 3d folding patterned minimizing configurations, and dependence on k. The convergence of the algorithm is dramatically slower, with a number of gradient iterations between 2000 and 5000.

Also shown on the plot are the values of the energy for θ = 0, π/3, 2π/3 and π, which are slightly different: a difference of the order of 2 • 10 -3 is observed (also decreasing with k), which shows that W kY is not Q(π/3) nor Q(2π/3) rightinvariant, as expected. Here is a closeup of Figure 11 We have seen that W 0 is not convex, cf. Corollary 5.5, and there is plenty of numerical evidence that neither is W kY . The question of whether or not W hom is convex remains unanswered. All the numerical tests we performed tend to point to a convexification when k → +∞, even on segments the endpoints of which differ by a rank-two matrix. Even for good candidates for being points of non convexity, such as the points corresponding to the local maxima in Figure 11, if we plot the energies on the tangent to the curve θ → FQ(θ ), we obtain the behavior shown below in Figure 12. The small irregularities in some of the curves are execution dependent, which means that they are presumably due to the random initial conditions and the presence of several local minima or almost minima. They tend to disappear as k increases and the general features of the set of curves are nonetheless always the same. It would thus not be a huge surprise if W hom turned out to be convex.

Figure 1 :

 1 Figure 1: •: type 1 nodes, •: type 2 nodes

3 2 1 Figure 2 :

 312 Figure 2: Hexagonal structure and the three different types of bars

ω ω 0 Figure 3 .

 03 Figure 3. Nodes in ω and boundary condition of place in ω 0 .

Figure 4 .

 4 Figure 4. Extending functions to ω. Triangles added in lighter dashed lines."Small" triangles of ω ε ∂ in gray, see below. Proposition 3.1 There exists an extension operator A ε from the set of continuous piecewise affine functions on T ε ∩ ωε to the set of continuous piecewise affine functions onT ε ∩ ω such that ∇(A ε ψ) L 2 (ω) ≤ C ∇ψ L 2 (ω ε )where C is independent of ε, and if ψ is globally affine on the union of boundary triangles, then A ε ψ is the same affine function on ω \ ωε .

Figure 5 :

 5 Figure 5: Full triangle, empty triangle, unit cell.

Figure 6 :

 6 Figure 6: Recovering coercivity on empty triangles.

kY W 0

 0 (y, gS + Dθ (y)) dy = kY W 0 (y, -g + Dθ (y)) dy = kY W 0 (y, g -Dθ (y)) dy, and thus W kY (gS) = W kY (g)for all k, which implies W hom (gS) = W hom (g).

FirstFigure 8 .

 8 Figure 8. Planar Cauchy-Born configuration, 84 iterations of gradient.

Figure 9 .

 9 Figure 9. Nonplanar, non Cauchy-Born configuration, 467 iterations of gradient.

Figure 10 .

 10 Figure 10. Plot of W kY (F) against k.

Figure 11 .

 11 Figure 11. Plot of W 50Y (FQ(θ )) against θ .

Figure 11bis .

 11bis Figure 11bis. Closeup of the previous plot.

Figure 12 .

 12 Figure 12. Plots of W kY on the tangent to θ → FQ(θ ) at θ = 0.14. From top to bottom curves: k = 3, 5, 10, 20, 30.

We now are in a position to complete the proof of Theorem 4.5. Proposition 4.16 For all ψ ∈ H 1 ω 0 (ω; R 3 ), we have (Γlim ε→0

Proof. Let ε be a sequence that tends to 0. From any subsequence ε , we extract a subsequence ε such that I ε 0 is Γ-convergent. By Proposition 4.13, we have

By Lemma 4.15, we also have

for all η > 0. Hence Γ-lim I ε 0 (ψ) = I 0 (ψ) and we conclude by uniqueness of the Γ-limit.

The following is an easy consequence of Theorem 4.5.

Corollary 4.17 The function W hom is quasiconvex.

Properties of the homogenized energy

First of all, it is quite clear that since the original energy is frame-indifferent, i.e., invariant by the left action of SO(3), so is the homogenized energy W hom .

We next turn to material symmetry considerations. The material symmetry of W 0 is fairly simple, see also [START_REF] Raoult | Elastic lattices: equilibrium, invariant laws and homogenization[END_REF]. Let R ∈ SO(2) be the rotation of angle 2π 3 and S ∈ SO(2) the rotation of angle π.

Proposition 5.1 The material symmetry group of W 0 contains the circular group C 6 of all rotations of angle an integer multiple of π 3 .

Proof. We have Rs = ts and Rt = -s. Therefore, for all g ∈ L (R 2 ; R 3 ),

if we set τ = g(s) + τ. Taking infima with respect to τ on the left and τ on the right, we obtain W 0 (y, gR) = W 0 (y, g).

Similarly, we have


IMPROVING CHAOTIC OPTIMIZATION ALGORITHM USING A NEW GLOBAL LOCALLY AVERAGED STRATEGY

Tayeb Hamaizia, René Lozi * † Abstract. Recently chaotic optimization algorithms as an emergent method of global optimization have attracted much attention in engineering applications. Their good performances have been emphasized [START_REF] Coelho | Tuning of PID controller for an automatic regulator voltage system using chaotic optimization approach[END_REF][START_REF] Coelho | Reliability-redundancy optimization by means of a chaotic differential evoluytion approach[END_REF][START_REF] Davendra | Chaos driven evolutionary algorithms for the task of PID control[END_REF]. In the frame of evolutionary algorithms, the use of chaotic sequences instead of random ones has been introduced by Caponetto and al. [START_REF] Caponetto | Chaotic Sequences to Improve the Performance of Evolutionary Algorithms[END_REF]. Since their original work, the literature on chaotic optimisation is flourishing. They are used in the scope of tuning method for determining the parameters of PID control for an automatic regulator voltage, or in order to solve economic dispatch problems, or also for engineering design optimization and in many others physical, economical and biological problems. Different chaotic mapping have been considered, combined with several working strategies. The assessments of the algorithms have been done with respect to numerous objective functions in 1, 2 or 3dimension. In this paper we present an improvement of the COLM (Chaotic Optimization based on Lozi Map) presented in [START_REF] Coelho | Tuning of PID controller for an automatic regulator voltage system using chaotic optimization approach[END_REF], which is based on a new global locally averaged strategy. The simulation results are done with a 2-D objective function possessing hundreds of local minima, in order to test this new method vs. the previous one in very tough conditions. We emphasize an improvement of the optimisation.

Improved COLM Method

Chaos theory (the term chaos was coined par Li and Yorke [START_REF] Li | Period three implies chaos[END_REF]) is recognized as very useful in many engineering applications. An essential feature of chaotic systems is sensitive dependence on initial condition, (i.e. small changes in the parameters or the starting values for the data lead to drastically different future behaviours). Details about analysis of chaotic behavior can be found in [START_REF] Li | Period three implies chaos[END_REF][START_REF] Lozi | Un attracteur étrange? du type attracteur de Hénon[END_REF][START_REF] Strogatz | Nonlinear dynamics and chaos[END_REF][START_REF] Parker | Practical numerical algorithms for chaotic system[END_REF][START_REF] Alligood | Chaos: an introduction to dynamical systems[END_REF]. The application of chaotic sequences can be an interesting alternative to provide the search diversity in an optimization procedure. Due to the non-repetition of chaos, it can carry out overall searches at higher speeds than stochastic ergodic searches that depend on probabilities. A novel chaotic approach is proposed in [START_REF] Coelho | Tuning of PID controller for an automatic regulator voltage system using chaotic optimization approach[END_REF] based on Lozi map [START_REF] Lozi | Un attracteur étrange? du type attracteur de Hénon[END_REF] which is piecewise linear simplification of the Hénon map [START_REF] Hénon | A two dimensional mapping with a strange attractor[END_REF] and it admits strange attractors. It is given by

y 1 (k) = 1 -a|y 1 (k -1)| + by(k -1) y(k) = y 1 (k -1) ( 1
)
where k is the iteration number. In this work, the values of y are normalized in the range [0,1] to each decision variable in 2dimensional space of optimization problem. This transformation is given by

z(k) = (y(k) -α) β -α . ( 2
(k), x 2 (k), x 3 (k) with f (x 2 (k)) < f (x 3 (k)) < f (x 1 (k)). (3)
The local search following global one starts from the best global result x 2 (k) (from (3)) and gives x 2 (k + 1). Instead the local-global search around x 1 (k), x 2 (k) and x 3 (k), leads to

x 1 (k + 1), x 2 (k + 1), x 3 (k + 1) which verify f (x 1 (k + 1)) < f (x 3 (k + 1)) < f (x 2 (k + 1)). ( 4
)
The local search following the local-global one starts now from the best globally averaged result x 1 (k + 1) (from(4)) and leads to x

f (x) < f (x 1 (k + 1)). ( 5
)
During the chaotic local search, the step λ is an important parameter in convergence behavior of optimization. Hence, two different values of λ are successively employed during the local search. We call this method ICOLM (improved COLM). Many unconstrained optimization problems with continuous variables can be formulated as the following functional optimization problem. Find X to minimize f (X),

X = [x 1 , x 2 , , xn] Subject to x i ∈ [L i , U i ].
Where f is the objective function, and X is the decision solution vector consisting of n variables x i ∈ R n bounded by lower (L i ) and upper limits (U i ).The ICOLM can be illustrated as follows:

Inputs: Mg: max number of iterations of chaotic Global search.

x i (k) = L i + z i (k).(U i -L i ) if f (X(k)) < f then X = X(k); f = f (x(k)) end if -Step 2-1: sub algorithm of chaotic local search: while j ≤ M gl1 do for i = 0 to n do if r ≤ 0.5 then x i (j) = xi + λ gl1 z i (j).|(U i -L i )| else x i (j) = xi -λ gl1 z i (j).|(U i -L i )| end if end for if f (X(j)) < f then X = X(j); f = f (x(j)) end if j = j + 1 end while -Step 2-2: sub algorithm of chaotic local search: while s ≤ M gl2 do for i = 0 to n do if r ≤ 0.5 then x i (s) = xi + λ gl2 z i (s).|(U i -L i )| else x i (s) = xi -λ gl2 z i (s).|(U i -L i )| end if end for if f (X(s)) < f then X = X(s); f = f (x(s)) end if s = s + 1 end while k = k + 1 end while -Step 3: algorithm of chaotic local search: while k ≤ Mg × (M gl1 + M gl2 ) + M l do for i = 0 to n do if r ≤ 0.5 then x i (k) = xi + λz i (k).|(U i -L i )| else x i (k) = xi -λz i (k).|(U i -L i )| end if end for if f (X(k)) < f then X = X(k); f = f (x(k)) end if k = k + 1 end while

A tough objective function

In order to test this new method vs. the previous one in very tough conditions the simulation results are done with the following 2-D objective function possessing hundreds of local minima: The function f which is very complex, has serval local maxima.

f = x 4 1 -7x 2 1 -3x 1 + x 4 2 -9x 2 2 -5x 2 + 11x 2 1 x 2 2 + 99sin(71x 1 ) + 137sin(97x 1 x 2 ) + 131sin(51x 2 ). (6)
We test ICOLM on the search domain:

-10 ≤ x i ≤ 10, i = 1, 2.
The essential feature of this benchmark function is that location of minima is not symmetric. In a forthcoming paper we will extend our numerical analysis in higher dimension with an extended benchmark suite [START_REF] Cong | An Improved Algorithm of Chaos Optimization[END_REF]. We display few of the results we have obtained showing the better optimization results obtained by this new methods.In each case study, 48 independent runs were made for each of both the COLM and ICOLM methods involving 48 different initial trial conditions y 1 (0), y(0) (parameters of Lozi map). not far from the value of f on the global minimum.

For all studied cases, the four configurations, numbered from IC1 to IC4 and C1 ti C4, that are used are presented in tab. 1. The locally averaged strategy of ICOLM is illustrated on Fig. 3

Conclusion

In every test, with the same computational cost, ICOLM gives better than COLM best values and Mean Best values but in one case. The presented study allows us to conclude that the proposed method is fast and converges to a good optimum. because we used a sampling mechanism to coordinate the research methods based on chaos theory, and we refined the final solution using a second method of local search. Further research is needed to gain more confidence and better understanding of the proposed methodology.

The proposed algorithm has to be evaluated for a large number of test functions in higher dimension.

1 0Figure 1 :

 11 Figure 1: density of iterated values of y(k) of equation (1) over the interval [0, 1] splitted in 100 boxes for 10,000,000,000 iterated values.

Figure 2 :

 2 Figure 2: Heuristics of the global locally-averaged strategy.

Figure 3 :

 3 Figure 3: Locally-averaged strategy of chaotic search. Results of every Step 2-2 for f

Figure 6 :

 6 Figure 6: magnification of Fig. 4.

Figure 7 :

 7 Figure 7: magnification of Fig. 6.

Table 1 :

 1 on which the result of every step 2-2 is plotted. The set of parameters values for every run on the benchmark suite defined in Sec. 2. 2.

		λ	λ M gl1 λ M gl2	M g	M l M gl1 M gl2
	IC1 0.001	0.04	0.01	6	50	2	2
	IC2 0.01	0.04	0.01	10	50	2	2
	IC3	0.1	0.04	0.01	10	50	2	2
	IC4	0.1	0.04	0.01	100	50	5	5
	C1 0.001			24	50
	C2	0.01			40	50
	C3	0.1			40	50
	C4	0.1			1000 50

Figure 4: plot of test function used in this study

Figure 5: Position of the minima in the search domain Best value Mean value Std.Dev

Table 2 :

 2 ICOLM

Table 3 :

 3 6. COLM

	Best value Mean value Std. Dev	x ȳ
	C1 -371.0150	-368.5212	11.1135	0.3105 0.2442
	C2 -358.4331	-352.9766	2.4424	3.7283 6.6158
	C3 -377.9280	-368.7777	8.8169	3.2738 6.1685
	C4 -382.7108	-379.7557	1.7817	-5.8930 2.9309