Nicolas Zilber

Philippe Ciuciu

Patrice Abry
email: 2patrice.abry@ens-lyon.fr

Virginie Van Wassenhove

MODULATION OF SCALE-FREE PROPERTIES OF BRAIN ACTIVITY IN MEG

Keywords: MEG, scale-free activity, scaling, multifractal analysis

The analysis of scale-free (i.e., 1/f power spectrum) brain activity has emerged in the last decade since it has been shown that low frequency fluctuations interact with oscillatory activity in electrophysiology, noticeably when exogenous factors (stimuli, task) are delivered to the human brain. However, there are some major difficulties in measuring scale-free activity in neuroimaging data: they are noisy, possibly nonstationary ... Here, we make use of multifractal analysis to better understand the biological meaning of scale-free activity recorded with Magnetoencephalography (MEG) data. On a cohort of 20 subjects, we demonstrate the presence of self-similarity on all sensors during rest and visually evoked activity. Also, we report significant multifractality on the norm of gradiometers. Finally, on the latter signals we show how self-similarity and multifractality are modulated between ongoing and evoked activity.

INTRODUCTION

Starting with the discovery of brain oscillations in electroencephalography (EEG) and later in magnetoencephalography (MEG),many studies have attempted to explain their mechanism and their possible functional role in behavior and cognition. Some frequency bands were particularly conspicuous and were given names such as δ (1-3Hz), θ (4-8Hz), α (8-12Hz), β (12-30Hz) and γ (>40Hz). On the other hand, very low frequency activity (< 1Hz) (also named "infraslow") did not attract as much attention due to the power-law distribution of the spectral density in the range of low frequencies. This property is very ubiquitous in dynamic systems (fully developed turbulence, internet traffic, earthquakes, stock market exchange, ...) [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF]. In the context of brain imaging, such activity was initially attributed to the intrinsic 1/f electronic device noise and was systematically removed by high-pass filtering or normalization of M/EEG data.

However, several authors have shown that ongoing activity, the major part of brain activity [START_REF] Bullock | Are the electroencephalograms mainly rhythmic? Assessment of periodicity in wide-band time series[END_REF][START_REF] Wj Freeman | Simulated power spectral density (PSD) of background electrocorticogram (ECoG)[END_REF], has scale-free dynamics -i.e. a 1/f power spectrum [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF][START_REF] Wj Freeman | Simulated power spectral density (PSD) of background electrocorticogram (ECoG)[END_REF][START_REF] Bj He | Scale-Free Properties of the Functional Magnetic Resonance Imaging Signal during Rest and Task[END_REF][START_REF] Ciuciu | Multifractal analysis of resting state networks in functional MRI[END_REF]. Since then, the study of scale-free properties has emerged as a new research topic in neuroscience and brain neuroimaging. A major finding is that scale-free activity recorded with EEG interacts with ongoing oscillatory activity [START_REF] Sadaghiani | The relation of ongoing brain activity, evoked neural responses, and cognition[END_REF]. Moreover, the surprising correlation between EEG microstates and fMRI resting-state networks (RSN) could be explained by the scale-free properties of these microstates [START_REF] Van De Ville | EEG microstate sequences in healthy humans at rest reveal scale-free dynamics[END_REF]. Other experiments have shed light on the modulation of scale-free activity by exogenous (drugs,task) and endogenous (brain state, pathology) factors [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF][START_REF] Buiatti | Feedback modulates the temporal scale-free dynamics of brain electrical activity in a hypothesis testing task[END_REF][START_REF] Suckling | Endogenous multifractal brain dynamics are modulated by age, cholinergic blockade and cognitive performance[END_REF].

Nonetheless, some major difficulties remain when assessing scalefree properties in neuroimaging data: they are noisy, multivariate, nonstationary, and possibly corrupted by physiological artifacts (eye movements, heartbeats...). In particular, the nonstationarity and the long range dependence yield a lack of statistical consistency when estimating scaling properties in the time or frequency domain using methods like Detrended Fluctuation Analysis (DFA) [START_REF] Linkenkaer-Hansen | Long-range temporal correlations and scaling behavior in human brain oscillations[END_REF]. To overcome this difficulty, one solution is to use wavelet-based analysis of scaling properties.

Wavelet analysis implies a major change of paradigm: frequency (and hence the frequency shift operator) is no longer considered as the relevant descriptor and is replaced instead by scale (and hence by a dilation, or scale change operator). This enables to correctly estimate scale invariance parameters even when data present various forms of non-stationarity. Moreover, some empirical results [START_REF] Ciuciu | Multifractal analysis of resting state networks in functional MRI[END_REF] suggest that ongoing brain activity may be better modeled by the class of multifractal (MF) processes. This class contains and extends that of self similar processes (e.g. fractional Brownian motion, fBm) thereby offering a more versatile and richer description of the data in terms of scaling properties and departures from Gaussian property.

