Couro Kane
email: fcouro.kane@polytechnique.fr

Marc Schoenauer
email: marc.schoenauerg@polytechnique.fr

Topological Optimum Design using Genetic Algorithms

Structural topology optimization is addressed through Genetic Algorithms: A set of designs is evolved following the Darwinian survival-of-ttest principle. The goal is to optimize the weight of the structure under displacement constraints.

This approach demonstrates high exibility, and breaks many limits of standard optimization algorithms, in spite of the heavy requirements in term of computational e ort: Alternate optimal solutions to the same problem can be found; Structures can be optimized with respect to multiple loadings; The prescribed loadings can be applied on the unknown boundary of the solution, rather than on the xed boundary of the design domain; Di erent materials as well as di erent mechanical models can be used, as witnessed by the rst results of Topological Optimum Design ever obtained in the large displacements model. But these results could not have been obtained without careful speci c handling of the speci c aspects of topological genetic optimization: First, speci c genetic operators (crossover, mutation) were introduced; Second, special attention was paid to the design of the objective function; The nonlinear geometrical e ects of the large displacement model lead to non viable solutions, unless some constraints are imposed on the stress eld.

Introduction

Since the seminal work of Holland 27] and the comprehensive study of Goldberg 21], Genetic Algorithms (GAs) have gradually been recognized as powerful stochastic optimization algorithms. More recently, the initial framework of xed length bitstrings has been widened to other search spaces [START_REF] Michalewicz | Evolutionary Algorithms for Constrained Parameter Optimization Problems[END_REF][START_REF] Martin | Synthesis of optical multi-layer systems using genetic algorithms[END_REF][START_REF] Fogel | A comparative study of a penalty function, a repair heuristic and stochastic operators with set covering problem[END_REF][START_REF] Sch | Evolution strategies for mixed-integer optimization of optical multilayer systems[END_REF], emphasizing the need for problem-speci c modi cations of the basic algorithms. The eld of Evolutionary Computation covers all alternate evolutionary algorithms 49, [START_REF] Dejong | Are genetic algorithms function optimizers ?[END_REF][START_REF] Anagnostou | A computational procedure for part design[END_REF].

The main interest of stochastic methods in Engineering Sciences is to break the limits of standard deterministic methods in many optimization problems: when the search space involve both discrete and continuous domains (e.g. for the optimal design of truss structures); when the objective function or the constraints lack regularity; or when the objective function admits a huge number of local optima. In counterpart, stochastic methods are computationally expensive: GAs for instance, are slower than classical optimization methods by about one or two orders of magnitude, when comparison is possible, i.e. when classical methods apply. This paper focuses on applying GAs to some well-studied problems in mechanical engineering, namely the structural topology optimization of cantilever plates. Section 2 brie y introduces the mechanical problem, and surveys previous works in structural topology optimization, discussing their advantages and limitations;. Section 3 presents the broad lines of Genetic Algorithms in the historical bitstring perspective, for the sake of completeness. Section 4 is devoted to the construction of the problem-speci c Genetic Algorithm. The chosen representation of structures is discussed. The standard genetic operators are then tailored for topology optimization. Finally, the design of the tness function is thoroughly discussed. Experimentations in the linear elasticity case are presented in section 5: After validation results on the standard cantilever plate problem, this problem is modi ed to highlight the ability of the GA-based approach to handle problems having multiple optimal or quasi-optimal solutions. Multi-loading results are then presented on the problem of optimizing the structure of a bicycle, and nally, the problem of optimizing the shape of an underwater dome, involving loading on the unknown boundary of the structure, is addressed and solved. Section 6 presents the rst results of nonlinear topological optimization, in the context of linear elasticity with large displacements. The nonlinear geometrical e ects clearly show the need to take into account the stress eld in the tness function to obtain realistic solutions.

The out breaking results on Topological Optimum Design presented in this paper have been obtained during the rst author's PhD work. Hence, the reader is referred to the PhD dissertation 31] for all technical details.

2 Mechanical Background

The problem

The general framework of this paper is the problem of nding the optimal shape of a structure (i.e. a repartition of material in a given design domain) such that the mechanical behavior of that structure meets some requirements (e.g. a bound on the maximal displacement under a prescribed loading). The optimality criterion is is here the weight of the structure, but it could involve other technological costs.

Throughout this paper, except in Section 6, the mechanical model will be the standard twodimensional plane stress linear model, and only linear elastic materials will be considered (see [START_REF] Chapman | Genetic algorithms as an approach to con guration and topology design[END_REF][START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with GAs[END_REF]). The e ects of gravity are neglected.

One of the well-known benchmark problem of Optimum Design is that of the cantilever plate: a rectangular plate is xed on the left vertical part of its boundary (both displacements are set to 0), and the loading is made of a single force applied on the middle of its right vertical boundary. Figure 1 shows the design domain for the 2 1 cantilever plate problem.

Related Works

The main trends in structural optimization can be sketched as follows.

A rst approach is that of domain variation 11] (also termed sensitivity analysis in Structural Mechanics). It consists in successive small variations of an initial design domain, and is based on the computation of the gradient of the objective function with respect to the domain. This approach has two major defects: rst, it requires a good initial guess, as it demonstrated unstable for large variations of the domain; second, it does not allow to modify the topology of the initial domain (e.g. add or remove holes).

