Interval multi-objective optimisation of structures using adaptive Kriging approximations

Computers and Structures 119 (2013) 68–84

Contents lists available at SciVerse ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Interval multi-objective optimisation of structures using adaptive Kriging approximations Fangyi Li a, Zhen Luo b,⇑, Jianhua Rong a, Nong Zhang b a b

School of Automotive and Mechanical Engineering, Changsha University of Science and Technology, Changsha 410114, China School of Electrical, Mechanical and Mechatronic Systems, University of Technology, Sydney, NSW 2007, Australia

a r t i c l e

i n f o

Article history: Received 13 April 2012 Accepted 29 December 2012 Available online 10 February 2013 Keywords: Multi-objective optimisation Interval number programming Uncertainty Approximation model

a b s t r a c t This paper proposes an interval uncertain multi-objective optimisation (IUMOO) method for structures with uncertain-but-bounded parameters. An adaptive Kriging model is established to improve the computational efficiency and numerical accuracy in the approximation of design functions. Latin Hypercube Design (LHD) is applied to achieve a set of sampling points both in the design and uncertain spaces for calibrating the Kriging surrogate model. The interval number programming method is used to transform the uncertain optimisation into a corresponding deterministic multi-objective optimisation. Typical numerical examples are used to demonstrate the effectiveness of the proposed methodology. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Multi-objective optimisation of structures has become increasingly important for a range of engineering problems [1,2]. However, most multi-objective optimisation methods (e.g. [3–5]) operate under the assumption that all parameters involved in the design are deterministic. Although deterministic multi-objective optimisation has been widely studied and applied to a number of engineering design problems, many real-world problems are too complex to be defined deterministically because of a lack of information. A variety of uncertainties are inherent in the loading, structural parameters, material properties, dimensional tolerances, boundary conditions, and dimensions due to the complexity associated with engineering problems [6,7]. Such deterministic assumptions may lead to unfeasible designs, as uncertainties can result in significant changes in the structural performance [8]. There is an increasing demand for including various uncertainties into multi-objective design optimisation problems to ensure structural safety, avoid structural failure, and avoid collapse under extreme working conditions. Several methods e.g. [9,10] can be applied to model the uncertainties in the multi-objective optimisation of structures, including the probabilistic method, fuzzy-set method, and convex models. Probabilistic methods have been widely applied to solve a number of uncertainty-type problems. For probabilistic methods, uncertain parameters are typically treated as random variables via the predefined probability distribution functions with the ⇑ Corresponding author. Tel.: +61 2 9514 2994; fax: +61 2 9514 2655. E-mail address: [email protected] (Z. Luo). 0045-7949/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2012.12.028

assumption of possessing sufficient information. For instance, Wang and Yang [11,12] studied the multi-objective structural optimisation of worst-case uncertain anti-optimisations, in which the fat Bezier curve is used to represent the geometry of a structural shape and topology and a local search strategy is integrated with a hybrid genetic algorithm. Limbourg and Kochs [13] proposed a multi-objective optimisation method for a reliability-based design using a feature modelling scheme and an evolutionary algorithm to obtain the Pareto front. Fieldsend and Everson [14] studied an uncertainty multi-objective optimisation method in which probabilistic domination contours were used to represent the estimated Pareto front. However, for most engineering problems, it is difficult to practically obtain sufficient information to determine the probability density functions with precision. Furthermore, Ben-Haim and Elishakoff [15] have shown that even small deviations from the real values may cause relatively large errors in the probability distributions of the feasible region of the design space. As a result, probabilistic methods may experience difficulty when used to predict performance in complex engineering designs. In addition, it is difficult to obtain appropriate membership functions to measure uncertainties in different optimisation problems while using fuzzy set based methods [16,17]. To overcome the limitations associated with probabilistic methods in structural optimisation, non-probabilistic models have been developed as useful supplements (e.g. [18–24]). Among these methods, bounded uncertainty parameters are typically represented using convex models [6,7], such as ellipsoid [21,22] and interval models [18–20,24,25]. In particular, the interval method has attracted significant attention for use in structural optimisation. The interval model measures the uncertainty of parameters

69

F. Li et al. / Computers and Structures 119 (2013) 68–84

without requiring complete information; it only requires a set interval rather than precise probabilistic distribution functions. The determination of the lower and upper bounds of a parameter interval is relatively easier compared to the identification of a precise probability distribution. The interval bounds are defined as the minimum and maximum responses of the objective and constraints with uncertainties. The interval method has been successfully incorporated into the design of structures involving uncertain but bounded parameters ([26–30,20]). However, multi-objective optimisation of structures with interval uncertainty has been seldom studied, though many problems with multiple design requirements are characterised by uncertainties. It is noted that IUMOO problems generally involve a nested double-loop procedure in the design, which makes the computational costs expensive for use in practical engineering problems. Approximation models (meta- or surrogate models) have been widely studied for use in design optimisation [31–33] to improve the computational efficiency in numerical approximations of multi-objective optimisations. This is because surrogate models can provide fast estimates for the objectives and constraints at new design points. Several surrogate models have been developed for implementing the approximation of design functions to improve computational efficiency in the original optimisation, e.g. the response surface method (RSM) [34–36], the radial basis function (RBF) [37,38], and the Kriging model [32,39,40]. It is noted that the Kriging model is more applicable to the optimisation of structures with many design variables [32,39]. Using a constant global model and localised Gaussian correction function, the Kriging model is able to yield a more accurate approximation for the general optimisation of structures. In most IUMOO problems, the surrogate model approximations are not updated during the optimisation e.g. [41], which results in poor numerical accuracy that degrades the optimal solution because the Pareto front is determined by the numerical accuracy of the surrogate model. It is therefore more reasonable to adaptively update the approximation models during optimisation because the accuracy of the approximation is of great importance in finding the Pareto optimums [42–44]. One objective of this work is to combine adaptive numerical techniques with surrogate models to develop an adaptive approximation scheme for IUMMO problems. To the best of the authors’ knowledge, we can find few publications discussing adaptive approximation models for IUMOO problems, even though the approximation in IUMMO problems is very important to the computational efficiency and numerical accuracy. For instance, Yang et al. [44] proposed an adaptive approximation method based on the GA algorithm by sequentially updating the Kriging models. Li et al. [45] studied multi-objective genetic methods with adaptive design and the Kriging model. This paper proposes a new method for multi-objective optimisation of structures involving uncertainty. The proposed IUMOO method focuses on multi-objective optimisation, interval uncertainty and surrogate models. The multi-objective optimisation technique is applied to describe multiple design requirements for engineering problems. The interval model is employed to measure the uncertainty in uncertain but bounded parameters without complete information. An adaptive surrogate model is developed to improve the computational efficiency of the IUMOO problems, in which the Kriging model with the LHD is used to generate accurate approximations of the uncertainty objectives and constraints. 2. Uncertain multi-objective optimisation Generally, an interval uncertain multi-objective optimisation (IUMOO) problem can be defined as

8 min : fi ðx; aÞ; i ¼ 1; 2; . . . ; k; > > > x 8 h i > > > > > g j ðx; aÞ 6 v Ij ¼ v Lj ; v Rj ; j ¼ 1; 2; . . . ; m; > > < > > > < a 2 aI ¼ ½aL ; aR ; > s:t: > h i > > > I L R > > > > > > ap 2 ap ¼ ap ; ap ; p ¼ 1; 2; . . . ; q; > > : : xl 6 x 6 xr ;

ð1Þ

where x is the n-dimension design vector and xl and xr are the left and right lateral constraints denoting the allowable minimum and maximum vectors of x, respectively. The variables fi (i = 1, 2, . . . , k) and gj (j = 1, 2, . . . , m) are the objective and constraint functions, respectively. The superscript ‘‘I’’ denotes the interval, and the subscripts ‘‘L’’ and ‘‘R’’ represent the lower and upper bounds of the interval, respectively. The variable a is the q-dimension uncertainty vector modelled by the interval vector aI = [aL, aR], with all its components being interval numbers. The variable aIp is an interval number, which is actually a bounded set of real numbers between the bounds aLp and aRp . The variable v Ij is the interval number of the jth constraint and m is the number of constraints. 2.1. Treatments of uncertain objective functions In this section, the order relation ‘‘6mw’’ [18–20] from interval mathematics is used to describe the decision maker’s preference as to the midpoint value m and the radius w of the interval number. For the minimisation problem, the objective of the interval function is to minimise m and w of the uncertainty 0objective function. Each objective in Eq. (1) can be transformed into a deterministic two-objective optimisation via ‘‘6mw’’ as

