FORM reliability analysis using a parallel evolutionary algorithm

Structural Safety 65 (2017) 84–99

Contents lists available at ScienceDirect

Structural Safety journal homepage: www.elsevier.com/locate/strusafe

FORM reliability analysis using a parallel evolutionary algorithm Dorival M. Pedroso School of Civil Engineering, The University of Queensland, St Lucia, QLD 4072, Australia

a r t i c l e

i n f o

Article history: Received 19 March 2016 Received in revised form 2 January 2017 Accepted 2 January 2017

Keywords: Nonlinear optimisation problems Finite element method Differential evolution Constraints handling Code reliability

a b s t r a c t This paper presents a parallel evolutionary algorithm to solve reliability problems with accuracy and repeatability of results. The last characteristic is usually overlooked; however, it is critical to the reliability of the calculation method itself. Note that evolutionary algorithms are stochastic processes and may not always generate identical results. The optimisation problem resulting from the first order reliability method is considered with an implicit state function that can include a call to a finite element analysis (FEA). A strategy to handle failures from the transformation of random variables or from the finite element call during the evolution process is explained in detail. Several benchmark tests are studied, including some involving bounded random variables that introduce strong non-linearities in the mapping to standard Gaussian space. In addition, the solutions of 2D and 3D frame problems using the finite element method illustrate the capabilities of the algorithm including the convenience of the algorithm in handling discrete limit state functions. Finally, the ability to obtain similar results after many runs is demonstrated. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The evaluation of reliability in engineering has indeed secured its place in the design and risk analysis of structures. This happened especially after the dissemination of the concept of riskbased design which has been adopted in a number of codes and standards; see e.g. discussion in [1,2]. This situation is a natural consequence of materials, loads and any other external agents being always uncertain, even within a very small margin. In other words, perfection in manufacturing or exactness in loads quantification is only achievable in a stochastic manner and hence uncertainty factors should accompany every specification. Among popular methods to assess reliability, the first and second order reliability methods (FORM and SORM) became good options due to their simplicity and efficiency. Other alternatives include Monte Carlo simulations that require a greater number of calculations. It is also known that FORM and SORM are able to produce reasonable results in terms of accuracy. Therefore, their acceptance and use have grown widely, particularly in the recent years. This is in part thanks to the improvement of numerical methods to approximate the solution of mechanical problems and the increase of computer power to allow faster and larger computations. This paper considers the time-invariant first order reliability method (FORM) only [3–7]. In this method, a constrained E-mail address: [email protected] http://dx.doi.org/10.1016/j.strusafe.2017.01.001 0167-4730/Ó 2017 Elsevier Ltd. All rights reserved.

optimisation problem is deduced for which several techniques have been developed along the years to obtain its solution. This leads to the quantification of the so-called reliability index b. Iterative solutions [8–13] and stochastic ones [14–21] are available but are still being proposed from time to time since challenges yet exist. For instance, derivatives cannot be easily computed in large complex problems involving other numerical methods such as the finite difference or finite element; these derivatives may not be even defined. Moreover, the satisfaction of constraints results in dealing with computer numeric precision which in turn poses difficulties to the solution algorithm. This paper thus presents an approach based on an evolutionary algorithm aiming at solving FORM with accuracy, efficiency and robustness. Evolutionary algorithms (EA) are based on a progressive update of an initial set of trial solutions (population) by recombining these solutions in a stochastic manner. The original algorithms were inspired by natural selection and evolution where the recombination is made by mimicking genetics; see e.g. [22–25]. There are actually several EAs and related analogies to solve optimisation problems. In addition to overcoming challenges with local optima, EAs have the advantageous capability of generating multiple solutions at the same time. This feature is particularly essential to multimodal problems. One subset of EAs with a promising applicability in reliability analysis is the so-called differential evolution (DE) [26–28]. DEs uses vector differences to recombine solutions instead of genetic analogs.

85

D.M. Pedroso / Structural Safety 65 (2017) 84–99

It was shown in [29] that the adaptation of differential evolution in the framework of an evolutionary algorithm based on tournaments results in a very efficient algorithm. This algorithm is able, for instance, to produce the same results after several runs; i.e. with a good repeatability behaviour. Note that this feature is not guaranteed by EAs because of their stochastic nature. In addition, the algorithm is able to obtain exact solutions for some particular multi-constrained and multi-objective optimisation problems. The algorithm and related computer code, named Goga, has been made available on the web at https://github.com/cpmech/goga as open (free) software. In this paper, the aforementioned evolutionary algorithm (Goga) is applied to the solution of the optimisation problem in FORM. In particular, the handling of constraints and the strategy to deal with failures during the computation of the objective value are described in detail. The failures may happen due to the transformation of random variables or due to the inability of a solver such as the finite element method to solve the underlying mechanical problem; for instance when the input (random) parameters are invalid. Several solutions to existent reliability problems having a reference solution are performed in order to investigate the repeatability characteristics of Goga in addition to assess its accuracy. Some challenging reliability problems [30] involving bounded random variables are analysed and 2D/3D frames are studied. It is shown that the EA is fairly robust in dealing with these non-linear problems and with discrete limit state functions. The paper is organised as follows. A very brief review of FORM is given in Section 2 followed by the description of the evolutionary algorithm and its specialisation to FORM in Section 3. In Section 4, the numerical examples are presented and in Section 5 some conclusions are drawn.

Fig. 1. Resistance-load paradigm in probabilistic analysis.

2. Probabilistic analysis and reliability Reliability is a measure of satisfactory performance of a system (or component) with uncertain characteristics; e.g. geometry, properties, boundary conditions, etc. Therefore, reliability is the opposite of probability of failure. In fact, the probability of failure is a more general quantity that can be readily related to risk when an evaluation of consequences is available. Failure in this context is a generic word that can mean mechanical breakage, inadequate serviceability, or infeasibility, for example. The component reliability is considered in the following. For details on system reliability, the reader is directed to the cited references; e.g. [7]. In a simplified manner, the probability of failure of a component described by two indicators R (resistance or supply) and L (load or demand) can be written as

pf ¼ PR ½failure ¼ PR½R < L

ð1Þ

The paradigm resistance-load is illustrated in Fig. 1 where the area under the intersection between the two curves is related to the probability of failure; hence safety cannot be simply expressed by the differences between the mean values as in deterministic approaches. The difference ðR  LÞ is known as the limit state function (LSF); clearly, negative values indicate failure. A more general LSF can be represented by gðx Þ where x is the vector of n random variables: x ¼ fx0 ; x1 ; . . . ; xn1 g (Fig. 2). The hyper-surface gðx Þ ¼ 0 thus delimits the regions of values leading to failure from the regions of values representing a safe situation. The probability of failure in this case can be calculated by means of

Fig. 2. Limit state function.

Z pf ¼

gðxÞ0 

f x ðx Þ dx

ð2Þ

where f x is the joint probability density function (PDF) of x . The above integral is often difficult to be obtained by analytical means; hence the FORM approximation procedure has been developed. In this method, first, a transformation T is introduced to convert the space of x into an uncorrelated normal space y; this can be symbol ised as

yðx Þ ¼ T ðx Þ 

ð3Þ

where the components of y are statistically independent normal  variables with zero means and unit standard deviations; see e.g. [31,32]. The cases of correlated and uncorrelated variables are implicitly represented by T . Further, the transformation considered in this paper is a direct mapping to the standard Gaussian space. For uncorrelated variables,

yi ¼ U1 ðF xi ðxi ÞÞ

ð4Þ

86

D.M. Pedroso / Structural Safety 65 (2017) 84–99

with UðxÞ being the standard normal cumulative distribution function and F xi being the cumulative distribution function of xi . For correlated variables,

yi ¼

X 1 1 Lij U ðF xj ðxj ÞÞ

ð5Þ

j

where Lij are the components of the lower-triangular matrix resulting from the Cholesky decomposition of the correlation matrix [31]. Hereafter, the following simplified notation is considered