A more recent approach to topology design, introduced in 9], is that of homogenization; it consists in dealing with a continuous density of material in 0; 1]. In the end of this deterministic optimization, the current density is forced toward value 1 or 0, that respectively stands for material or void. However, this approach requires the design of the homogenized operator, as thoroughly described in 2], and is insofar limited to the linear elasticity case. Moreover, it cannot address loadings that apply on the actual boundary of the shape to be determined (e.g. uniform pressure, as described in section 5.4), and hardly handles optimization for multiple loadings.

A possible approach to overcome these di culties of topological optimum design is to use stochastic optimization methods, such as simulated annealing 34] and genetic algorithms 27]. Both methods have been successfully applied to other problems of structural optimization: in the framework of discrete truss structures, for cross-section sizing [START_REF] Goldberg | Genetic algorithms with sharing for multi-modal function optimization[END_REF][START_REF] Grierson | Discrete optimal design using a genetic algorithm[END_REF][START_REF] Leriche | Optimization of laminate stacking sequence for buckling load maximization by genetic algorithms[END_REF][START_REF] Schoenauer | Evolutionary computation at the edge of feasibility[END_REF] among others, as well as for topological optimization 26, 24, 53] and for the optimization of composite materials 35].

More recently, some problems of structural components optimization have been addressed by stochastic methods: Simulated Annealing is used to nd the optimal shape of the cross-section of a beam, a simple problem on which interesting theoretical results are proved in [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF][START_REF] Fogel | Evolutionary Computation. Toward a New Philosophy of Machine Intelligence[END_REF]; and GAs are used to solve the cantilever problem presented in Section 2.1 in [START_REF] Holland | Adaptation in natural and arti cial systems[END_REF][START_REF] Chapman | Genetic algorithm-based structural topology design with compliance and topology simpli cation considerations[END_REF]. However, the work presented in this paper goes further than these latter works: The geometrical constraints are weakened, increasing the range of possible solutions (Section 4.5.1); The tness function is carefully designed, and a ne control of the mechanical behavior of the solution is thus possible (Section 4.5); Last, this paper is not limited to feasibility results on problems that can be solved by the homogenization method, as demonstrated by the breakthrough results presented in Section 5 and 6.

Genetic Algorithms

This section gives the broad lines of basic GAs; the reader is referred to 21] or 39] for further details.

Historical GAs

Given a search space E and a tness function F de ned from E onto IR + , GAs evolve a set of p individuals (points of E), termed population. This evolution crudely mimics the Darwinian evolution: according to the Darwinian survival-of-the-ttest principle, the ttest individuals, i.e. near-optimal points of tness function F will appear in population P i for some i. The basic step in GAs, called generation, is the transformation used for population P i to give birth to population P i+1 . This transformation involves four steps : Evaluation: the value of the tness of all individuals on the current population is computed. Note that this step involves p independent computations that can easily be parallelized.

Selection builds population P 0 i by copying elements of P i ; the number of copies of an individual increases with its tness, the total number of elements in P 0 i being same as in P i .

Crossover applies on population P 0 i to build population P i ". From two individuals x and y in P 0 i , crossover builds two o springs x 0 and y 0 with probability p c (p c usually varies from .2 to 1.). When considering a bitstring representation (E = f0; 1g N ), a crossover c can be represented as a bitstring itself, c = (c 1 ; : : : c N ):

x 1 :: x N y 1 :: y N ! x 1 0 :: x N 0 y 1 0 :: y N 0

with x i 0 = ( x i if c i = 1 y i if c i = 0 and y i 0 = ( y i if c i = 1 x i if c i = 0
Most authors only consider one-point crossovers, corresponding to masks (1::10::0), or two-point crossovers, corresponding to masks (1::10::01::1). The general case represented here is called uniform crossover 52].

Mutation applies on population P i " to build population P i+1 . Mutation transforms an individual x in P i " into an o spring x 0 ; when considering a bitstring representation, a mutation can similarly be represented by a bitstring m = (m 1 ; : : : m N ):

x 1 :: x N ! x 1 0 :: x N 0 with x i 0 =

(

1 x i if m i = 1 x i if m i = 0
The probability for m i to take value 1 is noted p m (p m usually varies from 10 2 to 10 4 ).

The representation issue

One of the main di culty arising when applying GAs to some optimization problem is the choice of the search space. Two di erent spaces have to be considered. While the tness function is de ned on the phenotype space, genetic operators usually apply in the genotype space, be it the space of bitstrings, as described above, or any other space, as in modern Evolutionary Computation works. The representation of a phenotype in the genotype space involves a mapping (or coding) which can cause a loss of information (it is generally not bijective, not isometric, : : :).

The choice of the genotype space is therefore related to a compromise between the simplicity of the coding and the possibility to design useful genetic operators in the genotype space. Whereas standard GAs emphasize the use of xed bitstrings as genotypes, other trends of Evolutionary Computation try to use the same space for both the genotypes and the phenotypes, the e ort being then put on the design of genetic operators. A thorough discussion on the topic can be found in 19]. Section 4 will provide an instance of such a situation.

Pros and Cons

When considered as function optimizers (which they are not in the rst place 18]), the main advantage of GAs is to be a zero-th order method: The only prerequisite is to be able to compute values of the objective function, and of possible constraints of the problem. Furthermore, no regularity is required, neither on the functions nor on the search space. GAs can hence be used for continuous parameter optimization, for totally discrete problems as well as for mixed integer-continuous optimization 38, 5], provided genetic operators are de ned on the search space. However, the design of good genetic operators is still a matter of experience, and a posteriori numerical experimentation remains the only possible validation, though some general guidelines have been stated 42]. Section 4.4 will give examples of speci c operators, justi ed by the performance of the resulting algorithm.

