Pareto Set Computation in Convex Multi-objective Design using Adaptive Response Surface Method (ARSM)

. Optimal design problems normally involve high dimensional design spaces and multiple objective functions. Depending on the complexity of the model, the time required to explore the design space could become excessive. This paper describes the calculation of the Pareto-optimal set based on adaptive surface methodology (ARSM) in order to reduce simulation times given a finite element analysis (FEA) simulation model. The Pareto-optimal strategy consists in the solution of a set of different single-objective problems. Each of this points is found via ARSM. The implementation of ARSM aims to use a few initial simulation points to approximate accurately the set of single-objective functions required. The methodology reduces significantly the number of points required to compute the efficient set compared to other strategies (e.g the exhaustive method), proving to reduce the simulation time of a computationally intensive model.


Introduction
Co mputational numeric simu lations of complex geometric models have become essential in the mechanical design process of the last decades.The implementation of such methods make it possible to solve problems in which the solution by an analytical formulation is not possible.
The field of optimal design has been greatly benefited by its imp lementation given the relation between model complexity and time processing improvement.In a typical design optimizat ion where a selected design space is constrained by design restrictions, it's desired to find the best solution according to design requirements.However, even with the imp lementation of co mputational aided calculations, an exhaustive explorat ion of the solution space leads to extended simulation t imes that are not practical, particu larly fo r large co mputational intensive models.
One of the proposed tools to solve this kind of problems is the response surface method (RSM ) [1].It consists of the combination of several techniques including mainly a design of experiments (DOE) and model fitting.The main purpose of RSM is to find an approximate form of the objective function, which allows the use of few sample po ints instead of making an exhaustive exploration.

Response Surface Method Review
According to Gosavi [2], RSM begins by choosing the points to be calculated fro m the design space.The selection of these points should be performed in order to minimize the nu mber of simulat ions with the best possible quality of the fitting model.Several classic DOE techniques have been widely developed to achieve this goal, including factorial designs, central composite designs, Lat in hypercube designs and low d iscrepancy sequences, just to name a few.
After the selection of the design variable space, an expected behaviour of the response represented by a meta-model is selected.In the most classic approach [1], this behaviour is assumed for simp licity of second order or even higher order degree polynomial, depending on the model's co mplexity.An appropriate fitting method for that meta-model is used to generate the response surface.Techniques such as linear regression, splines regression and neural networks are co mmonly used as fitting methods.Approximation quality can be proved with a statistical test.For least squares regression method, the parameter r 2 is normally used.Nevertheless, r 2 misleads to good results when the assumed polynomial is high order o r there are several design parameters.Further validation tests can be found in [3], [4] and [5].
In synthesis, the classical RSM method involves DOE, meta-model selection, fitting techniques and validation trough statistical methods.
Avoiding to restrict the response to a predefined set of points in the design space, new approaches have been developed in the last years that iteratively update the design space with the goal of achieving an optimu m.The update consists in finding the minimu m of the first response surface obtained by the typical techniques described before with an optimization algorithm, adding the design variables of the min imu m point to the initial design space and removing the point (or points, given a threshold) with the highest objective function value.The methodology was first proposed by Wang ([6], [7]) with the name of Adaptive Response Surface Methodology (ARSM).The method proposed by Wang involves the domain reduction of each design variable, which its time consuming given the min imization process for both upper and lower limit.Steenackers [8] proposed a revised method that does not involve the reduction of the design space, but rather a strategy called "pan and zoo m" (for more details refer to [8]).Co mpared to Wang's method, the pan and zoom method showed a reduced time per iteration as well as the total number of iterations required to converge to the min imu m.For further examp les using adaptive method, refer to [9] - [14].
The present study is based on the ASRM methods proposed by both Wang and Steenackers to explore the solution space of a finite element analysis (FEA) 2D simu lation.In addition to p revious works, this paper aims to apply the method to a mu ltivariable design space and extend it to a multiple objective optimization functions.
Particularly, the study proposes the implementation of the ARSM to a single objective function that contains the informat ion of a mult iple objective problem using weights.Given a defined number of weight co mbinations to be determined, each set of them requires an implementation of the algorithm.The final result is a Pareto-set plot efficiently obtained by the use of the ARSM, avoiding an exhaustive method which involves the simulat ion of the comp lete design space.The consequence of this is a significant reduction in simu lation times (without even taking into account the fact that for several design variables and objectives, the enumeration method fails to produce a result in reasonable time).

Methodology -ARSM
This paper imp lements an ARSM method jointly with the Pareto set computation of a convex mult i-objective design problem.The proposed approach is based on the formulat ion and solution of a set of single-objective problems, where the single objective function is obtained as weighted sum of the objective functions of the orig inal problem ( [18], [19]).Weights are assigned to each single function as follows: Where F is a linear comb ination of each objective function f i ; and w i is the given weight of each function, being the sum of weights equal to 1.The minimu m of this function is calculated using ARSM for each pair of weights.The aim is to obtain an approximation of the Pareto set using considerably less points than in the case of simp le nu merat ion of an exhaustive exp lored design space.
Particularly, the ARSM method is based on the "Pan and Zoom" method [8].It consists in the following steps: 1-The method begins by the selection of points from the design space.For this study, a Sobol Low discrepancy sequence ( [15], [16]) was used to explore unifo rmly the space.The Sobol set was configured to also include the boundary of the design space.2-The selected points are simulated to obtain the real value of the objective function.3-The obtained objective function values are fitted with a second order degree polynomial.4-The response surface is min imized using a local optimization algorithm.In this case, an algorithm fro m a co mmercial software was used (i.e.Mat lab), which corresponds to Nelder-Mead method [17].5-The calculated minimu m is added to the initial set of points which represents a promising search direction fro m the real object ive function minimu m.Simu ltaneously, the point with the highest objective function value is erased.This step reduces the design space which increase the quality of the response surface in the minimal region.6-With a user defined convergence criteria, the objective function value obtained is compared to the response surface value.If the difference achieves the criteria, the minimu m o f the surface is taken as the real objective function.In the opposite case, the method returns to step 2 and repeat the process until the criteria is met.
The aforementioned methodology is repeated for each pair of weights used to construct an approximation of the Pareto set.
An exhaustive explo ration was performed to determine the real Pareto set.The set was adjusted to a bi-exponential regression.The same fit was applied to the approximated Pareto set found by ARSM.
To compare the accuracy of the appro ximation to the complete design space Pareto front, the determination coefficient R 2 [4] was used: 2  (2) Where y i is the i-th value of the Pareto front approximation using ARSM, y mean is the mean value fro m the Pareto front approximation using ARSM, and f i is the i-th value from the complete design space Pareto