gðx Þ ¼ gðx ðyÞÞ ¼ gðyÞ 

ð6Þ



where the dependency of y on x and the inverse relation is to be  understood from the context. Second, the LSF is linearised via Taylor series after which only the first order terms are retained. Taking advantage of the rotational symmetry of the standard normal joint PDF for y and its  exponential decay, the problem of solving the integral in Eq. (2) is then replaced by the following constrained optimisation problem

minimisey 

qffiffiffiffiffiffiffiffiffi yy 



ð7Þ

subject to gðx ðyÞÞ ¼ 0 

The solution of this problem leads to the so-called design point yH . The Euclidean norm of yH corresponds to the reliability index 

b¼



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi yH  yH 

ð8Þ



because it is the smallest distance from the origin to the LSF in the transformed space. Note that in Eq. (7), the inverse transformation x ¼ T 1 ðyÞ is required. 

Most algorithms (e.g. [3,4,6,7]) directly operate on the transformed space of y in order to solve Eq. (7). Nonetheless, for an evo lutionary algorithm, it is more convenient to operate on the original space x . In this way, the inverse transformation is never needed. Furthermore, since the satisfaction of equality constraints is more difficult than the satisfaction of inequality constraints, especially in EAs, an inequality-constrained version of the optimisation problem is selected instead. This is possible as long as gðy ¼ 0 Þ > 0; i.e. the LSF is positive at the origin of the transformed  space y [33,34]; see also [10] although this condition is not explic itly mentioned therein. Therefore, in this paper, the FORM optimisation problem is re-written as

minimisex

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi yðx Þ  yðx Þ 



subjec to gðx Þ 6 0

3.2. Algorithm The set of trial solutions (population) is indicated by

S ¼ fx 0 ; x 1 ; . . . x Nsol1 g ð9Þ

Upon obtaining the solution x H and yH ¼ T ðx H Þ, the reliability  index can be computed by Eq. (8). Note that no gradients of g are required in the EA. Finally, the probability of failure can be approximated by means of

pf  UðbÞ

domly generate a population of trial solutions (candidates) that are mixed during a pseudo time loop (evolution or iterations) aiming at producing better candidates. An essential step in EAs is the selection of candidates to be mixed one with another aiming at generating new and better solutions; hence updating the population to the next generation. There are several methods to implement selection including the tournament selection where trial solutions participate in a tournament after being randomly picked, with no bias, from the pool of candidates. On the other hand, other conventional selection methods are based on ranking or some measure of fitness to pick solutions. These other methods usually cause problems such as the loss of diversity and hence convergence. It is noted that diversity maintenance is an essential characteristic to a good performance of an EA; see e.g. [35]. Several methods exist for mixing candidates’ data (recombination), including some that try to mimic genetics by employing crossover and mutation operators—these are known as genetic algorithms (GAs) [22–24]. Another class of recombination methods known as differential evolution (DE) is also available [26–28] and is chosen in this work due to a number of advantages such as directly handling of real numbers without the need for a coding algorithm (e.g. float to binary). In fact, the DE method is a solution technique on its own and has proven to be very efficient [28]. When combined with GA, the DE has also proven to be more efficient [36] than other recombinations methods such as the somewhat popular SBX method [37]. By combining DE with tournaments and the so-called crowding approach to niching [38] aiming at preserving diversity, an accurate, efficient and robust EA named Goga was designed as demonstrated in [29]. With regards to accuracy, several constrained and/or multiobjective optimisation problems were solved therein resulting in exact answers to some problems. With regards to repeatability characteristics, several runs of the same problem were carried out illustrating the good performance of the algorithm. Finally, the efficiency of the parallel code was also demonstrated where, for example, an ideal speed up is obtained when increasing the number of CPUs (groups) from 1 to 16. This code is then used here for reliability analysis and is described next.

ð10Þ

3. Evolutionary algorithm 3.1. Introduction Evolutionary algorithms (EAs) were originally inspired by how life evolves in nature; a stochastic process towards the best. The main application of EAs is hence in optimisation. Algorithms exist for problems with several constraints, several objectives or both. Codes have also been developed to handle multimodal functions and objectives based on binary, integer or real numbers. EAs ran-

ð11Þ

where x i are the trial solutions and there are N sol of them. One objective value (OVA) bi and one vector of out-of-range values (OOR) u i are assigned to each trial solution x i . Note that the use of OORs is a design strategy developed in order to handle constraints in a robust way [29]. The OOR vector u is introduced to effectively implement the constraints—equalities, inequalities or both. The OOR can indicate that constraints are violated or other computations failed; e.g. a call to the FE analysis has failed. For the reliability problem, the OOR will have two components. The first one (u0 ) works as a binary number (a flag) with a positive number indicating either that a failure on the transformation of any random variable occurred or the FE analysis has failed. The second component (u1 ) measures how much the inequality constraint in Eq. (9) has been violated. A ramp function is adopted to model the constraint in Eq. (9) and thus to define the second OOR component as u1 ¼ hgðx Þi where the ramp function (Macaulay brackets) is given by



hsi ¼

0

if s < 0

s

otherwise

ð12Þ

D.M. Pedroso / Structural Safety 65 (2017) 84–99

Note that, in this way, as long as the constraint is satisfied, the ‘‘indicator” u1 will be zero. Otherwise, the more the constraint is violated, the larger the indicator will be. In Goga, this indicator is used to rank trial solutions where the constraint satisfaction is a priority. Therefore, the unfeasible solutions will be incrementally trimmed out of the population. Later, the feasible solutions will start to evolve towards the best. The tournament selection and recombination process using DE are considered here. They are implemented starting with the random selection of pairs of trial solutions in a subset G of S. In this step, the DE is performed, some metrics regarding how close candidates are from each other are computed, and the competition among the closest candidates is executed. The definition of closest competitors is simply specified by an assignment problem. In the end, the subset is updated with new trial solutions. The use of a subset G is made to facilitate the development of a parallel algorithm. Nonetheless, the whole set S could be used instead in serial codes. The steps are presented in Fig. 3 and all minor details are given in [29]. The major steps are explained as follows with focus on the reliability problem. The DE operation is also indicated in Fig. 3. In this figure, a call to the finite element analysis, indicated in red, is shown as well. This call can be readily replaced by the LSF when solving simpler cases; e.g. when the LSF is explicitly defined.

87

The DE recombination operator is quite simple actually and is defined by

xnew ¼ x0;i þ F DE  ðx1;i  x2;i Þ i

ð13Þ

where F DE is a scaling factor that induces a displacement of the trial solution along a random direction and x i are three auxiliary solutions. In this work, F DE is a uniform random variable in ½0; 1Þ. The above expression is only applied if the ‘‘throw of a biased coin” with probability C DE results true; otherwise the components of a fourth vector x (the original one) are simply copied into the new vector x new . Note that, to create a new solution, four randomly selected solutions are required. It is also worth noting that the C DE parameter has some influence on the convergence properties but the value C DE ¼ 0:8 has proven to be a good choice for several experiments. In Fig. 3, the main loop is executed in Dtexc increments corresponding to the time increment for the exchange of solutions among the subsets of solutions G (groups). This increment and the number of groups (N cpu ) can be set to 1 if a serial run is desired. Until the final pseudo time tmax is not reached, the evolution process is continued within each group in parallel. For each group, the evolution is simulated until the current pseudo time reaches the time for exchange texc :

ff

Fig. 3. Evolutionary algorithm to solve FORM.

The update of one group G is as follows. At each time step t, pairs ðA; BÞ of trial solutions are randomly collected generating smaller sets with eight unique solutions in the set that were not yet selected. Each quadruple lead by A or B is then submitted to the DE operation resulting in a new pair ða; bÞ. After that, the objective function is called and the OVA and OOR values are calculated. A detailed description of the objective function call is given in Algorithm 1. Next, the new pair ða; bÞ is stored in a temporary structure (backup set). After all new solutions are created, a number of variables (metrics) are computed. These include for example the neighbouring distances that are used to avoid candidates overlapping each other and hence damaging the diversity of the entire population. The tournament then starts by looping over all pairs ðA; BÞ again. The matching of competitors is decided by the simple formula:

88

D.M. Pedroso / Structural Safety 65 (2017) 84–99

fA; ag and fB; bg if dAa þ dBb < dAb þ dBa fA; bg and fB; ag otherwise

ð14Þ

where dij is the neighbouring distance between trial solutions i and j. The winner of a competition is decided by observing the number of constraint violations, the minimum OVA, the minimum OOR, and the crowding distance; see details in [29]. The current population (group) is then updated with the winner solutions. At the exchange time, two steps are performed. First, tournaments between solutions in different groups that are organised in a cyclic permutation are carried out. Second, one solution is randomly exchanged between randomly selected groups. These two steps allow for some interaction between groups (CPUs) and the process helps with widening the search space.

4. Numerical solutions The performance of the resulting algorithm is investigated by solving a number of reliability problems, including benchmark tests extensively studied in the literature. The meaning of ‘‘good performance” here is connected to accuracy, efficiency and repeatability of results. The last behaviour is particularly of great importance since it indicates that the code can generate similar results after every run. Note that the minimum and maximum values of each random variable must be specified as well because the values generated by an evolutionary algorithm are in principle unbounded. In Goga, values outside the bounds are simply truncated. Here, the min/max values are selected in a way to capture most of the PDF; however very small values are avoided in order to minimise the number of failures in the transformation process. Furthermore, these values are chosen in a way to pose some challenges to the EA code; i.e. the widest feasible range is to be selected. To illustrate, Fig. 4 is drawn where the range 10 6 x 6 10 is selected for a standard normally distributed variable. Thus, the actual nonlinear optimisation problem is

minimisex 

ð15Þ

xmin 6 xj 6 xmax j j where the limits xmin (lower) and xmax (upper) of the random varij j ables are considered as well.

In this subsection, eighteen functions collected from the literature are considered. In each test/problem (P), the limit state function is explicitly given. The probability distributions of the random variables and the reference solutions corresponding to all problems are listed in Table 1. Each problem is described as follows. P1 and P2: The first two tests are defined by

1 gðx Þ ¼ 0:5 ðx0  x1 Þ2  pffiffiffi ðx0 þ x1 Þ þ 3 2

gðx Þ ¼ 1 þ

1 ðx0 þ x1 Þ2  4 ðx0  x1 Þ2 4

ð20Þ

P6 and P7 correspond to case 1 and case 2, respectively, of example 1 in [43]; see also problems 7 and 8 in [40]. Both problems are defined by

ð21Þ

The only difference between P6 and P7 is a small perturbation to the mean value of the second random variable as shown in Table 1. P8 corresponds to problem 7 in [41] and is classified as a highly nonlinear limit state function. It is defined by

ð22Þ

See also problem 9 in [40]. P9 corresponds to example 2 in [43] and problem 10 in [40]. It is defined by

gðx Þ ¼ x30 þ x31  67:5

ð23Þ

P10 corresponds to problem 2 in [41] and lead to multiple failure points; see also problem 11 in [40]. It is defined by

ð16Þ

gðx Þ ¼ x0 x1  146:14

ð17Þ

Note also that this problem poses some challenges due to the large difference of values between x0 and x1 (Table 1). P11 corresponds to example 2 in [42]; see also problem 12 in [40]. It is defined by

and correspond to a convex and concave limit state functions, respectively. They were initially presented in [39] (Example A.1.2); see also Problem 1 in [40]. P3 and P4 correspond to tests 6 and 8 of [41] and are defined by

gðx Þ ¼ 2  x1  0:1 x20 þ 0:06 x30

ð19Þ

respectively; see also Problems 2 and 3 of [40]. They are classified as a non-linear function with a saddle point and a highly non-linear function, respectively. P5 is a modified version of Example 1 in [42]; see also Problem 5 of [40]. The limit state function in this case is a simple quadratic expression defined by

gðx Þ ¼ 2:5  0:2357 ðx0  x1 Þ þ 0:0046 ðx0 þ x1  20Þ4

4.1. Benchmark tests

1 gðx Þ ¼ 0:1 ðx0  x1 Þ2  pffiffiffi ðx0 þ x1 Þ þ 2:5 2

gðx Þ ¼ 3  x1 þ 256 x40

gðx Þ ¼ x30 þ x31  18

bðx Þ

subjec to ui ðx Þ ¼ 0

Fig. 4. Standard normally distributed variable in 10 6 x  10.

ð18Þ

pffiffiffi 0:025 2 gðx Þ ¼ 2:2257  ðx0 þ x1  20Þ3 þ 0:2357 ðx0  x1 Þ 27

ð24Þ

ð25Þ

P12 corresponds to example 7.6 in [6]; see also problem 14 of [40]. The random variables follow now a lognormal distribution. The LSF is defined by

89

D.M. Pedroso / Structural Safety 65 (2017) 84–99 Table 1 Random variables in benchmark tests. test

var

l

r

Da

min

max

1

x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x0 x1 x2 x3 x4 x5 x0 x1 x2 x3 x4 x5 x0 x1 x2 x0

0 0 0 0 0 0 0 0 0 0 10 10 10 9:9 10 10 10 10 78064:4 0:0104 10 10 38 54 0 0 0 0 0 0 0 0 0 0 120 120 120 120 50 40 120 120 120 120 50 40 21:2 20 9:2

1 1 1 1 1 1 1 1 1 1 5 5 5 5 3 3 5 5 11709:7 0:00156 3 3 3:8 2:7 1 1 1 1 1 1 1 1 1 1 12 12 12 12 15 12 12 12 12 12 15 12 0:1 0:2 0:1

N N N N N N N N N N N N N N N N N N N N N N L L N N N N N N N N N N L L L L L L L L L L L L L L L N

10 10 10 10 10 10 10 10 10 10 20 20 20 20 10 10 20 20 15000 0 10 10 20 40 10 10 10 10 10 10 10 10 10 10 60 60 60 60 5 5 60 60 60 60 5 5 20:6 19 8:6

10 10 10 10 10 10 10 10 10 10 40 40 40 40 30 30 40 40 140000 0:02 30 30 60 70 10 10 10 10 10 10 10 10 10 10 180 180 180 180 130 110 180 180 180 180 130 110 21:8 21 9:8

1  106

4  107 0:0002

2 3 4 5 6 7 8 9 10 11 12 13

14

15

16

17

x1 x2 x0

18

a

2  107 0:0001 4

x1

2  107 0:0001

x2

4

5  106 2  105 1

N

5  106

G L

2  105 1

G

L

1  105 1 1  106 1  105 1

10 4:5  107 0:00022 10

N:Normal, L:Lognormal, G:Gumbel.

gðx Þ ¼ x0 x1  1140

ð26Þ

P13 corresponds to problem 3 in [41]; see also problem 6 in [40]. The number of random variables is now 10 and the limit state function is defined by 8 X gðx Þ ¼ 2:0  0:015 x2i  x9

ð27Þ

i¼0

P15 corresponds to example 2 in [32]; see also example 16 in [40]. It is defined in such a way to examine the effect of noise in the LSF. The LSF is defined by

gðx Þ ¼ x0 þ 2:0 x1 þ 2:0 x2 þ x3  5:0 x4  5:0 x5 þ

5 X 0:001 sinð1000 xi Þ

ð29Þ

i¼0

P14 corresponds to example 1 in [32]; see also example 1 in [11]. The LSF is a linear function with 6 variables defined according to

gðx Þ ¼ x0 þ 2:0 x1 þ 2:0 x2 þ x3  5:0 x4  5:0 x5

ð28Þ

P16 corresponds to example 5 in [44]; see also example 17 in [40]. It is defined by

90

D.M. Pedroso / Structural Safety 65 (2017) 84–99

