<figure><img src="image_1.png" /></figure>

HAL
open science

# A multiparametric strategy for the two step optimization of structural assemblies

Bruno Soulier, Pierre-Alain Boucard

► To cite this version:

Bruno Soulier, Pierre-Alain Boucard. A multiparametric strategy for the two step optimization of structural assemblies. Structural and Multidisciplinary Optimization, 2013, 47 (4), pp.539-553. 10.1007/s00158-012-0854-9 . hal-01376463

HAL Id: hal-01376463

https://hal.science/hal-01376463

Submitted on 20 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multiparametric strategy for the two step optimization of structural assemblies

B. Soulier · P.-A. Boucard

Received: date / Accepted: date

**Abstract** Generally speaking, the objective and constraint functions of a structural optimization problem are implicit with respect to the design variables; their evaluation requires finite element analyses which constitute the most expensive steps of the optimization algorithm. The work presented in this paper concerns the implementation of a two step optimization strategy which consists in optimizing first an empirical model (metamodel), then the full model. In the framework of multilevel model optimization, the computation costs are related, on the one hand, to the construction of global approximations and, on the other hand, to the optimization of the full model. Thus, many numerical simulations are required in order to perform a multilevel optimization. In this context, the objective of associating a multiparametric strategy based on the nonincremental LATIN method with the two step optimization process is to reduce these computation costs. The performance gains thus achieved will be illustrated through the optimization of structural assemblies involving contact with friction. The results obtained will show that the savings associated with the multiparametric procedure can reach a factor of 30.

**Keywords** two step optimization · Metamodel · LATIN method · Multiparametric strategy

## 1 Introduction

Optimum and robust product design relies massively on numerical simulations in which the phenomena being simulated are becoming increasingly detailed. In order to assist the design engineer in optimization tasks which, today, are becoming more and more complex due to their multilevel and multidisciplinary nature, numerical methods are being developed for the exploration of the design space through global approximations.

Many works relative to the global approximation concept have focused on the reduction of computation costs, which constitute a real stumbling block for optimization problems.

Barthelemy and Haftka (1993) identified various approximation categories, including *global* approximations defined over all or part of the design space and *local* approximations defined only in the vicinity of a point. They distinguished the *approximation of a function*, leading to the substitution of an explicit expression of the parameters for the objective function, which is usually implicit, from the *approximation of a problem*, which consists in replacing an optimization problem by a series of approximate problems which can be solved much more easily. Local approximation functions are usually based on Taylor series. First-order approximations, which consist in evaluating the functions and their first derivatives at all points, are used in many works, particularly by Braibant and Fleury (1985); Fleury and Braibant (1986) who proposed a dual approach to the resolution of convex problems. A comparison of first-order approximations with second-order approximations (which include the quadratic term of the Taylor series) was presented by Haftka (1988). Pritchard and Adelman (1990, 1991) proposed approximations based on differential equations.

<footnote>B. Soulier · P.A. Boucard
LMT-Cachan (ENS Cachan/CNRS/UPMC/PRES Univers-
Sud Paris), 61 avenue du Président Wilson, 94235 Cachan
Cedex, France
E-mail: soulier@lmt.ens-cachan.fr</footnote>
Then, the solution of the optimization problem can be carried out using classical algorithms, such as the simplex algorithm (Dantzig et al, 1955) for linear problems or the SLP (Sequential Linear Programming) algorithm for the resolution of nonlinear problems using Taylor series to linearize the objective and constraint functions (Zienkiewicz and Campbell, 1973; Pedersen, 1981).

The most efficient optimization methods rely on multilevel optimization concepts, which can be divided into three categories: parallel model optimization based on domain decomposition methods (El-Sayed and Hsiung, 1991; Umesha et al, 2005), multilevel parameter optimization, which consists in replacing an optimization problem by several subproblems (each with a reduced set of parameters) (Liu et al, 2004; Kravanja et al, 2003; Chen and Yang, 2005; Bendsoe, 1995), and multilevel model optimization, which introduces several modeling levels (Robinson and Keane, 1999; Keane and Petruzzeli, 2000; Engels et al, 2004). The proposed two-step hybrid optimization approach can be classified in the multilevel model optimization. A major difficult in direct optimization with the FE model is to find the global optimum. It is well known that the direct global optimization with the FE model increases the computation time significantly and a local optimization led us to local optima. To increase the chance of finding the global optimum the two step optimization method used here is based on successive global optimization of a metamodel, representing the first modeling level, and local optimization on the full mechanical model, which constitutes the second level. Cutting the computation costs in the whole optimization process requires not only the reduction of the computation costs inherent to the construction of the metamodel, which involves simulations of the full mechanical model, but also their reduction in the full model's optimization stage. In order to reduce these costs, we associate what we call a "multiparametric" strategy based on the LATIN method with a metamodel constructed by diffuse interpolation. The method is illustrated through examples of optimization of assemblies in contact with friction (Hilding et al, 2001; Li et al, 2005).

The first part of the paper puts multilevel modeling in the context of a full multidisciplinary optimization scheme in which, starting from a series of exact analyses, an explicit global approximation is developed in the form of a response surface. Besides the construction of continuous and differentiable global approximations over the whole design space, the metamodel should enable updating the approximation iteratively throughout the optimization procedure. Among many

kind of metamodels with these properties like polynomial regression (Box and Wilson, 1951), neural networks (Haykin, 1994), radial basis functions (Hardy, 1971), kriging (Cressie, 1990), this paper focuses on cumulative interpolation (Soulier et al, 2003).

The second part reviews the fundamentals of the multiparametric method based on the LATIN approach and emphasizes its advantages in the optimization context. In particular, this multiparametric method takes advantage of the capability, with the LATIN method, to reuse already available calculation results in order to improve computation times.

In the third part, these developments are applied to the optimization of the connection parameters between the different parts of an assembly.

## 2 The two step optimization strategy

Multilevel optimization is generally used when the problem to be solved is relatively complex and, therefore, divided into several steps (i.e. on several levels). Consequently, from one study to the next, the meaning of the term *multilevel optimization* can be quite different. Nevertheless, to summarize what can be found in the literature, one can point out the following two main definitions. In what is classically called multilevel parameter optimization (Liu et al, 2004; Kravanja et al, 2003; Chen and Yang, 2005; Bendsoe, 1995), one considers an initial model which remains unchanged during the optimization process, but is treated with varying levels of accuracy. In what we will now call multilevel model optimization, which is at the heart of the approach we are proposing, the initial model is allowed to evolve during the optimization. Two main types of enrichment can be found in the literature.

*Hierarchical* descriptions (Robinson and Keane, 1999; Keane and Petruzzeli, 2000) are based on different models defined successively by increasing the refinement of the model. Generally, the required information is carried only from Level *i* to Level *i* + 1, which enables one to forget the more primitive models (e.g. a more or less empirical analytical functional, a beam model, a linear 3D model, then a nonlinear 3D model) permanently. For example, this is typically the underlying idea in a global-local approach. Conversely, in *nested* descriptions (Engels et al, 2004), the basic model is always preserved and is actually enriched from one level to the next so the information can be carried back and forth among the different levels. This leads to a multiscale description of the problem which can take into account the various modeling levels (or scales) simultaneously and can be enriched during the optimization process.
## 2.1 The two step model optimization strategy

