Differential evolution approach to reliability-oriented optimal design

Probabilistic Engineering Mechanics 36 (2014) 72–80

Contents lists available at ScienceDirect

Probabilistic Engineering Mechanics journal homepage: www.elsevier.com/locate/probengmech

Differential evolution approach to reliability-oriented optimal design Sara Casciati Struttura Didattica Speciale di Architettura, University of Catania, Siracusa, Italy

art ic l e i nf o

a b s t r a c t

Article history: Received 31 January 2014 Accepted 5 March 2014 Available online 15 March 2014

The adoption of an evolutionary optimization approach, to tackle Reliability Based Design Optimization (RBDO) problems for which classical optimization is not feasible, can potentially lead to expand the use of RBDO in engineering applications. The herein proposed novel approach consists of coupling a differential evolution strategy with the one-level reliability method based on the Karush–Kuhn–Tucker (KKT) conditions. By selecting a few significant examples from the literature, convergence is found within a number of iterations compatible with the current capabilities of standard computational tools. These examples are chosen because lack of convergence is reported in literature when adopting classical gradient-based methods coupled with the reliability assessment. The performance and the efficiency of the algorithm are discussed, together with possible future improvements. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Reliability-based design optimization Moment based methods Nonlinear performance function Multiple limit states Differential evolution algorithm Heuristic optimization methods

1. Introduction Traditional design procedures of engineering systems are based on deterministic models and parameters. The variability of the loads, the material properties, the geometry, and the boundary conditions is included by means of modeling idealizations and simplifying assumptions, such as the consideration of mean or extreme values and the introduction of safety factors derived from past practice and experience. This approach is not able to realistically capture the influence of the uncertainties in the system parameters on the structural performance, which might result to be unsatisfactory. Furthermore, it cannot provide a quantitative measure of the system safety since it is unable to capture a mathematical relation with the risk assessments based on which decisional procedures are carried out [1]. Despite the development of a broad spectrum of different reliability-based design optimization (RBDO) methods in literature over the last 50 years (see, for example, [2]), their practical application to engineering problems is still rather limited in comparison with the deterministic design procedures. This is mainly due to the numerical challenges and computational cost often encountered when explicitly taking into account the effects of the unavoidable uncertainties in the structural performance. In Valdebenito and Schueller [3], the need to improve the current strategies for solving RBDO problems is stated, with particular emphasis on the aspects of numerical efficiency and robustness. The selection of a particular optimization algorithm can be crucial to solve a specific RBDO problem. Indeed, the numerical efficiency, accuracy and stability of the solving algorithm determine the adequacy of the method to practical engineering applications. When traditional gradient-based algorithms are adopted for the objective function minimization, the efficiency related to http://dx.doi.org/10.1016/j.probengmech.2014.03.001 0266-8920/& 2014 Elsevier Ltd. All rights reserved.

the moment based methods comes at the price of a limited field of applications, which does not include problems characterized by a large number of random variables and multiple failure criteria considering nonlinear performance functions. Hence, there is a need to re-address the topic in the case that a different solution strategy, such as the one based on evolutionary algorithms of differential type, is pursued. The herein proposed novel approach consists of coupling a differential evolution (DE) strategy [4] with a one-level reliability method based on the Karush–Kuhn–Tucker (KKT) conditions [5]. The adoption of a one-level approach based on the direct introduction of the KKT conditions in the formulation of the general cost minimization problem was proposed in Kuschel and Rackwitz [6] for a single limit state function. Mathematical proof that the solution point of the FORM optimization problem fulfills the KKT conditions is given. The existence of an optimal design point guarantees that the vectors of the Lagrangian multipliers associated to the KKT conditions also exists. Nevertheless, the dependence of the failure domain on the design parameters prevents one from mathematically proving any general convergence property of the method. This observation holds also for the two-steps RBDO approaches [7], where two nested optimization problems are formulated, with the inner loop dedicated to the reliability analysis and the outer loop to the cost optimization. Further adaptations were studied in order to handle problems with multiple independent failure modes and/or time-variant stationary loads [5]. The solution was sought by traditional gradientbased optimizers which required the computation of the gradient of the objective as well as the gradient of the constraints. In particular, taking the gradient of the second KKT optimality condition implicitly assumes the existence of the second order derivatives of the limit state functions. Such a requirement was identified as one of the

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

main drawbacks of the one-level approach, thus suggesting that a heuristic search algorithm would be particularly suited to overcome this difficulty. The application of the evolutionary strategies to RBDO problems has been studied mainly in association to a Monte Carlo simulation approach where the reliability calculations are carried out within a two-levels framework (see, for example, [8]). In the present work, the motivation of coupling a DE solution strategy to a one-level RBDO formulation is two-folds. (1) By reformulating the optimization problem so that the convergence is simultaneously sought in both the random variables space and the design parameters space, the reliability analysis is avoided thus potentially reducing the overall computational cost. (2) The DE approach can potentially broaden the reach of the RBDO method to those cases where multiple limit states and/ or strong nonlinearities prevent the traditional gradient-based optimizer from reaching convergence.