On the other hand, it is well-known that the main drawback of GAs is their slowness. In particular, when both a GA and some deterministic method (e.g. a gradient method) can be used on the same problem, the latter is usually some orders of magnitude faster 48]. Parallelization certainly brings a partial answer to that issue. However, it only distributes, and does not reduce, the CPU requirements needed for a successful GA-based optimization problem solving.

Moreover, GAs involve numerous user-de ned parameters that have to be tuned carefully to get the best out of the algorithm. Here again, that tuning can only rely on the experience of Prescribed loading(s)

F.E.M. Nevertheless, these drawbacks of GAs are negligible whenever they are the only way to a solution of the problem at hand, handling problems beyond the limits of other optimization methods. The following of the paper is devoted to such successful applications of GAs.

Tentative shape

Simulated behavior

Fitness of current shape

Genetic Optimum Design

This section will detail the implementation of GAs on the problem of Optimum Design presented in section 2.1, describing all necessary steps toward the successful implementation of GAs on this problem.

Tentative tness

As emphasized in section 3.3, GAs only require values of the objective function. Hence, the broad lines of the implementation of GAs for the problem of Optimum Design described in section 2.1 can be represented as in Figure 2: The mechanical behavior of a structure can be numerically simulated by some Finite Element Method. From that simulation, and from the weight of the structure, a function to optimize can be designed (handling the constraints on the mechanical behavior by using penalty functions, for instance). This function can be used as the tness function of a Genetic Algorithm. However, the main question remains the choice of a search space: the target space is that of partitions of the design domain into two subsets (material and void). However, this general space is far too large (for instance, for obvious mechanical reasons, only structures with continuous boundaries between material and void need to be considered). Unfortunately, there is no natural representation for partition of continuous domains.

A priori Discussion for Shape Representations

A theoretical framework has been developed by Ghaddar, Maday and Patera 20] in the same context of Structural Optimum Design, but on the problem of optimizing the cross-section of a beam submitted to a bending momentum. The search space is restricted to partitions with polygonal boundaries. Theoretical results of existence and uniqueness of a solution are proven, approximation spaces are introduced and corresponding approximation results are obtained. Though the objective function considered in this paper is quite di erent from the one in 20], the same class of search space will be used here.

However, a signi cant di erence between the objective functions in 20] and the one to be used here is that the topological Optimum Design problem requires some Finite Element Analyses (FEAs) on the direct problem to compute the tness of a possible solution (i.e. a given repartition of material in the design domain), as shown by Figure 2, and detailed in forthcoming section 4.5. It is well-known that meshing is a source of numerical errors 17]. Hence, when comparing two structures of di erent shapes, using a tness function based on the outputs of two FEAs performed on di erent meshes is bound to failure, at least when the actual di erences of behavior becomes smaller than the unavoidable numerical noise due to remeshing. The use of the same mesh for all FEAs, at least inside the same generation, is thus mandatory in order to obtain signi cant results.

The bitarray representation

Once the decision to use a xed mesh has been taken, and with even very little knowledge of GAs, the straightforward representation for a partition of the design domain is that of bitstrings: each element of the xed mesh belongs to either one of the subsets of the partition, which can be symbolically labeled 0 or 1. Figure 3 is an example of a \chromosome", together with the corresponding repartition of the material, for a regular mesh of the design domain into quadrangles.

At rst sight, the resulting representation can be viewed as a bitstring, and the Optimum Design problem seems to have been brought back into the historical framework of GAs 27,[START_REF] Ghaddar | Analysis of a part design procedure[END_REF], where the search space is the space of xed-length bitstrings. Hence, all previous works using Genetic Algorithms on Optimum Design problems did use that representation [START_REF] Holland | Adaptation in natural and arti cial systems[END_REF][START_REF] Chapman | Genetic algorithm-based structural topology design with compliance and topology simpli cation considerations[END_REF].

However, this simple approach is far from optimal, as will be shown in next section.

Bitarrays are not bitstrings

This section is devoted to the study of the crossover operators for the representation of shapes given in section 4.3. All crossover operators presented in this section will be symbolically presented on Figure 4: two o springs of a black parent and a grey parent are plotted. The color of the bits of the children only tells which parent this bit comes from, regardless of its actual value.

Geometrical bias of one-dimensional crossover operators

In order to evaluate the bias induced by handling a bitarray as a bitstring, the standard (bitstring) crossover operators are rst considered.