The optimization scheme we chose, shown in Figure 1, is based on the successive optimizations of a metamodel and of the full mechanical model. First of all an initial data set is selected with a random strategy or a design of experiment (DOE) strategy (Montgomery, 1997). Starting from this initial data set a metamodel is built and its parameters are estimated with cross-validation (Currin et al, 1991) and the R² criterion is used to assessing the quality of the metamodel. Then in the first step, a global optimization of the metamodel using a genetic algorithm (Goldberg, 1989) enables one to locate areas corresponding to extrema of the objective function. Different algorithms like genetic algorithm (GA) (Goldberg, 1989), particle swarm optimization (PSO) (Kennedy and Eberhart, 1995) or simulated annealing (SA) (Kirkpatrick et al, 1983) are often used in global optimization. GA is used in this study mainly because of its ease of implementation and its ability to solve optimization of nonlinear problems. Finally, from the optimum found at the previous step, a local optimization phase is carried out using the full mechanical model. In this second step optimization and simulation are decoupled and exchange information through data files. During the n iterations of the optimization procedure, n₂ simulations are performed in order to evaluate objective function and its derivatives. The local optimization in this second step is performed with a sequential quadratic programming (Han, 1977) algorithm (the fmincon function in matlab). This two step optimization strategy requires n₀ numerical simulations for the construction of the metamodel plus n₂ additional simulations during the second phase. These simulation phases are usually very expensive in terms of CPU time and constitute one of the main stumbling blocks of the method. The numerical simulation strategy we intend

to use is perfectly suitable for the treatment of these very costly phases. Indeed, the multiparametric strategy we are about to present enables the computation costs to be reduced both for the calculations needed to build the metamodel and for the optimization of the full model. Each calculation corresponds to a set of values of the variable parameters are introduced into the mechanical analysis. As an alternative to carrying out full computation for each set of parameters, multiparametric strategy is based on the capability of the LATIN method (Ladevèze, 1999) to re-use the solution to a given problem (for one set of parameters) in order to solve similar problems (for the other sets of parameters) (Boucard and Champaney, 2003). In this manner, except for the first computation, an approximation of the solution to the new set of parameters with a strong mechanical content is immediately available from the start in order to decrease the CPU time.

For the implementation of the two step model optimization strategy, we used the cumulative interpolation metamodel presented in the following section.

## 2.2 Construction of the metamodel by cumulative approximation

The metamodel presented in this paper was inspired by Smooth Particle Hydrodynamics (SPH) techniques (Gingold and Monaghan, 1977; Nayrolles et al, 1992), which were introduced in the 1980's in the field of astrophysics. The principle of these methods consists in approximating a function y using a weighted sum $\tilde{y}$ of exponential-type or spline-type shape functions denoted $\phi_i$. Rasmussen (1998) proposed a model called Cumulative Approximation (1) defined by a weighted sum:

$$ \tilde{y}_R(x) = \frac{\sum_{i} \phi_i(x) y_i}{\sum_{i} \phi_i(x)} \qquad (1) $$

where $y_i$ is the magnitude of the measured response y at point $x_i$ and $\phi_i(x)$ is the shape function associated with the observation points.

This approximation is not exact at the simulation points. Here, we use an interpolation model based on a cumulative approach. This cumulative interpolation model (Soulier et al, 2003), which can be called diffuse interpolation, is defined by:

$$ \tilde{y}(x) = \frac{\sum_{i} \phi_i(x) a_i}{\sum_{i} \phi_i(x)} \quad \text{with} \quad \phi_i(x) = e^{-k \|x - x_i\|^2} \qquad (2) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1 The two step optimization scheme</figcaption></figure>
where the $\phi_i$ are Gaussian-type shape functions. The coefficients $a_i$ are obtained through a collocation method and are solutions of the linear system:

$$ \left( \begin{array}{c} \frac{\phi_1(x_1)}{\sum_i \phi_i(x_1)} \\ \vdots \\ \frac{\phi_1(x_n)}{\sum_i \phi_i(x_n)} \end{array} \right) \left( \begin{array}{c} a_1 \\ \vdots \\ a_n \end{array} \right) = \left( \begin{array}{c} y(x_1) \\ \vdots \\ y(x_n) \end{array} \right) \quad (3) $$

The vector $(x_1, ..., x_n)$ represents the set of the simulation points at which function $y$ is evaluated, and $\phi_i(x)$ is the value of the shape function at point $x$.

The proposed metamodel generates continuous and differentiable response surfaces over the design space, which enables one to calculate their gradients explicitly. Moreover, the flexibility of this metamodel makes the updating of the interpolation surfaces during the optimization process easier, which enables one to improve the quality of the approximations locally.

Let us illustrate the construction of global approximation models using Rosenbrock's function. This analytical function, which is highly nonlinear and presents a very elongated, arc-shaped trough, is defined by:

$$ f(x_1, x_2) = 100(x_1 - x_2^2)^2 + (1 - x_1)^2 \quad (4) $$

The cumulative interpolation model enables one to improve the quality of the model incrementally by updating the model during the addition of new simulations. Figure 2 shows the model differentiated for an orthogonal factorial design of experiment (Montgomery, 1997), (Box and Hunter, 1957) with 3 levels per factor. A three-level design is a design of experiment (DOE) whose each factors are considered at 3 levels. These are (usually) referred to as low, intermediate and high levels. A factorial design of experiments is built by taking on all possible combinations of these levels across all such factors. One can verify that the approximation provides a good representation of the actual shape of Rosenbrock's function thanks to the additional data points introduced in the zone with the highest nonlinearities.

Some one-dimensional correlation functions classically used to build statistical models like kriging (Cressie, 1990) in the field of the analysis of computer experiments (Currin et al, 1991) (Sacks et al, 1989) are the following:

- linear correlation function