Problem formulation
In order to prove the scope of the imp lementation, a 2D cantilever beam was simu lated including an elastomeric interface embedded on a poly meric matrix instead of a common fixed support.The simu lated geometry is shown on Figure 1.For the sake of simp licity, all materials were considered within their linear elastic range.Two design 05003-p.2parameters where considered for each simu lation, where P1 is the initial thickness of the elastomeric interface and P2 is the gap between the polymeric mat rix and the beam.A plane stress formulation was implemented.Given the nature of the case, a 10 N load was applied at the beam's tip.The supports were fixed in the upper, lower and rear sides of the polymeric matrix.Finally, the elastomeric interface p re-stress generated by the variation of P1 and P2 was represented as a distributed load on the contact zone between the elastomeric interfaces and the upper/lower beam sides.The deflection will be reduced with a greater distributed load, given that the elastomeric interfaces are able to support a greater reacting mo mentu m generated by the load force at the beam's tip.In the other hand, a greater distributed load will increase the system's natural frequency due to an increment in the stiffness.

Simulation conditions
Fro m the FEA point of view, the mesh used for the polymeric matrix and the elastomeric interface was parameterized as function of both variables [P1, P2] to keep the element size according to the change of geometry.Th is relationship ratio between element size and parameter size was determined after a convergence study.An 8 node high order quadratic element was used in order to accurately capture high deformat ion behavior.Table 1 shows the summarized information for each simulation component.In terms of contacts, all the contacts between the polymeric matrix and the elastomeric interfaces were considered as bonded contacts.This resembles an adhesive union between the elements.All the contacts between the beam and the elastomeric interfaces were simu lated as an infin ite friction coefficient union.Given that the load is applied in the vertical direction (and the fact that the problem is defined with a plane stress formulat ion), the behavior of this contact allows to accurately study the deflection at the beam's t ip, and first resonant frequency.The vib ration mode corresponding to the first resonant frequency is shown in Figure 3.As explained in the section 3.2, the contacts between the elastomeric interfaces and the surface were modeled as an infinite friction coefficient union.Therefore the parameters [P1, P2] were constrained to be P1 > P2 in order to simu late only co mpression conditions for the elastomeric interfaces.
In addition to the co mpression, the ranges for [P1, P2] were assumed for a maximal 20% deformation of the elastomeric interface.In this way, the maximal distributed load applied by the pre-stress will be 615 N. Table 2 shows the limits defined for each variable

Results
Figure 4 shows an examp le of a normalized surface corresponding to the single-objective problem for a g iven pair of weights as a result of the second order regression.Figure 5 shows the mult i-objective p lot between deflection and frecuency considering the points of the completely simulated design space and those obtained by the ARSM for 10 pair of weights.A total of 418 points where used in the first case, and it took 55 evaluations of the simu lation to compute the pareto set in the second case.
Figure 6 shows the comparisson between the fitting of a full exp lored pareto front (solid line), and the ARSM approximations obtained with different in itial simulat ion points (n).As expected, a higher number of init ial points will produce a better quality Pareto front approximation.This observation is confirmed numerically by calculating R 2 determination coefficient named in section 2. Table 3 resumes the simu lations results for the three cases (exhaustive method, 10 in itial points and 18 initial points).There exist a tradeoff between the quality of the approximation and the initial nu mber of points.In any case, both ARSM appro ximations required only a fract ion of the computational effort used to compute the Pareto set.

Conclusions
The methodology shows a significant imp rovement in computation times compared to the simulat ion of the complete design space followed by the exhaustive 05003-p.4 method of simp le enumeration, thus saving time in a highly computational intensive model.In this case, the co mputation of the Pareto-optimal set for a g iven set of weights represented only 13% of the total simulat ion time requiered for the simu lation of the complete design space: 14 hours for exahustive method, 1 hour 49 minutes for 18 in itial ARSM points and 1 hour 24 minutes for 10 in itial ARSM points; with a marg inal compro mise in accuracy.In the 18 points case, the determination coefficient was 0.9939 with respect to the exhaustive Pareto-set, while in the 10 points case, it was 0.9756.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 201

Figure 2
Figure2presents the load condition of the simulat ion.Given the nature of the case, a 10 N load was applied at the beam's tip.The supports were fixed in the upper, lower and rear sides of the polymeric matrix.Finally, the elastomeric interface p re-stress generated by the variation of P1 and P2 was represented as a distributed load on the contact zone between the elastomeric interfaces and the upper/lower beam sides.

Figure 6 .
Figure 6.Pareto front comparison between ARSM approximation and exhaustive method.

Table 1 .
FE Simulation characteristics

Table 2 .
Variable simulation range

Table 3 .
Simulations results