One-point and two-points crossovers are geometrically biased, as they can only exchange horizontal bands of genetic material between the parents. Figures 4-a and4-b symbolically witness that phenomenon. A detailed Schema analysis can be found in [START_REF] Kane | Structural topology optimization in linear and nonlinear elasticity using genetic algorithms[END_REF][START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire[END_REF], emphasizing this bias. On the opposite, the uniform crossover (random mixture of both parents' genetic material not presented in Figure 4) does not su er from such a bias.

Nevertheless, two speci c two-dimensional operators are introduced in next section to address this issue.

Speci c two-dimensional crossover operators Diagonal crossover

The basic idea of diagonal crossover is to generalize the popular one-point crossover to the

The block crossover

First introduced in 28], the idea of the block crossover is to cut the whole two-dimensional domain by two horizontal lines and two vertical lines and to exchange some of the large blocks de ned by these lines. The values v 1 and v 2 (respectively h 1 and h 2 ) are chosen uniformly, and the 2 (or 3) blocks to be exchanged are selected randomly. Figure 3-d shows an example of 2-block crossover.

Experimental Comparison of Crossover Operators

These crossover operators have been experimentally compared on the problem of the optimization of the cross-section of a beam 20, 33]. Figure 5 shows the results obtained without any mutation, and with a crossover rate of 1.

The rst conclusion that can be drawn from Figure 5 is the ine ectiveness of the onedimensional crossover operators, signi cantly outperformed by both uniform and two-dimensional crossovers.

The second remark is that the uniform operator performs poorly in the early stage of the evolution, before catching up in the late generations. This fact is con rmed by all other experiments. This phenomenon is fairly general: in the beginning of the evolution, the uniform crossover equivalently disrupts emerging schemas, while all other crossovers (even one dimensional) preserve some large areas. These disruptive e ects do not occur in the end of evolution, since neither kind of crossover does disturb regions where convergence already occurred: If the bits of both parents at a given position are the same, crossing over has no e ect on those bits. In the meantime, the one-dimensional crossovers fail to bene cially exchange potentially good vertical parts (e.g. long vertical bars).

In all experiments of section 5 and 6, the block-3 crossover was used (it demonstrated slightly better performances than diagonal, block-2 and uniform crossovers) 31].

Mutation operators

The standard bit-ip mutation applies on bitarrays without inducing any geometric bias. Nevertheless, two other directions are explored regarding the mutation operator of bitarray shapes. The rst one is problem-independent, and uses statistics on the whole population to keep some genetic diversity. The second mutation operator, purposely devised for the problem of shape optimization, favors small modi cations of the boundary of the shape.

Population-based mutation

This mutation operator, rst de ned in 10] aims at preventing the premature convergence of the population. The probability of ipping a given bit is adjusted by considering all bits in the same position in the whole population. More precisely, the probability p i for a given bit i to be ipped depends on the mean value m i of this bit over the population: if this bit takes a uniform value (m i = 0 or m i = 1), the probability to mutate is set to a high value p max ; on the opposite, it is set to a low value p min if there is about the same proportion of 0s and 1s (m i = 0:5). The probability of ipping a bit is then a parabolic function of the mean value m i between these points. This operator thus imposes high values of mutation rate at positions that have already converged.

Boundary mutation

The underlying idea of the boundary mutation is that, for a given topology, i.e. a given number of "holes" in the shape, the optimal design for that topology can be found by slightly moving the boundary. The boundary mutation is de ned such that boundary bits, i.e. bits having one edge on the boundary of the connected component of the shape, are given higher probability (p B;max ) to be ipped than the other bits (p B;min ).

The population-based mutation performed slightly better than the standard bit-ip mutation (though many trials did no show signi cant di erences), and therefore was used in all experiments presented below. The boundary mutation will only be used in the end of the evolution process, to ne-tune the solution. More details can be found in 33, 31].

Fitness computation

A given individual is evaluated in a 2-step process. The rst step is a geometrical analysis of the repartition of the material in the design domain, followed by some Finite Element Analyses on the actual structure de ned by this partition.

Geometrical analysis

Some seed material is imposed at the point(s) where dynamical loading is applied. The connected component containing that seed is computed. Grid elements are connected if and only if they share an edge. Note that no seed material is prescribed on the part of the boundary where the plate is xed, as it is in 15]. The optimization process \chooses" where to hang the structure on the xed vertical boundary. This allows for a greater exibility in solving the optimization problem, as witnessed in the range of alternative solutions proposed for the modi ed cantilever plate problem (see Section 5.2). On the other hand, more structures are likely to be disconnected from that xed boundary, leading to a mechanically ill-posed problem: such structures are arbitrarily assigned zero tness, and are therefore eliminated by the next Darwinian selection.

Another di culty arises when considering the di erent connected components of the structure: Only material actually connected to both vertical boundaries takes part participate to the mechanical behavior of the structure. The disconnected parts have to be removed before the the FEM is called (see the example of Figure 3-b). But there are di erent possible ways to handle this di culty from a GA point of view:

The disconnected parts can be simply removed during the computation of the tness. However this introduces a fairly high level of degeneracy in the representation: Many genotypes correspond to the same phenotype, and hence have exactly the same tness. And, as pointed out in 43], this is not a desirable feature, as there is no preferred direction for the GA to go on these large plateaus of tness landscape.

The disconnected parts can be de nitely removed from both the phenotype and the genotype, in a Lamarckian-like way. This is amenable to the so-called repair technique used in genetic constraint handling 40]. But these disconnected parts possibly are valuable genetic material (like the dominated part of diploid chromosomes), at least in the beginning of evolution: This strategy generally leads rapidly into some local minima. Experimental results account for that statement.

