Past destructive earthquakes (e.g. the 1994 Northridge earthquake and the 1995 Kobe earthquake) have left a clear signature on the engineering community worldwide, changing thinking of structural engineers [1-2]. As such, after holding several workshops and conferences, an innovative approach namely Performance-Based Design (PBD) was presented by modern guidelines [3-5]. In principle, a structure designed using PBD approach should meet performance objectives in accordance with a set of specified reliabilities over its service life. This is aimed to reach structural design candidates associated with more predictable seismic behavior, quantifying and controlling the risk at an engineered acceptable level.
Both seismic demands and capacity parameters, that are inherently uncertain, are highly influential on the acceptable performance level of a structure. Furthermore, due to the fact that a structure on underlying soil is not rigid, soil-structure interaction (SSI) affects the responses of structures during an earthquake. Obviously, ignoring the SSI effects could lead to unrealistic structural responses and seismic demands. Hence, the effects of SSI should be considered in the seismic responses of structures . Therefore, soil type, material properties of the structure, and ground motion characteristics randomly affect the seismic structural responses.
Deterministic structural optimization without considering the uncertainties in design manufacturing and operating processes may lead to unreliable design resulting in inappropriate balance between cost and safety. A proper design procedure must reasonably account for the inherent uncertain nature of a structural system and associated external load . In structural optimization, non-deterministic performance of structures can be taken into account using robust design optimization (RDO)  and reliability-based design optimization (RBDO) . RDO aims to minimize variation of the objective function, but RBDO optimizes the structural cost under reliability of the constraints.
A few studies have been implemented in a structural optimization problem, where RBDO is incorporated into PBD concept. Foley
Nonlinear dynamic analysis of structures using finite element method requires much computational effort. This drawback may accentuate when the nonlinear dynamic structural responses are required in RBDO of the structure using the Monte-Carlo Simulation (MCS) method, the importance sampling technique and the response surface method. In order to obtain an acceptable confidence within probabilities of the order close to 10-4 - 10-6, the MCS method requires a large number of structural analyses. Based on Lagaros
In this chapter, RBDO of RC structures with considering SSI effects under time-history earthquake loading is presented in accordance with the PBD concept of SEAOC guidelines . In this work, a new discrete gravitational search algorithm (DGSA) and an efficient proposed meta-model were introduced for performing RBDO of RC structures . The objective function is the total cost of the structure while the constraints are treated as deterministic and probabilistic. The annual probability of non-performance for each performance level is considered as the probabilistic constraint in RBDO procedure. The new DGSA based on the fundamental concept of the standard GSA  is introduced for finding the optimal designs in the RBDO procedure. In DGSA, the position of each agent is presented in positive integer numbers. Also, the velocity of each agent is modified based on the particle swarm optimizer with passive congregation (PSOPC) which was proposed by He
Numerical examples show that the wavelet as a kernel function is much better than those of the common kinds as kernel function in WLS-SVM. The accuracy and generalization of WWLS-SVM is improved using GSA. Furthermore, numerical results demonstrate the efficiency and computational advantages of the proposed DGSA for RBDO of structures.
2. RBDO of RC structures
2.1. Formulation of optimization
Seismic design optimization of RC structures under time-history earthquake loads is an ongoing research topic and has received great attention among researchers [14, 27-33]. As such, RBDO of RC structures with the consideration of SSI effects was investigated in accordance with PBD concept of SEAOC guidelines  under seismic loading. This work incorporates the acceptable performance levels and the RBDO theory to compare the achieved annual probability of non-performance with target values for each performance level. The objective of the RBDO problem is to minimize the total cost whereas the deterministic and probabilistic constraints should not exceed a specified target.
The RBDO problem of RC structures can be formulated in the following form:
2.2. Life- cycle cost assessment of RC structure
The total cost,
where and are the unit cost coefficients of each material;
The repair cost refers to the cost of damage level from earthquake that may occur during the life of a structure. In this study, the overall damage index,
where is the probability density function for the index,
2.3. Constraint handling approach
A comprehensive overview of the most popular constraint handling approaches used in conjunction with meta-heuristic optimization methods was presented in the literature review by Coello Coello . In the present study, the external penalty function method as one of the most common forms of the penalty function in the structural optimization [15, 27-29, 36-39] is employed to transform constrained RBDO problem into unconstrained one as follows:
This formulation allows solutions with violated constraints, and the objective function is always greater than the non-violated one.
3. Reliability assessment of RC structure
In order to evaluate the system reliability index corresponding to each of the performance levels, RC structures should be assessed in the RBDO procedure [14, 22]. The system reliability index corresponding to each of the performance levels are estimated by MCS method. In the following subsections, the procedure of assessment of RC structures is explained.
3.1. Required database
In PBD approach, many uncertain variables influence the structural seismic responses. In the studies by Khatibinia
3.2. Limit state functions
The operational, life safety and collapse prevention levels have been defined as the performance levels. A performance level depends on some limit state functions. A limit state function, , is determined by capacity, , and demand, , as follows [13, 22]:
where is the limiting value for a seismic response at a given performance level, with mean value of and coefficient of variation .
The limit state functions and their probability distribution function (PDF) for the performance levels, according to SEAOC guidelines (2000), are shown in Table 1.
|Elastic roof displacement||Lognormal|
|Max. local damage index||Beta|
|Global damage index||Beta|
|Max. local damage index||Beta|
|Global damage index||Beta|
In Table 1, is the yielding horizontal displacement at top storey of the frame which is determined by a pushover analysis. Furthermore, the demand, , corresponding to each of seismic responses are defined using the mean, , and the standard deviation,
3.3. Annual probability of non-performance
The non-performance probability, , is considered as a function of the limit state functions in proportion to a specified performance level. Using the evaluation of the multiple integral over the failure domain, , is calculated as follows:
where is the joint probability density function of .
For each performance level, the total exceeding probability, , is considered as a series system. Determination of the total exceeding probability, , is based on integration of a multi-normal distribution function. In order to estimate the integral, the MCS method is used concurrently for all limit state functions in proportional to the performance levels listed in Table 1. Therefore, the seismic reliability corresponding to the each performance level is defined by an annual probability of non-performance. The annual probability of non-performance, , is computed using the occurrence of earthquakes as a Poisson process [13, 22]:
where can be expressed as an reliability index as shown in Eq. (11), using the standard normal cumulative distribution function, .
4. Seismic responses and SSI system
4.1. Seismic responses of SSI system
There are two main approaches for modeling and analyzing SSI systems, namely the direct method and the substructure method either in time domain or in frequency domain . Considering the discretized dynamic equations of structure and soil simultaneously, the direct method models the soil and structure together, and the responses of soil and structure are determined simultaneously by analyzing SSI system in each time step .
In the direct method, the discretization of nonlinear dynamic equations can be expressed in FEM framework as:
where and are mass, damping, and tangent stiffness matrices of SSI model, respectively; is the incremental vector of the relative displacements for SSI system between times and ; and is the vector of internal forces at time
Over the past two decades, the use of damage and energy concepts for the seismic performance evaluation and design of structures has attracted considerable attention among the researchers [30, 44-47]. These concepts can be simultaneously used through a combined damage index namely Park-Ang damage index. The index is taken into account as a combined index, defined as the linear combination of the maximum displacement and the hysteretic energy dissipation for a structural element. For this reason, the damage index  is one of the indices that have widely been used for damage assessment and damage-based design of RC structures [14, 22, 28-30, 45-46, 48]. As shown in Table 1, some limit states of the performance levels depend on the damage indices.
An improved version of the index namely modified Park-Ang damage index  is defined based on the cross-section deformation of structural elements as:
where is the maximum rotation during loading history; and are the ultimate and yield rotation, respectively; is the yield moment; and is the hysteretic energy dissipated in the same cross-section. Two connected indices, storey and overall damage indices, are computed using the weighting factors based on dissipated hysteretic energy at components and storey levels, respectively, as follows:
where is energy weighting factor; and is total absorbed energy by the component or storey
4.2. Finite element model of SSI system
OpenSEES , as an open-source computational software framework, is used for by simulation of SSI system, and performing nonlinear dynamic analyses of SSI system depicted in Fig. 1. Assuming materials of constant properties over its depth, soil encompasses different layers, and the foundation is considered as rigid strip footing. Beams and columns of structure are modeled using force-based nonlinear beam-column element with considering the spread plasticity along the element’s length. The integration along each element is based on Gauss-Lobatto quadrature rule. Also, the infinite boundaries of soil are modeled using the artificial boundaries (Fig. 1). The model of soil-structure system shown in Fig. 1 was successfully used by [14, 22, 31].
The Kent-Scott-Park model  is utilized for modeling the confined and unconfined concrete of cross-sections of structural elements. The constitutive parameters of this model are:
where , and are the ultimate compressive strain of confined and unconfined concrete, and the ultimate strain of longitudinal steel reinforcement in tensile stress, respectively; and are the volumetric ratio, and the yield stress of confining steel reinforcement, and the peak strength of confined concrete in compression, respectively.
In this study, the one-dimensional
Soil layers are modeled using isoperimetric four-node quadrilateral finite elements and assuming bilinear displacement interpolation. The plane strain condition is assumed for the soil domain with considering a constant soil thickness corresponding to the inter-frame distance. The material of the soil is modeled using a modified pressure-independent multi-yield-surface
One of the major problems in SSI system for infinite media has been the modeling of the domain boundaries. Infinite boundaries have to absorb all outgoing waves and reflect no waves back into the computational domain. In this study, the standard viscous boundary proposed by Lysmer and Kuhlemeyer  is used for this purpose. This boundary can be described by two series of dashpots oriented normal and tangential to the boundary of a finite element mesh (Fig. 1) as follows:
where and are the normal and shear damping, respectively; and are the mass density and Poisson ratio of soil, respectively;
The material damping matrix,
5. Artificial earthquakes
For RBDO of structures, it is then necessary to utilize accelerograms of compatible characteristics with a desired site. It is often difficult or impossible in some cases to choose a proper record for a site, since historically recorded accelerograms for a given site could be limited or scare. Hence, artificial earthquakes that are statistically influenced by desired properties of the given site are very useful for seismic design of structures. In this work, spectral representation method based on time domain procedure is used for the generation of synthetic ground motion records. The non-stationary ground motion is simulated using this method as :
where and are the non-stationary ground motion, the modulation function and the specific power spectral density function (PSDF), respectively. is the number of sine functions or frequencies included, between 0 and , and are the coefficient of variation and a standard normal variable that used in ordinates of PSDF, is frequency step, and are random phase angles with a uniform distribution between 0 and . In this work, the modulation function expressed in  is used:
where , and are specific times and the duration of the simulated record, and are constants. Also, the PSDF of the non-stationary ground motion suggested by Clough and Penzien  is considered as:
where is the constant PSDF of input white-noise random process; and are the characteristic ground frequency and the ground damping ratio; and are parameters for a high-pass filter to attenuate low frequency components. As listed in Table 2, the parameters for the generation of simulated ground motion are selected according to the values proposed by Möller et al. .
The PGA values are obtained corresponding to hazard curves and produced for a specific region. As shown in Table 3, in this work the hazard curves presented by Möller et al.  are used. An artificial earthquake generated based on Eq. (19) is shown in Fig. 4.
|50% in 50 years||73||0.27|
|10% in 50 years||475||0.6|
|5% in 50 years||975||0.8|
6. The new discrete gravitational search algorithm
Based on the work presented by Khatibinia
6.1. Gravitational search algorithm
GSA was introduced by Rashedi
To compute the acceleration of an agent, total forces from a set of heavier masses applied to it should be considered based on the law of gravity (Eq. (23)). Afterwards, the next velocity of an agent is calculated as a fraction of its current velocity added to its acceleration (Eq. (24)). Then, its next position could be calculated using Equation (25):
where , and present the acceleration, velocity and position of
6.2. The proposed discrete GSA
The binary GSA (BGSA) for solving discrete problem was developed by Rashedi
Based on Eq. (26), a large computer memory is needed for the position of agents in BGSA. Also, coding and encoding of the position of agents is a time consuming process. In order to overcome the shortcomings of BGSA, a new DGSA based on the fundamental concept of the standard GSA with passive congregation is presented herein. The passive congregation strategy as perturbations operator can transfer information among agents in the optimization procedure . Therefore, the search performance of DGSA can be improved using the passive congregation. To achieve this purpose, Khatibinia
where , and are the best previous position of the
In DGSA, the scalar corresponds to discrete values of the set
Therefore, the coding and encoding of the position of agents are omitted; and the position of agents is calculated as the integer value. The current position of agents may violate from the values of the set
where denotes the integral part function.
7. Approximating the structural seismic responses
MCS requires excessive computational cost for RBDO of structures in order to obtain an acceptable accuracy . Because of the drawback, Khatibinia
7.1. Weighted least squares support vector machines
WLS-SVM was introduced as excellent machine learning algorithms in large-scale problems by Suykens
WLS-SVM is described as the following optimization problem in primal weight space :
Subject to the following equality constraints:
where is a training data set; and represent the input and output data. Operator is a function which maps the input space into a higher dimensional space. The vector, , represents weight vector in primal weight space. The symbols, and , are the error variable and bias term, respectively. Using the optimization problem, Eq. (30), and the training data set, the WLS-SVM model could be expressed as:
It is impossible to indirectly compute from Eq. (30), for the structure of the function is generally unknown. Therefore, the dual problem shown in Eq. (30) is minimized by the Lagrange multiplier method as follows:
Based on the Karush-Kuhn-Tucker (KKT) conditions, by eliminating and the solution is given by the following set of linear equations:
According to Mercer’s condition, a kernel is selected, such that:
Consequently, the final WLS-SVM model for the prediction of functions becomes:
where is a robust estimation of the standard deviation for the error variables (); constants and are typically chosen as and . Here denotes the
In WLS-SVM, Gaussian radial basis function (RBF) is frequently used as the kernel function, and it is expressed as:
where is a positive real constant usually called the kernel width.
Based upon the Suykens
7.2. The new meta-model-based wavelet kernel
Wavelets as kernel function have been introduced and developed in ANNs and SVMs [60-63]. It has been shown that wavelet kernel functions are superior to other kernel functions in the training ANN and SVM. Accordingly, the kernel function of WLS-SVM is substituted with a specific kind of wavelet functions proposed by Khatibinia et al. . The meta-model based on wavelet kernel function is called WWLS-SVM. The cosine-Gaussian Morlet wavelet is used as the kernel function of WLS-SVM. The wavelet function is mathematically written as follows:
According to Zhang
Therefore, according to Eqs. (40) and (41), the wavelet kernel function of the cosine-Gaussian Morlet wavelet is given as follows:
The accuracy of WWLS-SVM prediction depends on the good selection of its parameters. Selecting appropriate values of these parameters is important for obtaining the excellent predicting performance. Hence, in this study, GSA is used to find the WWLS-SVM optimal parameter, , and the wavelet kernel parameters,
The WWLS-SVM training stage during GSA is performed according to the
The converged solution is affected by the setting value of parameters in GSA. In this study, the values are selected based on the general recommendations by Rashedi
8. Predicting failure probability of structures
In the RBDO procedure, nonlinear time-history analysis of SSI system is used and it may be failed regarding a number of random structures . In fact, a number of structures collapse and then lose their stability. Hence, these structures should be identified and eliminated from optimization process. For this purpose, a failure probability is considered as stability criterion. An efficient method is presented to train the failure probability with high performance . This efficient method is consisted of a modified adaptive neuro fuzzy inference system (ANFIS) with a hybrid of fuzzy c-means (FCM)  and fuzzy particle swarm optimization (FPSO) . To train the modified ANFIS, the input–output data are classified by a hybrid algorithm consisting of FCM-FPSO clustering. The optimum number of ANFIS fuzzy rules is determined by subtractive algorithm (SA).
8.1. Hybrid of FCM and FPSO for clustering
The FCM algorithm has been extensively studied and is known to converge to a local optimum in nonlinear problems. Moreover, the FPSO algorithm is robust method to increase the probability of achieving the global optimum in comparison with the FCM algorithm. The FCM algorithm is faster than the FPSO algorithm because it requires fewer function evaluations. This shortcoming of FPSO can be dealt with selecting an adequate initial swarm .
In this study, a hybrid clustering algorithm called FCM-FPSO is presented to use the merits of both FCM and FPSO algorithms and increase the procedure of convergence. In this way, the FCM algorithm finds an adequate initial swarm FPSO algorithm for commencing the FPSO. For this purpose, first, the FCM algorithm is utilized to find a preliminary optimization that shown by . This optimum solution is copied
8.2. Modified ANFIS
An ANFIS model depends on the number of ANFIS fuzzy rules and membership functions. In other words, creating an ANFIS model with a minimum number of fuzzy rules can eliminate a well-known drawback. Therefore, for overcoming of this drawback, Khatibinia
9. Numerical examples
In this work, two RC frame structures shown in Fig. 8 are selected as illustrative examples. Three layers of sand associated with material properties varying over its depth are considered as the soil under the frames. The depth of each soil layer and the entire width of soil domain are considered to be 10 m and 100 m, respectively. The soil is also assumed to have plane strain condition with a constant thickness of 5.0 m in proportion to the inter-frame distance. Inertia properties of the soil mesh are considered using lumped mass matrices modeling with soil mass density of 17 kN/m3 for all soil layers. The values of the dead and live loads are considered to be
For vertical continuity on the dimensions along the height of a column, the section database of columns is divided into three types in the height of RC frame. Hence, a database shown in Table 4 is generated. Similarly, the section database of beams is divided into three types in the height of RC frame. Distribution of beam dimensions along the height of the frame is shown in Table 4. The diameter of longitudinal bars for beams and columns is laid between 12 mm and 32 mm in the databases.
The initial cost is calculated for and Euros. To calculate the total expected cost of repair, first, the cumulative distribution is obtained using MCS and the proposed meta-model for the response
The concrete material parameters shown in Table 5 are considered for the cover of column cross-sections. The strain corresponding to the peak strength, , and the residual strength,
where and are the mean and the standard deviation of PGA, respectively. The values of these parameters are shown in Table 6.
The presented DGSA requires the user to specify several internal parameters that can affect convergence behavior at the search space. It is found that a population of 50 agents can be adequate. Higher values are not recommended, as this will increase significantly computation time in RBDO. In addition, different optimization runs are carried out for RBDO model in this study, so optimum designs are found by DGSA about 150 iterations. Due to the effect of decreasing gravity, the actual value of the gravitational constant,
9.1. Example 1: Six-storey RC frame
Six-storey RC frame is shown in Fig. 8(a). In the frame, the length of each bay and the height of stories are 5 m and 3 m, respectively. The members of the structure are divided into four groups for the columns C1, C2, C3, C4 and four groups B1, B2, B3 and B4 for the beams. The groups of structural elements are presented in Fig. 8(a).
9.1.1. Training and testing the meta-model
To predict the mean, , and the standard deviation, , of
In order to validate the performance and accuracy of the proposed meta-model, relative root-mean-squared error, i.e.
As given in Table 8, the proposed meta-model trained for the mean and the standard deviation of seismic responses has proper performance generality. Thus, the approximating performance of the meta-model based on WWLS-SVM and GSA is better than the WLS-SVM with RBF kernel in predictive ability and precision.
9.1.2. Results of RBDO
In this example, RBDO of the RC frame is performed using DGSA associated with WWLS-SVM-based MCS. In the reliability process, the reliability indices, , are estimated using WWLS-SVM-based MCS with 106 samples generated with the LHD method. The cross-section of beams and columns are selected from Types 2 and 3, which are shown in Table 4. The optimum designs of the RC frame are listed in Table 9. Furthermore, the optimal solutions of DGSA are also compared with those of BGSA in Table 9.
As shown in Table 9, the optimal solutions of DGSA are better than those of BGSA in terms of the total cost and the number of iterations. The minimum reliability index, , obtained corresponding to each performance level by DGSA and BGSA is shown in Table 9.
The convergence histories of the optimum objective function are shown in Fig. 9 for DGSA and BGSA models. As can be seen in Fig. 9, DGSA method is more efficient than BGSA method. Optimum designs are found by DGSA and BGSA in 4450 and 5900 required approximate analyses by the meta-model, respectively.
9.2. Example 2: Nine-storey RC frame
Nine-storey RC frame is shown in Fig. 8(b). In the frame, the length of each bay and the height of stories are 5 m and 3 m, respectively. The members of the structure are divided into six groups for the columns and six groups for the beams. The groups of structural elements are presented in Fig. 8(b).
9.2.1. Training and testing the meta-model
After training database using the presented WWLS-SVM optimal parameters of the meta-model associated with the mean and the standard deviation of seismic responses are shown in Table 10. Furthermore, the performance generality of the proposed meta-model and WLS-SVM is given in Table 10 in terms of
The results of Table 10 demonstrate that the meta-model is better than the WLS-SVM method in terms of performance generality. Therefore, the meta-model is reliably employed to predict the necessary responses during the RBDO process.
9.2.2. Results of RBDO
As the first example, in this example RBDO of the RC frame is performed using DGSA and BGSA associated with WWLS-SVM-based MCS. In this example, the cross-section of beams and columns are selected from Types 1, 2 and 3, which are shown in Table 4. The best optimum designs of the RC frame are listed in Table 11.
As revealed in Table 11, the optimal solutions of DGSA are better than those of BGSA in terms of the total cost and the number of iterations. The minimum reliability index, , obtained corresponding to each performance level by DGSA and BGSA is shown in Table 11.
The convergence histories of the optimum objective function are shown in Fig. 10 for DGSA and BGSA models. As can be seen in Fig. 10, DGSA method is more efficient than BGSA method. Optimum designs are found by DGSA and BGSA in 4150 and 6150 required approximate analyses by the meta-model, respectively.
In general, the optimum design of structures depends on a number of parameters that are inherently uncertain. Reliability-based design optimization (RBDO) has been employed as the only method that assesses the influence of uncertain parameters and balance both cost and safety of structures. To account for all necessary uncertain and random parameters in RBDO of RC structures and to achieve the realistic optimum design of RC structures, the uncertain material properties of soil and structure, as well as the characteristics of ground motions should be considered as random parameters. Furthermore, the realistic seismic responses of RC structures can be account by consideration of soil-structure interaction (SSI) effects. In this work, a new discrete gravitational search algorithm (DGSA) and a new meta-modeling framework were incorporated for RBDO of RC structures with Performance-Based Design (PBD) under seismic loading. The objective of the RBDO problem was to minimize the total cost whereas the deterministic constraints and the system reliability index corresponding to each of the performance levels should not exceed a specified target. Based on this study, the following conclusions can be drawn:
To reduce the computational effort and computational cost of the Monte-Carlo Simulation (MCS) method, a new meta-model based on a wavelet weighted least squares support vector machine (WWLS-SVM) and gravitational search algorithm (GSA) was utilized in the RBDO procedure. Therefore, the proposed meta-model, as a substitute for the nonlinear dynamic analysis of SSI system, can estimate the reliability index through MCS with a small computational cost.
The WWLS-SVM and kernel parameters were simultaneously optimized in the proposed meta-model in order to improve performance generality of WWLS-SVM. Numerical results of training and testing the meta-model indicated that performance generality of the meta-model was higher in comparison to WLS-SVM. Hence, the proposed meta-model can predict the nonlinear dynamic analysis of SSI system in terms of accuracy and flexibility.
The proposed DGSA was presented based on the standard GSA with passive congregation. The passive congregation strategy can be considered as perturbations operator in the optimization procedure. Therefore, the presented DGSA using the passive congregation can transfer information among agents avoiding local minima. Furthermore, the coding and encoding of the position of agents as a time consuming process is omitted in DGSA. To eliminate this drawback, the position of agents was calculated as the integer value. The optimum designs obtained by DGSA were compared with those produced by BGSA model. Numerical examples showed that the proposed DGSA can converge and reach the optimum design more quickly than the BGSA model.
Future extension of current research could include reducing the computations involved in the PBD by replacing MCS with the response surface method or the importance sampling technique. The constraints imposed on the objective function could be also treated as random quantities (see ).
Our special thanks go to Dr. Eysa Salajegheh (Distinguished Professor of Structural Engineering) in Department of Civil Engineering at Shahid Bahonar University of Kerman, Iran, for his cooperation in this research work.