Surrogate-Assisted Multiobjective Evolutionary Algorithms for Structural Shape and Sizing Optimisation

Thework in this paper proposes the hybridisation of the well-established strength Pareto evolutionary algorithm (SPEA2) and some commonly used surrogate models. The surrogate models are introduced to an evolutionary optimisation process to enhance the performance of the optimiserwhen solving design problemswith expensive function evaluation. Several surrogatemodels including quadratic function, radial basis function, neural network, and Kriging models are employed in combination with SPEA2 using real codes. The various hybrid optimisation strategies are implemented on eight simultaneous shape and sizing design problems of structures taking into account of structural weight, lateral bucking, natural frequency, and stress. Structural analysis is carried out by using a finite element procedure. The optimum results obtained are compared and discussed. The performance assessment is based on the hypervolume indicator. The performance of the surrogate models for estimating design constraints is investigated. It has been found that, by using a quadratic function surrogate model, the optimiser searching performance is greatly improved.


Introduction
Since they have been invented for several decades, the implementation of evolutionary optimisers on a wide variety of real-world design problems has successfully been made in numerous fields.The algorithms are attractive and popular as they can responsd to unavoidable disadvantages of classical mathematical programming.The evolutionary algorithms (EAs) can deal with almost all kinds of design problem as their search mechanisms, to some extent, rely on randomisation and need no function derivatives, for example, [1][2][3][4].Recently, single-objective evolutionary methods that have outstanding performance in several applications are realcode ant colony optimisation (ACOR) [5], covariance matrix adaptation evolution strategy (CMA-ES) [6], and differential evolution (DE) [7].The methods are robust and capable of reaching global optima.Most importantly, for multipleobjective design cases, they can be used to explore a Pareto optimal front of the problem within one simulation run.However, EAs have some unacceptable drawbacks that are a complete lack of consistency and low convergence rate.As a result, the obtained design results are considered near local optima for single objective cases and approximate Pareto fronts for multiobjective design [8].
With those aforementioned shortcomings, it is difficult or even impossible to employ the methods for solving the design problems with expensive function evaluation or limited number of function values available.Such obstructions can however be alleviated by introducing a surrogate model (SM) to an evolutionary optimisation procedure.With the use of such a model, an approximate function can be constructed, and inexpensive function prediction can be achieved.The hybridisation of EAs and SM for design optimisation has been studied since the last decade, and it has been found to greatly enhance the performance of EAs.Many researchers and engineers have focused their efforts on this research issue.The design problems can have single or multiple design objectives, for example, in [9][10][11][12][13][14] while mostly the demonstration problems are unconstrained or box-constrained.The hybrid EAs with SMs have been implemented on a wide variety of engineering application, for example, aerodynamics [15][16][17], heat and mass transfer [18][19][20], and structures [13,14,[21][22][23].Design problems requiring a surrogate-assisted optimisation may be roughly categorised to have 9 classes as given in [10].Surrogate-assisted EAs not only employ function approximation but also use some mathematical programming to enhance the searching performance [24,25].A comprehensive survey on the surrogate models used for enhancing EAs performance can be found in [26][27][28][29].
The surrogate models commonly employed include polynomial regression techniques [11,20,22,26,27,30], radial basis function interpolation (RBF) [9,11,12,16,18,27], artificial neural network (ANN) [9,15,23,26], support vector machines [31], moving least square method [32], and Gaussian process model (Kriging model) [10,16,[25][26][27].Recent alternative surrogate-based approaches are concerned with the approximation of field vectors (such as velocity and pressure fields for computational fluid dynamics and displacement and stress fields for structural analysis) rather than approximating objective and constraint functions directly.The proper orthogonal decomposition (POD) technique is employed for this task [33][34][35].With this concept, a designer can deal with more complicated design problems such as design optimisation when aeroelasticity is taken into consideration.Apart from that, there exist optimisers with the use of grid-based computing integrating nongeneration evolutionary algorithms, surrogate models, local search, and Lamarckian learning process together.This approach is termed asynchronous metamodel-based evolutionary algorithms [36,37].A technique used in the design of computer experiment phase is usually a Latin hypercube sampling technique.Recently, there have been several papers focusing on the development of optimum Latin hypercube sampling [38,39].Another interesting approach is an infill sampling technique [40].
The most challenging task for surrogate-based optimisation is to improve the convergence rate.To the authors' best knowledge, the work in the literature mostly demonstrated applying the design approaches to unconstrained or box-constrained optimisation which is simple to handle by using EAs.Nevertheless, there are usually many nonlinear constraints in the real-world design problems particularly in structural optimisation.The approximation of objective function, although being inaccurate, can lead to real optima (referred to as "Blessing of Uncertainty" in [11]).On the other hand, in estimating constraint functions, the approximate function needs to be as precise as possible so that the evolutionary search can land in the real feasible region.Inaccurate constraint function approximation can lead the optimiser to an undesirable design space.
In this paper, the hybridisation of the well-established multiobjective evolutionary algorithm (MOEA), namely, strength Pareto evolutionary algorithm [41] and some popular surrogate models is developed.The work is concerned with the studies of using several SMs for performance enhancement of multiobjective evolutionary algorithms in solving structural constrained optimisation.A number of surrogate models including quadratic function, radial basis function, neural network, and Kriging models are employed in combination with SPEA2 using real codes.The various hybrid optimisation strategies are implemented on eight simultaneous shape and sizing design problems of a structure having design functions as structural weight, lateral bucking, natural frequency, and stress.Structural analysis is carried out by using finite element analysis (FEA).The optimum results obtained are compared and discussed.The performance of SMs for estimating design constraints is investigated.It has been found that, by using a quadratic surrogate model, the searching performance of SPEA2 is greatly improved.

Hybridisation of SPEA2 and a Surrogate Model
There have been several ways to properly integrate SMs into EAs.The hybrid algorithm proposed in this paper is similar to what is presented in [18].The SPEA2 with real codes is employed as the main MOEA procedure.The surrogateassisted SPEA is illustrated in Figure 1.The computational steps in the figure can be separated into two parts as the main SPEA procedure and the SPEA subproblem that uses an approximation model.The main SPEA is slightly modified from the original version presented in [41].The real code recombination and mutation are similar to the work presented in [18].The procedure starts with an initial population, an (empty) Pareto archive, and some other predefined parameters, for example, an external archive size and recombination and mutation rates.The fitness of the individuals in the population is then assigned by performing actual function evaluation.The current population and their corresponding fitness values are brought to the surrogate environment as some of them are picked to build a surrogate model.Having such an approximation model, the SPEA is employed to solve the optimisation problem with approximate function evaluation.
The nondominated solutions from this stage are put back to the main SPEA procedure where their actual fitness values are determined.The current population, the external archive, and the nondominated solutions from performing a surrogate-based optimisation are put together while their nondominated solutions are determined.The external archive is then updated by replacing it with these nondominated solutions.In case that the number of the nondominated solutions exceeds the predefined archive size, some of them will be removed from the archive by means of the nearest neighbourhood method [41].Subsequently, the so-called binary tournament selection is performed to select some members in the newly updated archive to reproduce a set of offspring.The algorithm goes back to the step of fitness function evaluation and is repeated until the termination criterion is fulfilled.
Design solutions are selected from a current population to build a surrogate model based upon the idea that they should be evenly distributed throughout the design space.The proposed approach is operated on the domain of objective functions, which is somewhat similar to the adaptive grid algorithm, a Pareto archiving technique used in the Pareto archive evolution strategy [8].The procedure to select the solutions for surrogate modelling is illustrated in Figure 2. On the objective domain, all of the design points in the current population are plotted, and a proper grid is generated covering all the points.A solution that is the closest to the centre of each mesh is selected for function estimation.
Apply NNT if necessary Binary tournament selection

Surrogate Models
Let  = (x) are a function of a design vector x sized  × 1.Given a set of design solutions (or sampling points) X = [x 1 , . . ., x  ] and their corresponding function values Y = [ 1 , . . .,   ] selected during an evolutionary optimisation process, a surrogate or approximation model is constructed by means of curve fitting or interpolation.Approximation models used in this study are as follows.

Polynomial Regression.
The most commonly used polynomial surrogate model or a response surface model (RSM) is of the second-order polynomial or quadratic model, which can be expressed as where   for  = 0, . . ., ( + 1)( + 2)/2 are the polynomial coefficients to be determined.The coefficients can be found by using a regression or least square technique [26].

Kriging Model. A Kriging model (also known as a
Gaussian process model) used in this paper is the famous MATLAB toolbox named design and analysis of computer experiments (DACE) [42].The estimation of function can be thought of as the combination of global and local approximation models, that is, where (x) is a global regression model while (x) is a stochastic Gaussian process with zero mean and nonzero covariance representing a localised deviation.In this work, a linear function is used for a global model, which can be expressed as where The selected points for building a surrogate model for ,  = 1, . . .,  where  is the correlation function between any two of the  design points, and R is the symmetric correlation matrix size  ×  with the unity diagonal [26].The correlation function used herein is where   are the unknown correlation parameters to be determined by means of the maximum likelihood method.

Radial Basis Function Interpolation.
The radial basis function interpolation has been used in a wide range of applications such as integration between aerodynamic and finite element grids in aeroelastic analysis [43].The use of such a model for surrogate-assisted evolutionary optimisation is said to be commonplace [9,18,27].The approximate function can be written as where   are the coefficients to be determined, and  is a radial basis kernel.(Here it is set to be linear splines.)The coefficients can be found from the  sampling points as detailed in Srisomporn and Bureerat 2008 [18].

Neural Network. Artificial neural network (ANN) is
developed in response to the need to use complicated function expression rather than a simple quadratic or radial basis function.The model consists of interconnecting neurons or nodes that somewhat mimic the constructs of biological neurons as shown in Figure 3. ANN has been implemented  on a wide variety of real world applications.The use of ANN for surrogate-assisted optimisation can be found in the literature [9].The neural network model used in this paper is the feedforward multilayer perceptrons.The details of the network in this work are the learning algorithm = backpropagation and network training function = Levenberg-Marquardt backpropagation.
In this work, a simple heuristic algorithm is used to construct a network topology.The procedure starts with a network with two hidden layers while the number of nodes   on each layer is randomly chosen with the bound   ∈ [15,50].The transfer function for the hidden layer can be either hyperbolic tangent sigmoid or logarithmic sigmoid based on randomisation.If the termination criterion (training goal) is not met, randomly increase a number of nodes in the layers.If the number of nodes reaches the upper bound   = 50, add the third hidden layer to the network where a transfer function and node number are randomly chosen.In case that the stopping criterion is not fulfilled, increase the number of nodes.The algorithm repeats until the stopping criterion is met.

Testing Problems
Six testing multiobjective design problems are posed to design a torque arm [44] structure shown in Figure 4. Since it is a plate-like structure under inplane loading, design criteria to be taken into consideration are weight, stress, buckling, and natural frequency.The function calculation can be carried out by using FEA where shell elements are employed.The design variables and their bound constraints are given in Table 1 where other dimensions are set as   = 30 mm,   = 190 mm,  1 = 26 mm, and  2 = 27 mm.The six biobjective optimisation problems are termed F1, F2, F3, F4, F5, and F6, respectively, and can be expressed as: where  = structural weight,  = first mode natural frequency,  = lateral bucking factor (it is the ratio of critical load to applied load herein),  von = maximum von Mises stress on the structure, and  0 = 7.2688 rad/s.Furthermore, material properties are Young's modulus  = 60 × 10 6 N/cm 2 , Poisson's ratio ] = 0.3, yield stress   = 80000 N/cm 2 , density  = 0.00718 kg/cm 3 , and safety factor   = 1.5.
In fact, the F1, F2, and F3 problems are similar to the F4, F5, and F6 problems, respectively, except for the introduction of natural frequency constraints to the last three problems.This means that the performance of several surrogate models in estimating a natural frequency during an optimisation procedure can be examined.The bound constraints in Table 1 are not presented in the design problems as they can be dealt with in the SPEA procedure.It should be noted that the finite element analysis for these design studies is not time consuming as the work focuses on a comparative study of the several multiobjective optimisers.
Since the F1-F6 problems are small-scale, other two test problems are posed so as to investigate the capability of those surrogate models when dealing with large-scale problems.Two design problems of a simple wing-box displayed in Figure 5 are expressed as objective function: subject to where   is the th mode natural frequency of the wing-box.The wing-box is made up of Aluminium with Young modulus  = 70 GPa, Poisson's ratio ] = 0.34, density  = 2700 kg/m 3 , and yield strength   = 100 MPa.The structure consists of 2 spars, 4 ribs, and 4 skins (front, rear, top, and bottom skins) resulting in 52 wing segments as shown.The thickness of each segment is assigned to be a design variable; therefore, there are 52 design variables with the lower and upper bounds being 0.002 m and 0.010 m respectively.This can be classified as a large-scale sizing optimisation problem.The wing is subject to static loads on the front spar as shown in Figure 5. Shell elements are used for finite element analysis similarly to F1-F6.
The multiobjective optimisers used in this paper are named as follows: (i) SPEA-O the original SPEA2 without using a surrogate model, (ii) SPEA-Q SPEA using a quadratic response surface model, (iii) SPEA-R SPEA using the radial basis function interpolation, (iv) SPEA-N SPEA using a neural network model, (v) SPEA-K SPEA using a Kriging model.
Each evolutionary algorithm is employed to solve each problem ten runs with 30 generations (50 for F7-F8) and a population size of 20 (50 for F7-F8) for SPEA-O.For the algorithms using an approximation model, the number of generations and population size used in the main procedure is 15 (25 for F7-F8) and 20 (50 for F7-F8), respectively, while on each generation additional 20 (50 solutions for F7-F8) solutions are taken from optimising a surrogate subproblem.In the surrogate subproblem, the number of generations and population size is set to be 100 and 100, respectively.The external Pareto archive size is set to be 100.During the optimisation process, constraints are handled by using the approach presented in [45].It should be noted that there have been a number of techniques used to handle design constraints in evolutionary optimisation [46].We choose the aforementioned technique because it is simple and sufficiently effective.

Results and Discussion
The performance comparison of the five MOEAs for solving the eight design problems is based on the hypervolume indicator.Two types of performance assessment are conducted as detailed in Tables 2-4.In Table 2, each cell is the normalised mean value of 10 nondominated front hypervolumes obtained from using the optimiser on the corresponding column for solving the problem on the corresponding row.For each design problem, the best and worst mean values are normalised to be "1" and "0", respectively.The total normalised mean values given in the table are used for evaluating the overall performance of the surrogate-assisted SPEAs.From the results, the best algorithm for F1 is SPEA-R while the second best is SPEA-Q.SPEA-Q is the best for F2 while the second best is SPEA-N.The best optimisers for F3, F4, F5, and F6 are SPEA-K, SPEA-R, SPEA-K, and SPEA-K, respectively, whereas the second best for those 4 design problems is SPEA-Q.SPEA-O without using a surrogate model is the worst performer for F2, F5, and F6 while SPEA-N is the worst performer for F1, F3, and F4.The overall best performer Table 4: Comparison of ranking scores by t-test.
based on the evaluation is SPEA-Q whereas the second best is SPEA-K.For the large-scale problems, only SPEA-Q can surpass SPEA-O.SPEA-R is said to be equal to SPEA-O while SPEA-K gives that worst performance.Note that SPEA-N cannot be applied to these problems as it takes excessively long running time for network training.It can be said that only the quadratic response surface model is useful for the large-scale wing-box design.
The second performance assessment is based on the statistical t-test.For each design problem, to compare the 5 MOEAs, a performance matrix T sized 5 × 5 whose elements are full of "0" is generated.For the element   of the matrix, if the mean value of the 10 hypervolumes obtained from method  is significantly larger than that obtained from method  based on the statistical t-test at 95% confidence level, the value of   is set to be one.Table 3 shows the performance comparison of the problem F1.Having had the performance matrix, the 6th row on the table displays the total scores of each optimiser.The ranking is made in such a way that the best method is ranked as 1 while the worst is ranked as 5.The overall performance assessment based on the t-test is given in Table 4 where the last row on the table shows the total score of all the MOEAs.The lower the score the better the performer.Similar to the first assessment, the overall best method for F1-F6 based on this assessment is SPEA-Q while the close second best is SPEA-K.The worst method is SPEA-O without using a surrogate model.For the large-scale test problems, the only surrogateassisted method that is superior to the non-surrogate SPEA is SPEA-Q.Figures 6 and 7 show the search history as plots of hypervolume against the number of function evaluations of SPEA-O and SPEA-Q for solving F7 and F8, respectively.The values of front hypervolume can fluctuate due to the use of a Pareto archiving technique to remove some nondominated solutions from the Pareto archive.From the figures, it can be seen that SPEA-Q requires approximately half of the total number of function evaluations employed by SPEA-O to have equally good nondominated fronts.
The best approximate Pareto front of F1 obtained after numerous optimisation runs of the surrogate-assisted SPEAs is chosen and plotted in Figure 8 whereas the structures of some selected design points are shown in Figure 9. Figure 10 displays the best front of the F2 problem where some selected design solutions are illustrated in Figure 11.The best front and some selected optimal structures of the design problem F3 are displayed in Figures 12 and 13, respectively.Similarly, the best approximate fronts and some selected design solutions of F4, F5, and F6 are illustrated in Figures 14,15,16,17,18,and 19.Figures 20 and 21 show the best fronts obtained from using the various SPEAs for solving F7 and F8, respectively.It can be seen that the best fronts are from SPEA-Q, and SPEA-O.The SPEA-Q front does not totally dominate that of SPEA-O but they overlap each other while the SPEA-Q front has obviously, greater front extension.
According to the aforementioned performance assessments, it can be said that SPEA-R is the best algorithm in the cases of the F1 design problem while the close second best is SPEA-Q.When adding a natural frequency constraint to F1 becoming the F4 problem, the overall best and second best methods are still SPEA-R and SPEA-Q, respectively.For the F2 case and the F5 design problem which is F2 with additional natural frequency constraint, the top two performers are SPEA-Q and SPEA-K.The same can be said for the case of F3 while SPEA-K outperforms the others for the F6 design problem.The surrogate models, that is, the quadratic response surface model, the radial basis function interpolation, and the Kriging model, are said to enhance the searching performance of SPEA for the small-scale problems.The addition of constraints to the design problems can cause the searching performance of the hybrid algorithms.The SPEA using ANN is inferior to the other surrogate-assisted methods or even the original SPEA2 without using a surrogate model because it requires efficient network topology optimisation before being used to estimate functions while in this paper the simple heuristic approach is employed.More investigation on this issue would greatly improve the performance of SPEA-N.In cases of large-scale problems, only SPEA using a quadratic function is superior to the nonsurrogate SPEA.

Conclusions
The use of multiobjective evolutionary algorithm for the design of a torque arm leads to multiple optimum structures for decision making.The hybridisation of SPEA and several surrogate models is developed.The various surrogate models including the quadratic regression technique, the radial basis function interpolation, and the Kriging model can help improving the searching performance of the strength Pareto evolutionary algorithm.For a design problem having mass and natural frequency as objective functions and stress and buckling as constraints, SPEA using RBF is the best method.The SPEA using a quadratic regression is the best for the design case of F2.For the other multiobjective design problems, the best and the second best performers are SPEA using the Kriging model and the quadratic RSM, respectively.Nevertheless, when considering all of the design problems the overall top performer is SPEA2 using the quadratic regression while the close second best uses the Kriging model.Moreover, SPEA using a quadratic response surface model is the only strategy that is useful for solving large-scale problems.The use of ANN for approximating constrained optimisation problems is not effective as a proper network topology optimisation needs to be constructed before entering the surrogate subproblem.It is however possible to enhance the performance of the ANN-based optimiser if an efficient network topology optimisation is performed during the optimisation process.

P
k , f(P k ), and A k = { } Evaluate fitness values of

Figure 2 :
Figure 2: Selection of sampling points for constructing a surrogate model.