The tness can be modi ed to slightly favor the disappearing of such disconnected material: the Finite Element Method is run without the disconnected parts, and some penalty term is added, relative to the amount of disconnected material, in order to guide the algorithm toward better solutions. This latter possibility will be used throughout this paper. Next section will give the precise formulation of the tness (equation 1), including the penalization term to cope with the unconnected material.

The penalized tness

Once the component connecting the seed(s) and the xed boundary has been computed (if any), one FEA is performed for each one of the considered loading cases (generally one, except in multi-loading problems, see Section 5.3), using the same regular quadrangular mesh that supports the representation of the individuals. The FEM tool is detailed in 30]. Note the actual material boundary is used, in contrast with both 28] and 15] in which the FEA was done on the whole design domain, after assigning a very low Young modulus to void elements. Though the results of the analysis does not di er signi cantly from one method to the other, computing the actual boundary allows to take into account loading applied on this boundary (e.g. pressure, heat exchanges) as demonstrated in section 5.4.

Di erent criteria have been used in previous works on Optimum Design. The objective function for the homogenization method is based on the compliance of the structure, computed as the work of external loadings 2]. The main advantage of this objective function is to be di erentiable. However, it does not allow a precise control of the mechanical behavior of the solution. As GAs do not require di erentiable tness function, the compliance will not be considered here. It has however been considered in the GA framework [START_REF] Kane | Structural topology optimization in linear and nonlinear elasticity using genetic algorithms[END_REF][START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire[END_REF] in order to make more precise comparisons between the homogenization method and the GA-based method proposed here [START_REF] Kane | Algorithmes g en etiques et Optimisation topologique[END_REF][START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire[END_REF].

Yet another approach is to try to maximize the sti ness of the structure 15, 14]. Here again this does not allow a ne control of the desired behavior of the solution. Moreover, comparative experiments 31] demonstrate that using the sti ness as tness function generally leads to heavier solutions without much improvement of the mechanical properties.

Hence, the following of the paper will address the problem described in section 2.1 of the minimization of the weight of the structure subject to some upper limits on the maximal displacement of the structure when subject to prescribed loadings.

This problem is a constrained optimization problem. Many methods have been designed to constrained evolutionary optimization (see 47] for a survey of such methods). However, the standard penalization method remains the simplest one to use, as it only requires a modi cation of the objective function. The drawback lies in the di culty of adjusting the penalization parameters: Section 4.6 below will discuss that issue in in the context of this paper. Whatever the choice of penalization parameter, the formulation for the penalized tness is the following:

F = 1 A con + "A dis + (D Max D Lim ) +
(1) where A con and A dis are the area of the connected and disconnected material in the repartition (see section 4.5.1 above); D Lim the imposed limit value for the displacement; D Max is the maximal displacement of the structure when the prescribed loading is applied (computed using the FEM); " and are positive user-supplied penalty parameters (a + denotes the positive part of a).

Note that " is also a penalization parameter, but whose value does not in uence much on the performance of the algorithm: it was xed to 0:1 throughout the experiments described in the following sections.

Adjusting the penalization parameter

As quoted above, using penalization method to handle constraints raises the di culty of tuning the strength of that penalization, i.e. in the case of tness function 1, adjusting precisely the parameter .

The drawbacks of a xed value of are well-known (see 41] for a more thorough discussion):

A small value of may result in an optimal solution that violates the constraints; A large value of ensures that the constraints will be strictly met, but forbids exploration and short-cuts in infeasible regions that might be essential for the overall success of the algorithm, at least during the beginning of evolution. Hence a natural idea is to use for a dynamic schedule, starting from small values, and thus letting the population explore even infeasible regions, and gradually enforcing the penalization by increasing such that the constraints are nally met by the whole population. Two di erent solutions to achieve this goal have been tried on the problem of Optimum Design.

The rst method uses a prescribed geometrical evolution schedule for , and is termed exogenous scheme (similar dynamic penalty terms have been proposed in 29]). The value of is increased by some multiplicative factor every M generations. Hence the value of at generation i is given by

i = 0 i M ]
Typical values of M and are 10 and 1:001, while the initial value 0 is computed from the average weights and violations of the constraints in the initial population.

The second scheme is an adaptive scheme for , that uses the information contained in the current population to compute further values, in the line of the adaptive penalties proposed in 8, 50]. Some i is computed from the actual average weights and violations in population i, and is set to i if i is greater than the current value of , to ensure monotonicity of .

Comparative experiments between the xed scheme, the exogenous geometrical scheme and the adaptive scheme were performed on di erent instances of the cantilever problem. The rst strong conclusion is that the xed scheme is constantly outperformed by both dynamic schemes. A typical example of such a situation is shown on Figure 6-a. Note that the plots of both dynamic schemes are hardly distinguishable. But another typical result is shown on Figure 6-b, demonstrating a slight improvement of the adaptive scheme over the exogenous scheme. However, the latter was used in all experiments of next sections, as the being much simpler to tune.

Results in Linear Elasticity