gðx Þ ¼ 240758:1777 þ 10467:364 x0 þ 11410:63 x1 þ 3505:3015 x2  246:81 x0 x0  285:3275 x1 x1  195:46 x2 x2

ð30Þ

P17 and P18 correspond to case 1 and case 2, respectively, of example 4 in [43]; see also Problems 18 and 19 in [40]. Both problems are defined by

gðx Þ ¼ x0 x1  78:12 x2

ð31Þ

The only difference between P17 and P18 is the replacement of the normal distribution to a lognormal one for variables x1 and x2 in P18. Note that the ranges of variables are quite different as well. There are indeed challenges in some of these problems. For instance, as discussed in [41], conventional solution methods are not robust in the case of a complex limit-state, such as a highly non-linear failure function, the existence of multiple failure points or a combination of failure functions. Therefore, other authors have developed stochastic methods such as the radial-based importance sampling [41] to estimate the solution. To illustrate these challenges, the search spaces of Problems 2, 5, 11, 15, 16 and 18 are drawn in Fig. 5 alongside the initial trial candidates, the colour map of the b scalar field and the limit state function. Note that the evolutionary algorithm will at times select the opposite best solution, but will not suffer from problems such as getting stuck at a local minimum or numerical instabilities due to large differences between values such as in the problem 18 (Fig. 5f). As explained above, the minimum and maximum values of each random variable must be defined as well. These are listed in Table 1. To illustrate, the PDF of problems 14 and 18 are shown in Figs. 6 and 7, respectively. Each problem is solved 1000 times with the proposed algorithm. The input data are given in Table 2. The results with regards to the reliability index are statistically analysed and presented in Table 3. The following variables are considered:  N sol is the number of trial solutions; e.g.  10  N x with N x being the number of random variables  N cpu is the number of groups; it is chosen such that N sol is not greater than 40  tmax is the maximum number of iterations (evolution time)  Dtexc is the step for exchange of solutions  C DE is the differential evolution coefficient  N ev al is the total number of function evaluations; e.g. calls to Algorithm 1  T sys is the average computer time to solve one sample in a 4core Intel(R) i7-4770 CPU @ 3.40 GHz (Debian-Sid/GNU/Linux)  bref is the reference reliability index from other papers  bmin ; bav e and bmax are the minimum, average and maximum reliability indices, respectively, found in the 1000 samples  bdev is the standard deviation of b in the set of 1000 samples

Fig. 5. Search space for selected problems illustrating the possibility of multiple solutions and the existence of high non-linearity. The contour indicates the bðxÞ  scalar field. The limit state function is drawn in yellow. The initial candidates are indicated by black circles whereas the best solution is represented by a green star. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

As shown in Table 3, the results clearly indicate that the algorithm produces nearly the same results every time it is run. The standard deviation of b is very small for almost all benchmark tests. Exceptions are problems 11 and 16 with a deviation of the 102 order. Still, these results are quite reasonable. Furthermore, the reference results from the cited papers are achieved within a small range. Finally, with more iterations, the accuracy of results can be improved further. 4.2. Benchmark tests with bounded random variables This subsection presents an interesting benchmark problem involving a simply-supported beam with randomly located concen-

Fig. 6. Problem 14. Lognormal distribution of variable x0 .

91

D.M. Pedroso / Structural Safety 65 (2017) 84–99

dom variables. Specifically, the bounded random variables introduce strong non-linearities in the mapping to standard Gaussian space, making the search for design points more demanding [30]. The problem is illustrated in Fig. 8 where both the location (Di ) and the magnitude of loads (Pi ) are random variables represented by uniform distributions. The case of loads being represented by normal and lognormal distributions are studied as well for the sake of comparisons. Hereafter, this problem is named as the BA benchmark test or the BA problem. Three cases regarding the number of concentrated loads are considered in the BA tests: (1) a single concentrated load that can be randomly located along the whole span of the beam (tests BA01 to BA24); (2) four concentrated loads with positions varying within equally spaced ranges along the beam (tests BA25 and BA26); and (3) eight concentrated loads also with random positions within equally spaced ranges along the beam (tests BA27 and BA28). For the case with one concentrated load, hence with two random variables (x0 ¼ D1 and x1 ¼ P 1 ), the limit state function is

Fig. 7. Problem 18. Gumbel distribution of variable x2 .

trated loads. This problem was recently presented in [30] and poses some challenges to any FORM solver because of the bounded ran-

gðx Þ ¼ R  D1 

! D21 P1 L

ð32Þ

with

Table 2 Benchmark tests: input parameters. P

N sol

N cpu

tmax

Dt exc

C DE

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

20 20 20 20 20 20 20 20 20 20 20 20 80 60 60 30 30 30

1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

1 1 1 1 1 1 1 1 1 1 1 1 10 1 1 1 1 1

0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

R¼

pb L  e 4

ð33Þ

As in [30], the constants in the above equation are treated as non-dimensional and the following values are adopted: pb ¼ 1 and L ¼ 1. Note that there is a mistake in the values of e in [30]— the correct one is e=4. The cases with one concentrated load are additionally studied considering a correlation of the two random variables. In these cases, the off-diagonal terms of the correlation matrix equal q, and the inverse of the lower-triangular matrix from the Cholesky decomposition of the correlation matrix is

"

L1 ¼

1 0 1 pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi

#

ð34Þ

q

1q2

1q2

All combinations of parameters considered in the tests with one concentrated load are collected in Table 5, resulting in 24 tests. For the cases with n concentrated loads, the limit state function is

Table 3 Benchmark tests: results (N samples ¼ 1000). N ev al

T sys

bref

bmin

b av e

bmax

bdev

x0 6

: xref 0

x1

: xref 1

1

2020

4 ms

(2.5000)

2.5000

2.5000

2.5000

1:08  10

1.7678

(1.7677)

1.7678

(1.7677)

2

2020

4 ms

(1.6580)

1.6583

1.6583

1.6643

2:15  104

1.4716

(-0.7583)

0.7645

(1.4752)

3

2020

4 ms

(2.0000)

2.0000

2.0005

2.0728

4:68  103

0.00

(0.00)

2.00

(2.00)

4

2020

4 ms

(3.0000)

3.0000

3.0000

3.0000

7:72  107

0.00

(0.00)

3.00

(3.00)

5

2020

4 ms

(0.3536)

0.3536

0.3536

0.3545

3:09  105

0.25

(-0.25)

0.25

(0.25) (N/A)

6

2020

4 ms

(2.2401)

2.2401

2.2401

2.2401

1:19  106

2.0801

(N/A)

2.0801

7

2020

4 ms

(2.2260)

2.2260

2.2260

2.2260

1:03  106

2.0859

(N/A)

2.0742

(N/A)

8

2020

4 ms

(2.5000)

2.5000

2.5000

2.5003

8:24  106

15.3033

(N/A)

4.6966

(N/A)

9

2020

4 ms

(1.9003)

1.9144

1.9144

1.9144

2:08  106

3.2317

(N/A)

3.2316

(N/A)

10

2020

4 ms

(5.4280)

5.3333

5.3337

5.4021

3:27  103

18378.9

(N/A)

0.0

(N/A)

11

2020

4 ms

(2.2257)

2.2257

2.2292

2.4054

1:94  102

5.2785

(N/A)

14.7215

(N/A)

12

2020

4 ms

(5.2127)

5.2127

5.2127

5.2135

4:01  105

23.7530

(N/A)

47.9940

(N/A)

13

8080

33 ms

(2.0000)

2.0004

2.0024

2.0107

1:12  103

14

6060

45 ms

(2.3480)

2.3482

2.3490

2.3512

4:24  104

15

6060

46 ms

(2.3482)

2.3482

2.3490

2.3509

4:25  104

16

3030

9 ms

(0.8292)

0.8292

0.8330

0.8805

1:01  102

17

3030

9 ms

(3.3221)

3.3221

3.3221

3.3222

6:88  106

18

