Evolutionary identification of macro-mechanical models

This chapter illustrates the potential of genetic programming ( GP ) in the field of macro­ mechanical modeling, addressing the problene1 of identification of a mechanical model for a material. Two kinds of models are considered. One-dimensional dynamic models are rep­ resented via symbolic formulations termed ·rheological models, which are directly evolved by GP. Three-dimensional static models of byperelastic materials are expressed in terms of strain energy functions. A model is rated. based on the distance between the behavior predicted by the model, and the actual behavior of the material given by a set of me­ chanical experiments. The choice of GP is n:totivated by strong arguments, relying on the tree-structure of rheological models in the fi:rst case, and on the need for first and second order derivatives in the second case. Key issues are the exploration of viable individuals only, and the use of Gaussian mutations to optimize numerical constants.


Introduction
This chapter is concerned with applying G�P [Koza 1992] in the field of solid mechan ics. More precisely, the aim is to discover a model, or constitutive law, characterizing the mechanical behavior of a material from mechanical experiments.
The design of modern structures requires ever more detailed knowledge of the constitutive properties of materials. Suc1n knowledge is needed to predict through numerical simulations the behavior of the� structure under external loads. Reliable predictions allow for meeting the engineering requirements at a lower cost. A number of new materials (composite materials, polymers), have recently entered into wide usage and no accurate model of their behavior is so far available. The identification of a constitutive law for nLew materials therefore becomes a major challenge for mechanical science and industry.

State of the art
The design of an accurate law requires cwnsiderable insight in mechanics. When a new material is expected by experts to resemble a well known material, the model of the latter is adjusted. Otherwise, an ever stronger sense of mechanics is needed: Brand new models are elaborated through a time-consuming trial and error process, where the successive models guessed by mechanical engineers are checked against test experiments; these e:Kperiments may in tum suggest new mod els [Lemaitre and Chaboche 1985]. One c::an also start with a thorough analysis of the mechanical behavior of the material. at the microscopic scale; a macroscopic law is then derived from the microscopic model, for instance by homogenization (Sanchez-Palencia and Zaoui 1987]. However, such models often result in tremen dously time-consuming numerical simulations because of their complexity. Once the current law is found plausible, its numerical parameters are tuned by minimizing a distance between the observed behavior of the material and the behavior predicted from the law for a set of experimental conditions. This stage of identification is equivalent to parametric optimization (Gittus and Zarka 1985).
Much attention has been paid in designing experiments in order to adequately check a given law (Zarka & al. 1988). But searching for an accurate law remains a critical issue, as all above approaches fail in many cases: when the target material does not resemble any previous material; when the microscopic analysis does not allow a tractable model to be built ; when the (sometimes highly nonlinear) behavior of the material prevents the expert from fitting any appropriate model.

Goal
The goal of this chapter is to build a system able to provide a mechanical engineer with laws fitting data from available experiments. This amounts to non-parametric optimization, since it involves finding both the structure of the law, i.e., a set of partial differential equations, and the nun1erical coefficients in these equations. Two kinds of models are considered: • the one-dimensional dynamic aspects of the mechanical behavior are handled by means of rheological models (Persoz 1960). These models describe a material as an assembly of elementary elastic, plastic or viscous components. Only mixed rheolog ical assemblies composed of serial and/ or parallel branches are considered here. GP operates directly on such rheological models, represented by trees within an ad hoc set of functions and terminals. An interpreter of rheological trees then produces the equivalent system of partial differential equations.
• the three-dimensional static aspects of the mechanical behavior of hyperelastic materials are captured by a strain energy function [Green and Zerna 1968]. Energy functions are described via standard function and terminal sets; but not all trees represent valid energy functions and much attention is paid to exploring only vi able individuals, i.e., leading to a successful numerical resolution of the underlying system of partial differential equations.
In both cases, solving the system of partial differential equations gives the simu lated response of the material when external loads are applied. The fitness of a law is computed as the difference between tbe simulated and the actual responses, for a set of mechanical experiments. It is worth noting that an automatic building of such laws is beyond reach for all standard identification approaches. As far as we know, this chapter reports the first atten1pt to use GP to this end.