This section presents results obtained using the algorithm detailed in the preceding section in the framework of linear elasticity. Hence, as the rst results (Section 5.1) deal with the standard cantilever problem, they could have been obtained { much more rapidly { by the homogenization method. But after these validation results, the limits of homogenization methods are broken, and the presented results deal with multiple quasi-optimal solutions (Section 5.2), multi-loadings cases (Section 5.3) and nally address the problem of the underwater dome, involving loading on the unknown boundary (Section 5.4).

Validation Results

Figure 7-a is a typical result of a GA run for the 2 1 cantilever plate, discretized according to a 64 44 regular mesh. The population size for all runs is 125, and the number of generations arbitrarily xed to 2000. One run thus require about 150,000 FEM analyses, taking approximately 24 hours of a powerful HP workstation. The genetic operators are applied with the following probabilities: 0.6 for the block crossover (see Section 4.4.2), 0.1 for the population-based mutation, with p min and p Max set to 0.0001 to 0.01 (see Section 4.4.4). All these parameters were adjusted after exhaustive tests, and details can be found in [START_REF] Kane | Structural topology optimization in linear and nonlinear elasticity using genetic algorithms[END_REF][START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire[END_REF]. The result of Figure 7-a is to be compared to the result of homogenization of Figure 7-b. Both structures look alike, they have almost the same weight and same maximal displacement: these results validate the GA approach.

Moreover, the GA solution exhibits numerous small holes (of size one element) in the structure. This seems to indicate that the GAs somehow tries to approach the actual optimum, which is known to lie in the homogenized space 2]: The homogenized optimal solution can be viewed as the limit of a structure in which in nitely many holes of in nitely small size are drilled. Note however that this could have been avoided by imposing some lower bound on the thickness and the co-thickness of the structures, as in the theoretical study developed in 20].

Of course, complete validation of a stochastic method involves to check its robustness with respect to the random initialization. Di erent runs of the same program with di erent random have been performed: All lead to similar structures, but slight di erences could be observed in their appearances, though their weights and maximal displacements were very close. GAs are known to be able to nd quasi-optimal solutions on multi-modal problems [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF][START_REF] Lin | Genetic search strategies in large scale optimization[END_REF]. This latter remark suggested to consider a problem with known multiple solutions.

Multiple quasi-optimal solutions

The problem considered in this section is a modi cation of the standard cantilever plate problem to ensure the existence of multiple quasi-optimal solutions: the 1 4 cantilever plate is discretized according to a 10 40 mesh and both its left and bottom boundaries are xed (see Figure 8-a).

Depending on the height of the point where the loading is applied and on the constraint on the displacement, the same problem can have multiple solutions. A simple example of such situation is given in Figure 8 (b) and(c): if the loading is applied at height 10, and provided that the displacement constraint is large enough, both structures (b) and (c) are optimal solutions.

And indeed, as no material is prescribed on the xed boundary (see Section 4.5.1), the GA-based algorithm was able to nd both solutions.

Moreover, if the height of the loading point is xed, and the displacement constraint D Lim is gradually relaxed, di erent quasi-optimal solutions exist for some ranges of D Lim . The GA method was then able to nd out di erent multiple solutions, as demonstrated in Figure 9. Here again only the most signi cant results are shown, for di erent values of D Lim , while the loading is applied at height 15. The GA parameters are those of section 5.1, except for the population size (100) and the maximal number of generations (500 The gures below the structures are respectively D Lim (the limit on the displacement), D Max (the maximal displacement of the structure) and A (the area of the structure).

Multiple loadings

All problem addressed in the preceding sections considered constraints on the mechanical behavior of the structure in a single case of loading. But actual structures in real-world problem are generally subject to di erent loadings, and hence should be optimized while taking into account more than one loading case.

Modi ed tness

Taking into account this new situation is straightforward in the context of GA-based Optimum Design: for each structure, one Finite Element Analysis is performed for each of the loading cases, and the tness function de ned by equation ( 1) is modi ed to agglomerate the constraints corresponding to the L loading cases:

F = 1 A con + "A dis + P i=L i=0 i (D i Max D i Lim ) + (2)
where A con , A dis and " are as in equation 1; D i Lim is the limit value for the displacement for loading case i; D i Max is the maximal displacement of the structure when loading i is applied; i is the penalization parameter corresponding to the i th constraint; all i are adjusted as described in Section 4.6.

Of course, the computation of the tness of a single structure requires in that context L FEAs, and the computational time is therefore L times larger than when a single loading is considered.

An example

An example of such a situation is given by the structure of a bicycle: the forces applied to that structure are very di erent depending on the position of the rider, which in turn heavily depends on the slope of the road. For the sake of simplicity, only the following three di erent cases have been considered: on at landscape, the greatest force is applied on the saddle, on up-hill ground, the rider pushes hard on the pedals and pulls on the handlebars and on steep roads, he does not sit on the saddle any more, pushing extremely hard on the pedals.

Figure 10 shows the result of the three single-loading optimizatio problem, while Figure 11 is the solution of the Three-loadings optimization problem. The advantages of solving the multiloading problem appear clearly when comparing these two gures, as the resulting multi-loading solution di er signi cantly from any of the three solutions of the single-loading problems.

Another approach to the same problem could be to use multi-objective optimization methods. Speci c modi cations of GAs were made to handle multi-criterion optimization, taking advantage of the ability of GAs to nd mutiple optima [START_REF] Radcli | Fitness variance of formae and performance prediction[END_REF][START_REF] Smith | Genetic optimization using a penalty function[END_REF]. The aim of the resulting algorithm is to sample the Pareto set. To our knowledge, such approach was not tested yet on the multi-loading Topological Optimum Design problem.

Loading on the unknown boundary

All problems addressed up to this point dealt with xed loading, independent of the structure at hand: the loading(s) was (were) applied on the boundary of the design domain. There are situations, however, where some loading is applied on the actual boundary of the structure, which can be di erent from the boundary of the design domain. This is the case if one some uniform pressure is applied on the upper boundary of a civil engineering structure, or, in a di erent context, when optimizing the shape of heaters, which implies heat transfers through the unknown boundary of the solution 20].

This situation is of course intractable for the homogenization method, as the boundary of the target structure is not de ned precisely before the projection step (nor sometimes is it precisely de ned after that step either). From a GA perspective, this situation does not di er from the situation where all loading are applied on the boundary of the design domain. The boundary of any structure in the population is well-de ned, and any loading can be applied there onto. Hence no modi cation of the algorithm is required to handle this case.

Consider the problem of optimizing the shape of an underwater dome: a uniform pressure is applied on the upper boundary of all structures. Some material is imposed at both ends of the lower boundary of the design domain, as well as at some point above the middle, to ensure that a minimum height of the structure1 .

Figure 12 shows the resulting structures, for both a small pressure (12-a) and a large pressure (12-b). The latter structure is indeed a better solution than what was expected (a structure looking like a reinforced 12-a structure) as is certainly is lighter: This only shows that the material requirements that were imposed here were not the right ones: some void area should have been imposed in the middle of the design domain.

Nonlinear Geometrical E ects in Topology Optimization

This section considers standard plane stress problems in the context of large displacements. The material still obeys a linear law (the extension to any other constitutive law is straightforward), but the nonlinear geometric e ects due to the large displacement hypothesis are taken into account. A thorough description of the theoretical model together with the numerical algorithm can be found in 16]. Details on the numerical model and implementation used here can be found in 30].