3030

10 ms

(4.4527)

4.4282

4.4282

4.4285

1:62  105

see Table 4

92

D.M. Pedroso / Structural Safety 65 (2017) 84–99

Table 4 Benchmark tests: resulting x values. P

x0

x1

x2

x3

x4

x5

x6

x7

x8

x9

13 14 15 16 17 18

0.0109 117.422 117.338 21.198 3957896.1 10420508.8

0.0082 115.181 115.128 20.137 0.0 0.0

0.0044 115.300 115.191 9.246 4.5 8.7

0.0015 117.251 117.204

0.0193 83.515 83.715

0.0134 55.613 55.321

0.0195

0.0120

0.0104

2.0001

Table 5 Constants and random variables in BA Problems 1 to 24.

Fig. 8. BA problem: simply-supported beam subject to multiple randomly-located concentrated loads of random magnitude.

gðx Þ ¼ R 

! D2 Di  i Pi L

n=2 X i¼1

ð35Þ

where the value of R is either fixed to a given value or set to be equal to the critical resistance Rcrit , computed by means of

Rcrit ¼

2 pb L ðs1 þ s2 Þ n

ð36Þ

with

s1 ¼

n=4 X i¼1

2

i

2i n

! s2 ¼

n=2 X i¼n=4þ1

" ði  1Þ 

2 ði  1Þ n

2

# ð37Þ

Four additional tests are thus defined with four or eight concentrated loads. The combinations and parameters considered are collected in Table 6. Tests BA25 and BA27 employ a prescribed R value whereas tests BA26 and BA28 employ Eq. (36) as the R value, but with a small perturbation. The reason for the perturbation is explained next. An additional challenge arises when the limit state curve intersects the boundary of the space of the bounded random variables. In [30], this situation defines the ‘‘corner” problems, and it is noted that these problems cannot be solved by conventional methods because the convergence criterion cannot be reached. On the other hand, stochastic methods such as the evolutionary algorithm can handle the corner problems as long as there is a ‘‘small gap” between the curve and the limit of variables; see, e.g. 9–f. The gap between the limit state curve and the boundary of the random variables is controlled by the perturbation e. In the present work, the adopted perturbations for tests BA26 and BA28 are 104 and 102 , respectively (Table 6). These values were selected in a way to guarantee a minimum gap for the evolutionary algorithm. With smaller values, more trial solutions would have been required. For example, with e ¼ 108 , a huge number of function

a

Test

e

q

Var

l

r

Da

Min

Max

BA01

0:0125

0:0

BA02

0:0125

0:0

BA03

0:0125

0:0

BA04

0:0125

0:5

BA05

0:0125

0:5

BA06

0:0125

0:5

BA07

0:0125

0:5

BA08

0:0125

0:5

BA09

0:0125

0:5

BA10

108

0:0

D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1

– – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – –

U U U U U U U U U U U U U U U U U U U

0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0

1:0 1:0 0:6 1:0 0:5 1:0 1:0 1:0 0:6 1:0 0:5 1:0 1:0 1:0 0:6 1:0 0:5 1:0 0:5

BA11

108

0:5

P1 D1

– –

– –

U U

0:6 0:0

1:0 0:5

BA12

108

0:5

P1 D1

– –

– –

U U

0:6 0:0

1:0 0:5

BA13

0:0125

0:0

BA14

0:0125

0:0

BA15

0:0125

0:0

BA16

0:0125

0:0

BA17

0:0125

0:0

BA18

0:0125

0:0

BA19

0:0125

0:0

BA20

0:0125

0:0

BA21

0:0125

0:0

BA22

0:0125

0:0

BA23

0:0125

0:0

BA24

0:0125

0:0

P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1 D1 P1

– – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8 – 0:8

– – 0:04 – 0:04 – 0:04 – 0:16 – 0:16 – 0:16 – 0:04 – 0:04 – 0:04 – 0:16 – 0:16 – 0:16

U U N U N U N U N U N U N U L U L U L U L U L U L

0:6 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:0 0:0 0:0 0:0 0:0 0:0 0:6 0:0 0:6 0:0 0:6 0:0 0:0 0:0 0:0 0:0 0:0

1:0 1:0 1:0 0:6 1:0 0:5 1:0 1:0 1:6 0:6 1:6 0:5 1:6 1:0 1:0 0:6 1:0 0:5 1:0 1:0 1:6 0:6 1:6 0:5 1:6

N:Normal, L:Lognormal, U:Uniform

evaluations are required. This issue is equivalent to the one mentioned in [30] where no solver can handle the corner problems. All BA tests are now solved by the evolutionary algorithm times for investigating its reproducibility capabilities. The input parameters for the EA are listed in Table 7 and the results are listed

D.M. Pedroso / Structural Safety 65 (2017) 84–99

93

Table 6 Definitions and random variables in BA Problems 25–28.

a

Test

R

Var

Da

Min

Max

BA25

0:8

BA26

0:8750  104 ðRcrit   eÞ

D1 D2 D3 D4 P 1 to P 4 D1

U U U U U U

0:00 0:25 0:50 0:75 0:60 0:00

0:25 0:50 0:75 1:00 1:00 0:25

D2 D3 D4 P 1 to P 4 D1 D2 D3 D4 D5 D6 D7 D8 P 1 to P 8 D1

U U U U U U U U U U U U U U

0:25 0:50 0:75 0:60 0:000 0:125 0:250 0:375 0:500 0:625 0:750 0:875 0:600 0:000

0:50 0:75 1:00 1:00 0:125 0:250 0:375 0:500 0:625 0:750 0:875 1:000 1:000 0:125

D2 D3 D4 D5 D6 D7 D8 P 1 to P 8

U U U U U U U U

0:125 0:250 0:375 0:500 0:625 0:750 0:875 0:600

0:250 0:375 0:500 0:625 0:750 0:875 1:000 1:000

BA27

1:4

BA28

1:5625  102 ðRcrit   eÞ

U:Uniform

in Table 8. The reference reliability indices and random variables are also listed in Table 8. In addition, the solution (random variables) from problems BA25 to BA28 are listed in Table 9. It can be concluded that the EA works well for this challenging set of reliability analysis problems. 4.3. Reliability of frames using the finite element method The performance of the EA code is now investigated using 2D and 3D frame models solved by the finite element method [16]. The elements are standard Euler–Bernoulli linear elastic beams with the local system as illustrated in Fig. 10. The geometry and boundary conditions for the 2D problem are illustrated in Fig. 11 where the 12 storey frame is defined and five groups of beams are indicated. Each group corresponds to a different combination of cross-section properties and random parameters. The equivalent 3D problem is illustrated in Fig. 12. For the 2D problem, the nodal load P at the left-hand side is also a random variable. Thus, the set of random variables to be found comprises of six quantities as follows

x ¼ fAG1 ; AG2 ; AG3 ; AG4 ; AG5 ; Pg

ð38Þ

where Gi means group i. For the 3D problem, two cases are studied. Case A is equivalent to the 2D problem and has Py only with a Gumbel distribution. Case B has P x and P y , both with uniform distribution. The unknown vector for the 3D problem is, therefore,

  x ¼ AG1 ; AG2 ; AG3 ; AG4 ; AG5 ; Py ¼ P

ð39Þ

for Case B.

The random values are listed in Table 10. After finding the cross-sectional area A corresponding to the design point, the sectional moments of inertia in each group are computed by the following formula

I22 ¼ I11 ¼ a A2

ð40Þ

ð41Þ

where the a coefficients for groups from 1 to 5 are, respectively equal to 0:08333; 0:08333; 0:08333; 0:26670; 0:20. The Young modulus for all groups is deterministic and equal to E ¼ 2  107 [units of pressure]. For the 3D problem, the torsion constant is also required and is computed according to

J 00 ¼ c A2

for Case A, or

  x ¼ AG1 ; AG2 ; AG3 ; AG4 ; AG5 ; Px ; Py