$$ \left\{ \begin{array}{ll} r(d) = 1 - (1 - \rho) |d| & \text{if } |d| < \frac{1}{1-\rho} \\ r(d) = 0 & \text{otherwise} \end{array} \right. \quad (5) $$

<figure><img src="image_4.png" /><figcaption>Fig. 2 Comparison metamodel vs. actual surface</figcaption></figure>

- exponential correlation function

$$ r(d) = \rho^{|d|} \quad (6) $$

- Gaussian correlation function

$$ r(d) = \rho^{d^2} \quad (7) $$

where $d$ represents the one-dimensional distance and $\rho$ the correlation between two points.

These three types of correlation functions lead to different interpolation models. In the case of linear correlation functions, the resulting models are not differentiable at the measurement points due to $|d|$. We choose Gaussian correlation functions, which produce differentiable surfaces over the whole design space. Extension to more dimensions is obtained with the product correlation rule which consists to multiply one-dimensional correlation function. The shape functions are of the Gaussian type:

$$ \phi_i(x) = e^{-k \|x - x_i\|^2} \quad (8) $$

where $\|x - x_i\|$ denotes the Euclidian distance between points $x$ and $x_i$, and $k$ is a strictly positive standardization parameter.

Parameter $k$ controls the shape of the surrogate surface by modifying the size of the support of the shape functions (Figure 3). Small values of $k$ produce very smooth surfaces. Large values of $k$ produce strongly non linear surfaces.

<figure><img src="image_5.png" /><figcaption>Fig. 3 Influence of parameter $k$: left, $k = 2$; right, $k = 10$</figcaption></figure>
### 3 The multiparametric strategy using the LATIN method

Since the optimization process generates numerous calculations, it is crucial to introduce a suitable reanalysis strategy leading to reduced computation times when the parameters change. In order to do that, we use the LATIN iterative algorithm developed by Ladevèze (1999) along with a multiparametric strategy (Boucard and Champaney, 2003).

#### 3.1 Review of the LATIN method for contact problems with friction

The basic principle of the LATIN method consists in isolating difficulties by solving two groups of equations successively: the local (possibly nonlinear) equations and the linear (possibly global) equations. An assembly is a set of substructures (parts) which communicate through interfaces (Figure 4). Each interface is a mechanical entity in its own right with its unknowns and constitutive relations. In order to simplify the presentation, let us consider only two substructures $\Omega_E$ and $\Omega_{E'}$ connected by an interface $\Gamma^{EE'}$. The interface variables are two force fields ($\mathbf{f}^E, \mathbf{f}^{E'}$) and two dual velocity fields ($\dot{\mathbf{w}}^E, \dot{\mathbf{w}}^{E'}$) (Figure 4). By convention, ($\mathbf{f}^E, \mathbf{f}^{E'}$) are the actions of the interface on the substructures and ($\dot{\mathbf{w}}^E, \dot{\mathbf{w}}^{E'}$) are the velocities of the substructures viewed from the interface.

<figure><img src="image_6.png" /><figcaption>Fig. 4 Decomposition of an assembly; interface variables</figcaption></figure>

The displacement field at any point $M$ of $\Omega_E$ and at any time $t$ in $[0, T]$ is $\mathbf{u}^E(M, t)$ and $\mathcal{U}_{ad}^{[0,T]}$ is the associated space. Then, the problem to be solved in each substructure is to find the evolutions of the displacement field $\mathbf{u}^E(M, t)$ and stress field $\sigma^E(M, t)$ which verify:

- kinematic admissibility:

$$ \varepsilon = \varepsilon(\mathbf{u}^E); \quad \mathbf{u}^E \in \mathcal{U}_{ad}^{[0,T]} \qquad (9) $$

$$ \mathbf{u}^E(M, t)|_{\partial\Omega_E} = \mathbf{w}^E(M, t) \qquad (10) $$

- static admissibility: one consider $\dot{\mathbf{u}}^*$ a virtual velocity field and $\mathcal{U}$ the set of finite-energy velocity fields on $\Omega_E$ which vanish on $\partial\Omega_E$, $\forall \dot{\mathbf{u}}^* \in \mathcal{U}$

$$ \int_{\Omega_E} \mathrm{Tr}(\sigma^E \varepsilon(\dot{\mathbf{u}}^*)) d\Omega_E - \int_{\Omega_E} \mathbf{f}_d \cdot \dot{\mathbf{u}}^* d\Omega + \int_{\partial\Omega_E} \mathbf{f}^E \cdot \dot{\mathbf{u}}^* dS = 0 \qquad (11) $$

- elastic behavior: $\forall M \in \Omega_E, \forall t \in [0, T]$,

$$ \sigma^E(M, t) = \mathcal{K}_e \varepsilon(\mathbf{u}^E(M, t)) \qquad (12) $$

($\mathcal{K}_e$ being Hooke's operator)

The problem to be solved at the interface is to find the evolutions of the force fields ($\mathbf{f}^E(M, t), \mathbf{f}^{E'}(M, t)$) and velocity fields ($\dot{\mathbf{w}}^E(M, t), \dot{\mathbf{w}}^{E'}(M, t)$) which verify:

- equilibrium: $\forall M \in \Gamma^{EE'}$ and $\forall t \in [0, T]$,

$$ \mathbf{f}^E(M, t) + \mathbf{f}^{E'}(M, t) = 0 \qquad (13) $$

- the behavior described as a nonlinear evolution law $\mathcal{R}$ between the forces and the rate $\dot{\mathbf{w}}^{EE'} = \dot{\mathbf{w}}^{E'} - \dot{\mathbf{w}}^E$ of jump in velocity across the interface (typically: contact, friction, etc.) (Blanzé et al, 1995): $\forall M \in \Gamma^{EE'}$ and $\forall t \in [0, T]$,

$$ \mathbf{f}^E(M, t) = \mathcal{R}(\dot{\mathbf{w}}^{EE'}(M, t), t \in [0, T]) \qquad (14) $$

The solution $s$ is described a priori as a set of time-dependent fields relative to both the interface and the substructures. Here, the substructures have linear elastic behavior and the interior solution (i.e. the displacement $\mathbf{u}^E(M, t)$ and the stress $\sigma^E(M, t)$) can be easily calculated from the boundary quantities ($\dot{\mathbf{w}}^E(M, t)$ and $\mathbf{f}^E(M, t)$). The solution $s$ can be represented using only the force and velocity fields on both sides of the interface.

$$ s = \sum_E s^E, \quad s^E = \{\dot{\mathbf{w}}^E(M, t), \mathbf{f}^E(M, t)\}, $$

$$ \forall t \in [0, T] \qquad (15) $$

The equations are divided into two groups such that the substructures are considered to be elastic and all the nonlinearities are concentrated at the interface:
- The set $\mathcal{A}_d$ of the solutions $s^E$ which verify the linear equations relative to the substructures (Equations (10) to (12));

- The set $\Gamma$ of the solutions $s^E$ which verify the local (possibly nonlinear) equations relative to the interface (Equations (13) and (14)).

<figure><img src="image_7.png" /><figcaption>Fig. 5 Diagram of an iteration of the LATIN method</figcaption></figure>

Then, the solution of the problem is determined iteratively by seeking successive approximations $s$ which verify the two groups of equations alternatively, using search directions $E^+$ and $E^-$ (Figure 5). Thus, the two steps of the iterative algorithm are:

- local step: given $s_n \in \mathcal{A}_d$, find $\hat{s}$ such that:

$$ \hat{s} \in \Gamma \quad (\text{interfaces}) \quad (16) $$

$$ \hat{s} - s_n \in E^+ \quad (\text{search direction}) \quad (17) $$

- global step: given $\hat{s} \in \Gamma$, find $s_{n+1}$ such that:

$$ s_{n+1} \in \mathcal{A}_d \quad (\text{substructures}) \quad (18) $$

$$ s_{n+1} - \hat{s} \in E^- \quad (\text{search direction}) \quad (19) $$

Here, we use conjugate search directions which depend on a single scalar parameter $k_0$:

$$ \hat{s} - s_n \in E^+ \quad \equiv \quad (\hat{\mathbf{f}}^E - \mathbf{f}_n^E) = k_0(\hat{\mathbf{w}}^E - \dot{\mathbf{w}}_n^E) \quad (20) $$

$$ s_{n+1} - \hat{s} \in E^- \quad \equiv \quad (\mathbf{f}_{n+1}^E - \hat{\mathbf{f}}^E) = -k_0(\dot{\mathbf{w}}_{n+1}^E - \hat{\mathbf{w}}^E) \quad (21) $$

The solution of the problem is independent of the value of parameter $k_0$, which affects only the convergence rate of the algorithm. In the case of quasi-static calculations such as those we are concerned with, $k_0$ is given by:

$$ k_0 = \frac{E T}{L_c} \quad (22) $$

where $E$ is Young's modulus, $[0, T]$ the time interval being considered and $L_c$ the largest dimension of the structure.

## 3.2 The multiparametric strategy

The multiparametric approach (Boucard and Champaney, 2003, 2004) is based on the fact that the LATIN algorithm can be initialized using any solution which verifies the admissibility conditions. In the case of a parametric study, for a given set of parameters, one reinitializes the LATIN loop with the converged solution corresponding to another set of parameters (Figure 6). If a parameter changes only slightly, so does the global solution of the FE problem. Thus, thanks to the multiparametric strategy, convergence is achieved more rapidly and after a reduced number of iterations. Typically, the parameters involved here are prescribed

<figure><img src="image_8.png" /><figcaption>Fig. 6 Reinitialization of the LATIN method</figcaption></figure>

loads/displacements, gaps, friction parameters, preloads or tightening loads, i.e. all the parameters which describe the nonlinear behavior of the interface. The essential point which leads to an efficient strategy is that in the particular case of elastic substructures the interfaces play a fundamental role. Indeed, they enable one to initialize the calculation of the problem associated with a new set of parameters without having to save all the data for the substructures. Moreover, the information concerning the interfaces obtained from the reference calculation enables the method to be initialized using a first approximation with a strong mechanical content which is suitable for the problem being considered. Thus, if the solution of the reference problem is not too different from the solution of the new target problem, one can expect to achieve the latter at low cost. For a calculation associated with a new set of parameters, several other results (associated with other sets of parameters) may be available. Then, one can use any of these results for the initialization. In our study, we chose to initialize each new calculation with the converged solution obtained for the closest set of parameters, in the sense of the Euclidian distance in the parametric space.
## 4 Parametric study of contact with friction in an assembly

In this final part, the capabilities of the multiparametric strategy in the context of two step model optimization are illustrated using a technological example of a force-fit assembly.

The objective of our study was to optimize the performance of shaft-pinion connections achieved by cone clamping as a function of the binding parameters. A shrink disk (Figure 7) is a technological component which consists of a biconical inner ring fitted to the pinion and two external conical pressure flanges, one of which is threaded. A clamping load is applied between the external flanges through a series of screws distributed along the circumference. The tightening of the screws presses the conical surfaces against each other and generates radial forces which create the adhesion binding which is necessary to transmit the torque $M_a$ and axial load $F_a$ between the shaft and the pinion.

<figure><img src="image_9.png" /><figcaption>Fig. 7 A shrink disk</figcaption></figure>

Here, the optimization concerns the binding parameters, i.e. the gaps, tightening forces and friction coefficients between the different elements of the assembly. The objective functions studied were the transmissible torque $M_a$ and axial load $F_a$ between the pinion (or the hub) and the shaft. The strategy known as two step optimization consists in carrying out first a global optimization using the metamodel, then a local optimization based on that first optimum using the full FE model and a multiparametric strategy (see Figure 1).

The global optimization based on the metamodel was performed at low cost using a genetic algorithm (Goldberg, 1989) (the ga function in MATLAB) which generated $n_1$ evaluations of the objective function. This first step enabled us to localize the trough in which the global optimum of the objective function was to be sought. The metamodel was built in a preliminary phase using $n_0$ evaluations of the full mechanical model. The settings of the GA algorithm were specified as follows: a population size equal to 20, a crossover rate

of 0.8, a mutation rate of 0.2 and a number of generations equal to 100. The local optimization using the full mechanical model was performed with a sequential quadratic programming (Han, 1977) algorithm (the fmincon function in MATLAB) which generated $n_2$ evaluations of the objective function.

The total optimization time $t_{multilev+multipar}$ associated with this two step / multiparametric strategy (23) is the sum of:

- $t_0$, the time for generating the metamodel using the multiparametric strategy,

- $t_1$, the time for optimizing the metamodel,

- and $t_2$, the time for optimizing the full FE model using the multiparametric LATIN technique.

$$t_{multilev+multipar} = t_0 + t_1 + t_2 \quad (23)$$

We seek to compare this two step optimization method based on a multiparametric calculation strategy with the same optimization method without the multiparametric aspect. In order to compare the different calculation strategies, it is necessary to define the reference time associated with a calculation. If $t_{ref}$ denotes, for each study, the duration of the first calculation (which, because of the nonlinearities, is obviously a function of the associated set of parameters), one can use as the reference time the average of the different $t_{ref}$ obtained. Let $\overline{t_{ref}}$ designate this average value.

Thus, the total optimization time $t_{multilev}$ associated with the two step optimization strategy with no multiparametric strategy (24) is defined as the sum of:

- the time for generating the metamodel with no multiparametric strategy, which is considered to be the product of $n_0$, the number of sets of parameters used to build the metamodel, by $\overline{t_{ref}}$, the reference CPU time,

- $t_1$, the time for optimizing the metamodel,

- and the time for optimizing the full FE model with no multiparametric strategy, which is considered to be the product of $n_2$, the number of evaluations of the objective function using the local optimization algorithm, by $\overline{t_{ref}}$, the reference CPU time.

$$t_{multilev} = (n_0 + n_2)\overline{t_{ref}} + t_1 \quad (24)$$

The gain $G$ is defined as the ratio of the total optimization time associated with this two step / multiparametric strategy $t_{multilev+multipar}$ by the total optimization time $t_{multilev}$ associated with the two step optimization strategy without multiparametric strategy:

$$G = \frac{t_{multilev+multipar}}{t_{multilev}} = \frac{t_0 + t_1 + t_2}{(n_0 + n_2)\overline{t_{ref}} + t_1} \quad (25)$$
## 4.1 Optimization of the transmissible torque using an axisymmetric model

We studied the evolution of the transmissible torque $M_a$ at the interface between the shaft and the pinion as a function of five parameters (the friction coefficients and the gaps between the different parts in contact) using the axisymmetric model of Figure 8. The mesh is composed of 3,549 axisymmetric linear four-node elements (7,600 degrees of freedom). The calculations were carried out on a PC (Proc. Intel Xeon W5650(tm) 2.6 Ghz, RAM 12 GB) The constraints in terms of bounds

<figure><img src="image_10.png" /><figcaption>Fig. 8 The FE model</figcaption></figure>

on the parameters are given in Table 1. Different draws

<table><caption>Table 1 Bounds on the design parameters</caption><thead><tr><th>parameters</th><th>min</th><th>max</th></tr></thead><tbody><tr><td>μ<sub>shaft/gear</sub>, μ<sub>gear/inner</sub>, μ<sub>inner/outer</sub></td><td>0.05</td><td>0.5</td></tr><tr><td>j<sub>shaft/gear</sub> (μm)</td><td>0</td><td>49</td></tr><tr><td>j<sub>gear/inner-ring</sub> (μm)</td><td>36</td><td>106</td></tr></tbody></table>

of the sets of parameters were used in the metamodel construction preliminary phase:

- an orthogonal Design Of Experiments (DOE) combining 5 factors, either on two levels (8, 16 and 32 combinations),

- a random drawing (LHS) (McKay et al, 1979) of between 5 and 80 combinations in steps of 5 combinations.

Let $n_0$ denote the number of combinations in each draw. For each draw, a first calculation is carried out with the first setting of parameters. Then, the new $n_0-1$ computations are initialized by the solution to the previous one. If the parameters change slowly, the two solutions are close and only a few iterations are needed to reach convergence in the new calculation. If the parameters have very different values, more iterations are necessary to converge, but it is proved in (Ladevèze, 1999), that the LATIN strategy will always converge to the solution. Figure 9 shown that for the three orthogonal arrays sample the CPU time decreases with the iterations of the multiparametric strategy. Table 2

<figure><img src="image_11.png" /><figcaption>Fig. 9 CPU time with the multiparametric strategy in metamodel generation step</figcaption></figure>

summarizes, for each draw —of $n_0$ combinations—the number of iterations ($i_0$) of the LATIN method for the calculation of the $n_0$ combinations and the computation time per set of parameters. The time $t_{ref}$ of the first calculation and the minimum and average CPU times are listed separately for each set of parameters.
<table><caption>**Table 2** Construction of the metamodel</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th rowspan="2">i<sub>0</sub></th><th colspan="3">CPU time per parameter set (s)</th><th rowspan="2">total t<sub>0</sub>(s)</th></tr><tr><th>1st calc. t<sub>ref</sub>(s)</th><th>min. (s)</th><th>mean (s)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-style: italic;">orthogonal array</td></tr><tr><td>8</td><td>3,076</td><td>63</td><td>14</td><td>34</td><td>275</td></tr><tr><td>16</td><td>5,576</td><td>64</td><td>8</td><td>31</td><td>500</td></tr><tr><td>32</td><td>4,496</td><td>53</td><td>0.5</td><td>13</td><td>405</td></tr><tr><td colspan="6" style="text-align: center; font-style: italic;">Latin hypercub sampling</td></tr><tr><td>5</td><td>2,046</td><td>42</td><td>32</td><td>37</td><td>187</td></tr><tr><td>10</td><td>3,666</td><td>55</td><td>17</td><td>34</td><td>336</td></tr><tr><td>15</td><td>4,706</td><td>49</td><td>19</td><td>29</td><td>429</td></tr><tr><td>20</td><td>6,541</td><td>50</td><td>14</td><td>30</td><td>598</td></tr><tr><td>25</td><td>7,641</td><td>50</td><td>13</td><td>28</td><td>697</td></tr><tr><td>30</td><td>8,451</td><td>48</td><td>7</td><td>26</td><td>776</td></tr><tr><td>35</td><td>10,656</td><td>53</td><td>13</td><td>28</td><td>972</td></tr><tr><td>40</td><td>11,346</td><td>45</td><td>9</td><td>26</td><td>1,059</td></tr><tr><td>45</td><td>12,436</td><td>54</td><td>9</td><td>25</td><td>1,140</td></tr><tr><td>50</td><td>13,206</td><td>47</td><td>12</td><td>24</td><td>1,205</td></tr><tr><td>55</td><td>15,296</td><td>45</td><td>15</td><td>25</td><td>1,402</td></tr><tr><td>60</td><td>15,491</td><td>52</td><td>8</td><td>24</td><td>1,418</td></tr><tr><td>65</td><td>18,171</td><td>61</td><td>7</td><td>26</td><td>1,682</td></tr><tr><td>70</td><td>18,021</td><td>53</td><td>9</td><td>24</td><td>1,650</td></tr><tr><td>75</td><td>21,061</td><td>46</td><td>9</td><td>26</td><td>1,926</td></tr><tr><td>80</td><td>20,631</td><td>46</td><td>11</td><td>24</td><td>1,895</td></tr></tbody></table>

The optimization problem that we consider can be formulated as:

$$
\left\{
\begin{array}{l}
\arg\max_{(\mu_{./}, j_{./})} M_a(\mu_{./}, j_{./}) \\
\text{such that} \\
\quad 0.05 \le \mu_{shaft/gear} \le 0.5 \\
\quad 0.05 \le \mu_{gear/inner} \le 0.5 \\
\quad 0.05 \le \mu_{inner/outer} \le 0.5 \\
\quad 0\mu m \le j_{shaft/gear} \le 49\mu m \\
\quad 36\mu m \le j_{gear/inner-ring} \le 106\mu m \\
\quad \mu_{inner/outer} - \mu_{shaft/gear} \le 0 \\
\quad \mu_{inner/outer} - \mu_{gear/inner} \le 0
\end{array}
\right.
\qquad (26)
$$

The optimum values of the transmissible torque at the interface given by the two optimization steps are shown in Figure 10 for the only bounds constraints problem and in Figure 11 for the problem with linear inequality constraints. The optimum issued from the first step is an approximation given by the metamodel (Step 1), which can be greater than the optimum obtained in the second step using the complete mechanical model (Step 2). At the end of this study, the optimum found was 21,706 N.m, associated with the following optimum parameters: $\mu_{shaft/gear} = 0.5$, $\mu_{gear/inner} = \mu_{inner/outer} = 0.05$, $j_{shaft/gear} = 0\mu m$, $j_{gear/inner-ring} = 36\mu m$. The CPU times and the numbers of iterations which were necessary to carry out these two optimization steps are summarized in Table 3 and 4.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Comparison of the optima obtained from Step 1 and Step 2 with bound constraints</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11** Comparison of the optima obtained from Step 1 and Step 2 with bound and linear constraints</figcaption></figure>

Let us compare this two step optimization method based on a multiparametric calculation strategy with the same optimization method and no multiparametric aspect (Table (3) (4)). In order to do that, we consider that the reference time for a FE calculation is considered to be equal to $\bar{t}_{ref} = 51s$.

Gain G, defined as the ratio of $t_{multilev+multipar}$ (the computation time for our two step strategy coupled with the multiparametric aspect) to $t_{multilev}$ (the computation time for the two step calculation strategy with no multiparametric aspect), shows that in a two step strategy the addition of the multiparametric strategy reduced the CPU time by a factor of between approximately 4.4 and 14.2.

One can also observe in Figures 10 and 11 that approximately half of cases reached the same value of the maximum torque, i.e. 21,706 N.m. Nevertheless, there were cases where this value was not found. These are probably cases where the optimization using the meta-
<table><caption>**Table 3** Global optimization with bound constraints using the metamodel and local optimization using the complete model</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th colspan="2">STEP 1</th><th colspan="2">STEP 2</th><th colspan="2">GAIN</th></tr><tr><th>n<sub>1</sub></th><th>t<sub>1</sub>(s)</th><th>n<sub>2</sub></th><th>t<sub>2</sub>(s)</th><th>t<sub>multilev</sub></th><th>G</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-style: italic;">Orthogonal array</td></tr><tr><td>8</td><td>1,140</td><td>0.7</td><td>86</td><td>261</td><td>4,837</td><td>9.0</td></tr><tr><td>16</td><td>1,580</td><td>0.8</td><td>96</td><td>407</td><td>5,763</td><td>6.3</td></tr><tr><td>32</td><td>1,240</td><td>0.8</td><td>200</td><td>435</td><td>11,936</td><td>14.2</td></tr><tr><td colspan="7" style="text-align: center; font-style: italic;">LHS sampling</td></tr><tr><td>5</td><td>1,040</td><td>0.7</td><td>94</td><td>371</td><td>5,094</td><td>9.1</td></tr><tr><td>10</td><td>1,040</td><td>0.6</td><td>171</td><td>598</td><td>9,312</td><td>10.0</td></tr><tr><td>15</td><td>1,400</td><td>0.8</td><td>201</td><td>649</td><td>11,113</td><td>10.3</td></tr><tr><td>20</td><td>1,040</td><td>0.7</td><td>186</td><td>452</td><td>10,598</td><td>10.1</td></tr><tr><td>25</td><td>1,620</td><td>0.9</td><td>167</td><td>495</td><td>9,878</td><td>8.3</td></tr><tr><td>30</td><td>1,160</td><td>0.7</td><td>113</td><td>238</td><td>7,357</td><td>7.3</td></tr><tr><td>35</td><td>1,780</td><td>0.9</td><td>98</td><td>254</td><td>6,843</td><td>5.6</td></tr><tr><td>40</td><td>1,200</td><td>0.9</td><td>100</td><td>233</td><td>7,203</td><td>5.6</td></tr><tr><td>45</td><td>1,880</td><td>1.1</td><td>109</td><td>244</td><td>7,924</td><td>5.7</td></tr><tr><td>50</td><td>1,440</td><td>1.0</td><td>180</td><td>499</td><td>11,833</td><td>6.9</td></tr><tr><td>55</td><td>1,400</td><td>1.0</td><td>152</td><td>392</td><td>10,650</td><td>5.9</td></tr><tr><td>60</td><td>1,140</td><td>0.9</td><td>184</td><td>442</td><td>12,554</td><td>6.7</td></tr><tr><td>65</td><td>1,380</td><td>1.0</td><td>138</td><td>313</td><td>10,444</td><td>5.2</td></tr><tr><td>70</td><td>1,520</td><td>1.1</td><td>97</td><td>230</td><td>8,592</td><td>4.6</td></tr><tr><td>75</td><td>1,140</td><td>0.9</td><td>125</td><td>267</td><td>10,290</td><td>4.7</td></tr><tr><td>80</td><td>1,300</td><td>1.0</td><td>119</td><td>195</td><td>10,239</td><td>4.9</td></tr></tbody></table>

<table><caption>**Table 4** Global optimization with bound and linear constraints using the metamodel and local optimization using the complete model</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th colspan="2">STEP 1</th><th colspan="2">STEP 2</th><th colspan="2">GAIN</th></tr><tr><th>n<sub>1</sub></th><th>t<sub>1</sub>(s)</th><th>n<sub>2</sub></th><th>t<sub>2</sub>(s)</th><th>t<sub>multilev</sub></th><th>G</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-style: italic;">Orthogonal array</td></tr><tr><td>8</td><td>1,417</td><td>2.4</td><td>81</td><td>249</td><td>4,581</td><td>8.7</td></tr><tr><td>16</td><td>1,037</td><td>2.3</td><td>94</td><td>257</td><td>5,661</td><td>7.5</td></tr><tr><td>32</td><td>1,037</td><td>2.3</td><td>96</td><td>237</td><td>6,587</td><td>10.2</td></tr><tr><td colspan="7" style="text-align: center; font-style: italic;">LHS sampling</td></tr><tr><td>5</td><td>1,037</td><td>2.3</td><td>124</td><td>415</td><td>6,639</td><td>11.0</td></tr><tr><td>10</td><td>1,037</td><td>2.3</td><td>148</td><td>390</td><td>8,131</td><td>11.2</td></tr><tr><td>15</td><td>1,037</td><td>2.3</td><td>94</td><td>255</td><td>5,610</td><td>8.2</td></tr><tr><td>20</td><td>1,037</td><td>2.3</td><td>122</td><td>237</td><td>7,308</td><td>8.7</td></tr><tr><td>25</td><td>1,037</td><td>2.3</td><td>132</td><td>685</td><td>8,079</td><td>5.8</td></tr><tr><td>30</td><td>1,097</td><td>2.4</td><td>160</td><td>288</td><td>9,777</td><td>9.2</td></tr><tr><td>35</td><td>1,037</td><td>2.3</td><td>174</td><td>464</td><td>10,754</td><td>7.5</td></tr><tr><td>40</td><td>1,157</td><td>2.5</td><td>132</td><td>323</td><td>8,851</td><td>6.4</td></tr><tr><td>45</td><td>1,057</td><td>2.5</td><td>94</td><td>247</td><td>7,153</td><td>5.1</td></tr><tr><td>50</td><td>1,037</td><td>2.5</td><td>111</td><td>247</td><td>8,285</td><td>5.7</td></tr><tr><td>55</td><td>1,037</td><td>2.5</td><td>161</td><td>326</td><td>11,115</td><td>6.4</td></tr><tr><td>60</td><td>1,037</td><td>2.5</td><td>137</td><td>346</td><td>10,137</td><td>5.7</td></tr><tr><td>65</td><td>1,097</td><td>2.5</td><td>197</td><td>418</td><td>13,481</td><td>6.4</td></tr><tr><td>70</td><td>1,237</td><td>2.6</td><td>124</td><td>255</td><td>9,983</td><td>5.2</td></tr><tr><td>75</td><td>1,057</td><td>2.5</td><td>123</td><td>272</td><td>10,189</td><td>4.6</td></tr><tr><td>80</td><td>1,037</td><td>2.5</td><td>102</td><td>235</td><td>9,366</td><td>4.4</td></tr></tbody></table>

model led to a trough corresponding to a local minimum due to the GA algorithm. Therefore, the optimization using the mechanical model led to that local minimum. For example, this was the case for the 5 LHS draws and for several others.

As explain in the introduction, our two step optimization approach is used for increasing the chance of finding the global optimum. In order to verify this point, the two step strategy is compared with another reference case based on a multi-start of the local optimization of the FE model (Step-2) with the multiparametric strategy. First of all, a couple of initial starting points by LHS 10 samples and OA 8 samples is selected. Then Step-2 run only with those points and the best of obtained optima are assigned as a global optimum. One consider only the optimization problem with bound and linear constraints. We can observe that

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Multi-start of the direct local optimization of the FE model with the multiparametric strategy</figcaption></figure>

for the LHS draw (Figure 12) the best optimum find is equal to 16,816N for a cumulative CPU time equal to 2169s and the global optimum is not achieved. Figure 12 shows that the global optimum is reached with the OA 8 samples draw in 1820s. Note that in the two step optimization approach the global optimum is reached in 526s for the OA draw of 8 combinations and in 728s for the LHS draw of 10 combinations. Our approach shows that:

- the global optimum is in most cases obtained;

- the CPU time is always lower.

## 4.2 Optimization of the transmissible axial load using a three-dimensional model

The evolution of the transmissible axial load $F_a$ at the interface was studied as a function of six parameters (the tension loads $F_{scr}$ in the screws, the friction coefficients and the gaps between the different parts in contact). We used a volume model (Figure 13). The mesh is composed of 21,195 four-node tetrahedron elements (17,736 degrees of freedom). The calculations
were carried out on a PC (Proc. Intel Xeon W5650(tm) 2.6 Ghz, RAM 12 GB)

<figure><img src="image_15.png" /><figcaption>Fig. 13 The substructures used for the LATIN method</figcaption></figure>

The constraints in terms of bounds on the parameters are given in Table 1. The tension load in the screws varied between 1,500N and 2,500N. Different draws of the sets of parameters were used in the metamodel construction preliminary phase:

- an orthogonal Design Of Experiments (DOE) combining 6 factors on two levels (8, 16, 32 and 64 combinations),

- a random draw (LHS) of between 5 and 45 combinations in steps of 5 combinations.

Let $n_0$ denote the number of combinations in each draw. Table 5 summarizes, for each draw —of $n_0$ combinations—the number of iterations ($i_0$) of the FE solver for the calculation of the $n_0$ combinations and the computation time per set of parameters. The time $t_{ref}$ of the first calculation and the minimum and average CPU times are listed separately for each set of parameters.

The optimization problem that we consider can be formulated as:

<table><caption>Table 5 Construction of the metamodel</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th rowspan="2">i<sub>0</sub></th><th colspan="3">CPU time per parameter set (s)</th><th rowspan="2">total<br>t<sub>0</sub>(s)</th></tr><tr><th>1st calc.<br>t<sub>ref</sub>(s)</th><th>min.<br>(s)</th><th>mean<br>(s)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-style: italic;">orthogonal array</td></tr><tr><td>8</td><td>1,036</td><td>91</td><td>60</td><td>71</td><td>568</td></tr><tr><td>16</td><td>721</td><td>90</td><td>3</td><td>26</td><td>414</td></tr><tr><td>32</td><td>2,491</td><td>92</td><td>3</td><td>45</td><td>1,426</td></tr><tr><td>64</td><td>1,856</td><td>95</td><td>3</td><td>17</td><td>1,077</td></tr><tr><td colspan="6" style="text-align: center; font-style: italic;">Latin hypercub sampling</td></tr><tr><td>5</td><td>521</td><td>77</td><td>34</td><td>47</td><td>237</td></tr><tr><td>10</td><td>906</td><td>84</td><td>27</td><td>41</td><td>413</td></tr><tr><td>15</td><td>1,261</td><td>81</td><td>14</td><td>39</td><td>584</td></tr><tr><td>20</td><td>1,586</td><td>80</td><td>23</td><td>36</td><td>725</td></tr><tr><td>25</td><td>1,921</td><td>79</td><td>12</td><td>35</td><td>882</td></tr><tr><td>30</td><td>2,326</td><td>84</td><td>18</td><td>35</td><td>1,064</td></tr><tr><td>35</td><td>2,686</td><td>77</td><td>21</td><td>35</td><td>1,230</td></tr><tr><td>40</td><td>2,921</td><td>84</td><td>7</td><td>33</td><td>1,336</td></tr><tr><td>45</td><td>3,201</td><td>77</td><td>16</td><td>32</td><td>1,461</td></tr></tbody></table>

$$
\left\{
\begin{array}{l}
\displaystyle \arg \max_{(\mu_{./}, j_{./}, F_{scr})} F_a(\mu_{./}, j_{./}, F_{scr}) \\[1.5ex]
\text{such that} \\
\quad 0.05 \le \mu_{shaft/gear} \le 0.5 \\
\quad 0.05 \le \mu_{gear/inner} \le 0.5 \\
\quad 0.05 \le \mu_{inner/outer} \le 0.5 \\
\quad 0\mu m \le j_{shaft/gear} \le 49\mu m \\
\quad 36\mu m \le j_{gear/inner-ring} \le 106\mu m \\
\quad 1500N \le F_{scr} \le 2500N \\
\quad \mu_{inner/outer} - \mu_{shaft/gear} \le 0 \\
\quad \mu_{inner/outer} - \mu_{gear/inner} \le 0
\end{array}
\right.
\tag{27}
$$

The optimum values of the transmissible axial load at the interface given by the two optimization steps are shown in Figure 14 for the only bounds constraints problem and in Figure 15 for the problem with linear inequality constraints.

At the end of this study, the optimum found was 1,038kN, associated with the following optimum parameters:

$$- \mu_{shaft/gear} = 0.5$$

$$- \mu_{gear/inner} = \mu_{inner/outer} = 0.05$$

$$- F_{scr} = 2,500N$$

$$- j_{shaft/gear} = 0\mu m$$

$$- j_{gear/inner-ring} = 36\mu m$$

The CPU times and the numbers of iterations which were necessary to carry out these two optimization steps are summarized in Table 6 and 7. Overall, the total time $t_{multilev+multipar}$ for this two step optimization strategy varied between approximately 329s and 1,602s, and the variation—in the case of LHS draws— was approximately proportional to $n_0$, the number of draws. Let us
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Comparison of the optima obtained from Step 1 and Step 2 with bound constraints</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Comparison of the optima obtained from Step 1 and Step 2 with bound and linear constraints</figcaption></figure>

<table><caption>**Table 6** Global optimization with bound constraints using the metamodel and local optimization using the complete model</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th colspan="2">STEP 1</th><th colspan="2">STEP 2</th><th colspan="2">GAIN</th></tr><tr><th>n<sub>1</sub></th><th>t<sub>1</sub>(s)</th><th>n<sub>2</sub></th><th>t<sub>2</sub>(s)</th><th>t<sub>multilev</sub></th><th>G</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-style: italic;">Orthogonal array</td></tr><tr><td>8</td><td>1,420</td><td>0.8</td><td>68</td><td>76</td><td>6,387</td><td>9.9</td></tr><tr><td>16</td><td>1,500</td><td>0.8</td><td>195</td><td>239</td><td>17,731</td><td>27.1</td></tr><tr><td>32</td><td>1,180</td><td>0.8</td><td>88</td><td>87</td><td>10,084</td><td>6.7</td></tr><tr><td>64</td><td>1,040</td><td>0.9</td><td>61</td><td>52</td><td>10,505</td><td>9.3</td></tr><tr><td colspan="7" style="text-align: center; font-style: italic;">LHS sampling</td></tr><tr><td>5</td><td>1,120</td><td>0.7</td><td>102</td><td>91</td><td>8,992</td><td>27.4</td></tr><tr><td>10</td><td>1,140</td><td>0.7</td><td>45</td><td>72</td><td>4,622</td><td>9.5</td></tr><tr><td>15</td><td>1,560</td><td>0.9</td><td>57</td><td>49</td><td>6,051</td><td>9.5</td></tr><tr><td>20</td><td>1,120</td><td>0.7</td><td>45</td><td>66</td><td>5,463</td><td>6.9</td></tr><tr><td>25</td><td>1,520</td><td>0.9</td><td>92</td><td>81</td><td>9,832</td><td>10.2</td></tr><tr><td>30</td><td>1,040</td><td>0.7</td><td>55</td><td>49</td><td>7,143</td><td>6.4</td></tr><tr><td>35</td><td>1,260</td><td>0.8</td><td>73</td><td>87</td><td>9,076</td><td>6.9</td></tr><tr><td>40</td><td>1,720</td><td>1.1</td><td>73</td><td>86</td><td>9,496</td><td>6.7</td></tr><tr><td>45</td><td>1,220</td><td>0.9</td><td>70</td><td>69</td><td>9,664</td><td>6.3</td></tr></tbody></table>

<table><caption>**Table 7** Global optimization with bound and linear constraints using the metamodel and local optimization using the complete model</caption><thead><tr><th rowspan="2">n<sub>0</sub></th><th colspan="2">STEP 1</th><th colspan="2">STEP 2</th><th colspan="2">GAIN</th></tr><tr><th>n<sub>1</sub></th><th>t<sub>1</sub>(s)</th><th>n<sub>2</sub></th><th>t<sub>2</sub>(s)</th><th>t<sub>multilev</sub></th><th>G</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-style: italic;">Orthogonal array</td></tr><tr><td>8</td><td>1,069</td><td>2.3</td><td>183</td><td>178</td><td>16,052</td><td>21.4</td></tr><tr><td>16</td><td>1,029</td><td>2.3</td><td>204</td><td>235</td><td>18,489</td><td>28.4</td></tr><tr><td>32</td><td>1,169</td><td>2.4</td><td>131</td><td>174</td><td>13,699</td><td>8.6</td></tr><tr><td>64</td><td>1,029</td><td>2.5</td><td>45</td><td>66</td><td>9,162</td><td>8.0</td></tr><tr><td colspan="7" style="text-align: center; font-style: italic;">LHS sampling</td></tr><tr><td>5</td><td>1,029</td><td>2.3</td><td>166</td><td>202</td><td>14,371</td><td>32.6</td></tr><tr><td>10</td><td>1,029</td><td>2.3</td><td>45</td><td>66</td><td>4,624</td><td>9.6</td></tr><tr><td>15</td><td>1,189</td><td>2.4</td><td>175</td><td>224</td><td>15,968</td><td>19.7</td></tr><tr><td>20</td><td>1,149</td><td>2.4</td><td>45</td><td>66</td><td>5,464</td><td>6.9</td></tr><tr><td>25</td><td>1,029</td><td>2.3</td><td>176</td><td>181</td><td>16,892</td><td>15.8</td></tr><tr><td>30</td><td>1,029</td><td>2.4</td><td>239</td><td>252</td><td>22,606</td><td>17.1</td></tr><tr><td>35</td><td>1,029</td><td>2.4</td><td>45</td><td>66</td><td>6,725</td><td>5.2</td></tr><tr><td>40</td><td>1,029</td><td>2.5</td><td>199</td><td>245</td><td>20,085</td><td>12.7</td></tr><tr><td>45</td><td>1,029</td><td>2.5</td><td>123</td><td>130</td><td>14,119</td><td>8.9</td></tr></tbody></table>

compare this two step optimization method based on a multiparametric calculation strategy with the same optimization method and no multiparametric aspect. In order to do that, we consider that the reference time for a FE calculation is considered to be equal to $\overline{t_{ref}} = 84.0$s. Gain $G$, defined as the ratio of $t_{multilev}+multipar$ (the computation time for our two step strategy coupled with the multiparametric aspect) to $t_{multilev}$ (the computation time for the two step calculation strategy with no multiparametric aspect), shows that in a two step strategy the addition of the multiparametric strategy reduced the CPU time by a factor of between approximately 5 and 32.5.

## 5 Conclusion

We implemented a two step model optimization strategy based on a global optimization using a cumulative interpolation metamodel and a local optimization using the full FE model. What is unique to this work is the coupling between this two step strategy and a multiparametric resolution strategy, which helps reduce computation times during the numerous calls to the simulator, both in the metamodel construction step using a predefined sample and in the local optimization step for the evaluation of the objective function and its gradients. We compared the computation times obtained using our strategy with those obtained using a two step optimization strategy without multiparametric aspect. We showed that this two step/multiparametric association led to a significant reduction in computation times. The combination of the two strategies led to a gain of 4.4 to 14.2 for the axisymmetric example and 5 to 32.5 for the 3D example.
The extension to larger problems, in terms of both the number of parameters and the number of simulations, remains to be carried out. Another direction for future research would be the introduction of multiple objective functions in order to determine the Pareto fronts.

**Acknowledgements** This work was supported by the French National Research Agency (ANR) as part of the RNTL 2005 program: Multidisciplinary Optimization Project (OMD).

## References

Barthelemy JF, Haftka R (1993) Approximation concepts for optimum structural design - a review. Structural Optimization 5:129-144

Bendsoe MP (1995) Optimization of Structural Topology, Shape and Material. Springer-Verlag, Heidelberg

Blanzé C, Champaney L, Cognard J, P L (1995) A modular approach to structure assembly computations. application to contact problems. Engineering Computations 13(1):15-32

Boucard PA, Champaney L (2003) A suitable computational strategy for the parametric analysis of problems with multiple contact. International Journal for Numerical Methods in Engineering 57(9):1259-1282

Boucard PA, Champaney L (2004) Approche multirésolution pour l'étude paramétrique d'assemblages par contact et frottement. Revue Européenne des Eléments Finis 13(5/7):437-448

Boucard PA, Buytet S, Guidault PA (2007) Une stratégie multi-échelle pour l'étude paramétrique de détails géométriques au sein de structures en contacts multiples. Revue Européenne de Mécanique Numérique 16(8):1011-1036

Box GEP, Hunter JS (1957) Multi-factor experimental designs for exploring response surfaces. The Annals of Mathematical Statistics 28(1):195-241

Box GEP, Wilson KB (1951) On the experimental attainment of optimum conditions. Journal of the Royal Statistical Society Series B (Methodological) 13(1):1-45

Braibant V, Fleury C (1985) An approximation concepts approach to shape optimal design. Comp Meth Appl Mech Eng 53:119-148

Chen TY, Yang CM (2005) Multidisciplinary design optimization of mechanisms. Advances in engineering software 36(5):301-311

Cressie N (1990) The origins of kriging. Mathematical Geology 22:239-252, 10.1007/BF00889887

Currin C, Mitchell T, Morris M, Ylvisaker D (1991) Bayesian prediction of deterministic functions, with applications to the design and analysis of computer

experiments. Journal of the American Statistical Association 86(416):pp. 953-963

Dantzig G, Orden A, P W (1955) Generalized simplex method for minimizing a linear form under linear inequality restraints. Pacific Journal Math 5:183-195

El-Sayed ME, Hsiung CK (1991) Optimum structural design with parallel finite element analysis. Computers and Structures 40(6):1469-1474

Engels H, Becker W, Morris A (2004) Implementation of a multi-level optimisation methodology within the e-design of a blended wing body. Aerospace Science and Technology 8(2):145-153

Fleury C, Braibant V (1986) Structural optimization. a new dual method using mixed variables. International Journal for Numerical Methods in Engineering 23:409-428

Gingold R, Monaghan J (1977) Smooth particle hydrodynamics: theory and application to non-spherical stars. Royal Astronomical Society, Monthly Notices 181:375-389

Goldberg D (1989) Genetic Algorithms in Search, Optimization and Machine Learning. Addison Wesley Company Inc.

Haftka R (1988) First- and second-order constraint approximations in structural optimization. Comp Mech 3:89-104

Han SP (1977) A globally convergent method for nonlinear programming. Journal of Optimization Theory and Applications 22:297-309

Hardy RL (1971) Multiquadric equations of topography and other irregular surfaces. J Geophys Res 76(8):1905-1915

Haykin S (1994) Neural Networks: A Comprehensive Foundation, 1st edn. Prentice Hall PTR, Upper Sadle River, NJ, USA

Hilding D, Torstenfelt B, Klarbring A (2001) A computational methodology for shape optimization of structures in frictionless contact. Comput Methods Appl Mech Engrg 190:4043-4060

Keane AJ, Petruzzeli N (2000) Aircraft wing design using ga-based multi-level strategies. In: Proceedings 8th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Long Beach, USA, pp A00-40 171

Kennedy J, Eberhart R (1995) Particle swarm optimization. In: IEEE Int. Conf. Neural Networks, vol 4, pp 1942-1948

Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by simulated annealing. Science 220:671-680

Kravanja S, Sorsak A, Kravanja Z (2003) Efficient multilevel minlp strategies for solving large combinatorial problems in engineering. Optimization and engineering 4(1/2):97-151
Ladevèze P (1999) Nonlinear Computational Structural Mechanics - New Approaches and non-Incremental Methods of Calculation. Springer-Verlag, Berlin

Li W, Li Q, Steven GP, Xie YM (2005) An evolutionary shape optimization for elastic contact problems subject to multiple load cases. Computer Methods in Applied Mechanics and Engineering 194:3394-3415

Liu B, Haftka RT, Watson LT (2004) Global-local structural optimization using response surfaces of local optimization margins. Structural and Multidisciplinary Optimization 27(5):352-359

McKay MD, Beckman RJ, Conover WJ (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21(2):239-245

Montgomery D (1997) Design and analysis of experiments. Wiley

Nayrolles B, Touzot G, Villon P (1992) Generalizing the finite element method: Diffuse approximation and diffuse elements. Computational Mechanics 10(5):307-318

Pedersen P (1981) Optimization of distributed parameters structures, vol 1, Springer Verlag, Amsterdam, chap The integrated approach of FEM-SLP for solving problems of optimal design, pp 757-780

Pritchard J, Adelman H (1990) Differential equation based method for accurate approximations in optimization. In: Proc. AIAA/ASME/ASCE/AHS/ASC 31st Structures, Structural Dynamics and Materials Conf. (held in Long Beach, CA)

Pritchard J, Adelman H (1991) Differential equation based method for accurate modal approximations. AIAA Journal 29:484-486

Rasmussen J (1998) Nonlinear programming by cumulative approximation refinement. Structural and Multidisciplinary Optimization 15(1):1-7

Robinson GM, Keane AJ (1999) A case for multi-level optimisation in aeronautical design. Aeronautical Journal 103:481-485

Roux WJ, Stander N, R-T H (1998) Response surface approximations for structural optimization. International Journal for Numerical Methods in Engineering 42:517-534

Sacks J, Schiller SB, Welch WJ (1989) Designs for computer experiments. Technometrics 31(1):41-47

Sakata S, Ashida F, Zako M (2003) Structural optimization using kriging approximation. Computer Methods in Applied Mechanics and Engineering 192:923-939

Soulier B, Richard L, Hazet B, Braibant V (2003) Crashworthiness optimization using a surrogate approach by stochastic response surface. In: Gogu G, Coutellier D, Chedmail P, Ray P (eds) Recent Advances in Integrated Design and Manufacturing in Mechanical Engineering, Kluwer Academic Publishers, pp 159-168

Umesha PK, Venuraju MT, Hartmann D, Leimbach KR (2005) Optimal design of truss structures using parallel computing. Structural and Multidisciplinary Optimization 29:285-297

Wu ZM, Schaback R (1993) Local error estimates for radial basis function interpolation of scattered data. IMA J Num Anal 13(1):13-27

Zienkiewicz O, Campbell J (1973) Optimum structural design. New York John Wiley & Sons