8 min : ðmðfi ðx; aÞÞ; wðfi ðx; aÞÞÞ; > > > x >   < where; mðfi ðx; aÞÞ ¼ 12 fiL ðxÞ þ fiR ðxÞ ;   > > wðfi ðx; aÞÞ ¼ 12 fiR ðxÞ  fiL ðxÞ ; > > : where; i ¼ 1; 2; . . . ; k:

ð2Þ

The minimisation of the derivation leads to a decrease in the variance of the objective function due to the uncertainty, which indicates that the optimum value will render the uncertain objective function insensitive to fluctuations caused by the uncertainty. In the inner layer, for a specific x0 , the following optimisation problems must be solved

fiL ðx0 Þ ¼ minq fi ðx0 ; aÞ; a2C

and fiR ðx0 Þ ¼ max fi ðx0 ; aÞ; q a2C

n   where i = 1, 2, . . . , k and Cq ¼ aaLp 6 ap 6 aRp ;

ð3Þ

o p ¼ 1; 2; . . . ; q .

A linear combination is then used to transform Eq. (2) into the following assessment function fdi min : fdi ðx; aÞ ¼ ð1  bÞ½mðfi ðx; aÞÞ þ n=/ þ b½wðfi ðx; aÞÞ þ n=w; x

ð4Þ

where 0 6 b 6 1 is a weighting factor, in which different values are applied to represent different preferences to different objectives. The variable n is a parameter used to ensure the terms [m(f(x, a)) + n] and [w(f(x, a)) + n] are non-negative. The variables / and w are normalisation factors [19,20] mathematically determined from

/ ¼ minn ½mðfi ðx; aÞÞ þ n and w ¼ minn ½wðfi ðx; aÞÞ þ n; X2X

X2X

ð4aÞ

where x 2 Xn is an n-dimension decision vector and Xn is the set of its range. For practical designs, / and w can be determined according to the order of magnitude of each individual objective function.

70

F. Li et al. / Computers and Structures 119 (2013) 68–84

2.2. Treatments of uncertain constraints

3. Adaptive approximation models

In reference [20], a definition of the satisfactory degree of interval is presented to describe the degree that one interval is larger or smaller than another. In terms of the relations of interval numbers between C and D, the probability that P(CI 6 DI) P ki, ki 2 [0, 1] [20] is used to denote the possibility that the degree of C is smaller than that of D and to facilitate the treatment of the inequality constraints. Uncertain constraints are represented by a pre-determined satisfaction degree; the uncertain inequalities can be changed to deterministic constraints by satisfying a specified possible degree level. Therefore, according to the aforementioned satisfactory degree method, all of the uncertain constraints in Eq. (1) can be transformed into deterministic ones by

In this section, an adaptive approximation method is developed to solve the optimisation problem given in Eq. (7). In each step, the Kriging model [32,47] is used to construct accurate global approximations for the uncertain objectives and constraints to facilitate the IUMOO. A modified adaptive scheme is included to update the Kriging model until the approximation optimisation converges. The Kriging model is well documented in references [32,47]. A basic description of the Kriging method is as follows.

 8  I I > P C 6 D P kj ; j ¼ 1; 2; . . . ; m; > j j > > < h i I where; C j ¼ g Lj ðxÞ; g Rj ðxÞ ; > > h i > > : DIj ¼ v Lj ; v Rj ;

ð5Þ

g Lj ðx0 Þ ¼ minq g j ðx0 ; aÞ and g Rj ðx0 Þ ¼ max g j ðx0 ; aÞ: q a2C

ð6Þ

So far, the uncertain objective functions and constraints have been transformed into deterministic constraints in terms of the order relation of the interval and the satisfaction of the possibility degree, respectively. Then, the penalty function method (e.g. [19,20]) is used to change the deterministic constrained objective optimisation into an unconstrained optimisation using x

i ¼ 1; 2; . . . ; k ðxl 6 x 6 xr Þ;

ð7Þ

where f(X) = [f1(X), f2(X), . . . , fm(X)]T is the known regression function consisting of m functions, and b = [b1, b2, . . . , bm]T is the related regression coefficients. The term bf(X) globally models the design space via polynomial-based functions similar to the polynomial response surface. In most cases, this term is denoted by the constant b. The variable Z(X) is a stationary stochastic process with a zero mean, a variance of r2 and a non-zero covariance, which defines the localised deviations. ~ðXÞ for the response y(X) at an unThe estimated predicator, y tried X is formulated from

~ ~ þ r T ðXÞM1 ðY  EbÞ: ~ðXÞ ¼ b y

rT ðXÞ ¼ fGðX; X1 Þ; GðX; X2 Þ; . . . ; GðX; Xn ÞgT :

ð1  bÞ½mðfi ðx; aÞÞ þ n b½wðfi ðx; aÞÞ þ nÞ þ / w n h   i X þ r u P C Ij 6 DIj  kj ;

ð11Þ

~ 2 for the stationary stochastic process and the estimated variance r Z(X) is written as

ð7bÞ

r~ 2 ¼

  1 L 1 fi ðxÞ þ fiR ðxÞ and wðfi ðx; aÞÞ ¼ fiR ðxÞ  fiL ðxÞ ; 2 2 ð7cÞ

/ ¼ minðmðfi ðx; aÞÞ þ nÞ and w ¼ minðwðfi ðx; aÞÞ þ nÞ; h   i n h    io2 u P C Ij 6 DIj  kj ¼ max 0;  P C Ij 6 DIj  kj :

ð10Þ

ð7aÞ

i¼1

mðfi ðx; aÞÞ ¼

ð9Þ

Here, the maximum likelihood estimates (MLE) method [32,47] is used to determine the unknown correction parameters, which are used to define the correction function G(Xi, Xj). ~ are obtained via the least squares The estimated coefficients, b method [32,47] as

~ ¼ ðET M1 EÞ1 ðET M1 EÞ b

where

fpi ¼

ð8Þ

In the above equation, rT(X) is the correction vector of length n between an untried X and all the sample points, which are given as the following vector

2.3. Multi-objective optimisation formulation

min : fpi ðx; aÞ;

Assume that X = [X1, X2, . .. , Xn]T is the data collected from n T sampling points, where Xi ¼ x1i ; x2i ; . . . ; xpi (i = 1, . . . , n). The variable p is the number of design variables. The associated response is y(Xi). The Kriging model combines a global model and localised deviations to define the unknown function y(X) as

yðXÞ ¼ ðbfðXÞÞ þ ZðXÞ;

where PðC Ij 6 DIj Þðj ¼ 1; 2; . . . ; mÞ is the satisfactory degree of constraints and a fuzzy definition of the possibility that the interval C Ij is smaller than interval DIj . The function kj (j = 1, 2, . . . , m) is the pre-determined satisfactory degree level of the constraints. The variable C Ij is the interval of the jth constraint at x as a result of the uncertainty. For a specific x0 , the lower and upper bounds of the intervals g Lj ðxÞ and g Rj ðxÞ are solved, respectively, using an optimiser such as Sequential Quadratic Programming (SQP) [46] for a2C

3.1. Kriging model for approximations

ð7dÞ ð7eÞ

Here, Xn is the design space spanned by all the decision vectors x, and C q is the space of the vector a. The factor r is required for a larger constant, which can be determined in terms of numerical testing experience. Eq. (7) can be solved by traditional multi-objective optimisation methods (see Fig. 1). This is a typical double-loop nest optimisation, so the computational cost is extremely expensive.

o 1n ðY  EbÞT M1 ðY  EbÞ n

ð12Þ

for

Y ¼ ½yðX1 Þ; yðX2 Þ; . . . ; yðXn ÞT ;

ð13Þ

where Y is a column vector of length n, which contains the sample values of the response, M is the matrix consisting of the correlation function G(Xi, Xj), and E is a unit column vector of length n. For more details on determining the previous terms, please refer to [32,47]. The resulting Kriging model is used to develop the approximation models ~f ðx; aÞ and g~ðx; aÞ for the objective functions and constraints, both in the design space and the uncertainty space, respectively. During the approximation approach, the sampling data points used to construct the Kriging model can be obtained using the Latin Hypercube Design (LHD) [42] method.

71

F. Li et al. / Computers and Structures 119 (2013) 68–84

Outer layer multi-objective operator Design space The trial design vector

Uncertainty space Inner layer optimization

Inner layer optimization

Call Expensive high fidelity computational models

Call

Interval of the objective function

Intervals of the constraints

Assessment functions

Satisfactory degrees of the constraints

N Penalty functions

Stopping criterion Y Pareto optimal set

Fig. 1. Nesting optimisation based on the computational models.

and

3.2. Adaptive approximation model In the optimisation process, a series of approximations for the IUMOO problems are created. For instance, the optimisation model at the sth iterative step has the general form