In Section 2, we summarize spectral and multifractal analyses. In Section 3, we show the presence of scaling in MEG data at the sensor level and compare standard spectral estimators with the waveletleader based multifractal (WLMF) analysis [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. Additionally, we study the modulation of multifractal attributes between two cognitive states: i.) Rest (R), and ii.) visualization of Random Dots Kinematograms (RDK). Although no significant difference could be found when analyzing direction-specific gradio-or magneto-meters (see Section4.1), statistical differences are exhibited when analyzing the norm of gradiometers (see Section 4.2). A criterion is derived in Section 4.3 to explain this surprising difference and an interpretation is suggested. Conclusions are drawn in Section 5.

SCALING ANALYSIS

Spectral analysis

Scaling is classically analysed using standard spectrum estimation such as windowed-averaged periodograms (e.g. Welch's estimator) or detrended fluctuation analysis (DFA) [START_REF] Linkenkaer-Hansen | Long-range temporal correlations and scaling behavior in human brain oscillations[END_REF]. The information is supposed to be mainly contained in the power spectrum ΓX (f ) of the signal X for f → 0, which behaves as a power law:

Γ Welch (f ) = 1 n n k=1 ΓX (f, t k ) ≃ C|f | -γ , γ > 0. (1)
This implicitly requires the data to be stationary. Accordingly, fractional Gaussian noise (fGn) is used as the paradigm model for such analysis if γ ∈ (-1, 1). However, on real neuroimaging data, the measured exponent γ usually satisfies γ ∈ (1, 3) [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF], the data are thus better modelled by the non-stationary fBm process, whose in-crements follow a fGn. For both fBm and fGn, H takes values in the range of (0, 1), hence corresponding to γ = 2H + 1, for the former and to γ = 2H -1 for the latter. However, if analysed data follows a fGn with H > 0.5, the underlying process is said long range dependent, meaning that its autocorrelation function decreases very slowly, while if the analysed time series follows a fBm and H > 0.5, the corresponding process is self-similar, hence nonstationary. The fBm/fGn pair is intrinsically coupled to Gaussian processes. However, empirical results on brain imaging data [START_REF] Ciuciu | Log Wavelet Leaders Cumulant Based Multifractal Analysis of EVI fMRI Time Series: Evidence of Scaling in Ongoing and Evoked Brain Activity[END_REF] have shown that data often departs from Gaussian distribution, hence ruling out these models. Also, empirical evidence on brain electrical field [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF] indicates that the estimated γ takes values oscillating around the fBm/fGn γ = 1 limit, hence suggesting a constant switch between nonstationary selfsimilarity and stationary long memory, a very unsatisfactory situation for data modeling and interpretation. By using wavelet transforms, the nonstationarity issue can be overcome. Furthermore, the wavelet-based scaling properties (eg, Hurst parameter) are known to be more accurately estimated. Last but not least, departure from Gaussianity can be captured by analyzing statistical properties of the wavelet coefficients at different statistical orders ( = 2) and interpreting such results in the context of multifractal processes.

Multifractal Analysis

Let dX (j, k) denote the discrete wavelet transform coefficients of a signal X, where j denotes the analysis scale (a = 2 j ) and k the time position (t = 2 j k), computed from a Daubechies mother-wavelet ψ0(t). Scale invariance is now commonly and operationally defined as the power law behaviors of the (time average of the q-th power of the) dX (j, k), with respect to the analysis scale j, for a given (large) range of scales j ∈ [jm, jM ], jM > jm:

S d (j, q) ∆ = 1 nj n j k=1 |dX (j, k)| q ≃ cq2 jζ d (q) . (2)
S d (j, q) is the structure function at the scale j and order q. ζ d (q) refers to the scaling exponents. If X is monofractal (or self-similar), ζ d (q) is linear and can be simply characterized by one parameter (namely, the Hurst exponent H). However, if X is multifractal, ζ d (q) is a concave function and its description is more complex. It is worth noting that the only use of the power spectrum (i.e. S d (j, q) for q = 2) does not give any information on the multifractal properties. Ideally, we should assess ζ d (q) for all values of q (positive and negative). It has been shown recently that replacing wavelet coefficients dX (j, k) by the so-called wavelet Leaders LX (j, k) would improve the estimation of ζ d (q), especially for negative values of q.

Wavelet Leaders LX (j, k) are defined as the local suprema of wavelet coefficients within a local neighborhood, and over all finer scales: [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]: 2) still holds by replacing |dX (j, k)| q with LX (j, k) q and ζ d (q) with ζ L (q). In practice, estimating of ζ L (q) for all values of q turns out very difficult. We make use of a polynomial expansion of ζ L (q) instead:

LX (j, k) = sup λ j ′ ,k ′ ⊂3λ j,k |dX (j ′ , k ′ )|, where λ j,k = [k2 j , (k + 1)2 j ) and 3λ j,k = ∪ m{-1,0,1} λ j,k+m . Eq. (
ζ L (q) = p≥1 c L p q p /p ! = c L 1 q + c L 2 q 2 2 + c L 3 q 3 6 + • • • .
As detailed in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF], the c L p coefficients are estimated by linear regressions, where c L 1 reflects the self-similarity (or long memory) of the process, especially for monofractal processes, where c L 1 = H and ∀p ≥ 2 c L p = 0. For multifractal processes, we specifically have c L 2 < 0 to quantify the amount of singularities and c L 3 ∈ | for describing the assymetric nature of the singularity distribution. In this study, we have focused on c L 1 and c L 2 .

The WLBMF toolbox has been designed to robustly estimate these c L p parameters [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. it has been already used successfully to analyse multifractal properties of fMRI signals [START_REF] Ciuciu | Multifractal analysis of resting state networks in functional MRI[END_REF] and scale-free dynamics in EEG microstates [START_REF] Van De Ville | EEG microstate sequences in healthy humans at rest reveal scale-free dynamics[END_REF]. In the following, it allows us to accurately characterize scale-free properties of MEG time series and their modulation by external stimulus.

POWER LAW SPECTRUM ON MEG RECORDINGS

Data acquisition

Brain magnetic fields were recorded in a magnetically shielded room using the whole-head Elekta Neuromag MEG system (Elekta LTD, Helsinki) with 102 triple-sensor elements (two orthogonal planar gradiometers and one magnetometer) in upright position. Data were sampled at 2000Hz and pass-band filtered between 0.03Hz and 600Hz. Signal Space Separation (SSS) was carried out with Max-Filter (Elekta LTD) to remove external interferences. Occular and cardiac artifacts were removed using principal component analysis based on the recorded electrooculograms (EOG) and electrocardiograms (ECG).

Twenty healthy right-handed subjects (aged 19-27 years) with normal or corrected-to-normal vision participated in this experiment. Each participant provided an informed consent in accordance with the Declaration of Helsinki ( 2008) and the Ethics Committee on Human Research at NeuroSpin (Gif-sur-Yvette, France). Before each experiment, empty room MEG recordings were acquired for 5 minutes. The experiment consisted of two sessions: 1) 'Rest' (R) session: the subject was asked to keep eyes open and to stare at the black screen in front of him during 5 minutes. 2) 'RDK' session: The subject fixated a cross on the center of the screen, while a random dot kinematogram of 2.5 s was presented continuously 120 times (radius of the annulus: 4-15 deg, dot radius: 0.2 deg , dot flow: 16.7 dots per deg 2 .s, dot speed: 10 deg/s). Initially, all dots were moving incoherently (i.e., each dot moves in a different direction). In half of the trials, 95% of dots took a coherent direction (ie move in the same direction) after the stimulus onset. In the other half, dots remain incoherent. Moreover, at each frame, 5% of all dots were reassigned to new positions with new motion directions for incoherent dots.

Empirical evidence of long memory

For each subject, each session and each sensor, we computed Welch's power spectrum estimate derived in Eq. [START_REF] Bj He | The Temporal Structures and Functional Significance of Scale-free Brain Activity[END_REF]. A representative spectrum computed by averaging all radial gradiometers is shown in Fig. 1 in black and dark blue lines for the R and VM sessions, respectively. The presence of scaling or 1/f behaviour clearly appears as a linear slope in this log-log plot over the (0.1-3Hz) frequency range. In addition, this 1/f power spectral density is clearly different from that measured during the empty recording (light blue line in Fig. 1), thus ruling out the hypothesis that the 1/f is simply driven by the electronic device noise. Interestingly, the power of acquisition noise is upper bounded by that measured during brain activity recordings. Finally, αand β-band oscillations emerge during the presentation of visual motion stimuli and even more during the rest, while they do not appear in the empty recording.

As a comparison with the previous method, we also estimated the structure fonctions S d (j, q) for q = 2 only. Structure fonctions at this 2nd statistical order are indeed equivalent to the power spectrum while estimated using discrete wavelets. f and j are related one another by f = 3 respectively the R, VM and empty recordings. In contrast to Welch's periodogram Γ Welch (f ), the S d (j, 2) estimates are less sensitive to oscillation peaks. Discrete wavelets are thus more appropriate to analyse the 1/f behaviour in low frequencies.

This preliminary result confirmed the presence of scaling in the data, which cannot be attributed to the sensor noise. It also shows the advantage of using discrete wavelets for analyzing the 1/f spectrum.

