An innovative reliability-based design optimization method by combination of dual-stage adaptive kriging and genetic algorithm

Abstract

Purpose

This study aims to propose an efficient method for solving reliability-based design optimization (RBDO) problems.

Design/methodology/approach

In the proposed algorithm, genetic algorithm (GA) is employed to search the global optimal solution of design parameters satisfying the reliability and deterministic constraints. The Kriging model based on U learning function is used as a classification tool to accurately and efficiently judge whether an individual solution in GA belongs to feasible region.

Findings

Compared with existing methods, the proposed method has two major advantages. The first one is that the GA is employed to construct the optimization framework, which is helpful to search the global optimum solutions of the RBDO problems. The other one is that the use of Kriging model is helpful to improve the computational efficiency in solving the RBDO problems.

Originality/value

Since the boundaries are concerned in two Kriging models, the size of the training set for constructing the convergent Kriging model is small, and the corresponding efficiency is high.

Citation

Feng, K. and Lu, Z. (2022), "An innovative reliability-based design optimization method by combination of dual-stage adaptive kriging and genetic algorithm", Multidiscipline Modeling in Materials and Structures, Vol. 18 No. 4, pp. 562-581. https://doi.org/10.1108/MMMS-04-2022-0058

1. Introduction

Due to machining error, assembly error, mutative load, abrasion and friction, uncertainties extensively exist in the design, service and maintenance of the structure or system (Babak et al., 2022; Feng et al., 2019b; Li et al., 2021). Conventional design optimization (CDO) employs the safety factor to deal with these uncertainties, which may lead to a conservative design with great weight or large size. Besides, CDO cannot give a quantitative index about the safety of the designed structure or system. Thus, reliability-based design optimization (RBDO) is developed to overcome the shortcomings of CDO, where the uncertainties of inputs are described as randomness by their probability density functions (PDFs). RBDO is able to help structural designers balance cost and safety, therefore produce designs which not only are economical but also satisfying reliability constraints (Aoues and Chateauneuf, 2008; Yang and Hsieh, 2011). The RBDO problem can be generally formulated as,