8 > min : ~f i ðx; aÞ; i ¼ 1; 2; . . . ; k; > > x > > h i 8 > > > > > g~j ðx; aÞ 6 v Ij ¼ v Lj ; v Rj ; j ¼ 1; 2; . . . ; m; > < > > > < I L R > > s:t: a 2 a ¼ ½ah ; a  i > > > > > > > ap 2 aIp ¼ aLp ; aRp ; p ¼ 1; 2; . . . ; q > > > > > : : xl 6 x 6 xr ;

ð14Þ

where ~f i and g~j are approximations of the uncertain objectives and constraints using Kriging model, respectively. It can be seen that these variables are both explicit functions with respect to x and a. The estimated uncertain multi-objective optimisation can be transformed into a deterministic problem using

min : ~f pi ðx; aÞðxl 6 x 6 xr Þ; x

ð15Þ

where

~ ~ ~f p ¼ ð1  bÞðmðf i ðx; aÞÞ þ nÞ þ bðwðf i ðx; aÞÞ þ nÞ i / w n     X e I 6 DI  kj þr u P C j j i¼1

ð15aÞ

h i e I ¼ g~L ðxÞ; g~R ðxÞ ¼ ming~j ðx; aÞ; maxg~j ðx; aÞ ; C j j j q q a2C

a2C

ð15bÞ

e I is the interval of the where ~f pi is the Kriging penalty function and C j jth approximated constraint. Inside the penalty term, we obtain the Pareto set for the sth approximation optimisation problem by using Eq. (15). The flowchart in Fig. 2 shows the management of the interaction between the IUMOO and the Kriging model, in which IUMOO predicts the optimum point sets using the approximation models. From the set of optimum points, one or more points will be selected for re-analysis and further adapted to the approximation models of the real analysis according to the space-filling design criteria [48]. An adaptive approximation technique is then developed to manage the approximation models and improve the numerical accuracy of the Kriging model for the IUMOO problem. The major numerical procedures used in the adaptive Kriging approximation are shown in Fig. 3, and the relevant explanations are as follows: n (1) Eval ¼ ðx; aÞjxi;l 6 xi 6 xi;r ; aLp 6 ap 6 aRp ; i ¼ 1; 2; . . . ; n; p ¼ 1; 2; . . . ; qg. Eval is the sampling space, which is a combination of the current design and uncertainty spaces. In the other words, (x, a) corresponds to a sample point from the sampling space for the simulation, and the LHD is used to sample the points in Eval. The number of samples is defined as Ns

72

F. Li et al. / Computers and Structures 119 (2013) 68–84

Sampling

Approximation models

Interval Uncertain Multi-Objective Optimization (IUMOO)

Update approximation models including new points

Predicted point sets Evaluate Selected points

Select a few points

Expensive high fidelity computational models

Selector

Fig. 2. Interaction between the IUMOO and approximation models.

and the iterative step is s = 1. Costly, high fidelity computational models are then used to calculate the objective and constraint functions for all the sample points. (2) The response values are calculated for the all sample points by evaluating the true values of the objective and constraint functions. The approximations ~f ðx; aÞ and g~ðx; aÞ are constructed according to the Kriging model based on the current solution to Eval, where x and a are used as the input variables to the Kriging model. (3) The IUMOO problem is formulated based on the Kriging model of the objective and constraints using Eq. (14), and ðsÞ then solved for the estimated Pareto set P d0 . This should be performed by only including x(s) (Eq. (15)) through the approximation of the IUMOO problem. (4) The same approximation models from Step (2) are used for ðsÞ each Pareto optimal point x(s) in P d0 . The interval L ðsÞ ~R ðsÞ ~ ½f i ðx Þ; f i ðx Þ for the objective function of the approximation model is evaluated via the corresponding uncertain parameters aLf and aRf , which are obtained from   ~f L ðxðsÞ Þ ¼ ~f L xðsÞ ; aL and ~f Ri ðxðsÞ Þ ¼ ~f Ri ðxðsÞ ; aRf Þ ði ¼ 1; 2; . . . ; kÞ: i i f ð16Þ

Similarly, the interval of the estimated constraint function is obtained via aLg and aRg from     and g~Rj ðxðsÞ Þ ¼ g~Rj xðsÞ ; aRg ; j ¼ 1; 2; . . . ; m: g~Lj ðxðsÞ Þ ¼ g~Lj xðsÞ ; aLg ð17Þ ðsÞ P a0