Although spectral analysis provides restricted information on the scale invariance properties, it allows us to determine the scale range on which the WLBMF analysis is applied. In the following, we analyse the multifractal properties in the scale range [jm, jM ] = [START_REF] Suckling | Endogenous multifractal brain dynamics are modulated by age, cholinergic blockade and cognitive performance[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of within-subject fMRI time series[END_REF] (i.e. in the frequency range 0.1-3Hz). The log-scale diagram is actually linear in this part. In the same manner, Fig. 3 shows the topographies of the mean µ r 2,k,c -values. Multifractality (µ r 2,k,c < 0) is observed on the gradiometers located in the occipito-parietal regions.

log 2 Γ(f )[log 2 S L (j, 2)] in dB log 2 f [j]
µ VM 2,k,• µ R 2,k,• c = grad R c = grad T c = mag
To this end, we computed the t-statistics T r i,k,c = (µ r i,k,cδi)/σ r i,k,c , where the latter term denotes the group-level standard deviation. Rejecting H0 r (1,k,c) clearly amounts to localizing regions (k) and sensor type (c) eliciting self-similarity in session r while rejecting H0 r (2,k,c) enables to discriminate multifractality from selfsimilarity. Statistical results (not shown) demonstrate that H0 r (1,k,c) was rejected everywhere even after Bonferroni correction for multiple comparisons (p-values p corr < 10 -6 ). However, due to a large between-subject variability, H0 r (2,k,c) was not rejected after Bonferroni correction ( p uncorr ≈ 10 -2 ) for a False Positive Rate (FPR) of 5%. Paired t-tests were also computed to compare the R and VM sessions, where the null assumption consists in assuming the same mean values:

µ R i,k,c = µ VM i,k,c
The null hypotheses H0 R,VM (1,k,c) and H0 R,VM (2,k,c) were never rejected indicating that no significant difference can be exhibited between ongoing and evoked activity. However, in several MEG studies, the signal of interest is usually considered as a non-linear combination of both types of gradiometers, namely the ℓ2-norm of gradiometers. Therefore, we decided to analyse the norm of gradiometers too.

WLBMF analysis on the norm of gradiometers

As each pair of gradiometers is orthogonal, their norm is simply defined by grad = grad 2 R + grad 2 T . As already done, we estimated c L,r 1,k, . (s) and c L,r 2,k, . (s) for grad and then computed mean values µ r i,k, . (i = 1, 2). A very noticeable result (Fig. 4[left]) is the global reduction of self-similarity shown by the decrease of µ r 1,k, . in both sessions. Nonetheless, the self-similarity remains large enough to be statistically significant everywhere. We then compared the two sessions by computing ∆µ 1,k, . = µ R 1,k, . -µ VM 1,k, . (see Fig. 4[Mid.]) and testing the statistical significance of this difference using a paired t-test where the null hypothesis was H0 R,VM : µ R 1,k, .

µ VM 1,k, . . As shown in Fig. 4[right], significant differences (pcorr < 10 -2 ) emerges in the occipital area where ∆µ 1,k, . > 0. This finding is consistent with the literature dealing with the taskinduced modulation of scale-free properties [START_REF] Bj He | Scale-Free Properties of the Functional Magnetic Resonance Imaging Signal during Rest and Task[END_REF] observed in fMRI.

The multifractal behaviour of grad measured through negative values of µ r 2,k, . is emphasized in Fig. 5. As shown in Fig. 3[Mid.], ∆µ 2,k, . < 0 in the occipito-parietal (MT) area. This means that multifractality decreases in the target region (MT) of the delivered stimuli during evoked activity (cold spots in Fig. 5[Left]) in comparison with resting-state. The statistical paired t-test performed with H0 R,Vm : µ R 2,k, .

µ VM 1,k, . is significant in the MT region but not enough to survive to Bonferroni correction for FPR=5%.

µ VM 1,k, . µ R 1,k, . ∆µ 1,k, .
puncorr (in -log10)

Simulation with fGn

To understand the impact of the ℓ2 norm on the multifractal properties, we synthesized 10000 fGn processes X (H = .75) with the circulant embedded method [START_REF] Dietrich | Fast and exact simulation of stationary Gaussian processes through circulant embedding of covariance matrix[END_REF]. We considered a simplified problem where grad R = grad T , meaning that grad ∝ grad 2 R = |grad R |. Hence, our simulation amounts to estimating the multifractal properties of X and |X|. We systematically observed a decrease of long memory: ∆ c1 c1 |X| -c1 X < 0, which is consistent with our experimental results. Likewise, multifractality appeared in |X| whereas it is absent in X (fGn process). Given that the number of sign changes nsign in X is the only parameter that induces regularity changes in |X| (i.e. no addition of irregularity if X is always positive or negative), nsign was expected to be the main cause of these observations. This hypothesis is validated in Fig. 6 when plotting ∆ c1 with respect to the proportion of sign changes psign in X. ∆ c1 decreases almost systematically (i.e. with a low variance) as a function of psign. It goes down dramatically when psign > 0.3. Regarding the multifractality (i.e. ∆ c2), the observation is similar but less pronounced owing to the higher estimator variance (figure not shown).

∆ c1 psign Fig. 6. c1 variation when analyzing |X| as a function of psign. 10000 fGn processes synthesized (in blue). Fitted curve in red.

Therefore, the apparent contrast in the occipital area when analyzing the norm of gradiometers can be explained by the number of sign changes in the gradiometers or by their phase. Because gradiometers measure the spatial derivative of the magnetic fields in two orthogonal directions, a change of sign suggests a change of source orientation or even perhaps of the source itself.

DISCUSSION AND CONCLUSION

We have demonstrated the presence of long memory in MEG data over all sensors. A small amount of multifractality was observed on the gradiometers in the occipital and parietal scalp regions. Interestingly, only the norm of gradiometers exhibits a modulation of the multifractal properties between ongoing and evoked activity. Additionally, this modulation was very localized to an area probably sensitive to the nature of the stimuli (as assessed by additional experimental work in progress). It would be very interesting to change the sensory context (e.g. auditory stimuli) in order to check whether the modulation moves to the expected auditory regions in order to further show the specificity of our results.

The norm of gradiometers seems to capture more information. This is largely due to the nonlinear nature of the norm and the number of sign changes in the gradiometer signals. The contrast between R and RDK in the occipital area can be interpreted as a higher rate of source orientation changes in hMT+/V5 during the presentation of RDK. Importantly, we propose that this result extends to a more general concept: in any imaging modality (e.g. fMRI), the knowledge of any nonlinear transformation in the generative model [START_REF] Vincent | Spatially adaptive mixture modeling for analysis of within-subject fMRI time series[END_REF] of the data can be crucial to correctly interpret its multifractal properties.

Fig. 1 .

 1 Fig. 1. Log-log plot of spectrum estimates averaged across all radial gradiometers. Welch's periodograms are plotted in black (R), dark blue (VM) and light blue (Empty). Wavelet-based estimates are plotted in green (R), red (VM) and yellow (Empty).

4 .Fig. 2 .

 42 Fig. 2. Group-level µ r 1,k,c -topographies computed for all sensors k, all channel types c and the two sessions r = R,VM.

Fig. 3 .

 3 Fig. 3. Group-level µ r 2,k,c -topographies computed for all sensors k, all channel types c and the two sessions r = R,VM.

Fig. 4 .

 4 Fig. 4. Left: Group-level µ r 1,k,c -topographies computed for the norm of the gradiometers and the two sessions r = R,VM. Middle: mean difference ∆µ 1,k, . between the R and VM sessions. Right: Uncor. p-values resulting from the between-session paired t-test. µ VM 2,k, .

∆µ 2 ,Fig. 5 .

 25 Fig. 5. Left: Group-level µ r 2,k,c -topographies computed for the norm of the gradiometers and the two sessions r = R,VM. Middle: mean difference ∆µ 2,k, . between the R and VM sessions. Right: Uncor. p-values resulting from the between-session paired t-test.

fs 2 j , where fs is the sampling frequency. The resulting wavelet-based spectrum also named log-scale diagram (log 2 S L (j, 2) vs. log 2 2 j = j) averaged over all radial gradiometers is shown in Fig.1where green, red and yellow lines represent


J L Marcelin

Optimization of the boundary conditions by genetic algorithms

Keywords: optimization, boundary conditions, genetic algorithms

This work examines the possibility of using a stochastic method, called the genetic algorithm for the optimization of boundary conditions in finite elements calculations. The examples show that using genetic algorithms in order to optimize boundary conditions is an efficient way.

I. Introduction

In numerous application cases (e.g. taking workpieces in machining, vibrations of a mechanical structure), the optimization of boundary conditions (location and nature of the boundary conditions) can bring an interesting improvement of the studied structure's mechanical behavior. These structures are most often calculated with the finite element method, and within this framework, the calculation of the sensitivities with regard to the boundary conditions remains enough complicated due to the discrete nature of the problem, and for example, in shape optimization, the boundary conditions problem is even most often eluded. Besides, the deterministic methods of optimization, called gradient methods, need a reliable calculation of these sensitivities. Some other stochastic or probabilistic methods of optimization are currently in vogue, like the simulated annealing method or that of genetic algorithms, which main benefits are they assure convergence without the use of derivatives, and can be used with possibly discrete variable and non-derivable functions. The detractors of this method point up without reason the high number of calculations, especially in the case of an analysis method of finite element type.

The author has a great experience in genetic algorithms ( [START_REF] Marcelin | Genetic optimization of stiffened plates and shells[END_REF] to [6]). Few researches have been made on the boundary conditions optimization, and in the structural mechanics field, there are almost none. The [START_REF] Su | Determining the equivalent impedance boundary condition for corrugated coatings based on the genetic algorithm[END_REF] work, described below, deals with the optimization of the boundary conditions in electromagnetism. In [START_REF] Su | Determining the equivalent impedance boundary condition for corrugated coatings based on the genetic algorithm[END_REF], a methodology based on the genetic algorithm is proposed to determine the equivalent impedance boundary condition for corrugated material coating structures. We have find only two papers at our knowledge in the structural mechanics field [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF] and [START_REF] Wang | Application of genetic algorithm for the support location optimization of beams[END_REF], and only one [START_REF] Wang | Application of genetic algorithm for the support location optimization of beams[END_REF] with the use of genetic algorithms. In [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF], optimization of boundary conditions for maximum fundamental frequency of vibrating structures is done. In [START_REF] Wang | Application of genetic algorithm for the support location optimization of beams[END_REF] the use of genetic algorithms for the selection of optimal support locations of beams is presented. Both elastic and rigid supports are considered. The approach of adapting the genetic algorithms into the optimal design process is described. This approach is used to optimize locations of three supports for beams with three types of boundary conditions. This article relates to fixed geometry and is to show that the optimization of boundary conditions is feasible by combining a genetic algorithm and the finite element method (the optimization of boundary conditions in shape optimization will be dealt with in future research). Indeed, in the case of a fixed shape and because of the characteristics of the finite element method, the calculation volume can be considerably reduced. The main reasons that are to be explicated in this work are the following: the stiffness matrix is calculated and assembled once and for all; in the case of a structure for which some boundary conditions can be fixed and other can be variables (i.e. entering in the optimization framework), it would be possible to triangularize the stiffness matrix once and for all, and to take into account the variable boundary conditions thanks to a penalization process of the energetic functional, to be minimized by the boundary conditions. In such conditions, even is the number of analyses is still important, the calculation time will remain reasonable because the analyses won't be systematically complete. Various examples will aim at showing that the implemented process helps in optimizing the boundary conditions and is fairly efficient.

II. The methods used

II.1. Genetic algorithms

The genetic algorithm method has been used several times within the various problems of mechanics. These algorithms were found to be very efficient, as in the case of the damping maximization of composite beams or plates or as somewhat diverse issues. The interest of these algorithms has also been showed in the difficult case of the optimization of gears. The genetic algorithms are now well known and this article is not to introduce them in details nor generally.

Although it may seem so, the genetic algorithm method is not magic at all. It is part of the methods called "stochastic". The most famous of this kind of method is the already old simulated annealing. The main benefit of these methods is that they operate simultaneously on a sample of the solution space. The genetic method differs from simulated annealing due to the operators used to make this population sample evolve. The convergence is always ensured toward an extremum which is not necessarily the absolute extremum, but which is more likely to be absolute than if the conventional gradient method is used. Actually, a stochastic method explores more largely the solution space.

II.2. Optimization of the boundary conditions

This kind of optimization consists in combining a standard calculation program by finite elements (FE) (called thereafter analysis program) and the genetic algorithm.

The analyze program is a standard FE code. This code is simply to be called each time the genetic algorithm must estimate the cost function for a given chromosome. This is done for all the individuals of the population; consequently, for example, for 20 individuals and 30 generations, there will be 600 finite element half-analyses (the total stiffness being calculated once and for all), which is relatively low compared to the 2 20 possible solutions. On the opposite, for the various tests that are done, especially those introduced after, there was not necessary to implement a penalization strategy of the "total potential energy" functional by the imposed boundary conditions, because convergence was fast enough. The programmer work simply consists in drafting a "pre-analysis" program that can decode the chromosome in question and that can automatically modify the finite element data file accordingly, and then in creating a "post-analysis" program that can extract the cost function from the finite element result file. Both these programs, as the calling of the finite element code, are built in the genetic program that drives the process.

Choosing the coding and the objective function

The problem contains two difficulties: First, the implementation of a solution code in the form of a simple and efficient chromosome and then the development of an objective function. The most generally used code is simple and natural (it has variants that are to be set forth in the examples): It can use the often used code for the boundary conditions in finite element programs, 0 being a free freedom degree and 1 a fixed freedom degree. The various codes of the concerned nodes are arranged end to end in a chromosome that is made of n binary digits that correspond to the n degrees of freedom that can be fixed. When it comes to the objective function, it depends on the posed problem. The first two examples are static cases, where the aim is to minimize the maximum displacement, or to minimize a deformation or a stress; the third example is a dynamic case, and the objective is to maximize the first natural frequency; it is also possible to try to remove two resonance frequencies. A lot of other choices are possible, such as multi-objective functions or penalizing the objective function by limitations.

Obtained results

Before each use, the genetic algorithm asks the user to specify the values of the following parameters: the number of individuals contained in a population, the maximum number of generations, the chromosomes length, the crossover probability, the mutation probability.

It is clear that the algorithm gives best results when the chosen values for the first two parameters are high (within the limits of capacity of the used hardware). Practically speaking, the number of individuals contained in a population will be around 1 to 5 times the number of digits contained in a chromosome. However, the crossover and mutation probabilities are more difficult to choose. It has already been said that mutation is a far less frequent phenomenon than crossover; in [START_REF] Marcelin | Genetic optimization of gears[END_REF], it is recommended the following values: Pcrossover = 0,60 ; Pmutation =0,03

These recommenced values come from a numerical experimentation on numerous examples. In any event, the crossover probability must be clearly superior to the mutation probability because mutation is less frequent. For example, if any mutation is removed, the algorithm yet converges toward an extremum but it is unlikely to be the absolute extremum. Theoretically speaking, convergence is obtained when all the cost values of a population stabilize around a maximum value. Practically speaking, convergence is rather slow, with ebb and flow, due to the very nature of the algorithm. The user only has to stop the process when the maximum value of a population cost does not evolve anymore; he then manually selects the most interesting individual(s) of the final population to compare their

III. Examples

Test 1

The first very easy, static test is made with the axisymmetric workpiece (of CL axis and z symmetry) illustrated on figure 1 and aims at verifying and making the implementation of the used techniques reliable. The stiffness has been calculated once and for all but no penalization has been applied to the boundary conditions; since the calculations are fast enough for the tests, this procedure has never been implemented.

For this test, the chromosome is a 10-binary-digit string, the first 5 digits are the codes of the boundary conditions of the 5 nodes that can be locked following z, and the following 5 are the codes of the boundary conditions of the 5 nodes that can be locked following y; therefore, the chromosome 1011001000 corresponds to the boundary conditions applied to nodes 1, 3, 4, and 7.

There are 2 10 possibilities. The objective is to minimize the d displacement of the node to which forces are applied. Since the genetic algorithm actually seeks the maximum of an objective function, the chosen objective is to maximize the 1-d function. The interest of this test is that the optimal solution is known: It is of course the 1111111111 chromosome, but the test helps in validating the process and in estimating how many steps are necessary for the genetic algorithm calculation to get this solution. We take here 40 individuals per population. The number of individuals in a population is usually around 1 to 5 times the chromosome's size (here the number of digits). The maximum is reached in only 5 generations (for the crossover and mutation probabilities provided in the last part), which corresponds to 200 half-analyses or a bit less (because a solution that appears several times during the process is calculated once and for all) and which is low compared to the possible 2 10 combinations. Non-consistent convergence is characteristic of genetic algorithms because the best individual of each population may very unlikely be eliminated; besides, if the algorithm is forced to keep only the bests, the method is not probabilistic anymore and the algorithm may be more efficient or diverge in some cases. Besides, if optimization is launched again with the same parameters, the obtained convergence is not at all the same, because the process is totally random. This test illustrates a first industrial application in the taking of workpieces in machining, always with simple data that help in validating and checking the implemented strategy.

The quality of the workpiece depends on the deformations caused by the machining, because of the machining process itself or because of the holding process of the workpiece on its support. The chosen workpiece is of z axis revolution and is illustrated on figure 2. There are three surfaces to be machined, S1, S2, and S3. For this test, the calculations are only done for the S2 surface.

Contrary to the preceding test, where the 10 selected nodes could be locked, the clamping chuck can only be applied to one of the 8 possible nodes (nodes 9 to 16); the spindle stopper can be applied to one of the 8 nodes, numbered from 1 to 8. The test remains easy and calculates the genetic algorithm's behavior because the optimal solution can be forecasted and the number of possible solutions is limited, which would not be the case with a thinner mesh. The same type of code that in the previous test can be taken, that is to say that the first 8 chromosome digits concern the nodes 1 to 8, but the possible number of 1s in the algorithm is limited to 1 in this part of the chromosome; the 8 following digits are for the nodes 9 to 16, but any chromosome having a number of 1 greater than 1 in this part will be removed from the process. For example, 1000000001000000 is an acceptable chromosome. This code type has not been kept for this example because it leads to 16-digitlong chromosomes and assumes the genetic algorithm is modified. Another type of possible coding is to build a 2-decimal-long chromosome; for example, 29 means that the nodes 2 and 9 are subjected to boundary conditions; the first digit varies between 1 and 8, and the second one between 9 and 16. The chosen code in this example uses 6 binary digits, as in 100011. The decoding is done as followed (let's recall that the decoding program and that of modification of the finite element data file is to be designed by the user for each new example and must be placed immediately before the analysis): The first 3 digits give the code of the forced node 1 to 8, according to the following correspondence: 000 (node 1), 001 (node 2), 010 (node 3), 011 (node 4), 100 (node 5), 101 (node 6), 110 (node 7), 111 (node 8), and the following 3 digits provide the forced node 9 to 16 code, according to the same type of correspondence; therefore the example 100011 matches to the forced nodes 5 and 12. Of course, this example is still an easy test because only 64 combinations are possible and they can all be calculated to reach the problem optimum. The objective is to minimize the maximum equivalent deformation or the equivalent Von Mises stress that appears where forces are applied. The best solution found by the genetic algorithm is the combination of nodes 8 and 16, for which the Von Mises stress equals to 17.009 daN/mm 2 . This result is found after a dozen finite element calculations (and from the second generation for a 6-individual population). In contrast, the genetic algorithm can be instructed to find the less good solution: the program is launched again with the objective of maximizing the main Von Mises stress; and this less good solution is the combination of nodes 1 and 9; for which the Von Mises stress equals to 17.195 daN/mm 2 . The test remains easy because the mesh size is limited. It could be more complicated if the mesh was thinner and if the genetic algorithm was instructed to find a compromise solution that would be valid to machine the S1, S2 and S3 surfaces.