(1)Mind  C(d)s.t. {Pr{gk(X|d)≤0}≤Pfk∗k=1,2,...,phl(d)≤0l=1,2,...,qdL≤d≤dUd∈Rm

where X=[X1,X2,⋯,Xn]T is the n-dimensional random input vector, d=[d2,d2,⋯,dm]T is the m-dimensional design parameter vector which usually is the distribution parameter (mean, standard deviation, etc.) of the random input, C(d) is the objective function, Pr{gk(X|d)≤0}(k=1,2,...,p) represents the failure probability for the kth performance function, and Pfk∗ is the corresponding threshold of the failure probability, hl(d)(l=1,2,...,q) denotes the lth deterministic constraint, dL and dU are the lower and upper bounds of the design parameter vector d, respectively. According to Eq. (1), it can be easily observed that the double-layer nested process is involved in the RBDO problem, where the outer is the multi-parameter design optimization of the design parameter vector and the inner is the reliability analysis by considering the randomness of the input vector. Hence, more steps and much computational cost are needed in solving the RBDO problem compared with the CDO problem.

In recent decades, numerous tools for RBDO have been developed by different scholars from various engineering fields (Allen and Maute, 2005; Huyse et al., 2002; Gu et al., 2001), and they can be mainly classified into two groups based on the reliability analysis method involved (Marcos and Gerhart, 2010), i.e. the approximately analytical technique based methods and the simulation technique based methods. In approximately analytical technique based methods, the failure probability Pr{gk(X|d)≤0}    (k=1,2,...,p) defined in Eq. (1) is solved by the First Order and Second Moment (FOSM) method (Roger et al., 1999; Zhao and Ono, 1999) or improved versions of FOSM, thus the RBDO problem can be rewritten as,

(2)Mind  C(d)s.t.  {βk(d)≥βk∗k=1,2,...,phl(d)≤0l=1,2,...,qdL≤d≤dUd∈Rm

where βk(d)(k=1,2,...,p) is the kth reliability index of the corresponding performance function gk(X|d), and βk∗=−Φ−1(Pfk∗)(k=1,2,...,p) is the kth threshold of the reliability index in which Φ−1(⋅) is the inverse function of the cumulative distribution function for the standard normal distribution. The methods for settling the RBDO problem defined by Eq. (2) can be divided into three categories, i.e. double-loop methods, single-loop methods and decoupling methods. Double-loop method is the most direct approach of the approximately analytical technique based methods, where the outer loop is to optimize the design parameters and the inner loop is to evaluate the reliability index for the given set of design parameters (Nikolaidis and Burdisso, 1988; Tu et al., 2001). To reduce the computational cost of solving the RBDO problem, several so-called single-loop methods were proposed by different structural designers. For instance, Kuschel and Rackwitz (1997) constructed a single-layer RBDO model by employing the Karush-Kuhn-Tucker (KKT) conditions (Bonnans et al., 2003) and Lagrange multipliers, where the inner loop for estimating the reliability index is avoided. Other examples of the single-loop methods can be discovered in Kharmanda et al. (2002) and Mohsine et al. (2006), where an extended formulation was established for working out the RBDO problem, i.e. the generalized objective function was indicated as the product of the reliability indices and the original objective function. The basic idea of the decoupling methods is to extract the information from the reliability analysis period, which can be used at the optimization process so as to further improve computational efficiency. Up to our knowledge, the first decoupling method was proposed by Li and Yang (1994), where each reliability constraint is replaced by a linear approximation formula with the help of the Taylor expansion and the reliability based sensitivities, then the RBDO problem can be transformed to solve a series of deterministic optimization problems. Subsequently, various alternatives have been studied as well, such as the adaptive sequential linear programming algorithm proposed by Chan and co-workers (Chan et al., 2007), sequential approximate programming strategy algorithm introduced by Cheng et al. (2006), sequential optimization and reliability assessment algorithm developed by Du and Chen (2004), etc. Although the approximately analytical technique based methods have been developed rapidly in recent decades and some of them are very efficient, the basic framework of these methods are still based on FOSM method. Because FOSM method is only an approximate reliability analysis technique, it may bring large error in reliability analysis for the strong nonlinear performance function, so the optimization results obtained by the approximately analytical technique based methods might not satisfy the failure probability constraints defined in Eq. (1).

In the simulation technique based methods, the failure probabilities with respect to each performance function is estimated by the simulation technique, which might be computationally expensive, especially for analyzing large scale structures with complicated finite element models. Thus, the surrogate model (Liu et al., 2019; Zhou et al., 2019) has been widely used to evaluate the failure probabilities involved in the RBDO, which can dramatically reduce the total computational cost in the whole optimization process. Surrogate models can be applied to a particular RBDO problem with two different forms. In the first form, the surrogate model is employed to directly represent the true performance function, then the RBDO problem can be solved by replacing the true performance function with the corresponding surrogate model. A typical case of this class of surrogate models is the polynomial chaos technique (Xiu and Karniadakis, 2003), and Anup and Debraj (2016) used the polynomial chaos technique to solve the RBDO in aeroelastic stability problems. It should be pointed out that the accuracy of the optimization result obtained by this directly surrogate model based method would subject to several factors, such as the global precision of the surrogate model, sample size in evaluating the failure probability, the optimization algorithm selected and so on. In the second form, the surrogate model is used to approximately construct the failure probability function (Feng et al., 2019a) which is defined as a function of failure probability with respect to the design parameter vector. Next, by substituting the failure probability constraint to the proxy failure probability function, the RBDO problem can be translated into a deterministic optimization problem. Some specific cases of the second form can be referred to Hurtado (2004), Missoum et al. (2007) and Vincent et al. (2011). This type of method is very superior in theory, but obtaining the failure probability function with enough precision especially for multi-dimensional problem is not easy in practice. In addition, most of existing tools for solving RBDO problem use the gradient based optimization algorithms, such as sequential quadratic programing and interior–point method, which have quick astringency, but may converge to the local optimum.

This contribution expects to develop a novel algorithm for RBDO by using dual-stage adaptive Kriging and genetic algorithm, which has three advantages: 1) accurately measuring the reliability; 2) better global convergence; 3) high computational efficiency. In order to achieve this purpose, the simulation technique is used to estimate the failure probability of the structure so as to accurately measure the reliability, the genetic algorithm (GA) is employed to search the global optimum solution, and the adaptive Kriging model is successively employed to accurately and efficiently judge whether or not a specific sample of random input vector belongs to failure domain and a certain realization of design parameter vector pertains to feasible region.

The main content of this contribution is organized as follows. The outline of the proposed algorithm is presented in Section 2. The algorithm details and implementation are described in Section 3. Three numerical studies are employed to demonstrate the accuracy and efficiency of the proposed algorithm in Section 4. Concluding remarks are summarized in Section 5.

2. Algorithm outline

Generally, solving an optimization problem with simple objective function and complex constraint is harder than solving that with complex objective function and simple constraint, therefore the RBDO defined in Eq. (1) is primarily transformed into the following form,

(3)Mind  C(d)−CminCmax−Cmin+∑k=1pPRk⋅I[Pr{gk(X|d)≤0}−Pfk∗]+∑l=1qPDl⋅I[hl(d)]s.t. dL≤d≤dU d∈Rm

where Cmin and Cmax respectively denote the minimum and maximum values of the objective function C(d) by only considering the lower and upper bounds of the design parameter vector d, which can be easily obtained without any evaluation of the performance function, therefore C(d)−CminCmax−Cmin would be in the interval [0,1], and the units of C(d) can be eliminated as discussed in Fu et al. (2017), I[input] represents the indicator function, i.e.