Fig. 9. EA solutions to BA tests illustrating also some ‘‘corner” problems (BA10, BA26, BA28). The contour indicates the bðxÞ scalar field. The limit state function is  drawn in yellow. The initial candidates are indicated by black circles whereas the best solution is represented by a green star. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

ð42Þ

where the c coefficients for groups from 1 to 5, are, respectively equal to 0:1; 0:1; 0:1; 0:3, and 0:3. Furthermore, the shear modulus is set equal to G ¼ 7:5  106 [units of pressure] for all groups.

94

D.M. Pedroso / Structural Safety 65 (2017) 84–99

Table 7 Input parameters for BA problems. P

N sol

N cpu

tmax

Dt exc

C DE

BA01 BA02 BA03 BA04 BA05 BA06 BA07 BA08 BA09 BA10 BA11 BA12 BA13 BA14 BA15 BA16 BA17 BA18 BA19 BA20 BA21 BA22 BA23 BA24 BA25 BA26 BA27 BA28

20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 80 80 160 160

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 4 4

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 20 20 20 20

0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

Table 8 Results from BA problems (N samples ¼ 1000). P

N ev al

T sys

bref

bmin

bav e

bmax

bdev

x0 ¼ D

: xref 0

x1 ¼ P

: xref 1

7

BA01

2020

4 ms

(1.1500)

1.1503

1.1503

1.1504

4:55  10

0.500

(0.500)

0.950

(0.950)

BA02

2020

4 ms

(1.4100)

1.4257

1.4257

1.4257

6:57  107

0.465

(0.460)

0.955

(0.950)

BA03

2020

4 ms

(1.7800)

1.7792

1.7792

1.7792

1:17  106

0.436

(0.430)

0.966

(0.970)

BA04

2020

4 ms

(1.3000)

1.3189

1.3189

1.3190

6:89  107

0.513

(0.510)

0.951

(0.950)

BA05

2020

4 ms

(1.1900)

1.2061

1.2061

1.2061

4:05  107

0.479

(0.480)

0.952

(0.950)

BA06

2020

4 ms

(1.4600)

1.4605

1.4605

1.4605

2:60  107

0.441

(0.440)

0.963

(0.960)

BA07

2020

4 ms

(1.2900)

1.3189

1.3189

1.3190

3:74  107

0.487

(0.480)

0.951

(0.950)

BA08

2020

4 ms

(1.9400)

1.9764

1.9764

1.9765

6:56  107

0.457

(0.450)

0.957

(0.950)

BA09

2020

4 ms

(2.4600)

2.5077

2.5077

2.5081

1:51  105

0.433

(0.420)

0.967

(0.970)

BA10

2020

4 ms

(N/A)

6.4326

6.4328

6.4496

6:82  104

0.500

(0.500)

1.000

(1.000)

BA11

2020

4 ms

(N/A)

5.3996

5.3999

5.4716

2:38  103

0.500

(0.500)

1.000

(1.000)

BA12

2020

4 ms

(N/A)

9.0021

9.0025

9.0136

7:53  104

0.500

(0.500)

1.000

(1.000)

BA13

2020

4 ms

(3.7500)

3.7500

3.7506

4.3711

1:96  102

0.500

(0.500)

0.950

(0.950) (0.950)

BA14

2020

4 ms

(3.8700)

3.8660

3.8660

3.8660

1:27  106

0.492

(0.490)

0.950

BA15

2020

4 ms

(4.1400)

4.1377

4.1377

4.1392

4:74  105

0.468

(0.470)

0.954

(0.950)

BA16

2020

4 ms

(0.9400)

0.9375

0.9375

0.9375

5:56  107

0.500

(0.500)

0.950

(0.950)

BA17

2020

4 ms

(1.2000)

1.1987

1.1987

1.1988

4:74  106

0.438

(0.440)

0.965

(0.960)

BA18

2020

4 ms

(1.4500)

1.4534

1.4534

1.4535

5:61  106

0.401

(0.400)

0.989

(0.990)

BA19

2020

4 ms

(3.4600)

3.4641

3.4641

3.4642

1:08  106

0.500

(0.500)

0.950

(0.950)

BA20

2020

4 ms

(3.5900)

3.5876

3.5876

3.5877

2:60  106

0.490

(0.490)

0.950

(0.950)

BA21

2020

4 ms

(3.8700)

3.8557

3.8563

4.4949

2:02  102

0.464

(0.470)

0.955

(0.960)

BA22

2020

4 ms

(0.9700)

0.9668

0.9668

0.9668

1:27  108

0.500

(0.500)

0.950

(0.950)

BA23

2020

4 ms

(1.2100)

1.2105

1.2105

1.2106

1:33  106

0.433

(0.430)

0.967

(0.970)

BA24

2020

4 ms

(1.4400)

1.4410

1.4410

1.4411

6:00  107

0.393

(0.390)

0.996

(1.000)

BA25

16080

52 ms

(3.0800)

3.0883

3.0885

3.0890

1:03  104

BA26

16080

54 ms

(N/A)

9.4524

9.4596

9.4870

4:14  103

BA27

64160

236 ms

(3.6200)

3.6272

3.6277

3.6286

2:12  104

BA28

64160

240 ms

(N/A)

8.8255

8.8267

8.8286

4:39  104

The limit state for this example is defined as a serviceability control by monitoring the horizontal displacement (uxA or uyA ) of the highlighted node A in Figs. 11 or 12 such that it does not exceed

see Table 9

0.096 [units of length]. Therefore, the limit state function for the 2D and 3D-Case A (3D-A) problems is

gðx Þ ¼ 0:096  u A ðx Þ

ð43Þ

95

D.M. Pedroso / Structural Safety 65 (2017) 84–99 Table 9 Solutions and reference x values of some BA problems. P

x0

x1

x2

x3

x4

x5

x6

x7

x8

x9

x10

x11

x12

x13

x14

x15

BA25 Ref. BA26 Ref. BA27 Ref. BA28 Ref.

0.227 0.227 0.250 0.250 0.108 0.108 0.124 0.125

0.435 0.435 0.497 0.500 0.228 0.228 0.249 0.250

0.565 0.565 0.502 0.500 0.344 0.343 0.372 0.375

0.772 0.773 0.750 0.750 0.449 0.450 0.484 0.500

0.941 0.940 1.000 1.000 0.549 0.545 0.517 0.500

0.955 0.954 1.000 1.000 0.656 0.657 0.628 0.625

0.955 0.954 1.000 1.000 0.772 0.772 0.751 0.750

0.941 0.940 1.000 1.000 0.892 0.892 0.876 0.875

0.902 0.900 0.994 1.000

0.935 0.930 0.996 1.000

0.947 0.954 0.997 1.000

0.950 0.950 0.997 1.000

0.950 0.950 0.997 1.000

0.946 0.954 0.997 1.000

0.936 0.930 0.996 1.000

0.901 0.900 0.994 1.000

Fig. 10. Euler–Bernoulli beam element. Local system.

Fig. 12. 3D frame for reliability analysis. There are five groups as in the 2D problem. All lengths correspond to the lengths in the 2D frame. The lengths along x and y are equal to each other. All nodes indicated by an arrow are loaded. Node A is the control node for Case A. Py (Gumbel) is applied in Case A only whereas Case B uses Px (Uniform) and Py (Uniform).

Table 10 Random variables for 2D frame problem.

Fig. 11. 2D frame for reliability analysis. The horizontal displacement of the circled node A @ (18.5, 48.0) is monitored. There are five groups of beams; each group is identified by a different colour and tag from 1 to 5.

a

Problem

Group

Param

l

r

Da

Min

Max

2D, 3D 2D, 3D 2D, 3D 2D, 3D 2D, 3D 2D, 3D-A 3D-B 3D-B

1 2 3 4 5

A A A A A P Px Py

0.25 0.16 0.36 0.20 0.15 30.0 – –

0.025 0.016 0.036 0.020 0.015 7.5 – –

L L L L L G U U