Test 3

This test takes up the dynamic test offered in [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF] and helps in validating the implemented strategy, once again on an easy case. The chosen example is that of a square plate, measuring 30.5 cm with a thickness of 0.328 cm in deflection vibrations (Young's modulus 73.1 GPa, density 2,821 kg/m 3 ). This plate rests on 4 points that are located symmetrically on the diagonals (figure 3). The objective is to find the optimal location of the supports, maximizing the first fundamental frequency. Figure 3 In [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF], this problem is solved with a conventional gradient method, from a calculation of the frequencies' sensitivities with regard to the boundary conditions location. Since this is only the first symmetric mode of deflection, only a quarter of the plate is meshed. [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF] finds two equivalent optimal points (A and B on figure 3) that correspond respectively to frequencies of 169.46 Hz and 169.67 Hz. Actually, In [START_REF] Son | Optimization of boundary conditions for maximum fundamental frequency of vibrating structures[END_REF], it is only used a 36-element mesh for the whole plate and a study with thinner meshes has shown that the optimum is actually located between the A and B points. This test is often used in the literature; all the authors find that the optimal point is located between A and B.

Implementing a genetic algorithm strategy assumes that the support point location is coded under a chromosome on the main diagonal of the quarter of the plate. With the chosen mesh, that is 15 X 15 elements, there are only 16 possibilities that can all be calculated to get the reference solution that actually corresponds to the points 7 and 8, with frequencies around 205 Hz. The code of the 16 possible points is simply a binary one: 0000 corresponds to the node 1, 0001 to the node 2, 0010 to the node 3 and so forth until the node 16 (1111). Let's recall that the objective is to maximize the first frequency. For a 4-individual population, a 2individual population (0110, node 7 and 0111, node 8) is obtained after the thirtieth generation, which shows the genetic algorithm convergence, but its efficiency is more convincing with longer chromosomes (as in the example 1).

IV. Conclusion

This study has shown the efficiency of genetic algorithms in responding to the problem of the optimization of the boundary conditions in finite elements. This study is above all a feasibility study and will soon be complemented by industrial examples. The study can easily be spread to other fields than mechanics; for example, in thermal science, it would be easy to design chromosomes containing not only the information on the boundary condition type, but also that regarding the condition value to be optimized (flow value, heat transfer coefficient value); it could also be applied in fluid mechanics. The efficiency can still be improved in the case of important calculations (e.g. shape optimization), using neural networks to analyze the problem, instead of using a conventional finite element analysis. Actually, the use of neural networks to model mechanical structures appears to give good results [START_REF] Marcelin | Cognitive optimization of mechanical structures[END_REF]. It would be possible to make the learning of a neural network in parallel with the first generations that would be calculated by finite elements (that is using the results of the finite element analyses). Once the learning stage is over, the neural network would completely replace the calculations by finite elements. The calculations would therefore become much faster, and the genetic algorithm method, contrary to the deterministic methods, does not need extremely precise calculations of the objective function.

Translated by Amandine MARCELIN, on behalf of AMTrad'gram (www.amtradgram.com).

Figure 1 Test 2

 12 Figure 1 Test 2

First

  Author, Second Author, and Third Author Copyright © 2006 Praise Worthy Prize -All rights reserved International Review of Mechanical Engineering, Vol. xx, n. x

Figure 2

 2 Figure 2