(4)I[input]={0if  input≤01if  input>0

PRk(k=1,2,...,p) and PDl(l=1,2,...,q) respectively express the penalty factors which are used as a punishment when the design parameter vector does not satisfy the reliability constraint and the deterministic constraint, respectively.

The penalty function in Eq. (3) is used to identify whether a realization of design parameters is located in the feasible region. If the realization of design parameters is not located in the feasible region, a big value will be added to the objective function as a punishment through the penalty factor. Generally, a big value of penalty factor means a good effect of punishment. However, a big value of penalty factor also implies high nonlinearity of the generalized objective function containing the original objective function and penalty function, which increases the solving difficulty of the optimization problem. On the other hand, if the value of penalty factor is too small, the event of misjudging the state of the realization of deign parameters may occur, which will further lead to an inaccurate optimal result. As the original objective function is primarily normalized in the interval [0,1], the penalty factor is set to 2 in the proposed method and it is enough to punish the realization of design parameters not locating in the feasible region.

In order to acquire a global optimum solution, GA is employed to solve the optimization problem defined in Eq. (3) in this contribution. GA is a computational model inspired by the Darwin's biological evolution theory, which learns from the evolutionary law of the biological world (Mitchell, 1996). The evolution often begins with a population composed of randomly generated individuals, and then iterative process is gradually performed with the population in each iteration referred to as a generation. In each iteration, the fitness of every individual in current population is calculated, and the fitness function is chosen as one monotonic function of the optimization goal in general. The higher the fitness of individuals is, the greater possibility these individuals will be selected to the next generation. And each individual's chromosome is modified (combined the crossover and mutation with the genetic operators of natural genetics) to construct a new generation. The new generation representing the new solution set is then employed in the next iteration of the algorithm. Frequently, the stop condition of genetic algorithm can be classified into two types. The first one is that the number of generations reaches to the given maximum number, and the second one is that the fitness level of the population reaches to the threshold. Compared to the latter one, the former one is more conservative but simpler and more robust. Thus, the latter one is chosen as a more stable stop condition in the proposed method. Although this stop condition may lead to a waste of computational effort, this waste is quite small thanks to the inclusion of Kriging model.

Compared with the gradient based optimization algorithms, GA has two important advantages, the first one is that it directly operates on design parameters where the derivatives of the objective function with respect to design parameters are not needed, and the second one is that it has the capability of searching global optimal solution. Nevertheless, repeated fitness function calculation for complicated problems is usually the most limitation of GA. Fortunately, the occurrence of surrogate models dramatically reduces the computational cost of GA and expands its application range, especially for the complicated engineering problems.

For the optimization problem described in Eq. (3), the fitness function F(d) can be defined as the exponential function of the optimization goal, i.e.

(5)F(d)=exp(−{C(d)−CminCmax−Cmin+∑k=1pPRk⋅I[Pr{gk(X|d)≤0}−Pfk∗]+∑l=1qPDl⋅I[hl(d)]})

From Eq. (5), it can be observed that the smaller the objective function of the optimization problem defined in Eq. (3) is, the bigger the fitness of the individual is. In addition, it can be concluded that the fitness function F(d) is comprised of three components, i.e. the first part related to the original objective function C(d), C(d)−CminCmax−Cmin; the second part related to the reliability constraints ∑k=1pPRk⋅I[Pr{gk(X|d)≤0}−Pfk∗]; and the third part related to the deterministic constraints ∑l=1qPDl⋅I[hl(d)]. Generally, the most time-consuming work in reliability analysis or RBDO is the evaluation of the performance function, and the total number of the performance function evaluation is usually regarded as the total computational cost in such problems (Cheng et al., 2006; Feng et al., 2020; Wang et al., 2017; Yun et al., 2019). Thus, the computational cost in estimating F(d) only exists in the evaluation of its second part ∑k=1pPRk⋅I[Pr{gk(X|d)≤0}−Pfk∗], and for convenience of expression, the auxiliary function Lk(d)(k=1,2,...,p)  is introduced as,

(6)Lk(d)=Pr{gk(X|d)≤0}−Pfk∗,      (k=1,2,...,p) 

According to Eqs. (5) and (6), it can be seen that the most important and difficult issue in solving the RBDO by using GA is to judge whether or not every individual representing design parameter vector in each population is greater than zero, if it is, the indicator function defined in Eq. (4) is equal to 1, if not, that indicator function is equal to 0. In order to efficiently achieve this purpose, the surrogate model can be constructed to approximately replace the original auxiliary function Lk(d). It should be noticed that the surrogate model in this contribution is only used as a classification tool, i.e. the surrogate model judges whether the value of the auxiliary function Lk(d)  at a certain individual representing design parameter vector is larger than zero or not, while accurately evaluating the value of Lk(d)  is not necessary. Thus, the adaptive Kriging model based on U learning function (Echard et al., 2011) can be employed in this contribution, which is one of the most widely used classification tools, and more details of this technique will be introduced in section 3.1. In the first generation of GA, the rough Kriging model Lk(K)(d)(k=1,2,...,p)  is constructed by using a small part of the individuals in this generation, then this Kriging model will be continuously updated by orderly adding the new individual of the first generation according to the selecting criterion until the convergence criterion is satisfied. Next, in subsequent generations, the incipient Kriging model can employ the model obtained in the previous generation, then it is continuously updated so as to accurately distinguish the sign of the individuals in current generation.

In the last paragraph, the process of constructing and updating the Kriging model Lk(K)(d)  is explained in detail. In this process, the function value Lk(d)  at initial training individuals and updated training individuals should be accurately estimated. From Eq. (6), it can be observed that estimating the function value Lk(d)  is essentially a failure probability estimation problem. By using the Monte Carlo Simulation (MCS) and the Kriging model, Pr{gk(X|d)≤0}(k=1,2,...,p) can be estimated by,

(7)Prˆ{gk(X|d)≤0}=1−Prˆ{gk(X|d)>0}≈1−1NX ∑j=1NXI[gk(xj|d)](k=1,2,...,p)

where xj|d is the jth random sample vector when the design parameter vector is fixed at d, and NX denotes the total number of the random samples. Obviously, only the sign of the performance function at each sample, instead of its accurate value, is needed in estimating the failure probability. Besides, for different design parameter vector, the PDF of the random input vector X are different, but the structure of the performance function is the same. Hence, a single adaptive Kriging model based on U learning function is sufficient to estimate the failure probabilities under different design parameter vectors. The basic idea of using a single Kriging model to estimate various failure probabilities is that a Kriging model is firstly built for accurately recognizing the sign of the performance function with respect to the sample in the sample pool constructed based on the first individual (representing design parameter vector) of the first generation in GA. Next, for other individuals, the existing Kriging model is continuously updated in order to accurately distinguish the sign of the performance function with respect to the sample in the sample pool constructed according to the corresponding individual. The flowchart of the proposed algorithm for solving the RBDO problem is given in Figure 1.

3. Algorithm details

From section 2, it is clear that the main contents of the proposed algorithm include two parts; the first one is that GA is employed as the optimization tool to obtain the global optimal solution of the RBDO defined in Eq. (3), and the detailed steps of GA can be found in Leung and Wang (2001). The second one is that the adaptive Kriging model based on U learning function is successively used to accurately and efficiently judge whether a certain realization of design parameter vector pertains to feasible region and a specific sample of random input vector belongs to failure domain. In this section, this adaptive Kriging model is briefly introduced at first, then the detailed implementation of the proposed algorithm in solving RBDO is summarized in subsection 3.2.

3.1 An adaptive kriging model based on U learning function

In this subsection, the surrogated function is expressed by h(Θ), which can be the auxiliary function Lk(d)(k=1,2,...,p)  or the performance function gk(X|d)(k=1,2,...,p), and the sample pool of the input vector Θ is given as θS={θ1,θ2,⋯,θNθ}, in which Nθ denotes the total number of the samples in the sample pool.

The fundamental idea of Kriging model is that the function h(Θ) could be regarded as a realization of a stochastic field h(K)(Θ) which is introduced as,

(8)h(K)(Θ)=f(Θ)ξ+Z(Θ)

in which f(Θ)=[f1(Θ),f2(Θ),⋯,fNf(Θ)] represents the Nf-dimensional basis function vector, ξ=[ξ1,ξ2,⋯,ξNf]T denotes the Nf-dimensional regression coefficient vector, and Z(Θ) expresses a stationary Gaussian process with mean zero and the covariance between any two samples θj1 and θj2 in the sample pool θS is defined in Eq. (9),

(9)cov[Z(θj1),Z(θj2)]=σZ2RZ[Z(θj1),Z(θj2)]

in which σZ is the standard deviation, RZ is the correlation function that can determine the smoothness of the Kriging model and the commonly used Gaussian correlative function (Echard et al., 2011) is employed in this contribution. For any untrained sample θ in the sample pool θS, the Kriging model prediction is,

(10)h(K)(θ)∼N(μh(K)(θ),σh(K)2(θ))

where N(⋅) denotes normal distribution, and μh(K)(θ) and σh(K)(θ) stand for the mean and standard deviation of the prediction h(K)(θ) respectively.

When h(K)(θ)≥0, the probability of misidentifying the sign of h(θ) can be represented as,

(11)PI=Φ(0−|h(K)(θ)|σh(K)(θ))=Φ(−|h(K)(θ)|σh(K)(θ))

When h(K)(θ)<0, the probability of misidentifying the sign of h(θ) can be represented as,

(12)PII=1−Φ(0+|h(K)(θ)|σh(K)(θ))=Φ(−|h(K)(θ)|σh(K)(θ))

Eqs. (11) and (12) show that whatever the sign of h(K)(θ), the probability of misidentifying the sign of h(θ) can be denoted as,

(13)Pmis=Φ(−U(θ))

in which U(x) is known as the U learning function, and it is given by,

(14)U(θ)=|h(K)(θ)|σh(K)(θ)

According to Eqs. (13) and (14), it can be concluded that the smaller the value of U(θ) is, the higher the probability of misidentifying the sign of h(θ). Thus, the new training sample θnew could be chosen as the sample with the smallest value of U(θ), i.e.

(15)θnew=argminθ∈θSU(  θ) 

It is recommended that the process for iteratively updating the Kriging model can stop if U(θ)≥2 is valid for any sample in the sample pool  θS, which demonstrates that the probability of misjudging the sign of h(θ) is equivalent to Φ(−2)=0.0228 (Echard et al., 2011).

3.2 The implementation of the proposed algorithm

A detailed summary of the implementation of the proposed genetic algorithm and adaptive Kriging model based approach in solving RBDO is revealed in this subsection. It is demonstrated as follows.

• Step 1. Generate an initial population d(0)={d(0,1),d(0,2),⋯,d(0,Nd)} of the design parameter vector d where Nd denotes the number of individuals in each population, and initialize the generation number gen to 0, i.e. gen=0.

• Step 2. Randomly select NdT (NdT<<Nd) individuals from d(0), which are recorded as {d(0,s1),d(0,s2),⋯,d(0,sNdT)}. Estimate the value of Lk(d(0,s1))(k=1,2,...,p)  by the following steps.

  • Step 2.1. According to the design parameter vector d(0,s1), randomly generate NX samples of the input vector and construct the sample pool xS(0,s1), i.e. xS(0,s1)={x1(0,s1),x2(0,s1),⋯,xNX(0,s1)}.

  • Step 2.2. Randomly select NXT (NXT<<NX) training samples from the sample pool xS(0,s1), and compute their model outputs. Then, construct the initial Kriging model gk(K)(X)  by employing the selected samples.

  • Step 2.3. According to the Kriging model gk(K)(X) , evaluate the U learning function U(x) of all the samples in xS(0,s1). If U(x)≥2 holds for all these samples, go to Step 2.5; otherwise, go to Step 2.4.

  • Step 2.4. Select xnew=argminx∈xS(0,s1)U(x) as the new training sample and compute the corresponding model output, then update the Kriging model gk(K)(X)  by adding this new training sample. Next, go to Step 2.3.

  • Step 2.5. Based on the Kriging model gk(K)(X), identify all the failure samples in xS(0,s1), and the total number of these failure samples is recorded as NF(0,s1), then Lk(d(0,s1))  can be estimated by Lk(d(0,s1)) =NF(0,s1)NX−Pfk∗.

• Step 3. Estimate the value of Lk(d(0,r))(r=s2,...,NdT)  by the following steps.

  • Step 3.1. According to the design parameter vector d(0,r), randomly generate NX samples of the input vector and construct the sample pool xS(0,r)={x1(0,r),x2(0,r),⋯,xNX(0,r)}.

  • Step 3.2. Based on the current Kriging model gk(K)(X), evaluate the U learning function U(x) of all the samples in xS(0,r). If U(x)≥2 holds for all these samples, go to Step 3.4; otherwise, go to Step 3.3.

  • Step 3.3. Select xnew=argminx∈xS(0,r)U(  x)  as the new training sample and compute the corresponding model output, then update the Kriging model gk(K)(X)  by adding this new training sample. Next, go to Step 3.4.

  • Step 3.4. According to the Kriging model gk(K)(X), recognize all the failure samples in xS(0,r), and the total number of these failure samples is denoted as NF(0,r), subsequently Lk(d(0,r))  can be estimated by Lk(d(0,r))=NF(0,r)NX−Pfk∗.

• Step 4. Construct the initial Kriging model Lk(K)(d)(k=1,2,...,p)  based on the training individuals {d(0,s1),d(0,s2),⋯,d(0,sNdT)}.

• Step 5. According to the Kriging model  Lk(K)(d), evaluate the U learning function U(d) of all the individuals in d(0). If U(d)≥2 holds for all these individuals, go to Step 7; otherwise, go to Step 6.

• Step 6. Select dnew=argmind∈d(0)U(d) as the new training individual and compute the corresponding function value Lk(dnew) by the following steps,

  • Step 6.1. Based on the design parameter vector dnew, randomly generate NX samples of the input vector and construct the sample pool xSnew={x1new,x2new,⋯,xNXnew}.

  • Step 6.2. Based on the current Kriging model gk(K)(X), evaluate the U learning function U(x) of all the samples in xSnew. If U(x)≥2 holds for all these samples, go to Step 6.4; otherwise, go to Step 6.3.

  • Step 6.3. Select xnew=argminx∈xSnewU(x)  as the new training sample and compute the corresponding model output, then update the Kriging model gk(K)(X)  by adding this new training sample. Next, go to Step 6.4.

  • Step 6.4. According to the Kriging model gk(K)(X), distinguish all the failure samples in xSnew, and the total number of these failure samples is recorded as NFnew, then Lk(dnew)  can be estimated by Lk(dnew) =NFnewNX−Pfk∗.

• Step 7. Update the Kriging model  Lk(K)(d) by adding the new individual dnew.

• Step 8. Evaluate the fitness with respect to each individual in d(gen) by using the current Kriging model  Lk(K)(d).

• Step 9. Population Evolution. Firstly, each individual in current population is selected for crossover with given probability  pc. Secondly, each individual in current population may mutate with given probability  pm. Thirdly, among the individuals in  d(gen) and those generated by crossover and mutation, select the Nd elements with largest fitness to form the next generation. If gen>genmax (where genmax is the given maximum number of generation), go to Step 12; otherwise, gen=gen+1 and go to Step 10.

• Step 10. On the basis of current Kriging model  Lk(K)(d), evaluate the U learning function U(d) of all the individuals in d(gen). If U(d)≥2 holds for all these individuals, go to Step 12; otherwise, go to Step 11.

• Step 11. Select dnew=argmind∈d(gen)U(d) as the new training individual and compute the corresponding function value Lk(dnew) by repeating Step 6. Then, go to Step 7.

• Step 12. The individual in current population  d(gen) with largest fitness can be regarded as the optimum solution.

4. Numerical studies

In this section, the proposed genetic algorithm and adaptive Kriging model based approach is applied to three examples from the literature in order to compare the accuracy and efficiency of the proposed approach and existing methods. The parameters of GA and adaptive Kriging model used in the examples are listed in Table 1.

4.1 Standard RBDO test problem

In this subsection, a common RBDO test case used in Mohsen et al. (2014), Yang and Gu (2004) and Youn and Choi (2004). is employed, which can be expressed as,

(16)Mind  C(d)=d1+d2s.t.  {Pr{g1(X)=X12X220−1≤0}≤Φ(−3)=0.00135Pr{g2(X)=(X1+X2−5)230+(X1−X2−12)2120−1≤0}≤Φ(−3)=0.00135Pr{g3(X)=80(X12+8X2+5)2−1≤0}≤Φ(−3)=0.001350≤dj≤10                            j=1,2

where X1 and X2 follow normal distribution with means (d1,d2) and standard deviation 0.3. Mohsen et al. (2014) used five methods, i.e. reliability index approach (RIA), performance measure approach (PMA), single-loop approach (SLA), sequential optimization and reliability assessment (SORA) method, and the simulation-based method (SBM) to solve the RBDO problem defined in Eq. (16), and the corresponding results are established in Table 2, in which Pfmax denotes the maximum failure probability of all the performance functions under the design parameter vector.

According to Eq. (16), it can be computed that the maximum and minimum values of the objective function C(d) by only considering the lower and upper bounds of the design parameter vector d are Cmin=0 and Cmax=20. Thus, the RBDO problem defined in Eq. (16) can be rewritten as,

(17)Mind  C(d)20+∑k=132⋅I[Pr{gk(X)≤0}−0.00135]s.t.     0≤dj≤10      j=1,2

By employing the proposed genetic algorithm and adaptive Kriging model based approach to solve the problem defined in Eq. (17), it can be obtained that the final optimum design is {d1=3.4498,d2=3.2824}, and the convergence of the design parameter variables and cost is depicted in Figure 2. From Figure 2, it can be observed that the convergence of design parameter variables and cost is achieved when gen=102. Table 2 indicates that the final optimum design obtained by different methods is almost the same, but the reliability constraints are not satisfied in PMA, SLA and SORA methods. Thanks to the adaptive Kriging model, the optimum is reached in the proposed method by using only 39 model evaluations, which demonstrates that the proposed method is more efficient than other methods listed in Table 2.

4.2 Bracket structure design

The bracket structure taken from Chateauneuf and Aoues (2008) and Vincent et al. (2011) is shown in Figure 3. This bracket structure is loaded by an external force at its right tip and by its own weight on account of gravity. The two failure modes taken into consideration are described as follows.

1. The maximum bending moment σB in the horizontal beam (recorded as CD) appearing at point B should not surpass the yield strength σS, thus the first performance function can be written as,

(18)g1(σS,ωCD,t,L,P,ρ)=σS−σB

in which the maximum bending moment σB can be expressed as,

(19)σB=6MBωCDt2MB=PL3+ρgωCDtL218

1. The maximum internal force FAB in the scaffold beam (recorded as AB) should not surpass the critical buckling force Fbuckling, hence the second performance function can be written as,

(20)g2(σS,ωAB,ωCD,t,L,P,ρ)=Fbuckling−FAB

in which the critical buckling force Fbuckling and maximum internal force FAB can be expressed as,

(21)Fbuckling=π2EILAB2=9π2EtωAB3⁡sin2⁡θ48L2FAB=1cos⁡θ(3P2+3ρgωCDtL4)

where the eight parameters, i.e. σS, P, E, ρ, L, ωAB, ωCD and t are considered as the random input variables in this example, and theirs contribution types and parameters are established in Table 3.

As shown in Table 3, the means of ωAB, ωCD and t are regarded as three design parameters, and the expected weight of the structure is regarded as the objective function, which can be introduced as,

(22)c(ωAB,ωCD,t)=μρμtμL(439μωAB+μωCD)

thus, the RBDO problem can be formulated as,

(23)Mind  C(d)=μρμtμL(439μωAB+μωCD)s.t.  {Pr{g1(X)≤0}≤Φ(−2)=0.02275Pr{g2(X)≤0}≤Φ(−2)=0.0227550≤μωAB,μωCD,μt≤300 

Based on Eq. (23), it can be estimated that the maximum and minimum values of the objective function C(d) by only considering the lower and upper bounds of the design parameter vector d are Cmax=6259kg and Cmin=174kg. Therefore, the RBDO problem defined in Eq. (23) can be rewritten as,

(24)Mind  C(d)−1746085+∑k=122⋅I[Pr{gk(X)≤0}−0.02275]s.t.     50≤μωAB,μωCD,μt≤300

In Chateauneuf and Aoues (2008), four existing approaches, i.e. RIA, SORA, subset simulation (SS), subset simulation and kriging surrogate model (SS + AK), are used for comparison, and the corresponding results are listed in Table 4. By using the proposed approach to settle the problem defined in Eq. (24), the final optimum design can be obtained as {μωAB=54mm,μωCD=76mm,μt=295mm}, and the convergence of the design parameter variables and cost is drawn in Figure 4. Figure 4 shows that the convergence is achieved when gen=119. According to Table 4, it can be seen that the optimal solution obtained by all the methods except for the SS + AK approach satisfies the two reliability constraints, but that obtained by the proposed method has a smallest cost, which illustrates that the proposed method can get better result compared with the RIA, SORA and SS approaches because the optimal reliability level in this method is formulated in terms of the failure probability instead of the reliability index.

4.3 Passive vehicle suspension design

In this subsection, the passive vehicle suspension studied in Mohsen et al. (2014). is considered as the verification case, and its schematic is shown in Figure 5. The objective of this problem is to minimize the mean square value of the vertical vibration acceleration of the vehicle body that satisfies the following four constraints:

1. the road-holding ability of the vehicle should not less than a certain threshold (g1);

2. the rolling angle should not exceed a certain threshold (g2);

3. the suspension's dynamic displacement should not less than a certain threshold so as to avoid bumper hitting (g3);

4. the tire stiffness should not less than a certain threshold because the tire life is an increasing function with respect to the tire stiffness (g4).

The RBDO problem is defined by Eq. (25) where the mean of suspension stiffness c(kg/cm), tire stiffness ck(kg/cm) and damping coefficient k(kg/cm  s), i.e. μc , μck and μk are considered as three design parameter variables. Besides, c, ck and k  are regarded as normal random variables with a standard deviation 10 because of manufacturing variability. Other deterministic parameters are selected as: A=1(cm2/cycle  m), b0=0.27, V=10(m/s), M=3.2633(kg  s2/cm) and m=0.8158(kg  s2/cm).

(25)Mind  C(d)=Z¨2‾=(πAV/m2)(μckμk+(M+m)μc2μk−1)s.t.  {Pr{g1(X)=1−(πAVmb0g2μk)((μckM+m−μcM)2+μc2Mm+μckμk2mM2)≤0}≤0.15Pr{g2(X)=1−7.6394(4000(Mg)−1.5μc−1)−1≤0}≤0.15Pr{g3(X)=1−0.5(Mg)1/2(μk2μckμc−1(M+m)−1+μc)−1/2≤0}≤0.15Pr{g4(X)=1−((M+m)g)0.877μck−1≤0}≤0.15380≤μc≤490,1430≤μck≤1530,10≤μk≤60

Mohsen et al. (2014). employed RIA, PMA, SLA, SORA and SBM to deal with the RBDO problem defined in Eq. (25), and the results are listed in Table 5. According to Eq. (25), the maximum and minimum values of the objective function C(d) by only considering the lower and upper bounds of the design parameter vector d can be easily computed as Cmax=2.74×108kg and Cmin=5.35×108kg. Hence, the RBDO problem defined in Eq. (25) can be rewritten as,

(26)Mind  C(d)−2.74×1082.61×108+∑k=142⋅I[Pr{gk(X)≤0}−0.15]s.t.     380≤μc≤490,1430≤μck≤1530,10≤μk≤60

By employing the proposed method to solve the above RBDO problem, the final optimal solution can be obtained as {μc=401.0263kg/cm,μck=1451.8459kg/cm,μk=31.5289kg/cms} with the maximum failure probability of 0.150, and the convergence of the design parameter variables and cost is drawn in Figure 6. From Figure 6, it can be observed that the convergence is achieved when gen=100. Besides, Table 5 indicates that the proposed method is more excellent and efficient in solving this RBDO problem because it can get a better solution from the feasible unit with lower cost and smaller number of model evaluations. μcμckμk

5. Concluding remarks

Starting with the premise that some existing approximately analytical technique based methods are lack of precision and some of the simulation technique based methods are not affordable for numerous complicated engineering problems, the aim of this contribution is to develop an accurate and efficient algorithm for solving the RBDO problem. The main contributions of the proposed algorithm are summarized as follows.

1. The original RBDO problem with simple objective function and complex constraints is transformed into the equivalent problem with simple constraints and complex objective function, which is more convenient to be settled by GA.

2. When solving the equivalent RBDO problem by using GA, the optimal reliability level is formulated in terms of the failure probability instead of the reliability index.

3. Among the process of GA, two Kriging models are constructed and adaptively updated by employing the U learning function in order to drastically improve the computational efficiency of GA.

The results obtained by employing the proposed technique in three numerical examples are comparable to those acquired by using the existing algorithms. From the results, it can be concluded that the proposed technique can produce superior optimal solutions with small number of model evaluations. However, it should be pointed out that the computational cost of building a Kriging model may be quite expensive for high dimensional models, thus the efficiency of the proposed algorithm will be reduced to some extent in dealing with such problems.

Figures

Figure 1

The flowchart of the proposed algorithm for solving the RBDO problem

Figure 2

Convergence of the design parameter vector and cost in example 1

Figure 3

Bracket structure

Figure 4

Convergence of the design parameter vector and cost in example 2

Figure 5

Passive vehicle suspension

Figure 6

Convergence of the design parameter vector and cost in example 3

Table 1

The parameters of GA and adaptive Kriging model used in test examples

Parameter	Symbol	Value
Maximum number of generation	genmax	200
Number of individuals in each population	Nd	500
Crossover probability	pc	0.7
Mutation probability	pm	0.1
Number of samples of the input vector in each population	NX	105
Penalty factor	PRk or PDl	2

Table 2

Optimization results for standard RBDO test problem

Method	Cost	Pfmax	Model evaluations
RIA*	6.7257	0.001349	590
PMA*	6.7251	0.001363	612
SLA*	6.7556	0.001351	144
SORA*	6.7251	0.001363	360
SBM*	6.7457	0.001296	5000
Proposed	6.7322	0.001343	39

Source(s): *Cited from Mohsen et al. (2014)

Table 3

Distribution types and parameters of the random input variables in example 2

Variable	Distribution type	Mean	C.o.V
σS(MPa)	Lognormal	225	0.08
P(kN)	Gumbel	100	0.15
E(GPa)	Gumbel	200	0.08
ρ(kg/m3)	Weibull	7860	0.10
L(m)	Gaussian	5	0.05
ωAB(mm)	Gaussian	μωAB	0.05
ωCD(mm)	Gaussian	μωCD	0.05
t(mm)	Gaussian	μt	0.05

Note(s): where C.o.V. denotes the Coefficient of Variance

Table 4

Optimization results for bracket structure

Method	Optimum Design(mm)	Cost(kg)	Pfmax	Model evaluations
RIA*	μωAB=61μωCD=157μt=209	1675	0.02027	1340
SORA*	μωAB=61μωCD=157μt=209	1675	0.02027	2340
SS*	μωAB=58μωCD=119μt=241	1550	0.02068	107
SS + Kriging*	μωAB=59μωCD=135μt=226	1610	0.02386	250
Proposed	μωAB=54μωCD=76μt=295	1359	0.02237	241

Source(s): *Cited from Chateauneuf and Aoues (2008)

Table 5

Optimization results for passive vehicle suspension

Method	Optimum design	Cost	Pfmax	Model evaluations
RIA*	μc=400.9787μck=1451.7168μk=31.6646	314,348,853.14	0.150	146
PMA*	μc=400.9784μck=1451.7168μk=31.5546	314,348,715.78	0.150	91
SLA*	μc=400.9787μck=1451.7168μk=31.5522	314,339,509.58	0.150	392
SLA*	μc=400.9787μck=1451.7168μk=31.5546	314,348,851.01	0.150	736
SBM*	μc=400.8658μck=1458.8673μk=31.2778	314,319,856.51	0.151	15,000
Proposed	μc=401.0263ck=1451.8459μk=31.5289	314,294,959.36	0.150	72

Source(s): *Cited from Mohsen et al. (2014)