0.15 0.10 0.24 0.12 0.10 10.0 5.0 5.0

0.35 0.24 0.52 0.28 0.21 80.0 20.0 20.0

L:Lognormal, G:Gumbel, U:Uniform.

96

D.M. Pedroso / Structural Safety 65 (2017) 84–99

Table 11 EA parameters for the 2D/3D frame problems. P

N sol

N cpu

tmax

Dtexc

C DE

2D 3D-A 3D-B

60 60 70

2 2 2

100 100 100

10 10 10

0.8 0.8 0.8

Table 12 Solutions of the 2D/3D frame problems. P

N ev al

T sys

bref

bgoga

x0

x1

x2

x3

x4

x5

x6

2D 3D-A 3D-B

6060 6060 6870

2.59s 2m29s 2m43s

1.463 – –

1.4543 1.1510 0.5938

0.244 0.243 0.245

0.158 0.157 0.159

0.353 0.356 0.368

0.192 0.193 0.196

0.147 0.148 0.149

40.774 37.662 15.296

12.539

Fig. 13. Results from the 2D frame problem. Combination of random variables. The contour indicates the bðxÞ scalar field. The limit state function is drawn in yellow. The  initial candidates are indicated by black circles whereas the best solution is represented by a green star. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

D.M. Pedroso / Structural Safety 65 (2017) 84–99

97

Fig. 14. Results from the 3D-B frame problem. Combination of random variables. The contour indicates the bðxÞ scalar field. The limit state function is drawn in yellow. The  initial candidates are indicated by black circles whereas the best solution is represented by a green star. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For the 3D problem-Case B (3D-B), the LSF is given by

gðx Þ ¼ min

0:096  uAmax 0:096  uAmax 280  MAmax y x 22 ; ; ; 0:096 0:096 280 ! 260  MAmax 12  T Amax 11 00 ; 260 12

ð44Þ

where v Amax indicates the maximum absolute value v among all nodes or integration points (e.g. for momentum) in the mesh. Therein, M 22 and M 11 are bending momenta and T 00 is the torsional reaction. Note again that the evolutionary algorithm is ‘‘derivative-free” and hence the derivative of gðx Þ with respect to x is not required. This property is especially a great advantage to the more

complicate LSF used in the 3D problem-Case B. Note that some variables are computed at the nodes while others are available at the integration points. Moreover, the search for the maximum absolute value in the whole mesh, together with the use of the minimum function, makes a quite complex LSF. The three problems (2D, 3D-A, 3D-B) are solved by the EA with the results being collected in Table 12. For these problems, the Goga parameters are listed in Table 11. The reference reliability index for the 2D problem is taken from [16] and corresponds to the GA-ANN algorithm presented therein. We observe that the reliability index of the 3D-A problem is smaller than the one for the similar 2D problem. As for the 3D-B problem, the index is much smaller due to the introduction of other strict constraints. To further demonstrate the solutions obtained with Goga, the objective value (b) and the LSF, considering all combinations of

98

D.M. Pedroso / Structural Safety 65 (2017) 84–99

Fig. 15. Deformed frame from the 3D-A problem using the optimal random variables leading to the minimum reliability index. The displacements are magnified by 50 times.

random variables, are illustrated in Fig. 13, for the 2D problem, and in Fig. 14, for the 3D-B. Note that the LSF for the 3D-A problem is similar to the 2D one and hence it is not presented here. We observe that the LSF for the 3D-B problem is discrete and may cause problems to direct solvers. In addition, the last two random variables in x are bounded due to the uniform distribution. Therefore, this problem is somewhat challenging, but can satisfactorily be solved by Goga. Note that most trial solutions converge to the minimum b; however still satisfying the constraint. To further illustrate, a single finite element analysis with the critical x values listed in Table 12 is carried out. The resulting deformed frame from the 3D-A problem is illustrated in Fig. 15. The deformed frame from the 3D-B problem is illustrated in Fig. 16. In this problem (3D-B), the critical constraint is the torsional reaction T 00 . The values used to compute each component of the LSF defined in Eq. (44) are listed in Table 13. In this table, we observe that the displacements are well below the limits, for instance.

5. Conclusions The solution of the optimisation problem defined by the FORM approximation to reliability problems is obtained using an

Fig. 16. Deformed frame from the 3D-B problem using the optimal random variables leading to the minimum reliability index. The displacements are magnified by 100 times.

Table 13 Critical (optimal) results for the 3D problem-Case B. Variable

Allowed Max Absa

Current Max Absb

LSF Componentc

ux uy M 22 M 11 T 00

0.096 0.096 280.0 260.0 12.0

0.0750838 0.0662067 235.70998 234.27843 12.000357

0.21787705 0.31034649 0.15817861 0.09892908 2:979  105

a

Numeric constants in Eq. (44). Computed maximum absolute values among all nodes or integration points for use in Eq. (44). c Each one of the five input to the min function in Eq. (44). b

evolutionary algorithm (EA) in this paper. The code is named Goga and is able to solve large scale problems due to its parallel computing capabilities. The relative new method to combine trial solutions (population) named differential evolution (DE) is considered. This method is able to search for solutions in a more efficient way than previous genetic algorithm (GA) operators. As an essential step to maintain diversity of the population, tournaments between

D.M. Pedroso / Structural Safety 65 (2017) 84–99

randomly selected trial solutions are carried out. By maintaining diversity, the convergence of the algorithm is also improved. The code is designed with accuracy, efficiency and repeatability properties in mind. After running several benchmark tests and 2D/3D frame problems modelled by the FEM, it is observed that the algorithm successfully generates the same results with multiple runs; thus has very good repeatability characteristics. Not only this behaviour is convenient when solving generic FORM problems, but is also important to make the code reliable, noting that EAs are stochastic processes. The specific steps to implement the constraints and failure of finite element analysis or transformation of random variables are explained in detail and are the key to apply an EA to FORM. The reference solutions from different papers, including recent ones, are captured fairly well, and the case of bounded random variables is conveniently handled by the EA. Finally, the efficiency of the algorithm can certainly be improved if information of each specific problem is taken into account as done in other solution techniques. Note that there are many other methods around (e.g. HLRF, iHLRF and nHLRF, more), and the user has to try each one to see ‘‘which one works”. The proposed EA aims for avoiding this situation (e.g. to be a general purpose solver). Therefore, it was tested with the same properties for all experiments. By tweaking these properties, another gain of efficiency would be possible. Nonetheless, the robustness and generality of the proposed code were the main goal of this study; thus the properties were kept fixed. Moreover, in large scale simulations, the parallel code can produce results in the range of minutes, without any fiddling of parameters by the user (a great convenience). References [1] Galambos T, Ellingwood B, MacGregor J, Cornell C. Probability based load criteria: assessment of current practice. J Struct Div ASCE 1982;108 (ST5):959–77. [2] Ellingwood B, MacGregor J, Galambos T, Cornell C. Probability based load criteria: load factors and load combinations. J Struct Div ASCE 1982;108 (ST5):978–97. [3] Hasofer AM, Lind NC. An exact and invariant first order reliability format. J Eng Mech Div ASCE 1974;100:111–21. [4] Rackwitz R, Flessler B. Structural reliability under combined random load sequences. Comput Struct 1978;9(5):489–94. http://dx.doi.org/10.1016/00457949(78)90046-9. [5] Melchers RE. Structural reliability analysis and prediction. 2nd ed. Wiley; 1999, ISBN 978-0-471-98771-0. [6] Haldar A, Mahadevan S. Probability, reliability and statistical methods in engineering design. Wiley; 2000, ISBN 978-0-471-33119-3. [7] Der Kiureghian A. First- and second-order reliability methods. In: Nikolaidis E, Ghiocel D, Singhal S, editors. Engineering design reliability handbook. CRC Press; 2004, ISBN 978-0-8493-1180-2. p. 1–24. http://dx.doi.org/10.1201/ 9780203483930.ch14. Chap.10. [8] Liu PL, Der Kiureghian A. Optimization algorithms for structural reliability. Struct Saf 1991;9(3):161–77. http://dx.doi.org/10.1016/0167-4730(91)900417. [9] Liu P, Der Kiureghian A. Finite element reliability of geometrically nonlinear uncertain structures. J Eng Mech 1991;117(8):1806–25. http://dx.doi.org/ 10.1061/(ASCE)0733-9399(1991) 117:8(1806). [10] Rackwitz R. Reliability analysis—a review and some perspectives. Struct Saf 2001;23(4):365–95. http://dx.doi.org/10.1016/S0167-4730(02)00009-7. [11] Cheng J, Xiao RC. Serviceability reliability analysis of cable-stayed bridges. Struct Eng Mech 2005;20(6):609–30. http://dx.doi.org/10.12989/ sem.2005.20.6.609. [12] Der Kiureghian A, Haukaas T, Fujimura K. Structural reliability software at the University of California, Berkeley. Struct Saf 2006;28(1–2):44–67. http://dx. doi.org/10.1016/j.strusafe.2005.03.002. [13] Periçaro G, Santos S, Ribeiro A, Matioli L. HLRF–BFGS optimization algorithm for structural reliability. Appl Math Model 2015;39(7):2025–35. http://dx.doi. org/10.1016/j.apm.2014.10.024. [14] Cui L, Sheng D. Genetic algorithms in probabilistic finite element analysis of geotechnical problems. Comput Geotech 2005;32(8):555–63. http://dx.doi. org/10.1016/j.compgeo.2005.11.005. [15] Deng L, Ghosn M, Shao S. Development of a shredding genetic algorithm for structural reliability. Struct Saf 2005;27(2):113–31. http://dx.doi.org/10.1016/ j.strusafe.2004.06.002.