The rst experiments used the penalized tness function described by equation (1). The on the displacement). In the purely linear case, all such problems are of course equivalent.

A typical result of the GA based optimization in the large displacement context is shown in Figure 13-b, together with the result on the same problem in the pure linear model (Figure 13-a): one does not need to be a mechanical expert to see that such a structure is a disaster, from a mechanical point of view. And di erent values of F and D Lim did produce similar results. And a closer look at the maximal stress of the solutions Max con rmed the existence of a problem there: The value of Max in Figure 13-b are far too high.

First, the large displacement model can have di erent solutions for the same loading, some with higher stress than others. Second, the stress eld itself, on a rough domain like the structures obtained for the standard cantilever plate problem, and with such coarse discretization, presents some singularities.

As a matter of fact, Table 1 gives some idea of the nonlinear geometrical e ects in the case of two simple structures, the perfect ">" shape and the straight beam; both the displacement and the maximum stress are given for varying loads. This suggests to incorporate a constraint on the maximal stress in the tness function. Of course, some nice structures might be missed, like the perfect ">" shape for F = 910 3 . But, hopefully, another structure very similar to it can arise, with the same mechanical properties, and with the FEM model giving a reasonable solution among the possible solutions. Following these ideas, the tness becomes

F = 1 A con + "A dis + (D Max D Lim ) + + ( Max Lim ) + (3)
with the notations of equation 1, Lim being the constraint on the stress, Max the maximum of the stress on the structure, and some positive penalty parameter which has to be adjusted. Figure 14 shows the optimal designs obtained with that modi ed tness: Reasonable solutions are found, as long as the stress is imposed a strong constraint (computed from a purely linear numerical simulation).

Discussion and Perspectives

Breakthrough results have been presented in this paper, on simple problems of topological optimum design. They demonstrate the potentialities of the GA-based topological optimization. However, the computational cost remains the main limitation of the proposed method: using a middle-range HP workstation (HP-755), from 6 hours, for the simplest problem with a 10 20 mesh, up to 30-40 hours, for the three-loading 44 64 bicycle, were necessary for a single successful run. Of course, the parallelization of the tness computations is straightforward, and would result in a linear speed-up without modi cation of the underlying algorithm. But it remains that the overall computational requirements are very heavy.

Moreover, the accuracy of all results presented above derives from that of the Finite Element Analyses performed during the computations of the tnesses. And that accuracy is dictated by the size of the underlying mesh. All the results presented in this paper were obtained on rather F coarse meshes (from 10 20 to 44 64) whereas real-world problems and accurate analyses of the mechanical behavior of the shapes require much ner meshes. Increasing the size of the mesh would not only increase the cost of the tness computation, but also the size of the chromosomes. And increasing the size of the individuals would in turn require to increase both the size of the population and the number of generations to reach the same level of convergence (Cerf 12,[START_REF] Cerf | Une th eorie assymptotique des algorithmes g en etiques[END_REF] proved, in the bitstring case, the existence of a critical population size to ensure convergence in nite time; and this critical size increases linearly with the length of the bitstring). Hence, the bit-array representation poorly scales up when re ning the mesh | not to speak of handling 3-D shapes.

Therefore other representations have to be used in order to overcome this limitation, and dissociate the complexity of the representation and the accuracy of the evaluation: First results 45, 46] seem to con rm the interest of this direction of research.

Conclusion

The feasibility of GA-based optimal design had already been witnessed by 28] and 15]. This paper intends to emphasize on its exibility, and to demonstrate its potentialities to break some limits of up-to-date deterministic methods.