Thus, the uncertainty set in the uncertain vector a can be calcuðsÞ lated in terms of each design vector x(s) in the Pareto set P d0 based on the above approximation models of the objective and the constraints. ðsÞ ðsÞ ðsÞ (5) The set of design points P 0 is obtained from P d0 [ P a0 , which ðsÞ

ðsÞ

ðsÞ

includes P d0 and P a0 , where ðxðsÞ ; aLf Þ 2 P 0 ; ðsÞ

ðsÞ

ðsÞ

ðxðsÞ ; aRf Þ 2

P 0 ; ðxðsÞ ; aLg Þ 2 P 0 and ðxðsÞ ; aRg Þ 2 P 0 , respectively. ðsÞ (6) Using the design point set P 0 , the distances between each point of the design set and all the other samples are calculated to obtain the minimum distance to each point. After this operation is completed, the isolated points, nadd, the

points farthest from the minimum distance points, are solved. These isolated points are added to the Padd and ðsÞ P add  P 0 variables belonging to the design so that the uncertainty spaces can be obtained. In general, nadd is chosen to be two or three. The maximum distance, dmax, from the design point set is defined as    o     n  ðsÞ  ðsÞ   ðsÞ dmax ¼ max min P0 1  Evalj ; . . . ; P 0 i  Evalj  i 6 n P0 ; j 6 N s ; ð18Þ

  ðsÞ ðsÞ represents the number of P0 values and Ns is the where n P0 number of Eval samples. (7) According to the number of sample points in Padd, implement the high fidelity computational models from Step (6). This is performed by computing the responses of the actual objective ½fiL ðxðsÞ Þ; fiR ðxðsÞ Þ and the constraint ½g Lj ðxðsÞ Þ; g Rj ðxðsÞ Þ from the newly selected experimental points Padd via the computational model, where i = 1, 2, . . . , k and j = 1, 2, . . . , m. (8) Add Padd to the existing sample points in Eval = Eval [ Padd (e.g. the previous experiments Eval plus the new experiment Padd) so that the previous Kriging models of the objective functions ~f ðx; aÞ and the constraints g~ðx; aÞ are updated correspondingly. (9) Obtain a new Pareto set Pdnew, including the design vector x ~ðx; aÞ, and from the updated approximation of ~f ðx; aÞ and g then calculate the number of the Pareto set Pd0(s) and the Pareto set Pdnew, respectively. (10) Check the convergence of the solution; if the difference between the number of the previous Pareto set and the new one is smaller than the given criterion d, stop the calculation. Otherwise, set s = s + 1 and return to Step 3 to start a new iteration, and insert Pdnew into the initial generation of the IUMOO. 3.2.1. Further explanations of some key procedures are provided to assist in the understanding of the proposed adaptive Kriging approximation In Step (1), a basic requirement for an initial design is that the entire input space should be sampled uniformly so that all regions have an equal chance to be searched. Additionally, the number of points in the initial design is dependent on the dimensions of the

F. Li et al. / Computers and Structures 119 (2013) 68–84

73

Fig. 3. Flowchart of the present method.

input space; namely, the higher the dimension of the input space, the more sampling points need to be included in the initial design. In the proposed approach, it is assumed that the initial design generated by the LHD method provides sufficient information on the design and the response spaces. The aim of this step is to capture the major characteristics of the practical objectives and constraints. In Step (4), at each calculation, the precision in the approximation of the response bounds of the objective functions and constraints is gradually improved. However, the major focus is to improve the precision of the local critical regions, rather than the whole approximation space. The optimisation problem in Eqs. (16) and (17) is then solved (see Section 3.3). When the inner layer ðsÞ problem is solved, the P a0 for all the Pareto optimal points are obtained without additional computational cost (through high-fidelity computational models) as a result of using the same approximation models from Step (2).

In Step (6), to avoid the sample points clustering in one portion, a new experiment is selected according to the information from both the existing experiments and the Kriging model so that the Kriging model accuracy can be improved in the inaccurate regions of the predicted response space. Because the initial points may be far away from the Pareto frontier or insufficiently sample the design space, the initial Kriging model will be inaccurate. By adding new points to the existing samples to update the approximation model, the accuracy of the Kriging model is adaptively improved when the algorithm evolves. This procedure can help identify the regions that should be more heavily sampled, such as the uncertainty boundary and the design regions. However, for conventional methods based on surrogate models, the entire space will be sampled uniformly (same sample density) to construct the high quality approximation model so that all regions have an equal chance to be sampled. In other words, the samples are distributed uniformly in the entire hybrid

74

F. Li et al. / Computers and Structures 119 (2013) 68–84

space created by the LHD method, which lowers the accuracy of the Kriging model. In Step (7), the number of inputs to each approximation model is the sum of the number of design variables and uncertain parameters. For example, let x be an optimal solution to the IUMOO probðsÞ lem and a be the uncertainty set P a0 obtained from Eqs. (16) and (17). The value of [x, a] corresponds to a new online sample point. The actual objective and constraint functions are used to evaluate the responses for the new sample point. After all the new sample points are collected in the current iteration, they are combined with the sample points from previous iterations to update the approximation model for the objective and constraint functions. In Step (8), the approximation space remains unchanged, and the previous samples generated by the LHD are kept in the optimisation. The samples can be added at each iterative step, which includes the initial and previous samples. To avoid changing the sampling space, re-sampling is performed. Compared with most conventional methods, this method requires fewer sample points, and improves the efficiency of the optimisation. 3.3. Numerical solution of the approximation optimisation problem The solution process of the approximation optimisation is shown in Fig. 4, and includes a two-layer nesting optimisation. At each step, the Non-dominated Sorting Genetic Algorithm II (NSGA-II) [1,49] and Sequential Quadratic Programming (SQP) [46] are adopted in the outer and inner layers, respectively. Using NSGA-II [49] in the outer layer, a number of individuals are generated with each individual chromosome denoting a candi-

date decision vector x. The NSGA-II is a computationally fast, nondominated sorting scheme for multi-objective optimisation, including non-dominated sorting for fitness assignments, a fast non-dominated sorting technique, and a crowding distance to rank for selecting population fronts. The algorithm then applies the crossover and mutation operators to combine the current population with the next generation offspring. Finally, the best individuals, in terms of non-dominance and diversity are selected as the solutions. All individuals not dominated by any other individuals are assigned a front number of 1. All individuals dominated solely by individuals with the front number 1 are assigned the front number 2 and so on. Selection is made by a tournament between two individuals. The individual with the lowest front number is selected if the two individuals have different front numbers. The individual with the highest crowding distance is selected if they originate from the same front; i.e., a higher fitness is assigned to individuals located on a sparsely populated part of the front. A certain amount of N parents and N new individuals (offspring) are generated in each iteration, and both parents and offspring compete with each other for inclusion in the next iteration. For each x, the intervals of the objective function and the constraints can be calculated in the inner loop using SQP [46], in which x is specified as a constant and a is the decision vector. The SQP is a nonlinear programming method, which starts from a single searching point to find a solution with the gradient information. It outperforms most nonlinear programming methods in terms of efficiency, accuracy, and percentage of successful solutions over a large number of test problems. To compute each interval of the objective function and constraints, the outer NSGA-II calls the

Outer layer multi-objective GA (NSGA- ) Design space Design vector of the kth iterative step

Uncertainty space

Inner layer SQP

Inner layer SQP Call

Call

Approximation models of the uncertain objective function and constraints k+1

Interval of the objective functions

Intervals of the constraints Sampling points fromLHD

Assessment functions

Expensive high fidelity computational models

Satisfactory degrees of the constraints Current design space and uncertainty space Approximate penalty functions

N

Stopping criterion

Y Pareto optimal set Fig. 4. Approximation of the IUMOO based on the approximation models.

F. Li et al. / Computers and Structures 119 (2013) 68–84

75

inner SQP twice to compute the lower and upper bounds of the objective function and constraints, based on the Kriging model for each individual. Then, based on the intervals, the midpoint value and radius of the objective function, and the possibility degrees of the constraints are calculated to determine the estimated penalty function. Finally, the Pareto set is obtained.

the sample points found via Eq. (18). The design points at maximum distance are selected from the most closely located sample points and added to the sample point sets to update the approximation models.

3.4. The evaluation of design points set P0 and isolated points

In this section, we use two typical numerical test functions from [1,49,50] to show the effectiveness of the proposed method on the IUMOO problems, including numerical accuracy and computational efficiency. The readers may refer to the related papers for details [1,49,50]. The first test function (1) is expressed as

ðsÞ

The process of obtaining P0 and calculating the isolated points is shown in Fig. 5. The Pareto set Pd0(s), for which design vector x(s) is included and obtained by approximating the IUMOO problem. By specifying each design vector x(s) from Pd0(s), the approximations of the objective and the constraints are constructed with respect to the uncertainty vector a. The vector sets aLf and aRf , as well as aLg and aRg are obtained via the SQP to call the approximation models in the uncertain space. ðaLf ; aRf Þ and ðaLg ; aRg Þ are the uncertain vectors of the intervals of the objectives and constraints from Eqs. (16) and (17), respectively. In this way, the both uncertain set Pa0(s) consisting of ðaLf ; aRf Þ, ðsÞ

and ðaLg ; aRg Þ and the design point set P 0 , created from ðsÞ

½x

; aLf ;

ðsÞ

½x

; aRf ;

ðsÞ

½x

; aLg 

ðsÞ

and ½x

; aRg 

can be obtained. The sam-

pling space Eval is a combination space made up of the points [x, a]. The closest sample points can be found by computing the ðsÞ

Euclidean distance between the design points from P 0 and all

4. Numerical tests

8 min : ðf1 ðx; aÞ; f2 ðx; aÞÞ; > > > x > > > 2 2 > > f1 ðx; aÞ ¼ a1 ðx1 þ x2  7:5Þ þ a22 ðx2  x1 þ 3Þ =4; > > > 2 2 > 2 3 > f2 ðx; aÞ ¼ a1 ðx1  1Þ =4 þ a2 ðx2  4Þ =2; < 8 : > g 1 ðx; aÞ ¼ a21 ðx1  2Þ3 =2 þ a2 x2  2:5 6 ½0; 0:3; > > > > > > < > 2 > > s:t: g 2 ðx; aÞ ¼ a31 x2 þ a22 x1  3:85  8a22 ðx2  x1 þ 0:65Þ 6 ½0; 0:3; > > > > > > a1 2 ½0:9; 1:1; a2 2 ½0:9; 1:1; > > > > > : : 0 6 x1 6 5; 0 6 x2 6 3: ð19Þ The second test function (2) is a more complex three-dimensional test problem, defined as

ðsÞ

Fig. 5. Flowchart showing the process of obtaining P0 and calculating the isolated points.

76

F. Li et al. / Computers and Structures 119 (2013) 68–84

2.2

2.2 True Pareto point

2

2

1.8

1.8

1.6

1.6

1.4

1.4

fp2

fp2

Obtained Pareto point based on adaptive approximiation model

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

4

5

6

7

8

9

10

True Pareto point Obtained Pareto point based on adaptive approximiation model

0.2

11

4

5

6

7

8

9

fp1

fp1

(a) 30 sample points

(b) 40 sample points

10

11

2.2 2

True Pareto point Obtained Pareto point based on adaptive approximiation model

1.8 1.6

fp2

1.4 1.2 1 0.8 0.6 0.4 0.2

4

5

6

7

8

9

10

11

fp1

(c) 50 sample points Fig. 6. The Pareto points from test function 1 using different initial sample points: (a) 30 sample points, (b) 40 sample points, (c) 50 sample points.

8 min : ðf1 ðx; aÞ; f2 ðx; aÞÞ; > > x > > ! > > 3  2 > X > ai > > ; f1 ðx; aÞ ¼ 1  exp  xi  pffiffi3 > > > < i¼1 ! 3  2 X > > ; xi þ paiffiffi3 > > f2 ðx; aÞ ¼ 1  exp  > > i¼1 > >

> > > 4 6 x1 ; x2 ; x3 6 4; > > : s:t: a1 2 ½1:1; 1:3; a2 2 ½1:1; 1:3;

Table 1 Details of NSGA-II specific parameters used.

ð20Þ

a3 2 ½1:1; 1:3:

In these two test functions, the parameter n is specified to be 0, and the normalisation factors / and w are both chosen to be one. The possibility satisfaction level k for the constraints is specified as 0.8, and the penalty factor r is defined as 1000. The weighting factor b is set to be 0.5. The parameters for the NSGA-II optimiser are listed in Table 1. Because the NSGA-II is a well-studied optimisation algorithm, the description of the determination of the proper parameters for a NSGA-II operation is left to the reader [1,49]. However, for different designs, numerical test experience is required to choose an appropriate set of optimisation parameters. Two design points are selected to update the Kriging model with the objectives and constraints at each iteration. The convergence criterion d is set at 2%, which are different solution numbers for two adjacent Pareto sets.

GA parameter name

Value

Population size Number of generation Probability of crossover Distribution index for crossover Distribution index for mutation

50 50 0.9 20 20

Table 2 The results from the two test functions. Test functions

Initial sample points

1

30 40 50

2

50 60 70

Total number of points evaluations

Numerical results

82 72 58

112 112 108

Fig. 7(a) Fig. 7(b) Fig. 7(c)

202 188 196

252 248 266

Fig. 9(a) Fig. 9(b) Fig. 9(c)

Number of adding points evaluations

77

F. Li et al. / Computers and Structures 119 (2013) 68–84

3

2.2 Number of points evaluations = 42 Number of points evaluations = 54 Number of points evaluations = 78 Number of points evaluations = 102 Number of points evaluations = 112 (Final Result)

2.5

2 Number of points evaluations = 50 Number of points evaluations = 64 Number of points evaluations = 88 Number of points evaluations = 100 Number of points evaluations = 112 (Final result)

1.8 1.6

2

fp2

fp2

1.4 1.5

1.2 1

1

0.8 0.6

0.5

0.4 0

4

5

6

7

8

9

10

0.2

11

4

5

6

7

8

9

fp1

fp1

(a) 30 sample points

(b) 40 sample points

10

11

2.5 Number of points evaluations = 58 Number of points evaluations = 70 Number of points evaluations = 82 Number of points evaluations = 94 Number of points evaluations = 108 (Final result)

2

fp2

1.5

1

0.5

0

4

5

6

7

8

9

10

11

12

13

fp1

(c) 50 sample points Fig. 7. The convergence histories of the Pareto set in test function 1 using different sample points: (a) 30 sample points, (b) 40 sample points, (c) 50 sample points.

4.1. Optimisation convergence In the above numerical test functions, the objective functions and constraints are constructed by the Kriging models. The relevant optimisation results are given in Table 2. Test function 1 has two design variables, two uncertainty parameters, and two discontinuous regions in the Pareto front. Fig. 6(a)–(c) show the Pareto points of function 1 obtained from three different initial sample points: 30, 40 and 50, respectively. The obtained Pareto set is similar to the true Pareto set, as shown in Fig. 6. The convergence of the Pareto set with different initial sample points is shown in Fig. 7. The obtained Pareto set is determined to be acceptable. In the beginning, a poor quality Pareto set may be obtained, which shows that the initial approximation model has a low accuracy in the region marked with arrows. However, with the evolvement of the optimisation, the Kriging models are updated with respect to the design and uncertain spaces, until the obtained Pareto set converges through the addition of more sampling points. These results show the effectiveness of the propose method. Test function 2 has three design variables and three uncertain parameters. Fig. 8 displays the results of from test function 2, which has two higher order objective functions with three different initial sample points (50, 60 and 70). The estimated Pareto set of the present method is almost the same as the true Pareto set in Fig. 8. The convergence of the Pareto set for test function 2 with different initial sample points is shown in Fig. 9. The obtained

Pareto set is determined to be acceptable. Initially, the Pareto set is located in a small region of the entire design space. It should be noted that the poor approximation accuracy in the region, marked with arrows, does not affect the optimal result. During the optimisation, the Kriging models are updated by adding new points during each iteration, and the obtained Pareto set is sequentially convergent. In such a six dimensional hybrid space, it is difficult to construct accurate approximation models with a limited number of sampling points and a uniform distribution. However, a satisfactory result is obtained with the adaptive approximation method. 4.2. Optimisation efficiency In the two test functions, it is assumed that the objective functions in Eqs. (19) and (20) are computationally intensive simulation models, and the evaluation number for the objective function is our major concern. The Pareto set of accurate solutions is obtained from the interval number programming method from the actual models, in which NSGA-II and SQP are employed to solve the nested double-loop optimisation [19,20]. The maximum number of generations, and individuals per generation (the population size) involved in the outer NSGA-II are both specified as 50. Thus, a number of 2500 individuals are generated, and each individual represents a trial design vector. The SQP in the inner layer is called twice to compute the upper/lower bounds, and each call needs

78

F. Li et al. / Computers and Structures 119 (2013) 68–84

0.7

0.7 True Pareto point Obtained Pareto point based on adaptive approximiation model

0.6

True Pareto point Obtained Pareto point based on adaptive approximiation model

0.6 0.5

0.4

0.4

fd2

fd2

0.5

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0

0.7

0

0.1

0.2

0.3

0.4

0.5

fd1

fd1

(a) 50 sample points

(b) 60 sample points

0.6

0.7

0.7 True Pareto point Obtained Pareto point based on adaptive approximiation model

0.6 0.5

fd2

0.4 0.3 0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

fd1

(c) 70 sample points Fig. 8. The Pareto points from test function 2 using different initial sample points: (a) 50 sample points, (b) 60 sample points, (c) 70 sample points.

10 evaluations to solve for the objective function. Thus, the total evaluation number for the actual objective function is 2500  10  2 = 5.0  104. However, with the present IUMOO method, the total evaluation number for the actual objective function of test function 1 is 112, 112 and 108 for the different initial sample points, respectively, as shown in Table 3. The total evaluation number for the actual objective function of test function 2 is 252, 248 and 266 for the different initial sample points, respectively. It is noted that the additional points are used to increase the density of the number of local evaluation points to improve the approximation accuracy, and the total number of point evaluations is not related to the evaluation of the approximation models. As a result, the total number of point evaluations is many times less than the conventional nested optimisation for test functions 1 and 2, respectively. From Table 4, it is seen that the numerical accuracy of the proposed method is better than the conventional deterministic multi-objective optimisation using the RSM approximations [5], as the Kriging model and interval uncertainty are considered in this study. It is noted that the present IUMOO method still includes nesting the double-loop optimisations related to the interval number programming, which increases the computational costs due to the interval uncertainty. However, in contrast to the RSM method, the Kriging model can provide better approximation accuracy, which is of particularly importance for multi-objective optimisation with interval uncertainties taken from searching all the local

optimums in order to construct the Pareto front. Both the RSM and the Kriging model have advantages and disadvantages. However, the Kriging model is more preferable for general structural optimisation problems, as noted in the references [32,47]. 4.3. Analysis of the uncertainty level This section focuses on how different uncertainty levels (e.g. 20%, 30%, 40% and 50%) influence the optimisation results (all other parameters and variables remain constant). The Kriging model and the LHD method are employed to produce approximations of the uncertain objectives and the constraints. With the proposed method, the optimisation results of the Pareto solutions from the uncertainty levels are shown in Figs. 10 and 11. The deviations in the optimal results corresponding to the 20% and 30% uncertainty levels are close to the optimum denoted by ‘‘True Pareto point’’, which indicates better numerical accuracy. When the uncertainty level is changed to 40%, the deviation to the optimised front of the ‘‘True Pareto point’’ front becomes larger. For the 50% uncertainty level, only a portion of the Pareto solutions can be found. Therefore, the proposed method has better numerical accuracy for a reasonable uncertainty level. During uncertain optimisation, the surrogate model is established with a combination of the current design and uncertain spaces. For larger uncertainty levels, the numerical accuracy of the Kriging model will decrease, due to the inaccurate evaluation of the penalty

79

F. Li et al. / Computers and Structures 119 (2013) 68–84

0.6

0.7

0.5

0.6 0.5

0.4

fd2

fd2

0.4 0.3 Number of points evaluations = 126 Number of points evaluations = 146 Number of points evaluations = 170 Number of points evaluations = 210 Number of points evaluations = 252 (Final result)

0.2

0.3

0.1 0

Number of points evaluations = 116 Number of points evaluations = 152 Number of points evaluations = 172 Number of points evaluations = 218 Number of points evaluations = 248 (Final result)

0.2 0.1

0

0.1

0.2

0.3

0.4

0.5

0

0.6

0

0.1

0.2

0.3

0.4

fd1

fd1

(a) 50 sample points

(b) 60 sample points

0.5

0.6

0.6 0.5

fd2

0.4 0.3 Number of points evaluations = 118 Number of points evaluations = 158 Number of points evaluations = 198 Number of points evaluations = 250 Number of points evaluations = 266 (Final result)

0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

fd1

(c) 70 sample points Fig. 9. The convergence histories of the Pareto set in test function 2 using different sample points: (a) 50 sample points, (b) 60 sample points, (c) 70 sample points.

function caused by the inaccurate bounds of the objective and constraint functions. In this case, with respect to the IUMOO method, it is necessary to determine reasonable uncertainty levels corresponding to the design to avoid the degradation in the accuracy of the numerical approximations. 4.4. Effectiveness of the penalty function method This section focuses on the effectiveness of the penalty function included in the optimisation formulation. The results for the five typical points with the Kriging model for uncertainty levels of 10% and 20% are given in Tables 5 and 6, respectively. From these results, it is found that the satisfaction degrees of the constraint based on the Kriging model satisfy a possibility degree of 80%, which is close to the possibility degree of the original constraint without  using the Kriging model. The lower and upper bounds ~L1 ; g ~R1 of the constraint obtained using the approximation model g     are also close to the bounds of g L1 ; g R1 and g L2 ; g R2 . Furthermore, the penalty function value is close to the accurate one. These results show that the penalty function used in the proposed method is effective, and the solution is feasible for a given Pareto design.

spond to one Pareto front. Therefore, different weighting factors will lead to different Pareto fronts via the proposed method. The proposed IUMOO method leads to a set of Pareto fronts corresponding to different weighting factors. We further describe this concept based on test function 1. For the first test numerical case, the optimisation problem is given as

min : ~f pi ðx; aÞ; x

i ¼ 1; 2;

ð21Þ

where

~ ~ ~f p ¼ ð1  bÞðmðf i ðx; aÞÞ þ nÞ þ bðwðf i ðx; aÞÞ þ nÞ i / w n h   i X eI 6 D e I  kj : þr u P C j j

ð21aÞ

i¼1

If the weighting factor is set to b = 0.5 (as previously used), the optimisation problem changes to n h   i ~ ~ X ~f p ¼ ð1  0:5Þðmðf 1 ðx;aÞÞ þ nÞ þ 0:5ðwðf 1 ðx;aÞÞ þ nÞ þ r u P C eI 6 D e I  kj ; ð21bÞ j j 1 / w i¼1

4.5. Analysis of Pareto fronts

n h   i ~ ~ X ~f p ¼ ð1  0:5Þðmðf 2 ðx;aÞÞ þ nÞ þ 0:5ðwðf 2 ðx;aÞÞ þ nÞ þ r u P C eI 6 D e I  kj ; j j 2 / w i¼1

Compared to conventional multi-objective optimisation formulations using the linear combination method, the work described herein uses individual weighting factors (e.g. b = 0.5) that corre-

where ~f 1 ðx; aÞ and ~f 2 ðx; aÞ are the approximation functions of f1(x, a) and f2(x,a), which can be obtained using the Kriging model-based approximation.

ð21cÞ

80

F. Li et al. / Computers and Structures 119 (2013) 68–84

3

4

2.5

3.5

True Pareto point

True Pareto point Obtained Pareto point based on adaptive approximiation model

Obtained Pareto point based on adaptive approximiation model

3

2

fp2

fp2

2.5 1.5

2 1 1.5 0.5 0

1 4

6

8

10

12

14

0.5

4

6

8

12

14

fp1

(a) Uncertainty level 20%

(b) Uncertainty level 30%

5

16

7 True Pareto point Obtained Pareto point based on adaptive approximiation model

4.5

6 True Pareto point Obtained Pareto point based on adaptive approximiation model

4 5

3.5

4

fp2

3

fp2

10

fp1

2.5

3

2 2 1.5 1

1 0.5

6

8

10

12

14

16

18

0

20

6

8

10

12

14

16

18

fp1

fp1

(c) Uncertainty level 40%

(d) Uncertainty level 50%

20

22

Fig. 10. Pareto points with different uncertainty levels for test function 1. (a) Uncertainty level 20%. (b) Uncertainty level 30%. (c) Uncertainty level 40%. (d) Uncertainty level 50%.

Table 3 The results from the two different methods. Test functions

Table 4 Initialisation and results from the test functions (from [5]).

Total number of points evaluations IUMOO with approximations

1

112 112 108

2500  10  2 = 5.0  104

2

252 248 266

2500  10  2 = 5.0  104

Finally the proposed IUMOO method is applied to ~f p1 and ~f p2 for determining the Pareto front corresponding to the weighting factor b = 0.5. Using the proposed optimisation method, different weighting factors will result in different Pareto fronts, in contrast with the conventional multi-objective optimisation methods using the linear combination method. Test function 2, a more complex threedimensional problem, can be solved in a similar way.

5. Engineering applications This section focuses on a family of practical engineering applications, particularly vehicle crashworthiness designs. Thin-walled beams are typical structures used to absorb energy during vehicle

Test functions

M

Ns

x(0)

1

8

10



 Dxo1 =2; . . . ; Dxon =2

2

11

18



 Dxo1 =2; . . . ; Dxon =2

D(1)

Number of point evaluations

Dxo/4 Dxo/2 3Dxo/ 4

87 89 87

Dxo/4 Dxo/2 3Dxo/ 4

224 221 225

frontal and rear collisions. Therefore, it is important to investigate the design optimisation of thin-walled beam structures with regards to vehicle crashworthiness [51,52]. In general, during a crash, the requirements on a vehicle are to retain the shape of the passenger cabin without any obvious deformations in order to provide passengers with survival space and to absorb as much crash energy as possible in the frame to minimise the collision force transmitted to the passengers. The design problem for thin-walled beam described herein is taken from numerical cases in the literature [51,52]. The computational model used for this study is shown in Fig. 12. The initial velocity of the structure is 13.8 m/s and the impact duration time

81

F. Li et al. / Computers and Structures 119 (2013) 68–84

0.7

True Pareto point Obtained Pareto point based on adaptive approximiation model

True Pareto point

0.6

Obtained Pareto point based on adaptive approximiation model

0.6 0.5 0.5 0.4

fd2

fd2

0.4

0.3

0.3 0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0

0.7

0

0.1

0.2

0.3

0.4

0.5

0.6

fd1

fd1

(a) Uncertain uncertainty level 20%

(b) Uncertain uncertainty level 30%

0.7

0.7

0.7 True Pareto point Obtained Pareto point based on adaptive approximiation model

0.6

True Pareto point

0.6 0.5

0.4

0.4

fd2

fd2

0.5

Obtained Pareto point based on adaptive approximiation model

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0

0.7

0

0.1

0.2

0.3

0.4

0.5

0.6

fd1

fd1

(c) Uncertain uncertainty level 40 %

(d) Uncertain uncertainty level 50%

0.7

Fig. 11. Pareto points with different uncertainty levels for test function 2. (a) Uncertain uncertainty level 20%. (b) Uncertain uncertainty level 30%. (c) Uncertain uncertainty level 40 % (d) Uncertain uncertainty level 50%.

Table 5 Results for an uncertainty level of 10%. Pareto points

1

2

3

4

5

Design vector  L R g~ ; g~  1L 1R  g ;g  1L 1R  g~ ; g~  2L 2R  g2 ; g2 k metamodel k h i ~f L ; ~f R 1 1  L R f ;f h1 1 i ~f L ; ~f R 2 2  L R f2 ; f2 ~f p

(0.9400, 3.0) [0.5255, 0.3112]

(2.4476, 2.4446) [0.2640, 0.2256]

(1.8895, 2.4955) [0.2584, 0.2442]

(1.4666, 2.5605) [0.2879, 0.2550]

(3.0268, 1.9098) [0.3423, 0.2563]

[0.5206, 0.3176]

[0.2635, 0.2433]

[0.2549, 0.2445]

[0.2874, 0.2551]

[0.3428, 0.2556]

[71.438, 46.696]

[2.9109, 1.6087]

[15.0129, 9.2362]

[29.645, 18.964]

[1.4190, 0.2438]

[71.617, 46.686]

[3.1580, 1.3260]

[15.0146, 9.2201]

[29.649, 18.962]

[1.4193, 0.2433]

0.8074, 1.0000 0.8000, 1.0000 [16.591, 21.684]

0.8268, 1.0000 0.8053, 1.0000 [7.9394, 10.1976]

0.8022, 1.0000 0.8005, 1.0000 [11.3659, 14.6069]

0.8003, 1.0000 0.8001, 1.0000 [14.2490, 18.3371]

0.8171, 0.9404 0.8180, 0.9407 [6.6301, 8.2984]

[16.591, 21.686]

[7.9393, 10.1975]

[11.3659, 14.6069]

[14.249, 18.337]

[6.6323,8.3011]

[0.3659, 0.6683]

[1.3064, 2.2440]

[0.8012, 1.4511]

[0.7993, 1.4448]

[2.4242, 4.1499]

[0.3652, 0.6666]

[1.3062, 2.2440]

[0.7993, 1.4448]

[0.7994, 1.4448]

[2.4243, 4.1502]

10.8420

5.0988

7.3034

9.1687

4.1492

fp1 ~f p

10.8431 0.3341

5.0988 1.1220

7.3034 0.7256

9.1687 0.7224

4.1506 2.0749

fp2

0.3333

1.1220

0.7224

0.7224

2.0751

1

2

is 30 ms. Based on research in references [51,52], the typical responses of interest for vehicle crashworthiness include deformation, velocity, acceleration, intrusion, section forces, rigid wall forces, and absorbed energy. For engineering structures, we can find several different indices to evaluate the structural crashworthiness performance. However, these performance indices usually conflict, such as minimising the crash force transferred to the passengers through weak energy absorbing capabilities and the max-

imising of the material energy absorbing abilities to dissipate the rigid wall forces. Hence, crashworthiness designs are naturally characterised by multiple objectives. In this study, the maximum internal energy Ed and the average rigid wall force Fm are used as the two objective functions. The thickness t, height h, width w, and space d of two neighbouring spot welds are taken as design variables. The nominal values for the Young’s Modulus E and the yield stress rs are 210,000 MPa

82

F. Li et al. / Computers and Structures 119 (2013) 68–84

Table 6 Results for an uncertainty level of 10%. Pareto points

1

2

3

4

5

Design vector  L R g~ ; g~  1L 1R  g ;g  1L 1R  g~2 ; g~2  L R g2 ; g2 k metamodel k h i ~f L ; ~f R 1 1  L R f ;f h1 1 i ~f L ; ~f R 2 2  L R f2 ; f2 ~f

(2.980, 1.829) [0.7381, 0.3730]

(0.7839, 3.0) [1.3988, 0.5239]

(2.2376, 2.3641) [0.6047, 0.3429]

(1.8254, 2.3648) [0.6119, 0.3372]

(1.2820, 2.4980) [0.7682, 0.3800]

[0.7356, 0.3725]

[1.3948, 0.5245]

[0.6044, 0.3466]

[0.6120, 0.3360]

[0.7682, 0.3791]

[2.4860, 0.3730]

[95.8341, 0.5239]

[6.3807, 0.3429]

[16.3048, 0.3372]

[40.8055, 0.3800]

[2.2915, 0.3725]

[95.8144, 0.5245]

[6.3625, 0.3466]

[16.3067, 0.3360]

[40.8390, 0.3791]

0.7993, 0.8025 0.7992, 0.8134 [6.3360, 9.9207]

0.8055, 1 [0.8055, 1] [15.403, 26.365]

0.7964, 1 0.7933, 1 [8.2836, 13.5988]

0.8027, 1 0.8038, 1 [10.768, 17.655]

0.7997, 1 0.7640, 1 [13.9152, 23.055]

[6.3402, 9.9205]

[15.405, 26.366]

[8.2839, 13.5988]

[10.768, 17.659]

[13.915, 23.005]

[1.8339, 5.4835]

[0.2628, 0.8886]

[0.9300, 2.8637]

[0.7936, 2.5556]

[ 0.5905, 1.9777]

[1.8339, 5.4836]

[ 0.2635, 0.8808]

[0.9302, 2.8637]

[0.7936, 2.5556]

[0.5903, 1.9778]

4.960

13.1824

6.8121

8.8280

11.5028

fp1 ~f p 2 fp2

4.9609 2.7422

13.1828 0.4443

6.8443 1.4445

8.8280 1.2778

11.5028 0.9890

2.7424

0.4404

1.4767

1.2778

0.9889

p1

d

w

h

t

Spot-welding points Fig. 12. A closed-hat beam impacting the rigid wall and its cross section (mm).

Fig. 14. The Pareto optimal points of crashworthiness for the close-hat beam (left).

Fig. 13. Typical deformation behaviour seen in the finite element model.

and 525 MPa, respectively. Due to the uncertainties caused by manufacturing and measurement errors, E and rs are treated as uncertain parameters, and the uncertainty is ±10% of from their nominal values, defined as: E 2 [189,000 MPa, 231,000 MPa] and rs 2 [472.5 MPa, 577.5 MPa].

The FE model consists of 5760 shell elements. A 300 kg mass is attached to the free end of the beam during the crash analysis to supply the crushing energy. The FEM simulation is carried out using the explicit non-linear finite element code LS-DYNA. The typical deformation behaviour of the numerical model is displayed in Fig. 13. The parameters for the NSGA-II are identical to those in Table 1 based on numerical testing experience. The total number of points evaluated is 152, including the initial 60 sampling points and an additional 92 points. The weighting factor is set to b = 0.3 in this engineering problem for the proposed IUMOO method,

F. Li et al. / Computers and Structures 119 (2013) 68–84

83

Fig. 15. The Pareto optimal points with initial and final sample points (right). Fig. 17. The Pareto optimal points with initial and final sample points (right).

Fig. 16. The Pareto optimal points of crashworthiness for the close-hat beam (left).

which is different from conventional multi-objective optimisation problems using the linear combination method to find the Pareto front. In this research, b = 0.3 corresponds to one Pareto front, and different weighting factors will correspond to different Pareto fronts. Therefore, this parameter is actually changed in order to obtain all of the Pareto fronts. Based on the uncertain approximation optimisation model, the parameter n is specified as 0, and / and w are specified as 50 and 20 for objective 1 and 12,000 and 500 for objective 2, respectively. The uncertain multi-objective optimisation problem is formulated as

8 min : ðf1 ; f2 Þ; > > t;d;w;h > > > > > f 1 ðt; d; w; h; E; rs Þ ¼ Ed ; > > > > > f > 2 ðt; d; w; h; E; rs Þ ¼ F m ; > > 8 > > 0:8 mm 6 t 6 2:5 mm; < > > > > > 10 mm 6 d 6 40 mm; > > > > > > < > 80 mm 6 w 6 120 mm; > > > s:t: > > > 80 mm 6 h 6 120 mm; > > > > > > > > > > E 2 ½189; 000 MPa;231000 MPa; > > > > > : : rs 2 ½472:5 MPa;577:5 MPa:

tained with the present method for an uncertainty level 10% is given in Figs. 14 and 15. The penalty function value for energy absorption is set from 1.576 to 1.425, and the penalty function value for average force is set from 1.558 to 1.722. The numbers added from the FEM evaluations are 92. From Fig. 14, it is found that in the initial stage, a poor quality Pareto set may be obtained, which shows that the initial approximation model has a relatively low accuracy. Along with the optimisation process, the approximation models are updated with respect to the design space and uncertainty space by adding points so that a good Pareto set can be obtained via the adaptive approximation scheme. Figs. 16 and 17 denote the results obtained using an uncertainty level of 20%. Fig. 16 shows that a larger variation in the penalty function value compared to an uncertainty level of 10%, and Fig. 17 shows that the initial Pareto set is worse that the final Pareto set. A total of 82 points were added as a result of the FEM simulation. This shows that the uncertainty level has an obvious influence on the Pareto solutions. Along with the optimisation process, the approximation models are updated with respect to the design and uncertainty spaces by more adding points so that a finer Pareto set can be obtained from the adaptive scheme. The numerical results show that the proposed method leads to a satisfactory level of numerical accuracy. The present method can be regarded as an efficient and effective methodology for engineering problems (e.g. at least for the current engineering design) based on the FEM numerical model and simulation tools. In our ongoing research, some advanced engineering designs with more design and uncertain variables will be conducted to further demonstrate the effectiveness of the proposed approach.

6. Conclusions

ð22Þ

Initially, 60 sampling points are generated through the LHD for constructing the Kriging approximation models for the objective functions in the design and uncertainty spaces. Two points are selected for re-analysis in the IUMOO at each iteration. The Pareto set ob-

This paper proposes a systematic design optimisation method for structures with uncertain-but-bounded parameters. The interval method is used to describe the uncertainty and the Kriging model is applied to generate for the approximation model. The interval number programming method is utilised to transform each uncertain optimisation problem with a single objective function into a deterministic multi-objective optimisation problem. Using the Kriging model and the LHD method, an adaptive approximation scheme is established to estimate the objective functions and constraints and to improve numerical accuracy and computational efficiency. According to the space-filling design criteria, the approximation model is sequentially updated in the design and uncertain spaces until the optimisation converges.

84

F. Li et al. / Computers and Structures 119 (2013) 68–84

The results from typical numerical examples demonstrate that the proposed method can effectively search the Pareto frontier using estimated approximation models. Compared to most conventional methods, this method is able to retain an unchanged approximation space. The added sample points are used in subsequent steps to construct more accurate approximation models. In doing so, the cost of the computationally intensive re-evaluation procedure can be greatly reduced, as it is not required to re-evaluate the response values of all the sampled points for new approximations. Only a limited number of sample points are required, compared to most conventional approximation methods that require a uniformly distributed set of sample points. The proposed method has also been applied to a family of typical engineering applications, namely, the crashworthiness vehicle design, using a closed-hat beam with interval uncertain parameters. Acknowledgements The research is supported by the National Natural Science Foundation of China (NSFC) (11272067), Natural Science Foundation Collaborative Research Grant for Overseas, Hong Kong and Macau Scholars (51228801), Hunan Provincial Natural Science Foundation of China (12JJ3044), The Huxiang Scholar funding from Changsha University of Science and Technology, The Open Fund in the Key Laboratory of Manufacture and Test Techniques for Automobile Parts (Chongqing University of Technology) Ministry of Education (2010KLMT13), The Open Fund of Engineering Research Center of Catastrophic Prophylaxis and Treatment of Road & Traffic Safety (Changsha University of Science & Technology) (KFJ110303). References [1] Deb K. Multi-objective optimization using evolutionary algorithms. New York: Wiley; 2001. [2] Marler RT, Arora JS. Survey of multi-objective optimization methods for engineering. Struct Multidiscip Optim 2004;26:369–95. [3] Lagaros ND, Plevris V, Papadrakakis M. Multi-objective design optimization using cascade evolutionary computations. Comput Methods Appl Mech Eng 2005;194:3496–515. [4] Lin J, Luo Z, Tong L. A new multi-objective programming scheme for topology optimization of compliant mechanisms. Struct Multidiscip Optim 2010;40:241–55. [5] Liu GP, Han X, Jiang C. A novel multi-objective optimization method based on an approximation model management technique. Comput Methods Appl Mech Eng 2008;197:2719–31. [6] Ben-Haim Y. A non-probabilistic measure of reliability of linear systems based on expansion of convex models. Struct Saf 1995;17:91–109. [7] Elishakoff I. A non-probabilistic concept of reliability. Struct Saf 1994;14:227–45. [8] Schuëller GI, Jensen HA. Computational methods in optimization considering uncertainties-an overview. Comput Methods Appl Mech Eng 2008;198:2–13. [9] Maute K, Frangopol DM. Reliability-based design of MEMS mechanisms by topology optimization. Comput Struct 2003;81:813–24. [10] Beyer HG, Sendhoff B. Robust optimization – a comprehensive survey. Comput Methods Appl Mech Eng 2007;196:3190–218. [11] Wang NF, Yang YW. Structural design optimization subjected to uncertainty using fat Bezier curve. Comput Methods Appl Mech Eng 2009;199:210–9. [12] Wang NF, Yang YW. Design structures subjected to uncertainty using wide Bezier curve. IEEE Congr Evol Comput 2009:3141–8. [13] Limbourg P, Kochs HD. Multi-objective optimization of generalized reliability design problems using feature models – a concept for early design stages. Reliab Eng Syst Saf 2008;93:815–28. [14] Fieldsend JE, Everson RM. Multi-objective optimisation in the presence of uncertainty. IEEE Congr Evol Comput 2005;1:243–50. [15] Ben-Haim Y, Elishakoff I. Convex models of uncertainties in applied mechanics. Amsterdam: Elsevier Science Publisher; 1990. [16] Luo Z et al. Fuzzy tolerance multilevel approach for structural topology optimization. Comput Struct 2006;84:127–40. [17] Luo Z et al. A new hybrid fuzzy-goal programming scheme for multi-objective topological optimization of static and dynamic structures under multiple loading conditions. Struct Multidiscip Optim 2006;31:26–39. [18] Ishibuchi H, Tanaka H. Multiobjective programming in optimization of the interval objective function. Eur J Oper Res 1990;48:219–25.

[19] Jiang C, Han X, Liu GP. Uncertain optimization of composite laminated plates using a nonlinear interval number programming method. Comput Struct 2008;86:1696–703. [20] Jiang C, Han X, Liu GR. Optimization of structures with uncertain constraints based on convex model and satisfaction degree of interval. Comput Methods Appl Mech Eng 2007;196:4791–800. [21] Kang Z, Luo Y, Luo Z, Li A. Continuum topology optimization with nonprobabilistic reliability constraints based on multi-ellipsoid convex model. Struct Multidiscip Optim 2009;39:297–310. [22] Kang Z, Luo Y. Reliability-based structural optimization with probability and convex set hybrid models. Struct Multidiscip Optim 2010;42:89–102. [23] Möller B, Beer M. Engineering computation under uncertainty-capabilities of non-traditional models. Comput Struct 2008;86:1024–41. [24] Qiu Z, Elishakoff I. Anti-optimization of structures with large uncertain-butnon-random parameters via interval analysis. Comput Methods Appl Mech Eng 1998;152:361–72. [25] Gao W, Song C, Tin-Loi F. Probabilistic interval analysis for structures with uncertainty. Struct Saf 2010;32:191–9. [26] Lombardi M, Haftka RT. Anti-optimization technique for structural design under load uncertainties. Comput Methods Appl Mech Eng 1998;157:19–31. [27] Qiu Z, Wang X. Comparison of dynamic response of structures with uncertainbut-bounded parameters using non-probabilistic interval analysis method and probabilistic approach. Int J Solids Struct 2003;40:5423–39. [28] Rao SS, Cao L. Optimum design of mechanical systems involving interval parameters. J Mech Des 2002;124:465–72. [29] Gunawan S, Azarm S. A feasibility robust optimization method using sensitivity region concept. J Mech Des 2005;127:858–65. [30] Li F, Luo Z, Sun G, Zhang N. An uncertain multidisciplinary design optimization method using interval convex models. Eng Optim 2012. http://dx.doi.org/ 10.1080/0305215X.2012.690871. [31] Barthelemy JFM, Haftka RT. Approximation concepts for optimum structural design – a review. Struct Multidiscip Optim 1993;5:129–44. [32] Simpson TW, Mistree F. Kriging models for global approximation in simulation-based multidisciplinary design optimization. AIAA J 2001;39:2233–41. [33] Simpson TW et al. Metamodels for computer-based engineering design: survey and recommendations. Eng Comput 2001;17:129–50. [34] Yun Y, Yoon M, Nakayama H. Multi-objective optimization based on metamodeling by using support vector regression. Optim Eng 2009;10:167–81. [35] Sinha K. Reliability-based multiobjective optimization for automotive crashworthiness and occupant safety. Struct Multidiscip Optim 2007;33:255–68. [36] Liu W, Yang Y, Xing Z, Zhao L. Springback control of sheet metal forming based on the response-surface method and multi-objective genetic algorithm. Mater Sci Eng A 2009;499:325–8. [37] Fang H, Horstemeyer MF. Global response approximation with radial basis functions. Eng Optim 2006;38:407–24. [38] Fang H, Rais-Rohani M, Liu Z, Horstemeyer MF. A comparative study of metamodeling methods for multiobjective crashworthiness optimization. Comput Struct 2005;83:2121–36. [39] Myers RH, Montgomery DC. Response surface methodology: process and product optimization using designed experiments. New York: Wiley; 1995. [40] Wilson B, Cappelleri D, Simpson TW, Frecker M. Efficient Pareto frontier exploration using surrogate approximations. Optim Eng 2001;2:31–50. [41] Li F, Li G. Interval-based uncertain multi-objective optimization design of vehicle crashworthiness. CMC Comput Mater Continua 2010;15:199–220. [42] Morris MD, Mitchell TJ. Exploratory designs for computational experiments. J Stat Planning Infer 1995;43:381–402. [43] Shan S, Wang GG. An efficient Pareto set identification approach for multiobjective optimization on black-box functions. J Mech Des 2005;127:866–74. [44] Yang BS, Yeun YS, Ruy WS. Managing approximation models in multiobjective optimization. Struct Multidiscip Optim 2002;24:141–56. [45] Li G, Li M, Azarm S, Hashimi S, Ameri T, Qasas N. Improving multi-objective genetic algorithms with adaptive design of experiments and online metamodeling. Struct Multidiscip Optim 2009;37:447–61. [46] Boggs PT, Tolle JW. Sequential quadratic programming. Acta Numer 1995;4:1–51. [47] Sakata S, Ashida F, Zako M. Structural optimization using Kriging approximation. Comput Methods Appl Mech Eng 2003;192:923–39. [48] Johnson ME, Moore LM, Ylvisaker D. Minimax and maximin distance designs. J Stat Planning Infer 1990;26:131–48. [49] Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans Evol Comput 2002;6:182–97. [50] L Soares G et al. Interval robust multi-objective algorithm. Nonlinear Anal Theory Methods Appl 2009;71:e1818–25. [51] Chen W, Wierzbicki T. Relative merits of single-cell multi-cell and foam-filled thin-walled structures in energy absorption. Thin Wall Struct 2001;39:287–306. [52] Xiang Y, Wang Q, Fan Z, Fang H. Optimal crashworthiness design of a spotwelded thin-walled hat section. Finite Elem Anal Des 2006;42:846–55.