99

[16] Cheng J. Hybrid genetic algorithms for structural reliability analysis. Comput Struct 2007;85(19–20):1524–33. http://dx.doi.org/10.1016/j.compstruc.2007. 01.018. [17] Cheng J, Li QS. Reliability analysis of structures using artificial neural network based genetic algorithms. Comput Methods Appl Mech Eng 2008;197(45– 48):3742–50. http://dx.doi.org/10.1016/j.cma.2008.02.026. [18] Deb K, Gupta S, Daum D, Branke J, Mall A, Padmanabhan D. Reliability-based optimization using evolutionary algorithms. IEEE Trans Evol Comput 2009;13 (5):1054–74. http://dx.doi.org/10.1109/tevc.2009.2014361. [19] Cheng J. An artificial neural network based genetic algorithm for estimating the reliability of long span suspension bridges. Finite Elem Anal Des 2010;46 (8):658–67. http://dx.doi.org/10.1016/j.finel.2010.03.005. [20] Gomes HM, Awruch AM, Lopes PAM. Reliability based optimization of laminated composite structures using genetic algorithms and artificial neural networks. Struct Saf 2011;33(3):186–95. http://dx.doi.org/10.1016/j. strusafe.2011.03.001. [21] Luo X, Li X, Zhou J, Cheng T. A Kriging-based hybrid optimization algorithm for slope reliability analysis. Struct Saf 2012;34(1):401–6. http://dx.doi.org/ 10.1016/j.strusafe.2011.09.004. [22] Goldberg DE. Genetic algorithms in search, optimization and machine learning. 1st ed. Addison-Wesley; 1989, ISBN 0201157675. [23] Michalewicz Z. Genetic algorithms + data structures = evolution programs. Springer; 1996, ISBN 978-3-662-03315-9. http://dx.doi.org/ 10.1007/978-3-662-03315-9. [24] Deb K. Multi-objective optimization using evolutionary algorithms. Wiley; 2001, ISBN 978-0-471-87339-6. [25] Pedroso DM, Williams DJ. Automatic calibration of soil–water characteristic curves using genetic algorithms. Comput Geotech 2011;38(3):330–40. http:// dx.doi.org/10.1016/j.compgeo.2010.12.004. [26] Storn R, Price KV. Differential evolution–a simple and efficient adaptive scheme for global optimization over continuous spaces. Tech. Rep. TR-95-012; International Computer Science Institute; 1995. [27] Storn R, Price K. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;11 (4):341–59. http://dx.doi.org/10.1023/A:1008202821328. [28] Price K, Storn RM, Lampinen JA. Differential evolution: a practical approach to global optimization. Springer-Verlag 2005. http://dx.doi.org/10.1007/3-54031306-0. [29] Pedroso D.M. Parallel evolutionary algorithm in Go language for single and multi objective optimisation problems with constraints. Applied Soft Computing 2016; Submitted for publication. [30] Beck AT, Ávilada Jr SCR. Strategies for finding the design point under bounded random variables. Struct Saf 2016;58:79–93. http://dx.doi.org/10.1016/j. strusafe.2015.08.006. [31] Der Kiureghian A, Liu PL. Structural reliability under incomplete probability information. J Eng Mech 1986;112(1):85–104. http://dx.doi.org/10.1061/ (ASCE)0733-9399(1986) 112:1(85). [32] Der Kiureghian A, Lin H, Hwang S. Second-order reliability approximations. J Eng Mech 1987;113(8):1208–25. http://dx.doi.org/10.1061/(ASCE)0733-9399 (1987) 113:8(1208). [33] Sudret B, Der Kiureghian A. Stochastic finite element methods and reliability. Tech. Rep. UCB/SEMM-2000/08; Department of Civil and Environmental Engineering, University of California, Berkeley; 2000. [34] Haukaas T, Der Kiureghian A. Strategies for finding the design point in nonlinear finite element reliability analysis. Probab Eng Mech 2006;21(2):133–47. http://dx.doi.org/10.1016/j.probengmech.2005.07.005. [35] Mahfoud SW. Niching methods for genetic algorithms. Illinois Genetic Algorithms Laboratory; University of Illinois at Urbana-Champaign; 1995 [Ph.D. thesis]. [36] Tiwari S, Fadel G, Deb K. AMGA2: improving the performance of the archivebased micro-genetic algorithm for multi-objective optimization. Eng Optim 2011;43(4):377–401. http://dx.doi.org/10.1080/0305215X.2010.491549. [37] Deb K, Agrawal RB. Simulated binary crossover for continuous search space. Complex Syst 1995;9:115–48. [38] Mengshoel OJ, Goldberg DE. The crowding approach to niching in genetic algorithms. Evol Comput 2008;16(3):315–54. http://dx.doi.org/10.1162/ evco.2008.16.3.315. [39] Borri A, Speranzini E. Structural reliability analysis using a standard deterministic finite element code. Struct Saf 1997;19(4):361–82. http://dx. doi.org/10.1016/S0167-4730(97)00017-9. [40] Santos SR, Matioli LC, Beck AT. New optimization algorithms for structural reliability analysis. Comput Model Eng Sci 2012;83(1):23–56. http://dx.doi. org/10.3970/cmes.2012.083.023. [41] Grooteman F. Adaptive radial-based importance sampling method for structural reliability. Struct Saf 2008;30(6):533–42. http://dx.doi.org/ 10.1016/j.strusafe.2007.10.002. [42] Grandhi RV, Wang L. Higher-order failure probability calculation using nonlinear approximations. Comput Methods Appl Mech Eng 1999;168(1– 4):185–206. http://dx.doi.org/10.1016/S0045-7825(98)00140-6. [43] Santosh T, Saraf R, Ghosh A, Kushwaha H. Optimum step length selection rule in modified HL–RF method for structural reliability. Int J Pressure Vessels Piping 2006;83(10):742–8. http://dx.doi.org/10.1016/j.ijpvp.2006.07.004. [44] Mahadevan S, Shi P. Multiple linearization method for nonlinear reliability analysis. J Eng Mech 2001;127(11):1165–73. http://dx.doi.org/10.1061/(ASCE) 0733-9399(2001) 127:11(1165).