Multiple quasi-optimal solutions can be found, allowing to take into account some other criteria that could not be incorporated in some objective function (e.g. technological or esthetic criteria). Actually, the results in section 5.2 were obtained by successive runs of the GA starting from di erent random populations, thanks to the stochastic nature of the algorithm. However, it is possible by using some niching techniques (e.g. the sharing scheme of 22]), to locate in a single run more than one near-optimal solutions.

Results regarding the multi-loading optimization (Section 5.3 were the rst of their kind at the time they were obtained. Moreover, they required very few modi cations of the method handling a single loading. Since then, the homogenization method has been modi ed 1] and is now able to optimize structures for multiple loadings. However, deep modi cation of the numerical procedure was necessary: the modi ed homogenization algorithm relies on some local optimization performed in every element of the mesh, that was solved analytically in the single loading case, and which is now handled numerically, which dramatically increases its computational cost { though it remains faster than the GA-based algorithm. A fair comparison between stochastic and deterministic approaches should consider the development cost (e.g. the time necessary to work out the adjoint problem in inverse problem solving) together with the computational cost. Finally, even in the linear framework, stochastic optimization is the only method able to handle topological optimum design with loading on the unknown boundary of the target structure, as demonstrated in section 5.4. Moreover, from the GA point of view, it is a strict application of the original method { the only modi cation took place during the computation of the tness, as the external loading had to be computed for each structure.

Last but not least, in the context of large displacements, the experiments presented in section 6 acknowledge for the usual statement "the extension to other mechanical models is straightforward". Note however that the tness function had to be adjusted in order to take into account the maximal stress of the structure: there is no free lunch for nonlinearity.

However, these good results must not hide the main drawback of the method, namely its computational cost, which makes it highly unlikely to be applicable to real-world optimization problems. Nevertheless, and this is con rmed by preliminary results of on-going work, we are convinced that a change of representation can remedy to that weakness, allowing to treat di cult problems of real size. And the exibility of the method will then prove essential for some other yet unsolved problems of Topological Optimum Design.

Figure 1 :

 1 Figure 1: The standard 2 1 cantilever plate.

Figure 2 :

 2 Figure 2: A tentative tness for genetic optimum design

  bitarray genotype. b -The phenotype, or repartition of void/material.

Figure 3 :Figure 4 :

 34 Figure 3: Bitarray representation: GA applies on the bitarray genotype (a), while the tness is computed on the phenotype (b).

Figure 5 :

 5 Figure 5: Performances of crossover operators without mutation

Figure 6 :Figure 7 :

 67 Figure 6: Comparative results (averaged over 15 independent runs) of di erent schemes for the penalization parameter . (a) The xed scheme is outperformed by both dynamic schemes. (b)The adaptive scheme sometimes slightly supersedes the exogenous geometrical scheme.

  The problem Two optimal solutions Figure8: The 1 4 modi ed cantilever plate.

Figure 10 :

 10 The three loading cases for the optimization of the structure of a bicycle, together with the solution of the corresponding single-loading optimizatio problem. (a) Steady ground. (b) Heavy slope (no force on the saddle). (c) Sitting up-hill position.

Figure 11 :

 11 Figure 11: The multi-loading bicycle

  initial idea was to use di erent loadings F, with a xed ratio F=D Lim (D Lim is the constraint (a) (b) Figure 12: Results for the optimization of the underwater dome, with uniform pressure on the upper boundary. (a) For small pressure. (b) for large pressure.

Figure 13 :

 13 Figure 13: Optimal (disastrous) designs: Only displacement constraints were considered. F = 0:009 and D Lim = 0:02285.

  ). Optimal structures on problems with multiple solutions. F is applied at height 15.

	0.32 0.294 0.38	0.3927 0.33	0.4	0.3919 0.34	0.42 0.3415 0.26	7.66 0.15	7.66	6.743 0.19	46.6 46.64 0.1
	Figure 9:

Table 1 :

 1 The displacement is what was expected, but the stress does present weird values. Note that this phenomenon depends on the numerical model[START_REF] Chapman | Genetic algorithms as an approach to con guration and topology design[END_REF][START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with GAs[END_REF]. Nonlinear e ects on two reference shapes. D Max and Max are respectively the maximum displacement and the maximum stress.

	">" shape D Max Max 910 6 0.0002360 0.0005968 0.04161 0.03035 straight beam Load D Max Max 910 5 0.00235 0.0059 0.34499 0.28344 910 4 0.02286 0.05312 0.80763 1.272 910 3 1.06731 10.99 1.08208 1.11683 910 2 1.249 2.0763 1.6871 4.3787 910 1 2.9574 break 3.177 break

  Figure 14: Optimal designs for nonlinear elasticity with displacement and stress constraints.

	0.009 D Lim 0.22856 Lim 0.53	0.018 0.457 1.0622	0.09 2.2856 5.3
	D Max 0.2143 Max 0.550 Area 0.21	0.4504 0.9835 0.47	1.687 4.379 0.1

Otherwise, the optimal solution is the at structure along the lower boundary, which indeed has a small displacement. Note that, due to an error, one run was performed without this material requirement above the ground, and that an almost at solution was found by thealgorithm 31]

Acknowledgments

The authors are deeply indebted to Fran cois Jouve (CMAP) for the invaluable help provided all along this work, from patient theoretical explanations to modi cations of his Finite Element software to t our requirements.