( s:t: :

73

P f i ¼ Pr½Gi ðd; XÞ r 0 rP Tfi ; hj ðdÞ r 0;

i ¼ 1; …; m

j ¼ m þ 1; …; M

ð1Þ

where C(d) is the objective function (or cost function) which can represent either the initial cost of the structure upon construction or the expected value of its life-cycle cost, including the expected lifetime and the operating costs; hj(d) are the deterministic constraints and usually refer to the upper and lower bounds of the design variables; m is the number of limit states (or performance functions [8]), and M is the total number of constraints. In the above formulation, the key issue is represented by the computational effort required to evaluate the probabilistic constraints which restrict the feasible design region by bounding the probability, Pfi, of violating each i-th limit state, Gi(d, X), i ¼1,…,m, to not exceed the given admissible failure probability, PTfi. Indeed, in the time-invariant case, each failure probability, Pfi, i ¼1,…,m, is defined by the integral Z Z ⋯ f X ðxjdÞdx ð2Þ P f i ¼ Pr½Gi ðd; XÞ r 0 ¼ Gi ðd;XÞ r 0

Evidence for these two statements is sought by applying the proposed methodology to a few numerical examples for which the unfeasibility of classical optimization is openly stated in literature. For the selected numerical examples, convergence is found within a number of iterations compatible with the current capabilities of standard computational tools. Proving any general conclusion upon the convergence of the method is out of the scope of this work. Despite the insights derived from experimentation and the good convergence properties of the DE algorithm empirically observed by Storn and Price [4] from extensive testing under various conditions, a rigorous mathematical convergence proof still lacks and could be the object of future research. At present, heuristic tools, as the DE scheme adopted in this paper, come as solution providers without a sound mathematical foundation. For this reason their potential is often denied, or, alternatively, the methods are just used on a trial and error basis. Nevertheless, their capabilities in extending the reach of existing methods and permitting the feasibility of new approaches are of interest for many engineering fields that involve the solution of optimization problems. The paper is structured as follows. In Section 2, the governing relations of the general RBDO problem are briefly provided. In Section 3, the strategy of using the DE technique to overcome the limitations of the one-level KKT approach is proposed. A formulation of the optimization problem suitable to this application is developed. In Section 4, the performance of the proposed solution technique is discussed with reference to numerical examples which are characterized by high nonlinearities due to both the analytical form of the limit state functions and the non-normality of the random variables.

2. Governing relations Let X be the 1
N vector of arbitrarily distributed random variables and let d be the 1
D vector of the design variables, which can be either defined as independent deterministic variables or as parameters of the marginal probability distributions of the X's. As an example, one considers the failure modes of a portal frame [5]; from a probabilistic point of view, the loads and the plastic hinge moments are assumed as the random variables in vector X, whereas the mean values of the beam and columns plastic moments are the design parameters in vector d. The basic formulation of a typical RBDO problem takes the form: min : d

CðdÞ

where fX(x) is the joint probability density function of the random variables X whose realizations are denoted with x. It is worth noting that the joint probability density function may depend on d in the cases in which some of the parameters of the probability distributions (e.g., the mean values) are considered as design variables. As mentioned above, in the two-level approach, the RBDO problem is tackled by formulating two nested optimization problems, where the inner loop is dedicated to the reliability analysis and the outer loop concerns the cost optimization. Since the exact computation of the integral in Eq. (2) is unpractical, approximate methods are applied and they can be either based on stochastic simulations or moment methods. The first-order reliability method (FORM) often provides an adequate accuracy and it is widely applied for the inner loop optimization. After transforming the random variables X into the uncorrelated and standardized normal variables, U¼U(X), whose realizations are denoted by u, the reliability index is computed by solving a constrained optimization problem: min : U

jjUjj;

s:t: :

Gi ðUÞ ¼ 0

ð3Þ

Typically, the HL-RF numerical method is adopted to carry out the FORM reliability analysis, because it requires no line search. The solution uni , i¼ 1, …, m, represents the most probable failure point for the considered i-th limit state, and the reliability index βi is given by || u ni ||. According to the FORM approximation, the failure probability can then be estimated by Pfi  Φ( βi), where Φ(.) is the standard Gaussian distribution. Alternatively, the probabilistic constraint in Eq. (1) can be reformulated in terms of the reliability index, which must not be less than a threshold value, βTi , given by Φ  1(1 P Tf i ). Then, the discussed two-loops method takes the name of the Reliability Index Approach (RIA) as reported in literature (see, for example, [9]). The adoption of a two-levels strategy implies the repeated evaluation of the performance functions based on the mechanical model of the structure. Therefore, it becomes unpractical when finite element models of real structures with material and/or geometrical nonlinearities are considered. The introduction of an explicit model, such as the response surface, to approximate the implicit model representative of the outcomes of the finite element analyses can largely accelerate the inner reliability loop [7], but the computational cost is still high due to the evaluation of the response surface terms and the introduction of approximation errors. Recent studies targeted to the development of efficient reliability estimation approaches to be used in the framework of

74

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

reliability-based optimization focus on the introduction of importance sampling techniques to reduce the computational effort of stochastic simulations (see, for example, [10]).

3. Using the DE strategy to solve the one-level KKT approach 3.1. The one-level KKT approach The aim of the one-level approach is to solve the RBDO problem by a single loop procedure. The probabilistic constraints are estimated by directly introducing in Eq. (1) some reliabilityequivalent optimality conditions, such as the ones expressed by the KKT null equalities of the first order reliability method in the time invariant case, provided that the reliability task is formulated in the transformed standard space [5]. In detail, the explicit solution of Eq. (3) is not necessary if the following KKT conditions are fulfilled at the β-point 8 > < Gi ðd; Ui Þ ¼ 0 fUi gk jj∇u Gi ðd; Ui Þjj þ jjUi jjf∇u Gi ðd; Ui Þgk ¼ 0 ð4Þ > : uL r U r uU ; i ¼ 1; …; m; k ¼ 1; …; n  1 i i where the Ui are the independent standard normal vectors of size ni each, with i¼ 1,…, m, and uU and uL are the upper and lower bounds for these variables, respectively. This approach requires the explicit implementation of the probabilistic transformation between the original random X-space and the U-space, as well as the computation of the derivatives of the limit state functions.

3.2. The objective function formulation

3.3. The DE solving algorithm Evolutionary techniques are, in general, capable to reach the optimal solution with simplicity in theory and easiness of programming, even when a highly nonlinear and partly nondifferentiable objective function which exhibits many local minima is considered. Whereas the gradient methods operate on a single potential solution and look for some improvements in its neighborhood, the global heuristic optimization techniques, represented herein by the so called evolutionary methods, maintain large sets (populations) of potential solutions and apply to them some recombination and selection operators. Let NP be the number of parameter vectors contained in each population. In the considered numerical examples, NP is empirically set equal to 80–160. By applying the procedure to a population of NP equally important (N þ D)-dimensional vectors, instead of referring only to a single nominal parameter vector, the dependence of the convergence path on the initial guess is significantly reduced. The differential evolution (DE) approach [4] is herein adopted as a robust statistical, parallel direct search method for cost function minimization. Besides its good convergence properties, the attractiveness of the method also consists of requiring only few input parameters which remain fixed throughout the entire optimization procedure, as reported in Table 1. The method described in the following is one of the several variants of DE which, however, differ only slightly. Its working scheme is schematically illustrated in Fig. 1. 3.3.1. Initialization The first step of the DE algorithm is to create an arbitrary initial population given by U L L yð0Þ r ¼ ξr ðy  y Þ þ y ;

The formulation of a pertinent objective function is crucial to the design process and the multiple objectives need to be adequately weighted [11]. Indeed, in the proposed approach, all the optimization tasks are combined into a single objective function which is computed as a weighted sum: Fðd; XÞ ¼ CðdÞ þ λðKKTÞ F ðKKTÞ ðd; XÞ þ λðβTÞ F ðβTÞ ðXÞ

m ni  1

i¼1k¼1

m ni  1

þ ∑ ∑ jjUi ðXi Þjj f∇u Gi ½d; Ui ðXi Þgk

ð6Þ

i¼1k¼1

and qffiffiffiffiffiffiffiffiffiffiffi F ðβTÞ ¼ βT  minð UTi Ui Þ i

ð7Þ

Such an approach can be naturally extended to include also a robustness requirement as discussed in Casciati [12] and Casciati and Faravelli [13]. The weighting factors λ(KKT) and λ(βT) in Eq. (5) are used to define the importance associated with the different objectives and constraints, as well as to normalize different physical units. In this work, their values are manually tuned on a case specific basis, as specified at the end of the next sub-section. Based on the objective function given in Eq. (5), the RBDO problem in Eq. (1) can then be restated as follows: min : y

FðyÞ;

s:t: :

yL r y ryU

where ξr is a vector of size (N þD) containing pseudorandom values drawn from the standard uniform distribution on the open interval [0 1]. From the first population forward, the following generation (P þ1) is created on the basis of the current population (P) by applying the breeding algorithm in Fig. 1. 3.3.2. Mutation According to the mutation scheme of DE algorithms, a trial parameter vector is built by adding the weighted difference between two randomly chosen population vectors to a third random vector as follows:

F ðKKTÞ ¼ ∑ jGi ðd; Ui Þj þ ∑ ∑ fUi ðXi Þgk jj∇u Gi ½d; Ui ðXi Þjj i¼1

ð9Þ

ð5Þ

where m

r ¼ 1; 2; …; NP

ð8Þ

with y ¼ ½d X and with yU and yL the corresponding upper and lower bounds, respectively.

ðPÞ ðPÞ þ 1Þ vðP ¼ yr1 þ Λ UðyðPÞ r r2  yr3 Þ;

r ¼ 1; 2; …; NP

ð10Þ

with r1, r2, and r3 denoting randomly selected integers such that r 1 ; r 2 ; r 3 A f1; 2; …NPg  Z þ and r 1 a r 2 a r 3 a r. The real positive constant Λ represents a control parameter of the DE-scheme since ðPÞ it governs the amplification of the differential variation ðyðPÞ r2  yr3 Þ. It þ is empirically chosen such that: Λ A ½0:5; 1  ℜ . Setting Λ ¼ 0:5, as reported in Table 1, leads to satisfactory results in all the considered numerical examples. With respect to other evolutionary strategies, the mutation is not performed based on some separately defined probability Table 1 The assigned DE control parameters. DE – control parameters

Value

Initial population size, NP Mutation amplitude, Λ Cross-over constant, CR Tolerance for convergence, ε Maximum iterations number, Imax

80–160 0.5 0.8 10  4–10  3 3000–10 000

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

75

Fig. 1. Block diagram of the differential evolution strategy.

density functions (PDF), but it is solely derived from the positional information of three randomly chosen distinct vectors in the current population. This scheme provides for automatic selfadaptation and eliminates the need to adapt the standard deviations of a PDF. 3.3.3. Cross-over As it is commonly done in the evolutionary strategies, a discrete recombination or cross-over is then introduced to increase the diversity of the new parameter vectors. The cross-over is performed as follows: 8 þ 1Þ < fvðP gl ifl ¼ hsiðN þ DÞ ; hs þ1iðN þ DÞ ; …; hs þ T 1iðN þ DÞ r ðP þ 1Þ fwr gl ¼ : fyðPÞ otherwise r gl r ¼ 1; 2; …; NP;

l ¼ 1; 2; …; ðN þ DÞ

ð11Þ

In Eq. (11), the acute brackets denote the modulo function with modulus (NþD); the starting index s is a randomly chosen integer from the interval [0, NþD 1]; the integer T is drawn from the interval [1, NþD] with cumulative distribution function PrðT Z tÞ ¼ CRðt  1Þ . The real positive constant CR is another control parameter of the DE-scheme since it represents the cross-over probability. Setting CR ¼ 0:8, as reported in Table 1, leads to satisfactory results in all the considered numerical examples. The random selections for both s and T are again conducted for each þ 1Þ trial vector vðP . As a result of the cross-over operation, only a r certain sequence of the vector elements of wrðP þ 1Þ are identical to þ 1Þ the elements of vrðP þ 1Þ ; the other elements of wðP acquire the r ðPÞ original values of yr . 3.3.4. Constraints check þ 1Þ The vectors wðP are checked to satisfy the box-type conr straints of Eq. (8). 3.3.5. Selection The newly created vectors are compared with the members of the current population and they are chosen as their replacements

in the next generation only when they lead to an improvement of the objective function value toward the search direction. Hence, the population of the next generation is formed as follows: ( ðP þ 1Þ þ 1Þ wr ifFðwðP Þ o FðyrðPÞ Þ r þ 1Þ yðP ¼ ; r ¼ 1; 2; …; NP ð12Þ r otherwise yðPÞ r

3.3.6. Termination criterion The described procedure is repeated until either the maximum number of iteration Imax is met, or the distance of the new generation from the previous one falls below a given tolerance ε (see Table 1). Such a distance is measured as follows: max

r ¼ 1;::;NP

þ 1Þ jjyðP  yðPÞ r r jj

jjyrðPÞ jj

oε

ð13Þ

where ||.|| indicates the L2-norm. The performance of the described algorithm is assessed by Storn and Price [4] on a rich testbed with highly nonlinear, nondifferentiable and multimodal test functions of medium to high dimensionality. The gathered experience is used to derive empirical rules for setting the DE control parameters, thus permitting easiness of use which is one of the DE’s major assets. In particular, a large value (below unity) of the mutation control parameter, Λ, is shown to avoid that convergence is reached prematurely, and a large value of the cross-over control parameter, CR, often speeds convergence. Indeed, mutation leads fine tuning, cross-over introduces far potential solutions. The values Λ ¼ 0:5 and CR ¼ 0:8 are indicated as good initial choices by Storn and Price [4] and they correspond to a prudent (i.e., modest) shift from the current position, while higher values would move farther with the risk of missing nearby optimal points. The author’s strategy consists of fixing these two values for all the analyses, and varying the vectors population size, NP, and accordingly the maximum number of iterations, Imax, in order to determine the weighting factors λ(KKT) and λ(βT) in Eq. (5). First, one checks if a fast solution is possible by tuning the weighting factors so that convergence is met for a very small NP within a corresponding low number of iterations, Imax.

76

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

The extremant feature of the achieved solution is then investigated by slightly increasing the maximum number of iterations, Imax, and running a further analysis for an enlarged population and the previously identified values of the weighting factors. If the optimal solution is still possible and consistent with the previous one, then it is likely that the algorithm is robust and has not prematurely converged on a local minimum. In other words, after a preliminary assessment of the functions in (6) and (7) in a restricted set of the space of the design variables, the factors are selected in such a way that a perturbation on the design variables produces comparable differences on the two functions. It is worth noting that, when the formulation of the objective function in Eq. (5) is adopted, the optimal value is known and is null for both the weighted terms which express the safety requirement and the KKT conditions as defined by Eqs. (6) and (7), respectively. Furthermore, both terms assume positive values when the search is performed in the safe region of the design space. These observations can be usefully integrated into the solving algorithm to further improve its performance and convergence property as follows. (1) If the sign of either limit state function in Eq. (6) is negative at any point of the search domain along the solution process, then the value of the objective function evaluated at this point is augmented of an extra-cost value so high that such vector becomes undesirable and is disregarded at the next step of the procedure; in this manner, convergence on local minima in the unsafe domain is avoided. (2) If a significant difference between the optimal value of the objective function in Eq. (5) and the cost C(d) is detected, then the selected weighting factors lead to unsatisfactory results in terms of the unfulfillment of the associated constrains and need to be re-calibrated. This estimate is separately implemented at the end of each run during the above described tuning procedure of the weighting factors in Eq. (5). (3) If a high number of iterations is needed to reach convergence, it can be reduced by adopting the optional technique of progressively shrinking the search domain. This technique consists of defining a sub-domain smaller than the previous one around the current optimal point each time that a sequence of an assigned low number of iterations has been concluded. The aspects listed above are clarified and discussed through numerical examples in the next section. They eventually lead to a modified version of the original DE algorithm proposed by Storn and Price [4] which is targeted to the objective function formulated in Eq. (5). The resulting scheme is the one shown in Fig. 1.

4. Numerical examples As mentioned in the Introduction, the proposed approach coupling the DE optimizer with the one-level RBDO problem based on the KKT equality conditions of Eq. (4) is targeted to those situations where classical optimization techniques become unfeasible. Indeed, for such cases, the adoption of an evolutionary strategy is justified even if it results in an increased number of function evaluations, as long as the running time is still moderate when using standard computational tools. Therefore, rather than comparing the numerical performance of the proposed method with the one shown by other algorithms, in the present work the focus is placed on investigating its capability to broaden the reach of the classical RBDO concepts. Only few examples are reported in literature where a failure of convergence is openly stated by the authors. Based on this

consideration, two case-studies are selected among the mathematical problems proposed as benchmark studies in Aoues and Chateauneuf [2]. They are characterized by different degrees of nonlinearity in the objective function. First, only a single limit state and normally distributed random variables are considered. Then, the influence of multiple limit states and not normally distributed random variables is investigated. It is worth noting that, in the first example, the solution procedure based on the differential evolution technique becomes rather straightforward, thus completely eliminating the challenges encountered when using the traditional gradient based methods. Instead, the second example represents the core of this paper since it is characterized by significantly more numerical difficulties. Indeed, several counteractions need to be taken along the optimization process in order to avoid local minima and to limit the computational burden. For the sake of completeness, the conclusions drawn in Aoues and Chateauneuf [2] based on the results achieved within a traditional optimization framework that applies a sequential quadratic programming (SQP) algorithm to solve the RBDO problem of Eq. (1) are also summarized in the following paragraphs. The intent is to enable the reader to better understand and evaluate the advantages introduced by the proposed DE-KKT method. 4.1. Problem I: limit state nonlinearity 4.1.1. Problem statement A first problem is considered where a single nonlinear limit state function depends upon two design parameters of deterministic nature d ¼ ½ d1 d2  and two independent normally distributed random variables X ¼ ½ X 1 X 2  whose probabilistic definition is given in Table 2. An admissible probability of failure of 1% is specified (corresponding to βT ¼ 2.32), and the optimization problem of Eq. (1) is solved under the assumptions summarized in Fig. 2. Namely, two different formulations of a single limit state function are considered and denoted as Cases A and B, respectively. They correspond to an increasing level of nonlinearity, being the limit state function of Case B characterized by high curvatures in the neighborhood of the optimum.

4.1.2. Results from SQP-KKT and SQP-RIA The results in Aoues and Chateauneuf [2] yield to the observation that SQP-KKT is very sensitive to the choice of the initial point. When the initial point lies in the neighborhood of the optimum, SQP-KKT shows a high convergence rate. For the limit state 0 function of Case A and a starting point d ¼ ½ 2 2 , the convern gence is reached at the optimal point d ¼ ½ 5:65 5:65 , corren sponding to a minimum cost Cðd Þ ¼ 63:88 after performing a number of limit state evaluations equal to 9. For the limit state 0 function of Case B and the same starting point d ¼ ½ 2 2 , the n convergence is reached at the optimal point d ¼ ½ 1:35 1:35  n corresponding to a minimum cost Cðd Þ ¼ 3:67, after performing a number of limit state evaluations equal to 6. It is worth noting that when using a two-levels approach based on SQP-RIA to solve the same RBDO problems, a higher number of limit state evaluations equal to 31 and 14, respectively, is reported. Table 2 Probability distributions for Problem I. Random Variables

Distribution

Mean

Coefficient of Variation

X1 X2

Normal Normal

5 3

0.3 0.3

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

77

− Vector of design variables :

{

}

d∈ Sd := d∈ ℜ2 | 0 ≤ d ≤ 15 ⊂ ℜ2 − Cost minimization target : 2

2

C (d) = d1 + d2

− Probabilistic target : PfT = 1%; β T = 2.32 − Nonlinear limit state function (m = 1) : - Case A : G(d, X) =

1

2

d1d2 X2 − X1 5 - Case B : G(d, X) = d1d2 X2 − ln( X1)

Fig. 3. Problem II assignments.

Table 4 Probability distributions for Problem II.

Fig. 2. Problem I assignments.

Table 3 DE-KKT results for Problem I. Case

λ(KKT) ¼ λ(βT)

dn1 ¼ dn2

Cost

Iterations

A B B

200 100 200

5.62 1.35 1.35

63.1 3.65 3.65

544 389 451

Random variables

Distribution

Mean

Coefficient of variation

X1 X2

Normal/Gumbel Normal/Gumbel

d1 d2

0.12 0.12

additional computational burden with respect to the quadratic polynomial form defining the limit state function of case A. 4.2. Problem II: multiple limit states

When the initial point is far from the optimum, numerical instabilities are observed when using the SQP-KKT method. For 0 the limit state function of Case A and a starting point d ¼ ½ 12 12 , the convergence is still reached with an increased number of limit state evaluations equal to 13, which remains lower than the one reported for SQP-RIA (equal to 40). However, increased values of the cost functions are observed during the first iterations of SQP-KKT, due to inaccurate computation of the descent direction. For the limit state function of Case B and the 0 same starting point d ¼ ½ 12 12 , both SQP-KKT and SQP-RIA fail to converge.

4.1.3. Results from DE-KKT By adopting the differential evolution strategy, the sensitivity of the KKT approach to the initial guess is no longer observed. Indeed, the results in Table 3 are retrieved independently of the choice of the initial point, which can either be located in the vicinity of the optimal point, e.g. d0 ¼[5 5], or far away from it, e.g. d0 ¼[12 12], without affecting the convergence path of the algorithm to the optimal point. This result is particularly significant when compared to the lack of convergence reported for Case B in the benchmark study, when the SQP technique is applied to both the KKT and RIA approaches starting from the farthest initial point. The values of the two lambda factors in Eq. (5) are the results of few replicates of the whole analysis. At the end of each step, one quantifies the error in the associated constraints and, if necessary, increases the value of the factors to be used in the next analysis until the two associated sets of constraints are actually satisfied as equalities. Their final values are reported in Table 3 and they are satisfactorily assumed equal to each other for the considered casestudies. When comparing the number of iterations in Table 3, it is worth noting that this number depends on the values selected for the Lagrangian multipliers, which influence the convergence accuracy. Therefore, one cannot say that Case B is solved in a faster way than Case A, but only that the computational effort is similar. Hence, it can only be concluded that the increase of the limit state nonlinearity induced by the introduction of a logarithmic term in case B (see Fig. 2 and Table 3) does not result in any

4.2.1. Problem statement In the second considered problem, three nonlinear performance functions are defined by the analytical expressions given in Fig. 3. The design parameters d ¼ ½ μX1 μX2  represent the optimal mean values of the two random variables X ¼ ½ X 1 X 2  that appear in each of these functions. Different assumptions are made on the type of the probability distributions of the random variables, as reported in Table 4. Their coefficient of variation is fixed to 0.12. 4.2.2. Results from SQP-KKT and SQP-RIA The mathematical example defined in Fig. 3 was originally proposed by Youn and Choi [14] to investigate the influence of nonlinearities on the efficiency and robustness of the RBDO process. In particular, the dependence of the traditional twolevels approach using SQP-RIA on the distribution type of the random variables is emphasized. The performance failure surfaces Gi(d, Xi), i ¼1,…,m, are originally described in the X-space and they are themselves nonlinear. The reliability constraints in RIA, as well as the equivalent KKT optimality conditions in Eq. (4), are based on the performance failure surfaces Gi(d, Ui), i¼1,…,m in the U-space of independent standard normal variables. Hence, a preliminary transformation in the U-space, U i ¼U i (X i), i¼ 1,…,m, is required and it may introduce additional nonlinearity in the original failure surface depending on the given types of probability distributions. Most transformations required in RBDO are, indeed, highly nonlinear, except when the original random variables are normally distributed. In the case of non-normally distributed random variables, the nonlinearity of the RBDO problem therefore increases due to the transformations that take place several times, along the two-levels RIA optimization process, between the original random X-space, where the design optimization is carried out, and the random U-space of independent standard normal variables, where the reliability analysis is performed. In particular, the assumption of a Gumbel distribution results in a probabilistic constraint with a high skew, which may produce numerical instability in the RBDO process. When applying SQP-RIA for Gumbel distributions, a failure of convergence is reported and it is due to the numerical

78

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

Fig. 4. Graphical representation of Problem II with βT ¼2, under the assumptions that the variables are normally distributed and the optimal solution belongs to the quadrant bisector.

Table 5 DE-KKT results for Problem II. βTmin

Distribution

d1

d2

Cost

Iter.

4 3 2 4 3 2

Normal Normal Normal Gumbel Gumbel Gumbel

n.a. 3.6382 3.3239 3.3024 3.1419 3.0445

n.a. 3.9358 3.2910 3.9831 3.8179 3.4585

n.a. 7.574 6.6291 7.2855 6.9598 6.503

n.a. 843 930 6541 7634 8511

difficulties encountered when having either a nonlinear failure surface or a failure surface away from a design point. Instead, for normally distributed X-variables with standard deviation equal to 0.6 and a target reliability index equal to 2, the convergence to the optimal point dn ¼[3.61 3.66] (corresponding to a cost C(dn)¼7.27) is reached by SQP-RIA with a number of function evaluations equal to 53 when d0 ¼ [5 5] is selected as initial point. In Aoues and Chateauneuf [2], the numerical example is revisited by considering a coefficient of variation of the original random variables equal to 0.6 and three target reliability indices βT ¼2,3,4. Therefore, the results from the benchmark study are not directly comparable to the previous ones due to the influence of the different c.o.v. assumptions. When the random variables are normally distributed, the authors claim that all the considered RBDO methods with the same initial point d0 ¼[5 5] converge to nearly the same optimum which corresponds to a cost value of 6.2, 8.9, and 7.3 for the different reliability targets, respectively. Nevertheless, numerical difficulties and increased computational effort are reported when applying SQP-KKT to multiple limit states. When the assumption of Gumbel distributions increases the limit state nonlinearities due to the probabilistic transformations, both the classical approaches based on SQP-RIA and SQP-KKT fail to converge for all target reliability indices.

4.2.3. Results from DE-KKT The results in Table 5 are obtained by adopting the DE algorithm within the one-level KKT approach. It is worth noting that these results refer to the probabilistic assumptions in Table 4 (c.o.v.¼ 0.12), which coincide with the ones of the original problem formulation by Youn and Choi [14]. In the case of normally distributed random variables, a low c.o.v. value ensures that the mathematical problem stated in Fig. 3 is well-posed at least for the target reliability indices βT ¼ 2,3. A better understanding of this statement is achieved by plotting the three limit

state functions in the space of the design variables, as depicted in Figs. 4 for βT ¼2 and a solution vector xn characterized by equal entries. When the random variables are assumed as normally distributed, a circle with center in the design point un and radius equal to the reliability index βT in the standard normal space will transform in a circle centered in the optimal point xn with radius given by a function of βT which increases with the c.o.v. value. Hence, for the circle in the design space to reside within the safe region bounded by all the three limit state functions so that the multiple performance criteria defined in Fig. 3 are simultaneously satisfied, it is favorable to adopt a low c.o.v. value especially when high reliability targets are considered. Such a measure results, however, to be insufficient for the case of βT ¼4 (see Table 5), because such a high reliability target corresponds to a diameter of the circle which exceeds the linear dimension of the safe region measured along the quadrant bisector, even if a low c.o.v. value of 0.12 is assumed. In other words, for the highest value of the target reliability index, no feasible solution exists under the assumptions of normally distributed random variables given in Table 4. For the reasons specified in the following, it would be desirable that such an ill-posed problem could be autonomously detected by the solving algorithm without the need of a priori knowledge of the safe region. Such a desirable feature is achieved by performing a simple check of the fulfillment of the problem constraints. Indeed, the objective function in Eq. (5) is the sum of three terms and only when its value is nearly identical to the one assumed by the cost function C(d), one has evidence that the constraints are satisfied. If this is not the case, more suitable values of the Lagrangian multipliers in Eq. (5) should be tested. Eventually, if a significant difference between the objective function in Eq. (5) and the cost function C(d) remains, one can conclude that the problem does not possess any viable solution (ill-posed problem). Therefore, the results in Table 5 not only demonstrate the reliability offered by the differential evolution algorithm in achieving the optimal solution, but they also emphasize its ability to detect ill-posed problems. As it can be detected by observing the position of the circle in Fig. 4, the optimal point is constrained to belong to a region where all the three requirements imposed by the limit state functions are fulfilled. The first one is associated to large values of X1, the second with large values of X2, whereas the third constraint imposes low values for both the random variables. Due to these competing requirements on the coordinates of the optimal point, the search domain needs to span over the entire range (0, 10) defined for each variable, but this circumstance has non trivial consequences on the search process. Since the points belonging to the unsafe region are also examined as optimal candidates by the DE algorithm, there is a risk that the search process ends in some local minima located in such an unfeasible area. To avoid this abnormal end, the following counteraction was adopted. At any point in the design space which is a candidate to represent the mean of the random variables, one also computes the sign of each limit state function in Eq. (6). If one of these signs is negative, then a fictitious extra-cost is added to the objective function, thus making this point undesirable within the solution search. For the highest value of the target reliability index, βT ¼ 4, and under the normality assumption, two couples of Lagrangian multipliers were tried in Eq. (5) to represent situations where either the KKT constraint dominates or the safety constrain does, respectively: λðKKTÞ ¼ 100 000;

λðβTÞ ¼ 10 000

ð25Þ

λðKKTÞ ¼ 100 000;

λðβTÞ ¼ 100 000

ð26Þ

The first attempt results in a solution violating the safety constraint (i.e., the resulting beta circle has a radius lower than 4). The

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

second couple of Lagrangian multiplier leads to a circle with a big enough radius, but the KKT constraint is violated, thus indicating that the retrieved radius does not represent the distance from the mean to the zero-valued curves associated to the performance functions, but to other parametric curves where the performance functions (or at least one of them) assume negative values. Therefore, the unfeasibility of the solution is detected. The advantage of the proposed method with respect to the traditional one-level and two-levels RBDO approaches is evident by observing, from the results in Table 5, that convergence is possible even when the Gumbel assumption is introduced. In this case, the circle in the standard normal space is distorted following the transformation to the original design space and the mean point is located closer to the lower zero-valued curves associated to the performance functions G1(X) and G2(X). This situation enables to satisfy the requirement of the minimum admissible reliability index equal to 4 without violating the KKT conditions, when the Lagrangian multipliers in Eq. (5) are set as specified in Eq. (25). For lower values of the admissible reliability index, an increasing number of iterations is needed in order to reach convergence. This trend is confirmed by the results in Table 5 which are achieved under the normality assumption with βT ¼3, 2 and the Lagrangian multipliers given by Eqs. (26) and (25), respectively. The number of iterations is drastically reduced by adopting the numerical technique of a progressive shrinkage of the search domain. In other words, the search on the whole domain is stopped after few hundreds of iterations and a sub-domain smaller than the initial search domain is built around the current optimal point. As the search domain shrinks, a lower number of iterations is required to achieve the convergence. Once this is achieved, a further run of the algorithm with a narrow search domain limited to the surrounding of the optimal point can be useful to confirm that the procedure did not end on a local minimum. 4.3. Discussion Reliability assessment and reliability based design optimization are deeply influenced by the performance of the numerical optimization tool. In the last decade, many research efforts were devoted to build a sort of system architecture in which these general topics could be concisely organized. Recently, benchmarks have been formulated with the attempt of giving a systematic arrangement to the topic. However, most of the available conclusions hold only for the adopted optimization algorithm. The selection of a particular optimization algorithm can be crucial for solving a particular RBDO problem. The introduction of efficient strategies and approximation concepts, such as the ones related to the stochastic search techniques, plays a fundamental role in yielding challenging RBO problems tractable. Indeed, although the gradient-based optimization methods are in general rather efficient, there are cases in which they fail in determining the optimal solution, the estimation of the gradients becomes a difficult issue, and/or the process ends finding local optima only. On the contrary, stochastic search algorithms such as the evolutionary techniques of differential type are shown, in this paper, to be able to find the global optimum of specific RBDO problems. Moreover, this class algorithms presents the advantage of not relying on gradient information and can be easily driven away from local optima by adjusting the exploration of the search domain to the specifically considered shape of the objective function. Despite the fact that stochastic search algorithms usually require many more objective function evaluations than gradient based methods, the numerical effort associated with each evaluation is usually negligible. Furthermore, the advent of high performance computing and the application of parallel computing techniques

79

provide the capability of performing numerical tasks in reduced time. The results achieved by applying the DE technique to solve RBDO problems show that the required number of iterations is consistent with the speed of current elaborators.

5. Conclusions In this paper, it is shown that, for the selected numerical examples, the adoption of a differential evolution strategy is quite reliable in terms of finding the solution when the uncertainties are explicitly taken into account in a structural optimization problem. Whereas the application of differential evolution algorithms presents the advantage of requiring only few input parameters, the key issues are represented by the formulation of an objective function with a shape that facilitates the convergence process and the proper calibration of the Lagrangian multipliers associated to different objectives so that the process is driven toward a global minimum. In some cases, several counteractions need to be taken along the optimization process in order to avoid local minima and to limit the computational burden. Once these issues are taken into account, the procedure is quite straightforward. Such an advantageous feature may enable the improvement of the robustness of the strategies for solving RBDO problems, in the sense that they can become applicable to a wide spectrum of problems that can be found in engineering. In particular, the high computational capabilities of the DE approach are desirable when problems characterized by strong nonlinearities are considered. Two further improvements are considered as possible themes for future work: (1) the adoption of response surface approximations for the performance functions to avoid large finite element computations; (2) the fragmentation of the design space in sub-domains for subsequent optimum search in order to perform the equivalent of importance sampling in standard structural reliability, and, hence, to decrease the number of iterations to convergence. Of course in these design sub-domains the accuracy of an overall response surface approximation would result highly improved.

Acknowledgments This research was supported by grants from the University of Catania (PRA 2008). The author is grateful to the anonymous referees for their constructive remarks. References [1] Faber MH, Stewart MG. Risk assessment for civil engineering facilities: critical overview and discussion. Reliab Eng Syst Saf 2003;80:173–84. [2] Aoues Y, Chateauneuf A. Benchmark study for reliability-based design optimization. Struct Multidiscipl Optim 2010;41:277–94. [3] Valdebenito MA, Schueller GI. A survey on approaches for reliability based optimization. Struct Multidiscipl Optim 2010;42:645–63. [4] Storn R, Price K. Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;11:341–59. [5] Rackwitz R. Optimization and risk acceptability based on the life quality index. Struct Saf 2002;24:297–331. [6] Kuschel N, Rackwitz R. Two basic problems in reliability-based structural. optimization. Math Methods Oper Res 1997;46:309–33. [7] Gasser M, Schueller GI. Reliability-based optimization of structural systems. Math Methods Oper Res 1997;46(3):287–307. [8] Papadrakakis M, Lagaros ND, Plevris V. Design optimization of steel structures considering uncertainties. Eng Struct 2005;27:1408–18. [9] Tu J, Choi K, Park Y. Design potential method for robust system parameter design. AIAA J 2001;39(4):667–77.

80

S. Casciati / Probabilistic Engineering Mechanics 36 (2014) 72–80

[10] Beaurepaire P, Jensen HA, Schuëller GI, Valdebenito MA. Reliability-based optimization using bridge importance sampling. Probab Eng Mech 2013;34: 48–57. [11] Casciati S. Stiffness identification and damage localization via differential evolution algorithms. J Struct Control Health Monit 2008;15(3):436–49. [12] Casciati, S. Including structural robustness in reliability-oriented optimal design. In: George Deodatis, Pol D. Spanos (editors), Computational stochastic

mechanics; Proceedings of the CSM-5 fifth international conference, Rhodes, 21–23 June, 2006. Rotterdam: Millpress; 2007. p. 157–61. [13] Casciati S, Faravelli L. Building a robustness index. In: Michael Faber (editor) Action TU0601 robustness of structures; Proceedings of the 1st workshop, ETH Zurich, 4–5 February 2008. Zurich: Reprozentrale ETH Hönggerberg; 2008. p. 49–56. [14] Youn BD, Choi KK. An investigation of nonlinearity of reliability based design optimization approaches. J Mech Des 2004;126(3):403–11.