23.2
On e-dimensional rheological 1models This section addresses the building of general (visco-elasto-plastic) one-dimensional models using mechanical tests. We restrkt ourselves to finding the equation which links the stress function u(t) to the strain applied on this material £(t) and its time derivative i(t), of the form u(t) = .1"(£(t),€(t)).
The presentation will mainly focus on the genetic programm ing aspects of this work: First, the search space is describEd in terms of trees, based on a set of ad hoc nodes and terminals. An interpreter iis devised to compute the constitutive law corresponding to a given tree. Much care! is then spent in designing specific strate gies to address the optimization of the numerical constants and the exploration of the restricted search space.

The search sp ace
Rheological models allow one to describe 1most mechanical laws in the one-dimensio nal frame; they can be thought of as an assembly of elementary components, namely springs, sliders and dashpots, which resp, ectively symbolize elastic, plastic and vis cous behaviors. In particular, a rheological model, far from being a black box, provides a deep understanding of the constitutive properties of the current ma terial and can be rearranged by the expert. We restrict ourselves to rheological assemblies composed of serial and/ or parallel branches, which allow for describing most known materials. Figure 23.1a shows an example of such rheological model, which was proposed for polyethylene [Kic::henin 1992].
Such restricted rheological trees can be represented as tree structures built from two functions (respectively the serial and parallel assembly modes) and three terminals, (respectively springs, sliders and dashpots) (Figure 23.1b). In terms of GP, the function set includes the serial and parallel connectors, both of variable arity n > 2, involved in rheological models. The terminal set includes three elements, all involving a floating-point coefficient: a spring is characterized by its stiffness, a slider, by its stress threshold, and a dashpot, by its viscosity. A terminal thus (l ui(t)l < uf), it inhibits (makes inactive) all branches parallel to it and all branches in serial below it. The interpreter normally proceeds by recursively checking the rheological tree (given the underlying se1mantics of parallel and serial assemblies). This parsing is simplified via a syntactic restriction on the search space, or language bias; the children of a parallel (respectively seria� node are stipulated to be either terminals, or serial (resp. parallel) nodes.
Note that this language bias significantly reduces the amount of degeneracy of the representation, i.e. the number of genotypes (rheological tree) standing for a pheno type (mechanical behavior). Different aspects of degeneracy are discussed in the lit erature: degeneracy is blamed for misleading the search (Radcliffe and Surry 1994a); but special cases of degeneracy, such as induced by introns, are shown beneficial in limiting the disruptive effects of crossover [Levenick 1991], or allowing recessive information (not considered for fitness cmnputation) to be transmitted (Koza 1992]. Nevertheless, there is still room for deg�eneracy in the so-biased language: reduc tions similar to arithmetic simplification would be possible (e.g. spring k1 I I spring k2 =spring k1 + k2), and sliders may cause parts of the tree to become inactive, . .

The solver
The next step consists of computing the response of the model when external dy namical load is applied. This load is de8cribed as the sequence of strains applied at instants to, ... , tT, noted €ezp ( tj), j = 0, ... , T. The global stress O"o ( t ) attached to the root of the tree is incrementally computed for any t;: 1. Derivatives are expressed via finite differences: ei(tJ) = €i(ti]=.;;�:-d. 2. The loading history is taken into account: €o(t; ) = €ezp(t;) and io(tj) = iezp(tj)· 3. These equations, together with & R, a: re handled as a set of linear equations in the unknown O'i(tj ), €i(tj ), ei(tj ) and €i(i� j -1 ).
4. Solving this set of linear equations giv� es the state of the model at time ti, given its previous state at time t;-b the initiall state of the model being (0, 0, ... , 0).

5.
For each slider i, the internal stress O'i(tJ) is compared to the threshold uf; when ever a slider becomes saturated (lui(t; )I > uf ) or leaves saturation (lui(t; )I < u!), the system £ R is rebuilt by running the interpreter.
Complexity: The complexity of this reBolution amounts to T x 2(3N)3 13, where T is the number of time steps of the loading history and N is the size of the tree (the resolution of a linear system of size n being 2n3 13).
Error: This resolution proces8 involves two kind8 of error. First, the mechanical measures are known with a given precision at given instants only; and the discretiza tion of the loading history is beyond our control. Second our handling of derivatives induces numerical errors depending on the current model. We heuristically propose to estimate both kinds of errors for a git1en model, through the difference between the response computed from the whole loading history €ewp(t), t = t0, ••• , tT, and the response computed from an excerpt of the same loading history, including only one in two consecutive instants (a loading history typically includes several dozens or a few hundred points): Error = E!T t o luR(t, €e x p, to, t1, t2, ... , t T) -O'R(t, €ezp, to, t2, t4, ... , t T )I.
This estimate intends to capture both the imprecision inherent in the available data, and the error introduced by the solver. A run is considered successful if the fitness of the best model falls below this estimate of the unavoidable error.

Fitness computation
A fitness case is an experimental curve involving the loading history €ezp(t) and the observed stress O' ewp(t) , fort= t0, ••• , tT. A model R is evaluated through the difference between the observed stress O'e! wp(t ) and the stress uR(t, €ezp) computed from the model R according to the expeiimental loading history: The total fitness associated to a model is the sum of the errors on all fitness cases, i.e. on the available mechanical experiments (typically three).
As usual, most of the GP computationa l effort is spent in calculating the fitness.
And, since tournament selection (Goldberg and Deb 1991) is used, the evolutionary process is only dependent on the orderin; g of fitnesses (Andrews and Prager 1994]; this suggests several heuristics that could lessen the global computational cost, at least in the earlier stages of evolution: • Gradually increase the number of experiments to fit: evolution typically starts by considering a unique fitness case, and additional fitness cases are gradually taken into account when the current population meets some criterion, to be defined.
• Gradually increase the number of time-·steps taken into account: evolution starts by considering only the first K time-step8 of any loading history, and K is similarly incremented when the current population meets some criterion.
• Gradually increase the precision of loaeHng histories: only K evenly spaced time steps are first considered, with K increas1ing with the number of generations.
In all cases, the strategy consists of successively optimizing more and more pre cise approximations of the actual fitness.. A similar iterative strategy proved use ful for constrained optimization (Schoenauer and Xanthakis 1993]: first stages are concerned with sampling the feasible reg;ion. Only then, does evolution deal with the actual optimization problem. However, the successive optimization scheme addresses a basic need for constrained optimization, while here it aims at com putational savings. Moreover, our timing of the fitness sequence is more critical, because there is no guarantee as to the smooth convergence of the intermediate fitnesses toward the actual fitness. Failing to change the biased approximation of fitness frequently enough may result in overfi tting its bias; evolution then may get stuck in a local optimum.
The change of the fitness function em:p•irically takes place when the current best fitness falls below 75% of the initial best value for the current fitness function.

Evolution operators
The two standard genetic operators, cro ssover by random swapping of subtrees and mutation by random replacement of a subtree [Koza 1992), are modified for offsprings to meet the syntactic restrictions described in section 23.2.2.2.
Several factors in our problem indicat4ed that a third operator might be useful. Like numerical functions, rheological ass: emblies allow for values to combine (e.g. the series of two springs with respective stiffness kt and k2 behaves as a spring with stiffness ;;�t! ). The adjustment of numerical values could then be left to random mutation and crossover, as in (: Koza 1992]. However, this leads to a dra matic increase of the size of the trees alo: ng evolution; and the fitness computation increases as the cube of the size (sectiont 23.2.2.2). Such adjustment of constants is therefore much more expensive for rheological identification than for classical re gression problems. A specific operator was thus devised to address the optimization of floating-points values. This operator b• asically is a random hill-climber, involving an elementary evolution operator termed. surface mutation.
A surface mutation modifi es all floating-point terminals in a given model, via the addition of a Gaussian noise; the variance of the mutation is attached to the coeffi cient and recombined as in evolution stra1tegies [Schwefel 198 1). A surface mutation is successful if it results in an improved fitness. The random hill-climber repeatedly performs surface mutations, until a number 0 of consecutive surface mutations is found unsuccessful. The influence of parameter 0, termed stubbornness, on the dynamics of evolution is discussed in nex:t section.
Any crossover or (classical) mutation is: followed by a hill-climbing stage, in order to get the locally best version of the new incoming tree-structure. This strategy introduces into the field of GP the for1mal memetic approach developed in the frame of GA� in [Radcliffe and Surry 19gl4b]: the aim �till i� to confine the genetic population in the region of local optima with respect to the floating-point values.

First results
The aim of GP is to find rheological modc�ls fitting the set of available experiments. The GP environment itself was written from scratch in c++: Available GP packages (including GPQuick used in section 23.3.5) required too many modifications, and there is no need for efficient tree evaluation, as each tree needs to be interpreted only once per fitness computation (see 23.2.2). The typical time for a run is about 3 hours on a HP 710 workstation. Table 23.1 below summarizes the different parameters of this implementation. All results are averaged over 10 independent runs.
In order to ease further comparisons, the results presented here are concerned with identifying a known rheological mod.el, which represents a hand-crafted model of polyethylene [Kichenin 1992) (Figure 23.1). Parametric optimization was used to adjust the numerical coefficients ( k1 == 790.45, k2 = 150.20, k3 = 41.60, 'fJ = 6248.60, us= 7.25) from mechanical test: s. The fitness cases consider the actual be havior of this model in the same experim• ental conditions, simulated from a regular PDE solver.
The dynamics over the course of evolution are plotted as the best fitness averaged . . over 10 independent runs obtained for a number of fitness computations; the best of all generations is also indicated. The horizontal line gives as reference the average best result obtained by blind random search for the same total number of trials. Typical results are shown in Figure ��3.2, obtained for 200 individuals with a random hill climber of stubbornness 2, without approximations of the fitness. The rate of success is 20%; the best individual so far is shown in Figure 23.2b, with coefficients k1 = 998.892, k2 = 133.085, j� = 39.6647, rJ = 8698.78, us= 19.0409.
According to experts, the difference with the model of Figure 23.1b comes from the absence of creep in all considered experirnents. We are particularly interested in the ilnpact of the random hill-climber (23.2.3) and of the iterated fitness scheme (23.2.2�.3), on the dynamics of evolution.
A first key issue deals with the RHC stubbornness 0. Only small values of 0 were found profitable (see Figure 23.3a): a thorough optimization of coefficients in the first stages of evolution is observed to mislead the search (curve 0 = 3); but a limited optimization ( 0 = 2) leads to beltter results than standard evolution ( 0 = 1). However, the level of improvement decreases when the size of the population increases (e.g. 0 = 1 seems the best choke with 500 individuals; but better overall results were obtained with 200 than witb 500 individuals).
Another important issue is the effect of the gradual refinement of the fitnesses. Gradually increasing the number of experiments taken into account was found to be prejudicial: after several generations have considered one fitness case only, the population gets trapped in local optima. Gradually increasing the size of truncated fitness cases (involving the first K time steps in a loading history) was found simi larly prejudicial: considering only the beJ�nning of the experimental curves results in a population of mostly linear trees and leads to premature convergence. \ , . · · · · · · ' · · · · · · · · · · · · · . . . . �· · · · · · · · · · · · · · · · · · · · · · · · · · · -;; 100 NUJDber of Fitness Computatioos Last, gradually increasing the precision of the discretization of the fitness cases (i.e. increasing the number K of time st. eps) is found to be helpful, especially for small populations. Considering a given precision for too long is prejudicial; but starting with a low precision seems to glUide the population toward promising re gions. Figure 23.3b shows the dynamics obtai ned for the three-phase scheme: • The fitness is first computed with precision 1/3 (only one in three consecutive time steps is considered), until the current best fitness falls below the 75% of the initial best fitness, • The fitness is then computed with pr· ecision 1/2 (only one in two consecutive time-steps is considered), until again the current best fitness falls below 75% of the initial best fitness using this precision, • Last, the fitness is computed with full : precision.
The dynamics of such evolution needs f3ome corrections in order to perform com parisons. First, the difference in computa. tional time is accounted for in the number of fitness computations. Second, the best performance of a population is plotted as the actual fitness (with full precision) of the best individual with respect to the current fitness. This performance is not that of the best individual of the popu lation with respect to the actual fitness, and it occasionally decreases during first and second phases (Figure 23.3b).

Further work
A next step, from a mechanical point of view, is concerned with improving the readability of the result, through syntac: tic simplifications of the produced trees. Such simplifications could intervene as a new evolution operator: this could, as in the case of classical regression, limit the growth of trees as evolution goes on.
The most exciting developments for t; he GP community are connected to the structure of the fitness function. Rather than a sum, it could be viewed as a multi objective function by separately consideling the partial fitnesses derived from the fitness cases. This leads to characterize the set of optimal solutions in the sense of Pareto (solutions not dominated on at least one objective). In the same perspec tive, a new crossover operator could be designed for multi-objective optimization. For instance, the selection-seduction scheme [Ronald 1995] could be turned into a "reciprocal-seduction" scheme, where the! reciprocal seduction of any two individu als depends on their complementarity. We restrict the discussion in this sectio][l to hyperelastic isotropic materials, e.g. rubber. A material is said to be elastic if its deformation under external loads is reversible, and if its equilibrium state do� es not depend on its history.
The central problem in nonlinear thre«�dimensional elasticity [Green and Zerna 1968] consists of finding the equilibrium position of an elastic body that occupies a reference configuration n in the absentee of applied forces. When subjected to external loads, the deformed configuration is characterized by its deformation, a mapping <p : 0 � R3. The equation of equilibrium forms a boundary problem [Ciarlet 1988) in terms of cp, T(x, Vcp) (the first Fiola-Kirchhoff stress tensor of the real positive 3 x 3 matrix V <p), the density of the applied body forces per unit volume, the density of the applied surfac, e forces per unit area on some part of the boundary r of n, and a given deformation on the other part of r.
A material is said to be hyperelastic if 1there exists a potential law W, named the

The numerical model
For a given valid (see section 23.3.6.2) strain energy function W(t 1 , 1-2, t3), the behavior of a structure submitted to boundary conditions and a given load can be numerically computed using the finit· e element method (FEM). The structure is discretized into elements and the syst•em of PDEs is discretized, giving a finite system in terms of the values of the unknowns at some discrete points, i.e. the nodes of the mesh. The result is a disc1rete system of PDEs, that can be solved using iterative numerical methods. We consider the [0, 1]3 cube (all physieal measures are given in standard units), a regular cubic mesh is used. Each element has 27 nodes (9 nodes per face, and a node at center of the element), to allow an interpolation of degree 2. Figure 23.4a illustrates such a 2 x 2 x 2 mesh, and the resulting system has 3 x 125 unknown variables (3 displacements per node). The discrete system is solved by Newton's method, which implies the computation of the derivative matrix of the system and its inversion, for each nonlinear iteration. The complexity of one finite element analysis thus depends on the size of the discrete system (number of elements x number of degrees of freedom per element) and on the number of iterations to solve each nonlinear system (which can only be upper bounded). It does happen that Newton's method fails to converge, and such a situation can only be detected after the maximal number of iterations have been run. The complexity is maximal when no meaningful solution can be derived. The direct problem is: For each loading case, the corresponding response of the structure, in term of displacement (or stfnin or stress if needed} is given by a finite element analysis, provided the energy function is given.

The inverse problem
We are now interested in the inverse problem: Jilrom experimental responses (dis placement, strain or stress) of a structure: under known loads, find the strain energy function of the structure's material . A fundamental issue in function identification lies in the model chosen for the function to identify. H a parametric modE�l is known, the identification then reduces to parameter optimization. However, wh«en no specific model has been identified, a variety of general models can be used. The quasi-parametric models of polynomials, rational fractions or any other nonlinear regression models are a priori discarded. This is because the precision of the result depends on the number of parameters, which has to be fixed: too low, no accur�Lte solution exists in the search space, too high, the complexity of the parametric O]ptimization gets too large. Two other widely used non-parametrk identification tools are neural networks (R umelhart and McClelland 1986] and genetic programming (Koza 1992]. In both models, the complexity of the solution iEI identified together with the parameters. Note that both feed-forward and recurrent neural networks have proven to be robust identification tools, in situations where no other approach succeeded (e.g. [Yao 1993;Fadda and Schoenauer 1995]).
But a characteristic of the numerical modeling problem for mechanical structures described in section 23.3.2 is that the finite element analysis implies the computa tion of numerical values for first-and second-order derivatives of the strain energy function.
• If the chosen model only allows the co1mputation of values of the function itself, numerical derivation is the only way to obtain derivative values. But the instabil ity and inaccuracy of numerical derivation forbids its use to compute second order derivatives.
• Direct computation of derivative valw� is impossible for recurrent neural net works in which output values are given as a fixed point of an iterated process, and barely tractable for feed-forward neural networks.
Before addressing the identification problem, we describe in next section our GP environment, and discuss two important; issues: the adjustment of floating-point terminals and the influence of using derivatives in the fitness computation.

The GP environment
The basis of our GP package is A. Singleton's GPQuick (in C++), tailored to handle floating-point terminals and to coJmpute derivatives. Throughout the end of the chapter, the function set is made of the standard binary operations { +, -, *} and the terminal set consists of the variables for the problem, plus the floating-point random terminal 'R (to be discussed in next subsection).
The evolution scheme is a combination of steady state and generational: copies of the 30% of the population, chosen using a. tournament of size four, are altered using the genetic operators. The replacement of the resulting trees is then performed by repeated tournaments of size two among the initial population.
The genetic operators are crossover, standard mutations (change of function, in sertion of a random subtree or promotion of a subtree), plus Gaussian mutation to adjust the floating-point terminals. Rather high values for the mutation rates were found to be necessary, probably due to the small sizes of population used.

Handling floating-point tE�minals
In past GP research, no great attention was paid to the floating-point terminals. In most cases, as in [Koza 1992), discreti2:ed versions of floating-point terminals are used, and their adjustment to precise values is left to the combination of functions or to random mutations. In [Iba, de Garis and Sato 1994], the difficulty is tackled using a regression method, which is workable only on data fitting problems. We feel that, at least for numerical applications, floating-point terminals must be given a better chance to precisely adjust. The rnodifications to the original GPQuick are: Initialization: The precision of floating· -point terminals is the machine precision. They are chosen randomly in (10-3, 103]. Evolution: Two new mutation operators are allowed: • A mutation by addition of a random Gaussian variable, as first used in evolution strategies ( Schwefel 1981]. The average of the number of terminals of a given tree that undergo this Gaussian mutation is fixed (five in our experiments).
• A random hill-climber resembling that described in section 23.2.3, with the differ ence that this hill-climber is used as an alternative mutation operator rather than being systematically applied after any other operator.
Comparative results: Comparative �:periments were performed on a symbolic regression problem with two variables (: see Table 23.2), involving three different methods to handle floating-point terminalls: the standard evolution by GP crossover and mutation, excluding Gaussian mutations of floating-point terminals, was com pared to the two mutation operators described above.
The best results (Figure 23.5) are obtained using the Gaussian mutation operator. The gradient-like hill-climbing strategy pt erforms quite poorly indeed, especially in the early generations. It does end up a llittle higher than the Gaussian mutation, but these plots do not account for the computational overhead.
The other point is that the GP-style n1ethod performs quite well, only slightly worse than the Gaussian mutation. However, its major drawback is a tremendous increase of the average structural complexity of the trees. As in Section 23.2.3, this greatly increases the computational comJPlexity of one fitness computation. So all further experiments reported in this chapter use Gaussian mutation. The results of Figure 23.5 witness that GP, as well as GAs (G oldberg 1989] can suffer from too much determinism: optirnizing the floating-point terminals of raw tree structures does not help. On the opposite, when partial convergence has be gun, such technique can speed up the evolution, suggesting a mixed technique that could choose between both strategiE!S adaptively along generations as done in (S pears 1995(S pears , Julstrom 1995.

Using derivatives
The most specific aspect of this application, from a GP point of view, is the use of derivatives during the computation of tht e fitness. The computation of the derived trees are straightforward, but the comp1lexity of the derived trees is significantly larger than that of the initial tree.
It is clear from the function set usedl here ( { +, -, *}) that the raw size of a derivate tree is increased by a factor depe� nding on the complexity of the initial tree and on the number of times function * appears. Some trivial simplifications can easily be implemented, (e.g. �(Constant terminal) = 0). But our experiments suggest the increase for the size of the de� rivate tree compared to that of the initial tree ranges from around 1.5 times for sin1ple trees (in the early stages of the runs) to up to more than 3 times for the usually more complex trees appearing in the end of the runs.

Strain en erg y function identification
We now come back to the problem of identifying the strain energy function of some hyperelastic material from mechanical ex:periments (section 23.3.3). 23.3.6.1 Fitness computation The geometry of the structures is shown in Figure 23.4a, together with the very coarse mesh used through all experiments. The structure is fixed below, and the loading cases consist of forces applied on the whole upper surface, with different directions. The experimental displacem · ent fields are, at the moment, results of numerical simulations using a known str.ain energy function. The influence of the computational error thus vanishes, and the best solution is known.
The computation of one fitness case aJmounts to between 12 and 50 (number of Newton's nonlinear iterations) x 8 (number of elements of the mesh) x 27 (number of integration-points per element) compuLtations of the outputs of the 3 first-order and the 6 second-order derivate trees, plus one numerical solution of the 378 by 378 linear system and one comparison with the experimental displacements.
For an average number of nonlinear iterations of 25, with 12 loading cases (see section 23.3.6.3), the time consuming part, of a fitness evaluation lies in the execution of the GP trees (around 35000 per fitness computation). As these trees are derivate trees, their average complexity, compared to that of symbolic regression trees is around 3 for first-order trees and 9 for second-order trees (see section 23.3.5.2). So the 35000 derivate-tree evaluations actually amount to about 250000 standard tree evaluations. But the most important; point is that the maximal computational cost is reached when Newton's method filils to converge: such failure can only be detected after the maximal number of it•erations (50) has been reached. The next section addresses this difficulty by trying · to a priori detect such situations.

Constraints on the strain energy function
The very first attempts proved to be com: plete failures. The reason for that is fairly simple: only around 2% of random trees actually are valid strain energy functions. Using a random function as a strain energy function leads in 98% of the cases to non-convergent iterations (best case) or to unpredictable abortion (worst case) of the finite element analysis.
A closer look at hyperelastic materials theory shows that there are many neces sary conditions energy strain functions rrmst satisfy: • The limit when any of the three variables goes to infinity must be + oo; • The limit when the third variable t3 goes to 0 must be +oo; • When the structure is in its initial stat1e (<p(x) = x, i.e. \l<p = 1 and (t1, t2, t3) = {3, 3, 1)), the following no-loading conditi;on holds: aw +2aw� + aw = o.

0£ 1 8t2 8t3
The a priori restriction of the search space to match these conditions is impossi ble. But the first two limit conditions can be checked a posteriori: the price to pay is that of a few computations of W with at least one variable taking large values. During the initialization of the population, trees failing the test are discarded, and after crossover or mutation, offspring tha. t fail are assigned zero fitness.
Taking the third condition into accoun1t is achieved by mechanical means: Ogden materials [Ogden 1984) are defined by a strain energy function of the form W(t1, t2, £3) = W 9(4�t, t2, t3 ) -a log( t 3 ), for any function W9 and positive a. The behavior of almost all known hyperelastic materials can be approximated using some Ogden law. The idea is to use GP to identify W9, W being defined by the: above Ogden equation. Parameter a is computed a posteriori such that the no-loading condition holds.

The fitness cases
The "experimental" data are built usin, g the simplest form of Ogden materials, Mooney-Rivlin materials [Ciarlet 1988) where the W9 function to identify is linear: Three basic loading cases are used throughout the experiments. In each loading case, forces are applied to the upper boundary of the cube structure. Loading #1 is made of vertical forces. Figure 23.4b plots the displacement fields of the structure submitted to loading case #2 (oblique forces), inducing traction and compression in the structure. Loading case #3 resembles loading case #2, with different direction of forces inducing torsion.

First results
Experimental settings are summarized in Table 23.3. The first results dealt with Mooney-Rivlin materials with a= b = c = 1. Such a problem was far too easy: the exact solution was found every time in less than 5 generations by both GP and blind random search. So the Mooney-Rivlin energy function of the "experimental" data for all following experiments uses values 1.1, 1.2, 1.3 for a, b, c.
The choice of strength for the loading cases raises the next concern: random strain energy functions of Ogden materials, although satisfying the necessary conditions listed in section 23.3.6.2, can model matedals with weird behaviors, ending in stiff, highly nonlinear systems. And the Newton method used to solve the nonlinear system is not robust enough to still be able to converge. This implies these energy functions get null fitness, and thus will likE�ly be eliminated by the selection pressure. But unfortunately, they represent a large majority of the initial populations of any GP run. Therefore, it is essential to be a· ble to assign partial credit (as emphasized in [Kinnear 1994]) to these early solutions.
Our approach for overcoming this difficulty is that any viable strain energy func tion behaves smoothly around the initial position, i.e. for very small loads. But on the other hand, very small loadings give ri se to quasi-linear behaviors, whatever the energy function. Under very small loads:, any function whose linearization around the no-load point (3, 3, 1) resemble the initial Mooney-Rivlin function gets a low error on the displacements. For instance, for a force of intensity 10-3, the solutions found by GP, achieving very low error on the fitness cases, are useless to predict Therefore not significant loads are necessary in order to distinguish between dif ferent energy functions having the same !linearized behavior. Figure 23.6a presents the first successful results in this respect, obtained with forces intensity of 8 10-2 N. The maximal displacement of the structure is around 0.5.
The results are at the moment limited to four runs of GP (Figure 23.6a). They are to be compared to results of Figure 23i.6b, obtained by three runs of a stochastic iterated hill-climber (SIHC). During an �elementary hill-climbing step, an individ ual undergoes mutation, being replaced by any fitter offspring. It ends after a prescribed (20) number of unsuccessful rr.mtations. SIHC repeatedly performs such elementary steps, starting from randomLly chosen individuals (and one "run" is therefore an artificial sequential grouping; of some of these hill-climbing steps). For both GP and SIHC, only the number of finite element analyses actually performed are recorded. Non viable trees are simply discarded.
These curves clearly indicate that GP outperforms the blind local optimization approach in terms of error on the displacements. Moreover , the best three GP runs did find out the linear form of the W9 fu nction, the overall best coefficients being (1.278493549, 1, 1.527635828). We are aware that fu rther experiments are needed to actually demonstrate the full power of this approach.
23.3.6.5 Discussion and Further �To rk The number of fitness evaluations is quit•e limited in the results of section 23.3.6.4, due to the large computational comple�ity of the FEM. As quoted in section 23.3.6.1, this complexity depends not only on the mesh (chosen as coarse as possi ble) and on the structural complexity of the tree at hand, but also on the number of iterations of the nonlinear numerical method, which is greater for large than for small loads. To address both the problem of complexity related to load inten sity and the issue of partial credit to initial trees, we propose to use an iterative scheme, modifying the fitness cases during evolution, in a similar way as for the one-dimensional rheological model identification (section 23.2.2.3). On-going exper iments use fo rce strengths gradually increasing from 10-3 N for the first generation to 10-1 N.
Another on-going attempt to minimize the computational cost is to use symbolic computation to simplify the structural complexity of the derivative trees. Prelimi nary results in that direction suggest that a benefit can be expected. Moreover, a post-mortem optimization of the floating points terminals of the simplified version of the overall best tree would certainly irnprove the results.

Summary
This chapter demonstrates the basic feasibility in structural mechanics of a GP based approach for the identification of both one-dimensional rheological models and three-dimensional hyperelastic strai1n energy functions. Moreover, the use of GP for these problems is fully justified: The structure of the solutions is unknown, and no other non parametric technique can apply.
From the point of view of mechanical :modeling, this work opens many perspec tives regarding the identification of materials from their experimental behavior, as they are insofar restricted to heuristi·cs and specific studies. Another exciting perspective would be to ease the design of new materials with stipulated behaviors.
This study represents one of the few applications of GP to Applied Science (Oakley 1994]. Such cross-fertilization of Evolutionary Computation and Applied Science was made possible through a thorough collaboration between experts of both fields. The background knowledge in mechanics was essential in building the wrapper for the rheological model identification (section 23.2.2.2), and restricting the search space, in order for GP to explore only viable energy functions (section 23.3.6.2). The most limiting factor of this approach is the computational complexity. We propose to gradually refine the fitness (sections 23.2.4 and 23.3.6.5). Preliminary results show that such an iterative scheme can significantly reduce the time-to solution for the rheological model identification. However, more sophisticated switching strategies remain to be investi.gated in order to avoid population to get trapped in the local optima of intermediate fitnesses. Another critical issue deals with the precise adjustment of the floating-point values in numerical GP applica tions (sections 23.2.3 and 23.3.5.1).
As a conclusion, we feel that the unique ability for GP to explicitly deal with derivatives is one of its most appealing aspects for Engineering and Design problems. The presented work suggests that actual breakthroughs can result from a tight interaction between GPers and open-minLded researchers and engineers in Applied Science.