An approach for the global reliability based optimization of the size and shape of truss structures

– This paper presents a methodology to perform global reliability based design optimization (RBDO) of the size and shape of truss structures. This methodology is comprised by the use of a global constraint and the response surface method to deal with the reliability analysis together with the ﬁreﬂy algorithm (FA) to carry out the structural optimization. The former is responsible for the reduction of the computational cost required in the evaluation of the probabilistic constraints. The latter overcomes the issues related to the non-convexity and mixed-variables of the optimization problem. Two examples are analysed in order to show the eﬀectiveness of the methodology. In these examples, the FA is compared to other two well-known algorithms, the harmony search (HS) and genetic algorithm (GA). As a result, the FA reached the best performance in the examples analysed. All the optima found were checked using a classical ﬁrst order reliability method (FORM) approach, validating the results provided by the response surface method.


Introduction
The deterministic optimization of truss structures has been widely studied in the literature [1][2][3].The optimization of this kind of structure usually leads to mixed variables, i.e. discrete and continuous design variables, and non-convex problems in the presence of many local minima [4][5][6].In order to overcome these two difficulties, several metaheuristic algorithms have been employed, being the most widely applied the genetic algorithm (GA) [7][8][9][10][11][12].
However, the use of GA has presented some drawbacks, mainly related to the long computational time required when dealing with large computational models.Hence, alternative approaches have been developed in order to reduce the computational demand needed.Some methods that have been employed for truss optimization are: particle swarm optimization (PSO) [13,14] simulated annealing (SA) [15], ant colony optimization [16], harmony search [17,18], charged system search (CSS) [19], artificial bee colony [20], big bang-big a Corresponding author: rafael.holdorf@ufsc.brcrunch [21], taboo search [22] and cuckoo search algorithm (CS) [23].For a comprehensive review of the deterministic optimization of truss structures using metaheuristics, the reader is referred to Saka [24], Lamberti and Pappalettere [25] and the references therein.However, it is widely acknowledged that deterministic optimization is not robust with respect to the uncertainties which affect engineering design.In deterministic optimization, potential failure modes are converted in deterministic constraints and uncertainty is addressed indirectly by means of safety factors and conservative assumptions.This approach is inherited from design through design codes, but it is essentially non-optimal for two reasons.First, safety factors approaches tend to be conservative, since it must take into account a large number of cases with different characteristics.Besides, designs obtained using the safety factors approach are potentially less safe than required if the problem under study does not match exactly the assumptions made in order to obtain the safety factors used.
The main techniques used to take into account uncertainties in optimization using probabilistic approaches Article published by EDP Sciences are: robust optimization (RO) [26] and reliability based design optimization (RBDO) [27].The RO has as main goal the minimization of the variability of some parameters related to system response.For example, Calafiore and Dabbene [28] applied this concept to the field of truss structures.On the other hand, the main goal of the RBDO is to optimize the structure ensuring that its probability of failure is lower than a certain level, chosen a priori by the designer [29].RBDO of the topology or size of truss structures has been studied by Nakib [30] and Thampan and Krishnamoorty [31].The multi-objective optimization of truss structures considering uncertainties was performed by Greiner and Haleja [32].Only a few papers have dealt with the reliability based shape and topology optimization of truss structures [5,31,33,34].It is worth to highlight that none of the above mentioned studies employed a global optimization algorithm to carry out the RBDO of truss structures.However, global optimization algorithms are required due to the non-convexity of the problem at hand.It has been observed that the computational cost of the reliability analysis is one of the main issues in the application of RBDO for real problems.The evaluation of the probability of failure is traditionally made using first and second order reliability methods (FORM and SORM, respectively) [35].However, a reliability analysis problem must be solved for each limit state function of the problem, being most problems composed of several limit state functions (e.g. the stress constraint inside each bar of the structure defines a limit state function).Besides, FORM and SORM require the evaluation of the gradient of the limit state function, what is done using finite differences in most cases, e.g. the structural behaviour is analysed using a finite element model.Consequently, the evaluation of the probability of failure for a given design vector requires, in general, a large number of calls of some finite element code.This number is drastically increased in the context of RBDO since several design vectors must be checked in order to carry out optimization.One approach to overcome these difficulties is the utilization of an analytical approximated model as a surrogate for the reliability analysis, which can be obtained with the help of the response surface methodology [35,[37][38][39].By choosing an appropriate approximation, it is possible to represent the response of the system accurately.This reduces the computational cost of the optimization process because the reliability analysis problem is solved for the approximated analytical solution and consequently all the information needed (i.e.gradients) can be evaluated efficiently.Besides, it can be used as a black-box in a non-intrusive manner.
Computational cost of metaheuristic algorithm for RBDO is a critical issue because both optimization and reliability analysis tasks are usually expensive.However, the firefly algorithm (FA) [40] has proved to be more accurate and efficient than popular metaheuristic algorithms like GA and the PSO.For this reason, several researchers have focused their attention on solving optimization problems using FA in a growing number of papers [41][42][43][44][45], including the deterministic optimization of truss structures developed by Miguel et al. [46].However, its implementation in the field of structural optimization is still fairly recent and requires a substantial amount of further study [47].Moreover, at the best of the authors' knowledge, the FA has never been applied to RBDO problems.In view of this, we propose a framework for sizing and shape RBDO of truss structures.The firefly algorithm is utilized to perform design optimization while response surfaces serve to reduce the computational cost of reliability analysis.Since FA is utilized as optimizer, the proposed methodology attempts to perform global RBDO of truss structures.The validity of the present approach is proven by comparing optimization results with those obtained by RBDO frameworks employing harmony search and genetic algorithms.
The paper is organized as follows.Section 2 presents the formulation of the RBDO problem.The methodology used to pursue the reliability analysis using the response surface methodology is described in detail in Section 3. A description of the FA is given in Section 4. Test problems and optimization results are presented in Section 5. Finally, Section 6 presents the main conclusions of the study.

Formulation of the RBDO problem
We consider that the structures studied in this paper are subject to random variables of the following type: cross-sectional areas of the bars A, nodal coordinates C, load magnitudes F, yielding stress of the material in tension S t and compression S c and material elastic modulus E.
The algorithm performs size optimization of the truss by changing the mean value of the cross-sectional areas (a ∈ m ) of the m structural members.Shape optimization is made by modification of the mean value of the nodal coordinates of the nodes q considered as design variables (c ∈ q ).Notice that we employ capital letters for random variables and small letters for deterministic quantities, such as the mean value of the random variables.
The optimization procedure seeks the structure of minimum weight subject to stress, displacement and local buckling reliability constraints.For convenience of notation, the design variables a and c are grouped into the vector d = [a 1 , . . ., a m , c 1 , . . ., c q ] and the random variables are grouped into the vector X = [A, C, F, S t , S c , E].The optimization problem can then be posed as: where w is the structural weight, m is the number of bars, ρ is the specific weight of the bar material, is the length of each bar, P (. ..) is the probability of the event in parenthesis to occur, P f is the maximum allowable probability of failure for each reliability constraint, Ω is the set of available discrete cross-sectional areas, c min i and c max i are lower and upper bounds for the nodal coordinate c i , respectively, and G k are the local limit state functions given by where δ k and δ max k are the displacement and maximum allowable displacement at node i, respectively, S j is the stress of bar j, S t j and S c j are the yielding stresses in tension and compression of bar j, respectively, and S b j is the Euler buckling stress of bar j.

First order reliability method
The main idea of the FORM consists in replacing each limit state function G k by a tangent hyperplane in the most probable point of failure (MPP).The reliability of the system is then approximated using the reliability index, which is usually denoted by β.To evaluate the β of each constraint, it is usual to introduce a vector of normalized and statistically independent random variables U ∈ n and a transformation T so that U = T (X).The mapping T transforms every realization x of X in the physical space into a realization u in the normalized space.Notice that the limit state function can also be written in the normalized space as The reliability index β can be obtained from the following optimization problem in the normalized space: The solution of Equation (3) u * (d) is the MPP, which is defined as the realization of the random vector U that lies over the limit state surface and that is closer to the origin of the normalized space.The reliability index β is defined as the distance from the origin of the normalized space to the MPP.Finally, the failure probability can be approximated as where Φ is the standard Gaussian cumulated function and β t is the target reliability index.In this paper, the probability constraints of Equation ( 1) are approximated using Equation ( 4).For a full description of the FORM, the reader is referred to Melchers [35].
It is important to point out that the optimization problem of Equation ( 3) must be solved to evaluate the reliability index β k of each constraint.Thus, to check the feasibility of a given design vector d, q + m optimization problems must be solved.This procedure leads to very high computational costs, especially when a metaheuristic algorithm is employed for the structural optimization, since it usually requires the evaluation of at least thousands of designs in order to converge to the global optimum of the problem.To overcome this issue, we restate the reliability analysis using a global constraint and employ the adaptive response surface strategy to evaluate the resulting reliability index.

Statement of the reliability problem
Here, we take into account all the reliability constraints at once by using a global constraint that is obtained employing the maximum operator.In the case that all the constraints given by Equation ( 1) are respected, the following condition is also respected: max {G(d, X)} ≤ 0 ( 5 ) where we assume that the vector G is composed by all the limit state functions G k of the problem.Instead of evaluating each probabilistic constraint, we assume that the limit state function of the problem is the one given by Equation (5).That is, all limit state functions are considered at once by taking the maximum value among them.This new limit state function is called here global limit state function in order to put in evidence that we are considering a set of local limit state functions at once.This limit state function can be written as where G global is the global limit state function in the physical space.In the normalized space, we have and the resulting reliability analysis problem is stated as: Note that this may not be an efficient approach when the FORM is applied directly to the problem, since one would have difficulties in dealing with the gradient of H global , since it is not differentiable at some points because of the maximum operator.However, we take advantage here of the use of the adaptive response surface approach since it does not require the evaluation of the gradient of the limit state function, as discussed later.

Adaptive response surface approach
The main idea behind the adaptive response surface approach is to make some approximation for the limit state function and then solve the reliability analysis problem using this approximation [35,[37][38][39].The approximation is then iteratively refined and the procedure repeated until convergence is achieved.
Here, we use a linear response surface of the form where Hglobal is the approximation of H global , u i are the variables of the problem, n is the number of random variables and b i are coefficients to be determined.Once the closed form approximation for the limit state function from Equation ( 9) is built, the reliability analysis problem can be solved efficiently with FORM algorithms.According to Equation ( 9), we note that for n random variables the number of coefficients to be determined is n+1.Consequently, in order to uniquely define these coefficients, we must evaluate H global at least in n+1 different sampling points, denoted here as ũi (i = 1, 2,. . ., n+1).In order to generate the points ũi , we build the approximation centred at some chosen point ũ0 .The sampling points are then obtained as where Δ is a scalar and e i are canonical basis vectors with all its components equal to zero but component i, which is equal to 1.The coefficients a, b i from Equation ( 9) can be found by imposing the conditions (11) which results in a system of n + 1 linear equations with unknowns b 0 , b 1 , b 2 ,. . ., b n .
The adaptive response surface approach is made using the following procedure.At the first iteration, we choose some initial value for Δ as used in Equation (10) and denote it as Δ (0) .We then build an approximation H(0) global centred at the origin of the normalized space ũ(0) 0 = 0. We use this value for ũ(0) 0 since we lack further information about a better choice at this stage of the search.The reliability analysis problem is solved for H(0) global using FORM, thus giving the MPP at the first iteration ũ * (0) .
The value Δ is then reduced by the update rule where k is the iteration number and λ ∈ (0, 1).
With the updated value Δ (k+1) , another approximation H(k+1) global is built, but now centred at the MPP from the last iteration ũ * (k) .This problem is solved thus giving an updated MPP ũ * (k+1) .This procedure is repeated until some convergence criterion is met and the last MPP found u * is taken as the solution to the problem.The convergence criterion can be checked on the change of the design vector, on the change of the reliability index, or in both quantities.More details on this adaptive response surface approach are described by Torii and Lopez [39].

Firefly algorithm (FA)
The Firefly Algorithm (FA) is a very recent metaheuristic optimization algorithm developed by Yang [40] and is inspired by the flashing behaviour of fireflies.According to Yang [40], FA optimization has three idealised rules.
(a) All fireflies are unisex, so that one firefly is attracted to other fireflies regardless of their sex.(b) Attractiveness is proportional to brightness, so for any two flashing fireflies, the less bright firefly will move towards the brighter firefly.Both attractiveness and brightness decrease as the distance between fireflies increases.If there is no firefly brighter than a particular firefly, that firefly will move randomly.(c) The brightness of a firefly is affected or determined by the landscape of the objective function.
Based on these three rules, the basic steps of the FA can be summarised as the pseudo-code shown in Figure 1 [40].
There are two essential components to FA: the variation of light intensity and the formulation of attractiveness.The latter is assumed to be determined by the brightness of the firefly, which in turn is related to the objective function of the problem under study.
As light intensity and attractiveness decrease and the distance from the source increases, the variation of light intensity and attractiveness should be a monotonically decreasing function.For example, the light intensity can be: (13) in which the light absorption coefficient γ is a parameter of the FA and r ij is the distance between fireflies i and j at d i and d j , respectively, which can be defined as the Cartesian distance  attractiveness is proportional to the light intensity seen by other fireflies, it can be defined by: (14) in which η 0 is the attractiveness at r = 0. Finally, the probability of a firefly i being attracted to another, more attractive (brighter) firefly j is determined by: where t is the generation number, ε i is a random vector (e.g., the standard Gaussian random vector in which the mean is 0 and the standard deviation is (1) and α is the randomisation parameter.The first term on the right-hand side of Equation ( 15) represents the attraction between the fireflies and the second term is the random movement.In other words, Equation (15) shows that a firefly will be attracted to brighter or more attractive fireflies and also move randomly.Equation (15) indicates that the user must set parameters η 0 , γ, α and the distribution of ε i to apply the FA, and also shows that there are two limit cases when γ is small or large.
(a) If γ approaches zero, the attractive and brightness are constants, and consequently, a firefly can be seen by all other fireflies.In this case, the FA reverts to the PSO.(b) If γ approaches infinity, the attractiveness and brightness approach zero, and all fireflies are short-sighted or fly in a foggy environment, moving randomly.In this case, the FA reverts to the pure random search algorithm.
Hence, the FA generally corresponds to the situation falling between these two limit cases.The stopping criterion employed in this paper is the maximum number of iterations it max .
One crucial aspect in the application of metaheuristic algorithms is the constraint handling.Here, we propose the following multiplicative penalty scheme to handle the probability constraints: where J is the new objective function to be minimised and P t is directly proportional to the level of unfeasibility of the design d.This scheme employs different penalizations if the design fails in the deterministic case or not.First, before running the reliability analysis, the current design is checked in the deterministic case, i.e. the finite element code is run using the mean value of all the random variables of the problem.If the design fails in the deterministic case, the following penalization is employed in the construction of the objective function of the problem: in which μ is a positive constant and we may see that the penalization is also dependent on the generation number t.Thus, if we set μ relatively close to zero (e.g.μ = 0.1), at the first few generations the designs that violate the deterministic constraints are not highly penalized and it is possible to maintain their characteristics in the pool of designs.It is interesting to do so in order to keep in the pool low area bars, for instance.However, if the design is feasible in the deterministic check, the reliability analysis code is called in order to evaluate β global (d).If the reliability constraint is not fulfilled, P t is given by: , |( • )| stands for the absolute value.Finally, if the reliability constraint is fulfilled J (d) = w (d) .The constraints on the bounds of the nodal coordinates are addressed by a coding approach.These bounds are imposed by not sampling unfeasible designs in the computer code.Thus, the optimization problem that the FA solves in this paper is given by: The next section contains a numerical analysis to demonstrate the effectiveness of the proposed methodology, comprised by the response surface method and the FA, for solving the reliability based shape and size optimization of trusses.

Test problems and optimization results
In this section, we adapted two benchmark examples that have been analysed in the deterministic optimization of truss structures.Due to the stochastic nature of the FA, the final result can vary depending on the seed used for the random number generation.Yet there is no established statistical benchmark criterion in the literature to evaluate the performance of metaheuristics in size and shape optimization of trusses, especially in the case of RBDO.With the goal of providing a statistical basis for further comparison, this paper presents the results of over 50 runs for each example.In this way, the average values and coefficients of variation are presented along with the optimal results.The problems are presented in increasing order of complexity.The parameters of the FA were adopted from Fadel Miguel et al. 0 and were kept constant in all examples: μ = 0.1, η 0 = 1, γ = 1, α = 0.5 and a value of ε i that follows a uniform distribution between -0.5 and 0.5.The parameters employed in the adaptive response surface approach were λ = 0.25 and Δ (0) = 1 (i.e., one time the standard deviation of the random variables).These values were obtained after some investigation by the authors.The FA is compared in this section to a HS [18] and a GA (MATLAB optimization toolbox).In order to make the comparison fair among the optimization algorithms, the total number of calls of the objective function is kept constant for all of them in each example.

Planar 10-bar truss problem
The first test problem is the weight minimization of the planar 10-bar truss structure shown in Figure 2. The structure has 10 elements connected by 6 nodes and is subject to two load cases, i.e. the truss is subject to either load F 1 or F 2 .The total length and total height of the ground structure (see Fig. 2) are 720 cm and 360 cm, respectively.The weight density of the bars is 78.5 kN.m −3 .The moment of inertia of each truss element (used for checking buckling constraints) is related to the cross-sectional area by the following equation: where a is the element cross-sectional area.The maximum allowable nodal displacement is δ max = 2.0 cm.All random variables included in the optimization problem are listed in Table 1.For convenience, load uncertainty was modeled by multiplying the concentrated loads F 1 and F 2 by uncertainty load factors.
Since nodal coordinates are continuous while crosssectional areas must be selected from a discrete set of 38 available values, this test case is a mixed-variable optimization problem including integer and continuous variables.
In the present study, n = 20 fireflies and it max = 3000 optimization iterations were considered: 60 000 reliability analyses were hence carried out.Details of the best design weighing 4.602 kN are given in Table 2.The feasibility of the optimized design was checked using a standard FORM approach which converged to the same reliability index.The corresponding optimized layout is shown in Figure 3.
Table 2 compares the optimized design obtained by using FA as metaheuristic optimizer with those obtained  by using HS and GA.Remarkably, the present algorithm designed the lightest structure with the smallest standard deviation.The mean optimized weight was 4.929 kN with 2.4% standard deviation.Figure 4 shows a typical convergence history obtained in the optimization process.
It should be noted that individual evaluation of each reliability constraint stated in Equation (1) using a standard FORM procedure may lead to have convergence issues and usually entails high computational cost especially in the case of very conservative designs [49].In the present study, we tried to solve the 10-bar problem with a standard FORM approach but several designs diverged because of the above mentioned reason and computational cost was unaffordable.This gives further evidence of the validity of the proposed approach.

Planar 21-bar tower problem
The second test problem considered in this study is the optimization of the planar 21-bar tower shown in Figure 5.The structure includes 21 elements connected by 10 nodes.The ground structure has base width of 2 m and total height of 4 m (see Fig. 5).The truss is made of steel, with mass density of 8002 kg.m −3 .The structure must be designed to withstand three independent loading conditions depending on concentrated forces F 1 = 100 kN and F 2 = 25 kN.The maximum allowable nodal displacement is δ max = 1 cm.
The moment of inertia of each truss element (used for checking buckling constraints) is related to the crosssectional area by the following equation All random variables included in the optimization problem have lognormal distributions (see Tab. 3).For convenience, load uncertainty was modelled by multiplying the concentrated loads F 1 and F 2 by uncertainty load factors.Simultaneous layout and sizing optimization of the planar tower was performed so to reach the target optimum reliability index β t = 3.0.The optimization variables chosen for this test problem are the mean values (a)    The mean values of element cross-sectional areas can be selected from the following set of discrete values Ω = (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360,370,380) cm 2 .Nodal coordinates taken as optimization variables can vary by + or -0.8 m with respect to the initial layout of the ground structure shown in Figure 5.The origin of the X-Y coordinate system coincides with node 1. 60 000 reliability analyses (resulting from n = 20 fireflies and it max = 3000 optimization iterations) were carried out also for this test problem.Details of the best design weighing 103.6 kN are given in Table 4.The feasibility of the optimized design was checked using a standard FORM approach which converged to the same reliability index.The corresponding optimized layout is shown in Figure 6.
Table 4 compares the optimum designs obtained by using FA, HS and GA as metaheuristic engines.The RBDO algorithm developed in this study again designed the lightest structure with the smallest statistical dispersion on optimized weight.The mean optimized weight was 110.4 kN with coefficient of variation equal to 5.6%.
Figure 7 shows a typical convergence history obtained in the optimization process.
Similar to the 10-bar truss problem, the standard FORM approach was not efficient also in this test problem: in fact, several designs diverged and computational cost was unaffordable.This confirms the validity of the proposed RBDO framework.

Conclusions
This study presented a novel methodology to perform global reliability based design optimization of truss structures including sizing and shape variables.The methodology combined global constraint formulation and response surface approximation to deal with reliability analysis and a firefly algorithm to carry out global optimization.In this way, it was possible to reduce the computational cost entailed by the evaluation of probabilistic constraints as well as to deal with design space non-convexity.
The proposed RBDO framework was successfully tested in two weight minimization problems of truss structures.The FA metaheuristic search engine outperformed other popular algorithms like harmony search and genetic algorithms designing lighter structures with lower statistical dispersion on optimized weight.Since there are no general benchmark criteria to compare relative merits of metaheuristic algorithms in sizing and shape RBDO optimization of truss structures, performance of FA, HS and GA was statistically evaluated in terms of average optimized weight and coefficient of variation of weight obtained over 50 independent optimization runs.
Optimized designs were also checked using a classical FORM approach, thus validating results provided by the response surface approximation method developed in this study.The implementation of the proposed approach not only allowed the computational cost of global RBDO to be drastically reduced (in particular, reliability analysis), but also overcame the typical limitation of FORM poor convergence.However, all other FORM drawbacks are not eliminated by the present formulation: for example, the accuracy achieved in the approximation of the failure probability for non-linear limit state functions and nonnormal random variables.The coupling of more precise reliability methods with global optimization algorithms is very important and should be investigated in future research.

Fig. 3 .
Fig. 3. Best result for the size and shape optimization of the planar 10-bar truss problem.

Fig. 4 .
Fig. 4. Convergence history for the size and shape optimization of the planar 10-bar truss problem.

Fig. 7 .
Fig. 7. Convergence history for the size and shape optimization of the 21-bar tower problem.
Because a firefly's

Table 1 .
Random variables of planar 10-bar truss problem.

Table 2 .
Optimum size and shape solution for the planar 10bar truss problem.

Table 3 .
Random variables of 21-bar tower problem.
of cross-sectional areas of truss elements and horizontal coordinates of free nodes.Because of structural symmetry, the number of nodal coordinates included as design variables was reduced to 4 (i.e.x 2 , x 3 , x 4 and x 5 ) and the number of areas was reduced to 13 (see Tab. 4).Therefore, this test case is a mixed integer-continuous optimization problem with 17 optimization variables (i.e. 13 sizing and 4 layout variables).

Table 4 .
Best design of the 21-bar tower problem.