**Seismic Reliability-Based Design Optimization of Reinforced Concrete Structures Including Soil-Structure Interaction Effects**

Mohsen Khatibinia, Sadjad Gharehbaghi and Abbas Moustafa

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/59641

## **1. Introduction**

Past destructive earthquakes (e.g. the 1994 Northridge earthquake and the 1995 Kobe earth‐ quake) have left a clear signature on the engineering community worldwide, changing thinking of structural engineers [1-2]. As such, after holding several workshops and confer‐ ences, 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 [6]. 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 manu‐ facturing 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 [7]. In structural optimization, non-deterministic performance of structures can be taken into account using

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

robust design optimization (RDO) [8] and reliability-based design optimization (RBDO) [9]. 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 *et al*. [10] proposed a state-of-the-art model code and a PBD methodology. The methodology was applied to multiple-objective optimization prob‐ lems for single storey and multi-storey structural frameworks with fully and partially restrained connections. Lagaros *et al.* [11] introduced a tool for the solution of realistic structural optimization problems that incorporate PBD under seismic loading into RBDO. The tool consisted of two distinctive methodologies based on artificial neural networks (ANNs). Fragiadakis *et al*. [12] presented optimum seismic design of reinforced (RC) structures considering the reliability constraints and PBD, where RBDO was implemented by the evolution strategies (ES). In the study by Moller *et al.* [13], the concept of PBD with RBDO was implemented by consideration of the uncertainties in the structural demands and capacities in order to evaluate reliability associated with each of the required performance levels. Khatibinia *et al*. [14] introduced RBDO of RC structures including SSI effects. In their study, the uncertainty of the structural demand was in terms of the random properties of the structure, underlying soil and the uncertain characteristics of artificially generated earthquakes. Also, the structural capacity associated with each of the required performance level in the concept of PBD was treated as an uncertain quantity.

Nonlinear dynamic analysis of structures using finite element method requires much compu‐ tational effort. This drawback may accentuate when the nonlinear dynamic structural re‐ sponses 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 *et al*. [11], for such results, the large number of analyses ranges from 1.6×10<sup>6</sup> to 1.6×10<sup>9</sup> is required to achieve 95% likelihood for actual probability. To eliminate such drawbacks, utilizing soft computing-based models (e.g. artificial neural networks (ANNs)) is of crucial importance. These models can efficiently approximate structural responses and limit state functions in reliability analysis and optimization [11, 13, 15, 16]. Recently, as another model, support vector machines (SVMs) due to their simplicity, ease of implementation, and good performance have been developed to represent actual limit state functions. A hybrid of SVM and genetic algorithm, called (SVM-GA), was proposed to forecast reliability in engine systems [17]. The results showed the feasibility of SVM-GA in the reliability prediction compared with those of the ANNs and the autoregressive integrated moving average model. Zhiwei and Guangchen [18] investigated the capability of least squares support vector machine-based MCS (LSSVM-MCS) rather than SVM-MCS in reliability analysis. Based on the results of their study, LSSVM-MCS was more accurate and required less computational effort in comparison with SVM-MCS. Tan *et al.* [19] stipulated that there is no difference between the performance of the SVM-based response surface method (SVM-RSM) and the radial basis function neural network-based response surface method (RBFN-RSM). As a comparative study, Moura *et al*. [20] assessed the SVM effectiveness in forecasting time-to-failure and reliability of engineered components based on time series data. The efficacy of SVM with respect to other learning methods was shown. In the context of structural reliability assessment, a sampling method based on the adaptive Markov chain simulation and support vector density estimation was also developed by Dai *et al*. [21]. In their study, the application of SVM was proposed as a density estimator for structural reliability analysis.

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 [3]. In this work, a new discrete gravitational search algorithm (DGSA) and an efficient proposed meta-model were introduced for performing RBDO of RC structures [22]. 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 [23] 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 *et al*. [24]. The modifications can improve the global exploration ability of DGSA and overcome the shortcomings of the Binary GSA (BGSA) model introduced by Rashedi *et al*. [25]. A meta-model-based MCS method is also presented herein to take into account the probabilistic constraint in conjunction with the nonlinear finite element analysis (FEM) of SSI system. Due to the fact that the computational cost of MCS for structural reliability analysis is high, the meta-model is proposed to predict the structural seismic responses of SSI system, and significantly reduce the computational effort of the RBDO process. The meta-model is a combination of weighted least squares support vector machine (WLS-SVM) [26] and Morlet wavelet kernel function, which is called WWLS-SVM [14, 22, 27]. The selection of WWLS-SVM parameters efficiently affects the prediction accuracy of WWLS-SVM. Hence, the parameters of WWLS-SVM and wavelet kernel are assigned by using the standard gravitational search algorithm (GSA).

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**

robust design optimization (RDO) [8] and reliability-based design optimization (RBDO) [9]. RDO aims to minimize variation of the objective function, but RBDO optimizes the structural

268 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

A few studies have been implemented in a structural optimization problem, where RBDO is incorporated into PBD concept. Foley *et al*. [10] proposed a state-of-the-art model code and a PBD methodology. The methodology was applied to multiple-objective optimization prob‐ lems for single storey and multi-storey structural frameworks with fully and partially restrained connections. Lagaros *et al.* [11] introduced a tool for the solution of realistic structural optimization problems that incorporate PBD under seismic loading into RBDO. The tool consisted of two distinctive methodologies based on artificial neural networks (ANNs). Fragiadakis *et al*. [12] presented optimum seismic design of reinforced (RC) structures considering the reliability constraints and PBD, where RBDO was implemented by the evolution strategies (ES). In the study by Moller *et al.* [13], the concept of PBD with RBDO was implemented by consideration of the uncertainties in the structural demands and capacities in order to evaluate reliability associated with each of the required performance levels. Khatibinia *et al*. [14] introduced RBDO of RC structures including SSI effects. In their study, the uncertainty of the structural demand was in terms of the random properties of the structure, underlying soil and the uncertain characteristics of artificially generated earthquakes. Also, the structural capacity associated with each of the required performance level in the concept

Nonlinear dynamic analysis of structures using finite element method requires much compu‐ tational effort. This drawback may accentuate when the nonlinear dynamic structural re‐ sponses 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 *et al*. [11], for such

likelihood for actual probability. To eliminate such drawbacks, utilizing soft computing-based models (e.g. artificial neural networks (ANNs)) is of crucial importance. These models can efficiently approximate structural responses and limit state functions in reliability analysis and optimization [11, 13, 15, 16]. Recently, as another model, support vector machines (SVMs) due to their simplicity, ease of implementation, and good performance have been developed to represent actual limit state functions. A hybrid of SVM and genetic algorithm, called (SVM-GA), was proposed to forecast reliability in engine systems [17]. The results showed the feasibility of SVM-GA in the reliability prediction compared with those of the ANNs and the autoregressive integrated moving average model. Zhiwei and Guangchen [18] investigated the capability of least squares support vector machine-based MCS (LSSVM-MCS) rather than SVM-MCS in reliability analysis. Based on the results of their study, LSSVM-MCS was more accurate and required less computational effort in comparison with SVM-MCS. Tan *et al.* [19] stipulated that there is no difference between the performance of the SVM-based response surface method (SVM-RSM) and the radial basis function neural network-based response surface method (RBFN-RSM). As a comparative study, Moura *et al*. [20] assessed the SVM effectiveness in forecasting time-to-failure and reliability of engineered components based on time series data. The efficacy of SVM with respect to other learning methods was shown. In

to 1.6×10<sup>9</sup>

is required to achieve 95%

cost under reliability of the constraints.

of PBD was treated as an uncertain quantity.

results, the large number of analyses ranges from 1.6×10<sup>6</sup>

#### **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 [3] under seismic loading. This work incorporates the acceptable performance levels and the RBDO theory to compare the achieved annual proba‐ bility 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:

$$\begin{aligned} \text{Minimize } & \quad \mathsf{C}\_{\text{TOT}}\\ \text{Subject to } & \quad \mathcal{g}\_{i} \Big( \mathbf{X}, \overline{\mathbf{X}}\_{\text{rand}} \big) \leq 0 & \quad i = 1, 2, \dots, N\_{\text{g}}\\ & \quad \mathcal{P}\_{\text{annual}} \Big[ \mathsf{G}\_{/} \Big( \mathbf{X}, \overline{\mathbf{X}}\_{\text{rand}} \big) \geq 0 \Big] \geq \mathcal{P}\_{\text{j}}^{\text{Tayst}} & \quad j = 1, 2, 3\\ & \quad \mathbf{X} = \{ \mathbf{x}\_{1}, \mathbf{x}\_{2}, \dots, \mathbf{x}\_{k}, \dots, \mathbf{x}\_{n} \} \in \mathbb{R}^{d} \end{aligned} \tag{1}$$

where *CTOT* is the total cost; *βannual* is the system reliability index corresponding to the *j*th performance level (performance function), *Gj* (*<sup>X</sup>* , *<sup>X</sup>***¯** *rand* ) ; *β<sup>j</sup> Target* is the prescribed target value of the reliability index; *Ng* is the number of deterministic constraints, *gi* (*X* ) ; A given set of discrete values, *X*, is expressed by *Rd* and design variables, *xk,* can take values only from this set. The vector *X***¯** *rand* represents the random variables.

#### **2.2. Life- cycle cost assessment of RC structure**

The total cost, *CTOT*, of a structure is the initial structural cost for a new structure construction and the repair cost from an earthquake and different levels of damage that may occur during the life of structure. This cost can be expressed as a function of the design vector *X* and the time *t* as follows:

$$\mathbb{C}\_{\text{TOT}}\left(\mathbf{X}\_{\prime}, \overline{\mathbf{X}}\_{\prime \text{und}}, t\right) = \mathbb{C}\_{\text{IC}}\left(\mathbf{X}\_{\prime}, \overline{\mathbf{X}}\_{\prime \text{und}}\right) + \mathbb{C}\_{\text{RC}}\left(\mathbf{X}\_{\prime}, \overline{\mathbf{X}}\_{\prime \text{und}}, t\right) \tag{2}$$

where *CIC* is the initial cost of structure; *CRC* is the present value of the repair cost. In this study, the initial cost is considered as the sum of the total cost of concrete, *CC*, and the total cost of reinforcing steel, *CS*, is given [14]:

$$\mathbf{C}\_{\rm IC} = \mathbf{C}\_{\rm C} + \mathbf{C}\_{\rm S} = \sum\_{i=1}^{\rm Ne} \mathbf{w}\_{\rm C} b\_{i} h\_{i} L\_{\gamma} + \sum\_{i=1}^{\rm Ne} \mathbf{w}\_{\rm S} A\_{\rm Si} L\_{i} \tag{3}$$

where *wC* and *wS* are the unit cost coefficients of each material; *bi* , *hi,* and *Li* are the section dimensions of *i*th element and its length; and *ASi* is the total reinforcement in the section of element. *Ne* is the number of the elements of the structure.

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, *DIoverall*, is considered as an indicator of structural damage. The total expected cost of repair based on the overall damage index is expressed as follows [13]:

$$\mathbf{C}\_{\rm RC} = \int\_{0}^{1} \mathbf{C}\_{\rm RC} \Big|\_{\rm DI\_{overall}} \cdot f\_{\rm DI\_{overall}} \left(\mathbf{D}I\_{\rm overall}\right) \cdot d\left(\mathbf{D}I\_{\rm overall}\right) \tag{4}$$

where *<sup>f</sup> DI overall* is the probability density function for the index, *DIoverall*. *CRC* <sup>|</sup> *DI overall* is the expected present cost which is defined as [34]:

Seismic Reliability-Based Design Optimization of Reinforced Concrete Structures Including Soil-Structure… http://dx.doi.org/10.5772/59641 271

$$\left. \mathcal{C}\_{RC} \right|\_{DI\_{overall}} = \int\_0^\infty \mathcal{C}\_{f0} \left( D I\_{overall} \right) f(t) dt = \mathcal{C}\_f \left( D I\_{overall} \right) \frac{\nu}{r + \nu} \tag{5}$$

$$\mathcal{C}\_{f}\left(DI\_{\text{overall}}\right) = \begin{cases} k\,\mathcal{C}\_{0}\left(\frac{DI\_{\text{overall}}}{0.60}\right) & \text{if} \quad DI\_{\text{overall}} \le 0.60\\ & k\,\mathcal{C}\_{0} & \text{if} \quad DI\_{\text{overall}} > 0.60 \end{cases} \tag{6}$$

where *υ* and *r* are the mean occurrence of earthquakes and the discount rate, respectively. *C0* is the complete replacement cost, with *k* =1.20 assumed to be a factor to account for demolition and clearing. In this study, the value of 0.04 is assumed for the discount rate, *r*.

#### **2.3. Constraint handling approach**

The RBDO problem of RC structures can be formulated in the following form:

270 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

( )

*X X*

*annual j rand j*

*xx x x R*

{ , ,..., ,..., }

of the reliability index; *Ng* is the number of deterministic constraints, *gi*

*rand* represents the random variables.

where *wC* and *wS* are the unit cost coefficients of each material; *bi*

*RC* 0 *RC DI DI overall overall*

expected present cost which is defined as [34]:

element. *Ne* is the number of the elements of the structure.

= Î

Subject to , 0 1,2,...,

*k n*

, 0 1,2,3

*d*

 b

*Target*

*rand* ) ; *β<sup>j</sup>*

*C tC C t TOT rand IC rand RC rand* (*XX XX XX* ,, , ,, ) = ( ) + ( ) (2)

=+= + å å (3)

(1)

*Target* is the prescribed target value

(*X* ) ; A given set of

, *hi,* and *Li* are the section

is the

is the system reliability index corresponding to the *j*th

*i rand g*

£ =

*g i N*

(*<sup>X</sup>* , *<sup>X</sup>***¯**

discrete values, *X*, is expressed by *Rd* and design variables, *xk,* can take values only from this

The total cost, *CTOT*, of a structure is the initial structural cost for a new structure construction and the repair cost from an earthquake and different levels of damage that may occur during the life of structure. This cost can be expressed as a function of the design vector *X* and the

where *CIC* is the initial cost of structure; *CRC* is the present value of the repair cost. In this study, the initial cost is considered as the sum of the total cost of concrete, *CC*, and the total cost of

> 1 1 *Ne Ne IC C S C i i i S Si i i i C C C w bhL wA L* = =

dimensions of *i*th element and its length; and *ASi* is the total reinforcement in the section of

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, *DIoverall*, is considered as an indicator of structural damage. The total expected cost of repair based on the overall damage index is

( ) ( ) <sup>1</sup>

is the probability density function for the index, *DIoverall*. *CRC* <sup>|</sup> *DI overall*

*C C f DI d DI overall overall* =× × ò (4)

*G j*

é ù ³³ = ë û

( )

*X X*

*TOT*

*C*

*X*

performance level (performance function), *Gj*

**2.2. Life- cycle cost assessment of RC structure**

b

Minimize

where *CTOT* is the total cost; *βannual*

reinforcing steel, *CS*, is given [14]:

expressed as follows [13]:

where *<sup>f</sup> DI overall*

set. The vector *X***¯**

time *t* as follows:

1 2

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 [35]. 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:

$$\text{fit}\left(\mathbf{X}, \overline{\mathbf{X}}\_{\text{rand}}\right) = \mathbb{C}\_{\text{TOT}}\left(1 + r\_p PF\right) \tag{7}$$

where *fit*(*X* ), *PF* and *rp* are the modified function, the penalty function, and an adjusting coefficient, respectively. The penalty function based on the violation of normalized constraints [36] is defined as the sum of all active constraints violations as indicated:

$$PF = \sum\_{i} \left[ \max\left(g\_{i}, 0\right)\right]^{2} + \sum\_{i} \left[ \max\left(1 - \frac{\beta\_{\text{annual}, j}}{\beta\_{j}^{Twyt}}, 0\right)\right]^{2} \tag{8}$$

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 *et al*. [14, 22], material properties of concrete, steel and soil, as well as earthquakes are considered as intervening uncertain variables. Because of a few historical records of earthquake for a selected site, selection of a proper ground motion record for a site is often difficult, even impossible in some cases. To overcome this problem, artificial earth‐ quakes, statistically influenced by desired properties of a selected site, are utilized in seismic design of structures. The spectral representation method based on time domain procedure can be used for generation of artificial earthquakes [40]. The proper parameters for generation of artificial earthquakes are selected according by values proposed by Möller et al. [13]. In order to perform RBDO of RC structures using the proposed meta-model-based MCS, samples are generated randomly, and are used to train and test the meta-model. Inputs of the meta-model include the random combinations of intervening variables. In order to generate the database, seven random Peak Ground Acceleration (PGA) values are chosen. The PGA values are equal to 260, 350, 400, 550, 650, 700 and 800 (cm/sec2 ). Accordingly, by using the Latin Hypercube Design (LHD) sampling method [41], and considering the 120 combinations based on the intervening variables for each PGA value, the total number 840 combinations are generated. Then, corresponding to the each PGA of each 840 combinations, five artificial earthquakes, as sub-combinations, with random phase angles are generated corresponding to each PGA value. For each of the 840 combinations, in the first phase, the structure is analyzed subjected to the combination of gravity loads according to ACI code [42]. Steel reinforcement ratios of longi‐ tudinal bars of structural elements' cross-sections, *ρ*, shall be satisfied in accordance with ACI code [42]. Furthermore, when choosing the steel reinforcement ratios, it is verified that they are sufficient to provide adequate strength against the combination of gravity loads. Based on ACI code [42], the strong column-weak beam concept shall be satisfied for every joint of designed structures for earthquake loads. In the second phase, nonlinear dynamic analysis of SSI system is performed for each sub-combination and seismic responses of SSI system are obtained. Using nonlinear dynamic analysis of SSI system, maximum roof displacement, *umax*, maximum inter-storey drift, *DRmax*, maximum local damage index, *DILmax*, and overall damage index, *DIoverall*, are considered as the structural seismic responses. After that, the mean, *R*¯ *i* , and the standard deviation, *σ<sup>R</sup> <sup>i</sup>* , of *i*th seismic response, *Ri* , corresponding to each of the 840 combinations subjected to five artificial earthquakes, are achieved. In this work, the modified Park-Ang damage index [43] as one of the most acceptable indices in seismic damage analysis of structures is used for calculating *DILmax* and *DIoverall*.

#### **3.2. Limit state functions**

The operational, life safety and collapse prevention levels have been defined as the perform‐ ance levels. A performance level depends on some limit state functions. A limit state function, *<sup>G</sup>*(*X***¯**), is determined by capacity, *RLIM* , and demand, *R*(*X***¯**), as follows [13, 22]:

$$\mathbf{G}\left(\tilde{\mathbf{X}}\right) = \mathbb{R}\_{\text{LIM}}\left(\mu\_{\prime}\delta\right) - \mathbb{R}\left(\tilde{\mathbf{X}}\right) \tag{9}$$

where *RLIM* is the limiting value for a seismic response *R*(*X***¯**) 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.


**Table 1.** The limit sate functions of the performance levels [14, 22]

**3.1. Required database**

to 260, 350, 400, 550, 650, 700 and 800 (cm/sec2

the standard deviation, *σ<sup>R</sup> <sup>i</sup>*

**3.2. Limit state functions**

of structures is used for calculating *DILmax* and *DIoverall*.

In PBD approach, many uncertain variables influence the structural seismic responses. In the studies by Khatibinia *et al*. [14, 22], material properties of concrete, steel and soil, as well as earthquakes are considered as intervening uncertain variables. Because of a few historical records of earthquake for a selected site, selection of a proper ground motion record for a site is often difficult, even impossible in some cases. To overcome this problem, artificial earth‐ quakes, statistically influenced by desired properties of a selected site, are utilized in seismic design of structures. The spectral representation method based on time domain procedure can be used for generation of artificial earthquakes [40]. The proper parameters for generation of artificial earthquakes are selected according by values proposed by Möller et al. [13]. In order to perform RBDO of RC structures using the proposed meta-model-based MCS, samples are generated randomly, and are used to train and test the meta-model. Inputs of the meta-model include the random combinations of intervening variables. In order to generate the database, seven random Peak Ground Acceleration (PGA) values are chosen. The PGA values are equal

272 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

Design (LHD) sampling method [41], and considering the 120 combinations based on the intervening variables for each PGA value, the total number 840 combinations are generated. Then, corresponding to the each PGA of each 840 combinations, five artificial earthquakes, as sub-combinations, with random phase angles are generated corresponding to each PGA value. For each of the 840 combinations, in the first phase, the structure is analyzed subjected to the combination of gravity loads according to ACI code [42]. Steel reinforcement ratios of longi‐ tudinal bars of structural elements' cross-sections, *ρ*, shall be satisfied in accordance with ACI code [42]. Furthermore, when choosing the steel reinforcement ratios, it is verified that they are sufficient to provide adequate strength against the combination of gravity loads. Based on ACI code [42], the strong column-weak beam concept shall be satisfied for every joint of designed structures for earthquake loads. In the second phase, nonlinear dynamic analysis of SSI system is performed for each sub-combination and seismic responses of SSI system are obtained. Using nonlinear dynamic analysis of SSI system, maximum roof displacement, *umax*, maximum inter-storey drift, *DRmax*, maximum local damage index, *DILmax*, and overall damage index, *DIoverall*, are considered as the structural seismic responses. After that, the mean, *R*¯

, of *i*th seismic response, *Ri*

combinations subjected to five artificial earthquakes, are achieved. In this work, the modified Park-Ang damage index [43] as one of the most acceptable indices in seismic damage analysis

The operational, life safety and collapse prevention levels have been defined as the perform‐ ance levels. A performance level depends on some limit state functions. A limit state function,

*<sup>G</sup>*(*X***¯**), is determined by capacity, *RLIM* , and demand, *R*(*X***¯**), as follows [13, 22]:

*GR R* (*X X* ) = - *LIM* (m d

). Accordingly, by using the Latin Hypercube

*i* , and

, corresponding to each of the 840

, ) ( ) (9)

In Table 1, *u*¯ *<sup>y</sup>* is the yielding horizontal displacement at top storey of the frame which is determined by a pushover analysis. Furthermore, the demand, *R*(*X***¯**), corresponding to each of seismic responses are defined using the mean, *R*¯, and the standard deviation, *<sup>σ</sup><sup>R</sup>* .

#### **3.3. Annual probability of non-performance**

The non-performance probability, *Pf* , 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, *G*(*X***¯**)≤0, *Pf* is calculated as follows:

$$P\_f = \int\_{\mathcal{L}(\overline{X}) \approx 0} \int \dots \int f\_{\overline{X}} \left( \overline{X} \right) d\overline{X} \tag{10}$$

where *<sup>f</sup> <sup>X</sup>***¯**(*X***¯**) is the joint probability density function of *X***¯**.

For each performance level, the total exceeding probability, *P f <sup>E</sup>*, is considered as a series system. Determination of the total exceeding probability, *P f <sup>E</sup>*, 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, *P f annual* , is computed using the occurrence of earthquakes as a Poisson process [13, 22]:

$$Pf\_{\text{annual}} = 1 - \exp\left(-\nu \, Pf\_{\text{E}}\right) \to \, \left. \mathcal{J}\_{\text{annual}} = -\Phi^{-1}\left(Pf\_{\text{annual}}\right) \tag{11}$$

where *βannual* 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 [6]. 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 [22].

In the direct method, the discretization of nonlinear dynamic equations can be expressed in FEM framework as:

$$\mathbf{M}\,\Delta\ddot{\boldsymbol{u}} + \mathbf{C}\,\Delta\dot{\boldsymbol{u}} + \mathbf{K}\_{\mathrm{T}}\,\Delta\mathbf{u} = -\boldsymbol{m}\_{\boldsymbol{x}}\,\ddot{\boldsymbol{u}}\_{\boldsymbol{\varrho},\boldsymbol{x}}\left(t + \Delta t\right) - \mathbf{F}\left(t\right)\tag{12}$$

where *M* , *C* and *KT* are mass, damping, and tangent stiffness matrices of SSI model, respec‐ tively; *Δu* is the incremental vector of the relative displacements for SSI system between times *t* and *t* + *Δt* ; and *F* (*t*) is the vector of internal forces at time *t*. The term *u***¨** *<sup>g</sup>*, *<sup>x</sup>*(*t* + *Δt*) is the freefield component of acceleration in *x* direction. The column matrix, *mx*, is the directional mass values of the structure only.

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 [44] 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 [43] is defined based on the cross-section deformation of structural elements as:

$$DI = \frac{\theta\_w - \theta\_y}{\theta\_u - \theta\_y} + \frac{\beta}{\theta\_u M\_y} \int dE\_H \tag{13}$$

where *θm* is the maximum rotation during loading history; *θ<sup>u</sup>* and *θy* are the ultimate and yield rotation, respectively; *My* is the yield moment; and *∫dEH* 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:

$$\text{IDI}\_{\text{stay}} = \sum\_{i=1}^{\text{Ne}} \text{DI}\_{i} \cdot \text{\mathcal{A}}\_{i, \text{component}} \; ; \quad \text{\mathcal{A}}\_{i, \text{component}} = \left\lceil \text{E}\_{i} / \sum \text{E}\_{i} \right\rceil\_{\text{component}} \tag{14}$$

$$\text{IDI}\_{\text{overall}} = \sum\_{i=1}^{\text{Vs}} \text{DI}\_{i,\text{stray}} \cdot \text{\mathcal{A}}\_{i,\text{stray}}; \quad \text{\mathcal{A}}\_{i,\text{stray}} = \left[ \text{E}\_i / \sum \text{E}\_i \right]\_{\text{stray}} \tag{15}$$

where *λ<sup>i</sup>* is energy weighting factor; and *Ei* is total absorbed energy by the component or storey *i*. *Ne* and *Ns* are the number of structural elements and stories, respectively.

#### **4.2. Finite element model of SSI system**

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

274 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

( ) ( ) <sup>1</sup> 1 exp *annual <sup>E</sup> annual annual Pf Pf Pf* - = - - ® = -F

 b

where *βannual* can be expressed as an reliability index as shown in Eq. (11), using the standard

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 [6]. 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

In the direct method, the discretization of nonlinear dynamic equations can be expressed in

where *M* , *C* and *KT* are mass, damping, and tangent stiffness matrices of SSI model, respec‐ tively; *Δu* is the incremental vector of the relative displacements for SSI system between times

field component of acceleration in *x* direction. The column matrix, *mx*, is the directional mass

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 [44] 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

*M u C u K u mu F* D + D + D =- +D - && & T , *x gx* && (*tt t* ) ( ) (12)

u

determined simultaneously by analyzing SSI system in each time step [22].

*t* and *t* + *Δt* ; and *F* (*t*) is the vector of internal forces at time *t*. The term *u***¨**

, is computed using the occurrence of earthquakes as a Poisson

(11)

*<sup>g</sup>*, *<sup>x</sup>*(*t* + *Δt*) is the free-

non-performance, *P f annual*

normal cumulative distribution function, *Φ* ( ⋅ ).

**4. Seismic responses and SSI system**

**4.1. Seismic responses of SSI system**

process [13, 22]:

FEM framework as:

values of the structure only.

the damage indices.

OpenSEES [49], as an open-source computational software framework, is used for by simula‐ tion 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 [50] is utilized for modeling the confined and unconfined concrete of cross-sections of structural elements. The constitutive parameters of this model are: *fc*=concrete peak strength in compression, *fu*=residual strength, *ε*<sup>0</sup> =strain at peak strength, and *εu*=ultimate compressive strain (Fig. 2(a)). The material of cover and core concrete used in the cross-section are modeled as unconfined and confined, respectively. The constitutive model of confined concrete developed by Saatcioglu and Razvi [51], are used. Furthermore, to determine the ultimate compressive strain of confined concrete, the relationship introduced by Paulay and Priestley [52] is utilized as:

$$
\mathfrak{s}\_{u\circ} = \mathfrak{s}\_{u\circ} + 1.4 \frac{\rho\_v \, f\_{yh} \, \mathfrak{s}\_{us}}{f\_{\alpha}} \tag{16}
$$

where *ε<sup>u</sup> <sup>c</sup>*, *ε<sup>u</sup> <sup>o</sup>* and *ε<sup>u</sup> <sup>s</sup>* are the ultimate compressive strain of confined and unconfined concrete, and the ultimate strain of longitudinal steel reinforcement in tensile stress, respectively; *<sup>ρ</sup><sup>v</sup>* , *<sup>f</sup> yh* and *<sup>f</sup> cc* 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 *J*2 plasticity model with linear hardening is utilized for modeling the constitutive behavior of the steel reinforcement. The material parameters defining *J*2 plasticity model are: *fy*=yield strength, *H*=hardening modulus and *E*=Young's modulus (Fig. 2(b)).

**Figure 1.** Direct method configuration for modelling of SSI system [22, 31]

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-

Seismic Reliability-Based Design Optimization of Reinforced Concrete Structures Including Soil-Structure… http://dx.doi.org/10.5772/59641 277

**Figure 2.** Material constitutive models; (a) Concrete, (b) Steel [22]

determine the ultimate compressive strain of confined concrete, the relationship introduced

1.4 *v yh u s*

r

where *ε<sup>u</sup> <sup>c</sup>*, *ε<sup>u</sup> <sup>o</sup>* and *ε<sup>u</sup> <sup>s</sup>* are the ultimate compressive strain of confined and unconfined concrete, and the ultimate strain of longitudinal steel reinforcement in tensile stress, respectively; *<sup>ρ</sup><sup>v</sup>* , *<sup>f</sup> yh* and *<sup>f</sup> cc* are the volumetric ratio, and the yield stress of confining steel reinforcement,

In this study, the one-dimensional *J*2 plasticity model with linear hardening is utilized for modeling the constitutive behavior of the steel reinforcement. The material parameters defining *J*2 plasticity model are: *fy*=yield strength, *H*=hardening modulus and *E*=Young's

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-

*f f*

*cc*

 e

(16)

*uc uo*

and the peak strength of confined concrete in compression, respectively.

**Figure 1.** Direct method configuration for modelling of SSI system [22, 31]

e e

= +

276 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

by Paulay and Priestley [52] is utilized as:

modulus (Fig. 2(b)).

yield-surface *J*2 plasticity model [53]. As shown in Fig. 3, this nonlinear model of soil material is described by a shear stress-strain backbone curve. The detailed description of the parameters of the shear stress-strain backbone curve can be found in [53].

**Figure 3.** Yield surfaces of multi-yield-surface *J2* plasticity model; (a) Octahedral shear stress-strain, (b) Von Mises multi-yield surfaces [22, 53].

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 [54] 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:

$$C\_u = a \, \rho \, V\_p \quad ; \quad V\_p = \sqrt{\frac{2G\left(1 - v\right)}{\rho\left(1 - 2v\right)}}\tag{17}$$

$$\mathbf{C}\_s = \mathbf{b} \,\rho \,\mathbf{V}\_s \quad ; \quad \mathbf{V}\_s = \sqrt{\frac{\mathbf{G}}{\rho}} \tag{18}$$

where *Cn* and *Cs* are the normal and shear damping, respectively; *ρ* and *v* are the mass density and Poisson ratio of soil, respectively; *a* and *b* are dimensionless parameters to be determined, and *Vp* and *Vs* are dilatational and shear wave velocity of propagation, respectively. Standard viscous boundary with the normal and shear damping is modeled using the Zero Length element. Moreover, the parameters of soil layers are consisted of *Gi* as low strain shear modulus, *Bi* as bulk modulus, and *τ<sup>i</sup>* that is shear strength, with *i=*1,2*,...,n* representing the soil layers' number. Other parameters of soil material depend on *Bi* , *Gi* , and *V <sup>p</sup>*, *<sup>i</sup>* .

The material damping matrix, *C*, of the SSI system is assembled by its corresponding damping matrices of structure and soil and considering the Rayleigh damping model. The factors of proportionality for damping matrices of structure and soil are computed based on 5% and 10% viscous damping respectively for structure and soil. The P-Δ effects are regarded in nonlinear time-history analyses. The accelerated Newton algorithm based on Krylov subspaces [55] is utilized for solving nonlinear equations of SSI system equilibrium.

#### **5. Artificial earthquakes**

For RBDO of structures, it is then necessary to utilize accelerograms of compatible character‐ istics 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 [13]:

$$a(t) = I\_m\left(t\right) \sum\_{n=1}^{NER} \left(4S\_{kT}\left(n\Delta f\right)\right) \left[1 + \delta\_s R\_N\right] \Big| \Delta f \quad \sin\left(2\pi n\Delta f t + \theta\_n\right) \tag{19}$$

where *a*(*t*) , *Im*(*t*) and *SKT* (.) are the non-stationary ground motion, the modulation function and the specific power spectral density function (PSDF), respectively. *NFR* is the number of sine functions or frequencies included, between 0 and *f* max, *δS* and *RN* are the coefficient of variation and a standard normal variable that used in ordinates of PSDF, *Δf* is frequency step, and *θ<sup>n</sup>* are random phase angles with a uniform distribution between 0 and 2*π*. In this work, the modulation function expressed in [13] is used:

Seismic Reliability-Based Design Optimization of Reinforced Concrete Structures Including Soil-Structure… http://dx.doi.org/10.5772/59641 279

$$I\_m\left(t\right) = \begin{cases} \left(t \,/\, T\_1\right)^d & 0 \le t \le T\_1 \\ 1 & T\_1 \le t \le T\_2 \\ e^{-c\left(t - T\_2\right)} & T\_2 \le t \le T \end{cases} \tag{20}$$

where *T*1, *T*2 and *T* are specific times and the duration of the simulated record, *d* and *c* are constants. Also, the PSDF of the non-stationary ground motion suggested by Clough and Penzien [56] is considered as:

; *ss s <sup>G</sup> C bV V* = =

278 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

where *Cn* and *Cs* are the normal and shear damping, respectively; *ρ* and *v* are the mass density and Poisson ratio of soil, respectively; *a* and *b* are dimensionless parameters to be determined, and *Vp* and *Vs* are dilatational and shear wave velocity of propagation, respectively. Standard viscous boundary with the normal and shear damping is modeled using the Zero Length element. Moreover, the parameters of soil layers are consisted of *Gi* as low strain shear

The material damping matrix, *C*, of the SSI system is assembled by its corresponding damping matrices of structure and soil and considering the Rayleigh damping model. The factors of proportionality for damping matrices of structure and soil are computed based on 5% and 10% viscous damping respectively for structure and soil. The P-Δ effects are regarded in nonlinear time-history analyses. The accelerated Newton algorithm based on Krylov subspaces [55] is

For RBDO of structures, it is then necessary to utilize accelerograms of compatible character‐ istics 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 [13]:

*a t I t S n f R f n ft*

() () ( ) { } ( ) 1 2

<sup>=</sup> D + D D+ é ù å ë û d

4 1 sin 2

*m KT S N n*

where *a*(*t*) , *Im*(*t*) and *SKT* (.) are the non-stationary ground motion, the modulation function and the specific power spectral density function (PSDF), respectively. *NFR* is the number of sine functions or frequencies included, between 0 and *f* max, *δS* and *RN* are the coefficient of variation and a standard normal variable that used in ordinates of PSDF, *Δf* is frequency step, and *θ<sup>n</sup>* are random phase angles with a uniform distribution between 0 and 2*π*. In this work,

pq

(19)

r

that is shear strength, with *i=*1,2*,...,n* representing the soil

, and *V <sup>p</sup>*, *<sup>i</sup>*

.

(18)

r

layers' number. Other parameters of soil material depend on *Bi* , *Gi*

utilized for solving nonlinear equations of SSI system equilibrium.

1

=

*n*

the modulation function expressed in [13] is used:

*NFR*

modulus, *Bi*

as bulk modulus, and *τ<sup>i</sup>*

**5. Artificial earthquakes**

$$S\_{\rm{AT}}\left(f\right) = S\_0 \left[\frac{1 + 4\,\xi\_g^2 \left(f \,\, \left/f\_g\right)^2}{\left(1 - \left(f \,\, \left/f\_g\right)^2\right)^2 + 4\,\xi\_g^2 \left(f \,\, \left/f\_g\right)^2\right)}\right] \left[\frac{\left(f \,\, \left/f\_f\right)^4}{\left(1 - \left(f \,\, \, \, \, \, \, f\right)^2\right)^2 + 4\,\xi\_f^2 \left(f \,\, \, \, f\_f\right)^2}\right] \tag{21}$$

where *S*<sup>0</sup> is the constant PSDF of input white-noise random process; *f <sup>g</sup>* and *ξg* are the charac‐ teristic ground frequency and the ground damping ratio; *f <sup>f</sup>* and *ξ<sup>f</sup>* are parameters for a highpass 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. [13].


**Table 2.** Parameters for the generation of simulated ground motion (\* *aG* = max(*a*(*t*)))

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. [13] are used. An artificial earthquake generated based on Eq. (19) is shown in Fig. 4.


**Table 3.** The PGA values of seismic hazard levels

**Figure 4.** An artificial earthquake with PGA=0.8g [22]

#### **6. The new discrete gravitational search algorithm**

Based on the work presented by Khatibinia *et al*. [14], a new discrete gravitational search algorithm (DGSA) based on the standard gravitational search algorithm (GSA) is utilized to find the optimum designs in the RBDO procedure.

#### **6.1. Gravitational search algorithm**

GSA was introduced by Rashedi *et al*. [23] as a new stochastic population based search algorithm based on the law of gravity and mass interactions. In GSA, each agent of the population represents a potential solution of the optimization problem. The *i*th agent in *t*th iteration is associated with a position vector, *X<sup>i</sup>* (*t*)={ *x <sup>i</sup>* <sup>1</sup> , ... , *<sup>x</sup> <sup>i</sup> <sup>d</sup>* , ... , *<sup>x</sup> <sup>i</sup> <sup>D</sup>*}, and a velocity vector, *V<sup>i</sup>* (*t*)={ *v <sup>i</sup>* 1 , ... , *v <sup>i</sup> <sup>d</sup>* , ... , *<sup>v</sup> <sup>i</sup> <sup>D</sup>*}. *D* is dimension of the solution space. Based on [23], the mass of each agent is calculated after computing the current population fitness for a minimi‐ zation problem as follows:

Seismic Reliability-Based Design Optimization of Reinforced Concrete Structures Including Soil-Structure… http://dx.doi.org/10.5772/59641 281

$$M\_i\left(t\right) = \frac{fit\_i\left(t\right) - worst\left(t\right)}{\sum\_{j=1}^{N} \left(fit\_j\left(t\right) - worst\left(t\right)\right)}\tag{22}$$

where *N*, *Mi* (*t*) and *fiti* (*t*) represent the population size, the mass and the fitness value of agent *i* at *t*th iteration, respectively; and *worst*(*t*) is defined as the worst fitness of all agents.

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):

$$a\_i^d(t) = \sum\_{j \text{ackbest}, \, j \neq i} r \, and\_j G(t) \frac{M\_j(t)}{R\_{i/j}(t) + \varepsilon} \left( x\_j^d(t) - x\_i^d(t) \right) \tag{23}$$

$$v v\_i^d \left(t + 1\right) = r \; \text{and} \; v\_i^d \left(t\right) + a\_i^d \left(t\right) \tag{24}$$

$$
\propto \mathbf{x}\_i^d \left( t + \mathbf{1} \right) = \mathbf{x}\_i^d \left( t \right) + \upsilon\_i^d \left( t + \mathbf{1} \right) \tag{25}
$$

where *ai <sup>d</sup>* , *vi <sup>d</sup>* and *xi <sup>d</sup>* present the acceleration, velocity and position of *i*th agent in dimension *d*, respectively. *r andi* and *r andj* are two uniformly distributed random numbers in the interval [0, 1]; *ε* is a small value; and *Ri*, *<sup>j</sup>* (*t*) is the Euclidean distance between agent *i* and *j*. *kbest* is the set of first *k* agents with the best fitness value and biggest mass, which is a function of time, initialized to *k*<sup>0</sup> at the beginning and decreased with time. Here, *k*<sup>0</sup> is set to *N* and is decreased linearly to 1.0. *G* is a decreasing function of time.

#### **6.2. The proposed discrete GSA**

**Probability of exceedance Recurrence interval PGA (g)** 50% in 50 years 73 0.27 10% in 50 years 475 0.6 5% in 50 years 975 0.8

280 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

0 2 4 6 8 10 12 14 16 18 20 Time (sec)

Based on the work presented by Khatibinia *et al*. [14], a new discrete gravitational search algorithm (DGSA) based on the standard gravitational search algorithm (GSA) is utilized to

GSA was introduced by Rashedi *et al*. [23] as a new stochastic population based search algorithm based on the law of gravity and mass interactions. In GSA, each agent of the population represents a potential solution of the optimization problem. The *i*th agent in *t*th

mass of each agent is calculated after computing the current population fitness for a minimi‐

(*t*)={ *x <sup>i</sup>*

<sup>1</sup> , ... , *<sup>x</sup> <sup>i</sup>*

*<sup>D</sup>*}. *D* is dimension of the solution space. Based on [23], the

*<sup>d</sup>* , ... , *<sup>x</sup> <sup>i</sup>*

*<sup>D</sup>*}, and a velocity

**Table 3.** The PGA values of seismic hazard levels

**Figure 4.** An artificial earthquake with PGA=0.8g [22]

**6. The new discrete gravitational search algorithm**

find the optimum designs in the RBDO procedure.

iteration is associated with a position vector, *X<sup>i</sup>*

*<sup>d</sup>* , ... , *<sup>v</sup> <sup>i</sup>*

**6.1. Gravitational search algorithm**


vector, *V<sup>i</sup>*

(*t*)={ *v <sup>i</sup>* 1 , ... , *v <sup>i</sup>*

zation problem as follows:

Ground acceleration (m/sec^2)

The binary GSA (BGSA) for solving discrete problem was developed by Rashedi *et al*. [25]. In the BGSA model, the positions of agents are indicated by one or zero. Hence, Eq. (24) can be used without any changes for updating the velocity of agents. Because of terms of probability, the velocity must be converted into interval [0, 1] by a transfer function, *S*(*v <sup>i</sup> <sup>d</sup>* (*t*)) [25]. For updating the position of the agents, Eq. (25) is re-defined as follows [25]:

$$\mathbf{x}\_{i}^{d}\left(t+1\right) = \begin{cases} \operatorname{complement}\left(\mathbf{x}\_{i}^{d}\left(t\right)\right) & \text{for} \quad r \text{ and} < \mathbf{S}\left(\mathbf{v}\_{i}^{d}\left(t\right)\right) \\\qquad \mathbf{x}\_{i}^{d}\left(t\right) & \text{for} \quad r \text{ and} \ge \mathbf{S}\left(\mathbf{v}\_{i}^{d}\left(t\right)\right) \end{cases} \tag{26}$$

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 [24]. Therefore, the search performance of DGSA can be improved using the passive congregation. To achieve this purpose, Khatibinia *et al*. [14] modified the velocity of agents based on the particle swarm optimizer with passive congregation (PSOPC) which is proposed by He *at al*. [24]. The modified velocity is expressed as follows:

$$\begin{aligned} \boldsymbol{\upsilon}\_{i}^{d}\left(\boldsymbol{t}+\mathbf{1}\right) &= \boldsymbol{r}\_{i,1}\boldsymbol{\upsilon}\_{i}^{d}\left(\boldsymbol{t}\right) + \boldsymbol{a}\_{i}^{d}\left(\boldsymbol{t}\right) + \boldsymbol{r}\_{i,2}\left(\boldsymbol{p}\_{\text{best},i}^{d}\left(\boldsymbol{t}\right) - \boldsymbol{\mathsf{x}}\_{i}^{d}\left(\boldsymbol{t}\right)\right) + \\ &\boldsymbol{r}\_{i,3}\left(\boldsymbol{p}\_{\text{s}}^{d}\left(\boldsymbol{t}\right) - \boldsymbol{\mathsf{x}}\_{i}^{d}\left(\boldsymbol{t}\right)\right) + \boldsymbol{r}\_{i,4}\left(\boldsymbol{p}\_{i}^{d}\left(\boldsymbol{t}\right) - \boldsymbol{\mathsf{x}}\_{i}^{d}\left(\boldsymbol{t}\right)\right) \end{aligned} \tag{27}$$

where *pbest <sup>d</sup>* , *pg <sup>d</sup>* and *pi <sup>d</sup>* are the best previous position of the *i*th agent, the best previous position among all the agents and a randomly selected agent from the population, respectively. *ri*,1 and *ri*,2 are the random number in the interval [0, 1], respectively; *ri*,3 is the uniform random number in the interval (0, 1).

In DGSA, the scalar *x <sup>i</sup> <sup>d</sup>* ∈ {1, 2, ... , *n*} corresponds to discrete values of the set *A*={*A1*, *A2*,...,*An*}. Hence, the position of agents is updated by the following equation instead of Eq. (26):

$$\mathbf{x}\_{i}^{d}\left(t+1\right) = \text{INT}\left(\mathbf{x}\_{i}^{d}\left(t\right) + \upsilon\_{i}^{d}\left(t\right)\right) \tag{28}$$

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 *A*. To avoid this problem, the current position of particles is limited as:

$$\mathbf{x}\_{i}^{d}(t+1) = \begin{cases} P\_{\boldsymbol{\mathcal{S}}}^{d} & \text{if } \operatorname{ ft}\left(p\_{\boldsymbol{\mathcal{S}}}^{d}(t)\right) = \operatorname{ft}\left(p\_{\boldsymbol{\mathcal{S}}}^{d}(t-1)\right) \\\ P\_{i,\text{best}}^{d}(t) & \text{if } \operatorname{ ft}\left(p\_{\boldsymbol{\mathcal{S}}}^{d}(t)\right) \neq \operatorname{ft}\left(p\_{\boldsymbol{\mathcal{S}}}^{d}(t-1)\right) \\\ \operatorname{INT}\left(\mathbf{x}\_{i}^{d}(t) + \boldsymbol{v}\_{i}^{d}(t)\right) & \text{otherwise} \end{cases} \tag{29}$$

where *INT* ( ⋅ ) 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 [11]. Because of the drawback, Khatibinia *et al*. [22] proposed a metamodel that predicted the mean and the standard deviation of the structural seismic responses in the RBDO process based on MCS. The meta-model consists of weighted least squares support vector machine (WLS-SVM) and wavelet kernel function, which is called WWLS-SVM [14, 22, 27].

#### **7.1. Weighted least squares support vector machines**

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 [24]. Therefore, the search performance of DGSA can be improved using the passive congregation. To achieve this purpose, Khatibinia *et al*. [14] modified the velocity of agents based on the particle swarm optimizer with passive congregation (PSOPC) which is proposed

282 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

( ) () () () () ( )

Hence, the position of agents is updated by the following equation instead of Eq. (26):

( ) () () 1 ( ) *<sup>d</sup> d d*

( () ())

*INT x t v t*

*d d d d i i best g g d d i i*

*x t P t fit p t fit p t*

,

<sup>ï</sup> <sup>+</sup> <sup>î</sup>

**7. Approximating the structural seismic responses**

ï

where *INT* ( ⋅ ) denotes the integral part function.

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 *A*. To avoid this problem, the current position of particles is limited as:

> *d d d g g g*

MCS requires excessive computational cost for RBDO of structures in order to obtain an acceptable accuracy [11]. Because of the drawback, Khatibinia *et al*. [22] proposed a meta-

<sup>ì</sup> = - <sup>ï</sup> <sup>ï</sup> + = <sup>í</sup> ¹ -

*P fit p t fit p t*

1 ( ) if 1

,1 ,2 ,

*d dd d d i i i i i best i i*

*vt rvt at r p t xt*

,3 ,4

( () ()) ( ( ))


*d d dd ig i ii i*

among all the agents and a randomly selected agent from the population, respectively. *ri*,1 and *ri*,2 are the random number in the interval [0, 1], respectively; *ri*,3 is the uniform random number

*r p t x t r pt x t* += + + - +

( )

*<sup>d</sup>* are the best previous position of the *i*th agent, the best previous position

*<sup>d</sup>* ∈ {1, 2, ... , *n*} corresponds to discrete values of the set *A*={*A1*, *A2*,...,*An*}.

*<sup>i</sup> i i x t INT x t v t* += + (28)

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

if 1

otherwise

(27)

(29)

by He *at al*. [24]. The modified velocity is expressed as follows:

1

where *pbest*

*<sup>d</sup>* , *pg*

in the interval (0, 1).

In DGSA, the scalar *x <sup>i</sup>*

*<sup>d</sup>* and *pi*

( )

WLS-SVM was introduced as excellent machine learning algorithms in large-scale problems by Suykens *et al*. [26]. In fact, assigning weights to SVM as well as to the least squares version of SVM (LS-SVM) is resulted in more robust and precise prediction of functions [26].

WLS-SVM is described as the following optimization problem in primal weight space [26]:

$$\min J(\boldsymbol{w}, \boldsymbol{e}) = \frac{1}{2} \left\| \boldsymbol{w} \right\|^2 + \frac{1}{2} \gamma \sum\_{i=1}^n \overline{v}\_i e\_i^2 \tag{30}$$

Subject to the following equality constraints:

$$\mathbf{x}\_{i} = \mathbf{w}^{T} \boldsymbol{\phi} \mathbf{(x}\_{i}) + \mathbf{b} + \mathbf{e}\_{i}, \qquad \mathbf{i} = \mathbf{1}, \mathbf{2}, \dots, \mathbf{n} \tag{31}$$

where {*x<sup>i</sup>* , *yi* }*i*=1 *<sup>n</sup>* is a training data set; *xi* <sup>∈</sup>*<sup>R</sup> <sup>n</sup>* and *yi* <sup>∈</sup>*R* represent the input and output data. Operator *φ*(.): *R <sup>n</sup>* →*R <sup>d</sup>* is a function which maps the input space into a higher dimensional space. The vector, *<sup>w</sup>* <sup>∈</sup>*<sup>R</sup> <sup>d</sup>* , represents weight vector in primal weight space. The symbols, *ei* <sup>∈</sup>*<sup>R</sup>* and *b*∈*R*, 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:

$$
\partial \rho(\mathbf{x}) = \mathbf{w}^T \rho(\mathbf{x}) + b \tag{32}
$$

It is impossible to indirectly compute *w* from Eq. (30), for the structure of the function *φ*(*x*) is generally unknown. Therefore, the dual problem shown in Eq. (30) is minimized by the Lagrange multiplier method as follows:

$$L\left(w, b, e; \mathbf{x}\right) = f\left(w, e\right) - \sum\_{i=1}^{n} \alpha\_i \left(w^T \boldsymbol{\phi}\left(\mathbf{x}\_i\right) + b + e\_i - y\_i\right) \tag{33}$$

Based on the Karush-Kuhn-Tucker (KKT) conditions, by eliminating *w* and *e* the solution is given by the following set of linear equations:

$$
\begin{bmatrix}
\boldsymbol{\Omega} + \mathbf{V}\_{\boldsymbol{\gamma}} & \mathbf{1}\_{\boldsymbol{u}}^{T} \\
\mathbf{1}\_{\boldsymbol{u}} & \mathbf{0}
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{\alpha} \\
\boldsymbol{b}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{y} \\
\boldsymbol{0}
\end{bmatrix} \tag{34}
$$

where

$$\begin{aligned} \mathbf{V}\_{\boldsymbol{\gamma}} &= \text{diag}\left\{ \mathbf{1}/\boldsymbol{\gamma}\,\bar{\boldsymbol{v}}\_{1\boldsymbol{\gamma}} \,\dots \,\,\mathbf{1}/\boldsymbol{\gamma}\,\bar{\boldsymbol{v}}\_{n} \right\} \; ; \; \boldsymbol{\Omega}\_{i,j} = \{ \boldsymbol{\varphi}(\boldsymbol{x}\_{i}), \,\boldsymbol{\varphi}(\boldsymbol{x}\_{j}) \}\_{\mathcal{H}} \; \; i, \; j = \mathbf{1}, \; \ldots, n \\ \mathbf{y} &= \begin{bmatrix} \mathbf{y}\_{1\boldsymbol{\gamma}} \,\dots \,\,\boldsymbol{v}\_{n} \end{bmatrix} \; ; \; \mathbf{1}\_{n}^{T} &= \begin{bmatrix} \mathbf{1} \; \wedge \; \ldots \; \,\mathbf{1} \end{bmatrix} \; ; \; \mathbf{a} = \begin{bmatrix} a\_{1\boldsymbol{\gamma}} \; \ldots \; \, a\_{n} \end{bmatrix} \end{aligned} \tag{35}$$

According to Mercer's condition, a kernel *K*(. , .) is selected, such that:

$$K\left(\mathbf{x}, \overline{\mathbf{x}}\right) = \left<\boldsymbol{\varphi}(\mathbf{x}), \boldsymbol{\varphi}(\overline{\mathbf{x}})\right>\_{H} \tag{36}$$

Consequently, the final WLS-SVM model for the prediction of functions becomes:

$$\log\left(\mathbf{x}\right) = \sum\_{i=1}^{n} \alpha\_{i} K\left(\mathbf{x}\_{i}, \mathbf{x}\right) + b \tag{37}$$

Weight *v*¯ *<sup>k</sup>* is estimated as follows [57-58]:

$$\overline{\upsilon}\_{k} = \begin{cases} 1 & \text{if } \left| \boldsymbol{e}\_{i} / \hat{\boldsymbol{s}} \right| \leq \boldsymbol{c}\_{1} \\ \boldsymbol{c}\_{2} - \left| \boldsymbol{e}\_{i} / \hat{\boldsymbol{s}} \right| & \text{if } \left| \boldsymbol{c}\_{1} < \left| \boldsymbol{e}\_{i} / \hat{\boldsymbol{s}} \right| \leq \boldsymbol{c}\_{2} \\ & \text{10}^{-4} & \text{otherwise} \end{cases} \tag{38}$$

where *s* ^ is a robust estimation of the standard deviation for the error variables (*ei* <sup>=</sup>*ai* / *Dii* −1 ); constants *c*1 and *c*2 are typically chosen as *c*<sup>1</sup> =2.5 and *c*<sup>2</sup> =3. Here *Dii* <sup>−</sup><sup>1</sup> denotes the *i*th primal diagonal element of inverse of matrix *D*, which is the matrix on the left-hand of the system of the linear Eq. (34) [59]. After that, weights *v*¯ *<sup>k</sup>* are determined, and the model (Eq. (31)) is obtained following the solving WLS-SVM problem (Eq. (34)).

In WLS-SVM, Gaussian radial basis function (RBF) is frequently used as the kernel function, and it is expressed as:

$$K(\mathbf{x}, \overline{\mathbf{x}}) = \exp\left(-\frac{\left\|\mathbf{x} - \overline{\mathbf{x}}\right\|^2}{\sigma^2}\right) \tag{39}$$

where *σ* <sup>2</sup> is a positive real constant usually called the kernel width.

Based upon the Suykens *et al*. [26], the WLS-SVM model based on RBF kernel function for predicting the output data is implemented using the following procedure:

**Step 1.** Assign training data {*x<sup>k</sup>* , *yk* }*<sup>k</sup>* =1 *Ntot*, set *N=Ntot*.

**Step 2.** Find an optimum (*γ*, *σ*) combination on the all range of *Ntot* training data by 10-fold cross-validation, then solve linear system (Eq. (34)), and give the model (Eq. (37)).

**Step 3.** Sort the values | *α* |.

where

where *s*

where *σ* <sup>2</sup>

and it is expressed as:

*V<sup>γ</sup>* =*diag*{1 / *γ v*¯ <sup>1</sup>, ..., 1 / *γ v*¯ *<sup>n</sup>*} ; *Ωi*, *<sup>j</sup>* = *φ*(*x<sup>i</sup>*

According to Mercer's condition, a kernel *K*(. , .) is selected, such that:

( ) () , , ( ) *<sup>H</sup> <sup>K</sup> xx x x* <sup>=</sup> j

Consequently, the final WLS-SVM model for the prediction of functions becomes:

284 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

*n*

*i y Kb* = *x xx* = + åa

2

ï î

obtained following the solving WLS-SVM problem (Eq. (34)).

2 1 4

*i k i*

*c es*

*c c*

constants *c*1 and *c*2 are typically chosen as *c*<sup>1</sup> =2.5 and *c*<sup>2</sup> =3. Here *Dii*

( )

*x x*

<sup>2</sup> *K* , exp

is a positive real constant usually called the kernel width.

predicting the output data is implemented using the following procedure:

10

() ( ) 1

1 / ˆ

*v if c e s c*

<sup>ï</sup> - <sup>=</sup> <sup>í</sup> < £ - <sup>ï</sup>

<sup>ì</sup> £ <sup>ï</sup>

/ <sup>ˆ</sup> / <sup>ˆ</sup>

*otherwise* -

diagonal element of inverse of matrix *D*, which is the matrix on the left-hand of the system of the linear Eq. (34) [59]. After that, weights *v*¯ *<sup>k</sup>* are determined, and the model (Eq. (31)) is

In WLS-SVM, Gaussian radial basis function (RBF) is frequently used as the kernel function,

æ ö - = -ç ÷ ç ÷ è ø

Based upon the Suykens *et al*. [26], the WLS-SVM model based on RBF kernel function for

*x x*

s

^ is a robust estimation of the standard deviation for the error variables (*ei* <sup>=</sup>*ai* / *Dii*

*i*

*if e s c*

*i i*

 j

,

1

1 2

2

*y* = *y*1, ..., *yn <sup>T</sup>* ; 1*<sup>n</sup>*

Weight *v*¯ *<sup>k</sup>* is estimated as follows [57-58]:

), *φ*(*x <sup>j</sup>*

*<sup>T</sup>* <sup>=</sup> 1, . . . , 1 ; *<sup>α</sup>* <sup>=</sup> *<sup>α</sup>*1, ..., *<sup>α</sup><sup>n</sup>*

) *<sup>H</sup> i* , *j* =1, . . . , *n*

(36)

(37)

(35)

(38)

−1 );

(39)

<sup>−</sup><sup>1</sup> denotes the *i*th primal

**Step 4.** Remove a small number of *M* points (typically 5% of the *N* points) which has the smallest values in the sorted | *α* |.

**Step 5.** Retain *N-M* points and set *N=N-M*.

**Step 6.** Go to 2 and retrain on the reduced training set.

#### **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. [22]. 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:

$$\Psi(t) = \frac{1}{\sqrt{a}} \cos \left( a \partial\_0 \left( \frac{t-b}{a} \right) \right) \exp \left( -0.5 \left( \frac{t-b}{a} \right)^2 \right) \tag{40}$$

where *a* and *b* are the scale factor and the translation factor, respectively.

According to Zhang *et al*. [64], the translation-invariant wavelet kernels can be explained as follows:

$$K\left(\mathbf{x}, \overline{\mathbf{x}}\right) = \prod\_{i=1}^{n} \Psi\left(\frac{\mathbf{x}\_{i} - \overline{\mathbf{x}}\_{i}}{a}\right) \tag{41}$$

where *n* is the number of samples; *x* and *x***¯** ∈*R <sup>n</sup>*

Therefore, according to Eqs. (40) and (41), the wavelet kernel function of the cosine-Gaussian Morlet wavelet is given as follows:

$$K(\mathbf{x}, \overline{\mathbf{x}}) = \prod\_{i=1}^{n} \frac{1}{\sqrt{a}} \cos \left( a o\_0 \left( \frac{\mathbf{x}\_i - \overline{\mathbf{x}}\_i}{a} \right) \right) \exp \left( -0.5 \frac{\left\| \mathbf{x}\_i - \overline{\mathbf{x}}\_i \right\|^2}{a^2} \right) \tag{42}$$

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, *a* and *ω*0. To achieve this purpose, a mean absolute percentage error (*MAPE*) is used to evaluate the performance of the WWLS-SVM model as follows:

$$MAPE = \frac{1}{n} \sum\_{i=1}^{n} 100 \times \left| \frac{y\_i - \overline{y}\_i}{y\_i} \right| \tag{43}$$

where *y* and *y*¯ are the actual value and the predicted value, respectively.

The WWLS-SVM training stage during GSA is performed according to the *k*-fold crossvalidation (CV) [26]. Consequently, an optimization problem based on *MAPE* of the *k*-fold CV (*MAPEk-fold CV*) is expressed as:

$$\begin{array}{ll}\text{Minimize} & \text{MAPE}\_{k-foldCV} \\ \text{Subject to} & \gamma\_{\text{min}} \le \gamma \le \gamma\_{\text{max}} \\ & a\_{\text{min}} \le a \le a\_{\text{max}} \\ & a\_{0\text{min}} \le a\_0 \le a\_{0\text{max}} \end{array} \tag{44}$$

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 *et al*. (2009). The flowchart of the meta-model based on WWLS-SVM [22] and GSA [14, 23] is shown in Fig. 5.

**Figure 5.** The prediction meta-model based on WWLS-SVM and GSA [14]

#### **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 [65]. 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 [65]. This efficient method is consisted of a modified adaptive neuro fuzzy inference system (ANFIS) with a hybrid of fuzzy c-means (FCM) [66] and fuzzy particle swarm optimization (FPSO) [67]. 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**

1 <sup>1</sup> <sup>100</sup> *n*

where *y* and *y*¯ are the actual value and the predicted value, respectively.

Minimize Subjectto

**Figure 5.** The prediction meta-model based on WWLS-SVM and GSA [14]

**8. Predicting failure probability of structures**

(*MAPEk-fold CV*) is expressed as:

*y y MAPE n y* <sup>=</sup>

286 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

*i i*

The WWLS-SVM training stage during GSA is performed according to the *k*-fold crossvalidation (CV) [26]. Consequently, an optimization problem based on *MAPE* of the *k*-fold CV

> min max min max 0 min 0 0 max

 gg

 ww

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 *et al*. (2009). The flowchart of the meta-model based on WWLS-SVM [22] and GSA [14, 23] is shown in Fig. 5.

In the RBDO procedure, nonlinear time-history analysis of SSI system is used and it may be failed regarding a number of random structures [65]. 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 [65]. This efficient method is consisted of a modified adaptive neuro fuzzy inference system (ANFIS) with a hybrid of fuzzy c-means (FCM) [66] and fuzzy particle swarm optimization (FPSO) [67].

*a aa*

g

w

*MAPEk foldCV*


*i i*


(44)

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 [65].

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 *XFCM* . This optimum solution is copied *NFCM* times to create the some part of the initial swarm FPSO. Other particles of the initial swarm, i.e. X*rnd* , *<sup>j</sup>* ( *j* =1, 2, ..., *NPSO* −NFCM), are selected randomly to complete the initial swarm. Then, the FPSO algorithm is used using this initial swarm. The algorithm flow of the FCM–FPSO strategy is shown in Fig. 6.

**Figure 6.** The algorithm flow of the FCM–FPSO for clustering [65]

### **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 *et al*. [65] proposed a modified ANFIS to predict the probability of failure. In this model, the number of clusters, the cluster centers and membership grades are considered as parameters which optimized by subtractive algorithm (SA) and the hybrid of FCM-FPSO and used in FIS for tuning ANFIS. The algorithm flow of the proposed model is shown in Fig. 7. The proposed method is executed in the following steps [65]:

**Step 1.** SA finds the optimum number of the clusters (*nc)*.

**Step 2.** The hybrid FCM-FPSO algorithm partitions training data to *nc* clusters and determines membership grades each of clusters. This parameters is used for optimizing the center of rules and membership functions for the input and output data.

**Step 3.** The FIS structure with a minimum number of fuzzy rules and membership functions is generated by using SA and the hybrid FCM-FPSO algorithm. The FIS uses Gaussian function and linear function for membership function of input and output, respectively. These param‐ eters are tuned for the ANFIS.

**Step 4.** The ANFIS is employed for training data. The ANFIS uses a hybrid learning algorithm to identify parameters of Sugeno-type fuzzy inference systems. a combination of the leastsquares method and the back-propagation gradient descent method are applied for training FIS membership functions.

**Figure 7.** The algorithm flow of the proposed method [65]

## **9. Numerical examples**

**8.2. Modified ANFIS**

method is executed in the following steps [65]:

eters are tuned for the ANFIS.

FIS membership functions.

**Figure 7.** The algorithm flow of the proposed method [65]

**Step 1.** SA finds the optimum number of the clusters (*nc)*.

and membership functions for the input and output data.

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 *et al*. [65] proposed a modified ANFIS to predict the probability of failure. In this model, the number of clusters, the cluster centers and membership grades are considered as parameters which optimized by subtractive algorithm (SA) and the hybrid of FCM-FPSO and used in FIS for tuning ANFIS. The algorithm flow of the proposed model is shown in Fig. 7. The proposed

288 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

**Step 2.** The hybrid FCM-FPSO algorithm partitions training data to *nc* clusters and determines membership grades each of clusters. This parameters is used for optimizing the center of rules

**Step 3.** The FIS structure with a minimum number of fuzzy rules and membership functions is generated by using SA and the hybrid FCM-FPSO algorithm. The FIS uses Gaussian function and linear function for membership function of input and output, respectively. These param‐

**Step 4.** The ANFIS is employed for training data. The ANFIS uses a hybrid learning algorithm to identify parameters of Sugeno-type fuzzy inference systems. a combination of the leastsquares method and the back-propagation gradient descent method are applied for training

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 *DL*=5.884 N/mm2 (600 kg/m2 ) and *LL*=1.961 N/mm2 (200 kg/m2 ), respectively.

**Figure 8.** Two illustrated RC frame structures [22]

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 *wC* =60 and *wS* =700 Euros. To calculate the total expected cost of repair, first, the cumulative distribution is obtained using MCS and the proposed metamodel for the response *DIoverall* adjusting a Beta distribution. Then, the density function is assigned by the derivative of the cumulative distribution. The target values of reliability indices corresponding to the three performance levels (Table 1) considered in RBDO of RC structures are equal to 1.276, 2.326 and 2.697, respectively [52]. For RBDO of RC structures, the probability density function (PDF), the mean and standard deviation (SD) values for each random constitutive parameter are listed in Table 5.

The concrete material parameters shown in Table 5 are considered for the cover of column cross-sections. The strain corresponding to the peak strength, *ε*<sup>0</sup> *<sup>o</sup>*, and the residual strength, *fu*o, for unconfined concrete are selected as 0.002 and 0.0, respectively. Shear-wave velocity, *Vs,* and friction angle, *φ*, are considered as random parameters of the soil layers. The other parameters of soil layers depend on their shear-wave velocities. Thus, in the process, first, the shear-wave velocities of soil layers are randomly selected; then, the other parameters are computed based on the shear-wave velocity. The PGA value, *ag*, and the central frequency, *fg*, for soil filter are considered as random variables in the RBDO process. The PGA value is also taken into account as follows [13, 22]:

$$
\overline{a}\_{\mathfrak{g}} = \overline{a}\_{\mathfrak{g}} \left( \mathbf{1} + \sigma\_{\overline{\mathfrak{g}}\_{\mathfrak{g}}} \right) \tag{45}
$$


where *a*¯ *<sup>g</sup>* and *σa*¯ *<sup>g</sup>* are the mean and the standard deviation of PGA, respectively. The values of these parameters are shown in Table 6.

**Table 4.** The section database of columns and beams. (a. The unit of sections is cm)

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, *G*(*t*), depends on the actual age of the universe. In this study, *G*(*t*) is considered as a linear decreasing function [23] in DGSA. Since the initial value of the gravitational constant is found to affect the optimization results significantly, the fixed value of *G0* = 50 is utilized in this study. In order to consider the stochastic nature of the optimization process, ten independent optimization runs are per‐ formed and the best solution is considered as the final results.

#### **9.1. Example 1: Six-storey RC frame**

model for the response *DIoverall* adjusting a Beta distribution. Then, the density function is assigned by the derivative of the cumulative distribution. The target values of reliability indices corresponding to the three performance levels (Table 1) considered in RBDO of RC structures are equal to 1.276, 2.326 and 2.697, respectively [52]. For RBDO of RC structures, the probability density function (PDF), the mean and standard deviation (SD) values for each random

290 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

The concrete material parameters shown in Table 5 are considered for the cover of column cross-sections. The strain corresponding to the peak strength, *ε*<sup>0</sup> *<sup>o</sup>*, and the residual strength, *fu*o, for unconfined concrete are selected as 0.002 and 0.0, respectively. Shear-wave velocity, *Vs,* and friction angle, *φ*, are considered as random parameters of the soil layers. The other parameters of soil layers depend on their shear-wave velocities. Thus, in the process, first, the shear-wave velocities of soil layers are randomly selected; then, the other parameters are computed based on the shear-wave velocity. The PGA value, *ag*, and the central frequency, *fg*, for soil filter are considered as random variables in the RBDO process. The PGA value is also

> ( ) 1 *g g ag a a* = +s

where *a*¯ *<sup>g</sup>* and *σa*¯ *<sup>g</sup>* are the mean and the standard deviation of PGA, respectively. The values of

1 65 × 65 55 × 55 45 × 45 55 × 45 50 × 45 45 × 40 2 60 × 60 50 × 50 40 × 40 55 × 40 50 × 40 45 × 35 3 55 × 55 50 × 50 40 × 40 55 × 35 50 × 35 45 × 30 4 55 × 55 45 × 45 40 × 40 55 × 30 50 × 30 45 × 25 5 55 × 55 45 × 45 35 × 35 50 × 40 45 × 40 40 × 35 6 50 × 50 45 × 45 40 × 40 50 × 35 45 × 35 40 × 30 7 55 × 55 45 × 45 35 × 35 50 × 30 45 × 30 40 × 30 8 - - - 50 × 30 45 × 30 35 × 30

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, *G*(*t*), depends on the actual

**Table 4.** The section database of columns and beams. (a. The unit of sections is cm)

**Column (h × h)a Beam (h × b)a Type 1 Type 2 Type 3 Type 1 Type 2 Type 3**

(45)

constitutive parameter are listed in Table 5.

taken into account as follows [13, 22]:

these parameters are shown in Table 6.

**Number**

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).


**Table 5.** The marginal probability distribution, mean and standard deviation of materials


**Table 6.** Properties of the random variables for generation of the artificial earthquakes

#### *9.1.1. Training and testing the meta-model*

To predict the mean, *R*¯ *i* , and the standard deviation, *σR i*, of *i*th seismic response during RBDO of the frame, the proposed meta-models are trained based on the generated database. The WWLS-SVM training during GSA is performed according to five-fold cross-validation. The lower and upper bonds of the parameters required in the optimization process are selected as *γ* ∈ 1.0, 500 , *a*∈ 0.5, 5.0 and *ω*0∈ 1.0, 10 , respectively. Therefore, the training optimal parameters of the meta-model associated with the mean and the standard deviation of seismic responses are shown in Table 7.

In order to validate the performance and accuracy of the proposed meta-model, relative rootmean-squared error, i.e. *RRMSE*, and *R2* as the absolute fraction of variance, during testing the meta-model and WLS-SVM, are defined using the following equations:

$$RRMSE = \sqrt{\frac{n\sum\_{i=1}^{n} \left(y\_i - \overline{y}\_i\right)^2}{\left(n-1\right)\sum\_{i=1}^{n} y\_i^2}}\tag{46}$$

$$\mathcal{R}^2 = 1 - \left(\frac{\sum\_{i=1}^n (y\_i - \overline{y}\_i)^2}{\sum\_{i=1}^n \overline{y}\_i^2}\right) \tag{47}$$


**Table 7.** Optimal parameters of the meta-model for training the mean and the standard deviation of the seismic responses

The smaller *RRMSE* and *MAPE* and the larger *R2* are indicative of better performance gener‐ ality. The comparison of WWLS-SVM and WLS-SVM, with respect to *MAPE*, *RRMSE*, and *R2* is shown in Table 8.


**Table 8.** Performance associated with the mean and the standard deviation of the seismic responses.

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*

*9.1.1. Training and testing the meta-model*

responses are shown in Table 7.

mean-squared error, i.e. *RRMSE*, and *R2*

*i*

, and the standard deviation, *σR i*, of *i*th seismic response during RBDO

as the absolute fraction of variance, during testing the

(46)

(47)

of the frame, the proposed meta-models are trained based on the generated database. The WWLS-SVM training during GSA is performed according to five-fold cross-validation. The lower and upper bonds of the parameters required in the optimization process are selected as *γ* ∈ 1.0, 500 , *a*∈ 0.5, 5.0 and *ω*0∈ 1.0, 10 , respectively. Therefore, the training optimal parameters of the meta-model associated with the mean and the standard deviation of seismic

292 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

In order to validate the performance and accuracy of the proposed meta-model, relative root-

( )


1

*n yy*

*i*

*n*

å

*i n i i*

=

=

2 1

1

*R*

=

*n*

å

( )


*i i*

*n i i*

å

=

2

**Mean Standard deviation**

*γ a ω*<sup>0</sup> *γ a ω*<sup>0</sup>

*n y*

2 1

*y*

=

*umax* 375.73 2.892 4.274 284.62 1.647 4.048

*DRmax* 386.04 3.804 6.347 314.42 2.104 4.702

*DILmax* 400.48 3.615 5.895 373.67 2.548 3.692

*DIoverall* 365.37 4.052 6.329 308.38 1.947 4.082

The smaller *RRMSE* and *MAPE* and the larger *R2* are indicative of better performance gener‐ ality. The comparison of WWLS-SVM and WLS-SVM, with respect to *MAPE*, *RRMSE*, and *R2*

**Table 7.** Optimal parameters of the meta-model for training the mean and the standard deviation of the seismic

å

( )

*y y*

æ ö ç ÷ - = - ç ÷ ç ÷ ç ÷ è ø

*i i*

2

2 1 1

meta-model and WLS-SVM, are defined using the following equations:

*RRMSE*

To predict the mean, *R*¯

**Parameter**

responses

is shown in Table 8.

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, *βannual*, 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, *βannual*, 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.


**Table 9.** Optimum designs obtained by DGSA and BGSA

**Figure 9.** Convergence histories of the best solution of DGSA and BGSA

### **9.2. Example 2: Nine-storey RC frame**

**Element groups no.**

*βannual*

*βannual*

*βannual*

4200

5200

6200

7200

Total cost (Euro)

8200

9200

**Table 9.** Optimum designs obtained by DGSA and BGSA

**Figure 9.** Convergence histories of the best solution of DGSA and BGSA

**DGSA BGSA Cross-section** *ρ* **(***%***) Cross-section** *ρ* **(***%***)**

C1 55 × 55 2.67 55 × 55 3.23 C2 55 × 55 2.64 55 × 55 3.20 C3 45 × 45 2.56 45 × 45 2.60 C4 45 × 45 2.33 45 × 45 2.48 B1 50 × 45 1.98 50 × 45 2.38 B2 50 × 40 2.13 50 × 40 2.30 B3 45 × 40 1.81 45 × 40 1.88 B4 45 × 35 1.69 45 × 35 1.78

294 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

*CIC* (Euro) 3448 3542 *CRC* (Euro) 1080 1045 *CTOT*(Euro) 4528 4587 Iterations 89 118

*Operational* 1.4435 1.6901

*Life safety* 2.4885 2.7342

*Collapse* 2.7986 3.1384

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Iteration

DGSA BGSA 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 metamodel 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 *MAPE, RRMSE* and *R2* .


**Table 10.** Performance associated with the mean and the standard deviation of seismic responses.

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, *βannual*, 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 approx‐ imate analyses by the meta-model, respectively.


**Table 11.** Optimum designs obtained by DGSA and BGSA

**Figure 10.** Convergence histories of the best solution of DGSA and BGSA

## **10. Conclusions**

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, *βannual*, obtained

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 approx‐

> C1 65 × 65 2.49 65 × 65 2.90 C2 60 × 60 2.50 60 × 60 3.08 C3 55 × 55 2.53 55 × 55 2.87 C4 50 × 50 2.49 50 × 50 2.79 C5 45 × 45 2.27 45 × 45 2.68 C6 40 × 40 2.29 40 × 40 2.89 B1 55 × 45 1.88 55 × 45 2.18 B2 55 × 40 1.80 55 × 40 1.98 B3 50 × 45 1.82 50 × 45 2.01 B4 50 × 40 1.87 50 × 40 1.98 B5 45 × 40 1.68 45 × 40 1.79 B6 45 × 35 1.66 45 × 35 1.85

*CIC* (Euro) 5571 5736 *CRC* (Euro) 1004 987 *CTOT*(Euro) 6575 6723 Iterations 83 123

*Operational* 1.4985 1.7101

*Life safety* 2.4953 2.7996

*Collapse* 2.8514 3.1837

**DGSA BGSA Cross-section** *ρ* **(***%***) Cross-section** *ρ* **(***%***)**

corresponding to each performance level by DGSA and BGSA is shown in Table 11.

296 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

imate analyses by the meta-model, respectively.

**Element groups No.**

*βannual*

*βannual*

*βannual*

**Table 11.** Optimum designs obtained by DGSA and BGSA

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 compu‐ tational cost.


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 [68]).

## **Acknowledgements**

Our special thanks go to Dr. Eysa Salajegheh (Distinguished Professor of Structural Engineer‐ ing) in Department of Civil Engineering at Shahid Bahonar University of Kerman, Iran, for his cooperation in this research work.

## **Author details**

Mohsen Khatibinia1 , Sadjad Gharehbaghi2 and Abbas Moustafa3\*

\*Address all correspondence to: abbas.moustafa@yahoo.com

1 Department of Civil Engineering, University of Birjand, Birjand, Iran

2 Department of Civil Engineering, Behbahan Khatam Alanbia University of Technology, Behbahan, Iran

3 Department of Civil Engineering, Minia University, Minia, Egypt

## **References**

analysis of SSI system, can estimate the reliability index through MCS with a small compu‐

**•** 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 metamodel was higher in comparison to WLS-SVM. Hence, the proposed meta-model can predict

**•** 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

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

Our special thanks go to Dr. Eysa Salajegheh (Distinguished Professor of Structural Engineer‐ ing) in Department of Civil Engineering at Shahid Bahonar University of Kerman, Iran, for his

2 Department of Civil Engineering, Behbahan Khatam Alanbia University of Technology,

and Abbas Moustafa3\*

the nonlinear dynamic analysis of SSI system in terms of accuracy and flexibility.

298 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

tational cost.

quantities (see [68]).

**Author details**

Mohsen Khatibinia1

Behbahan, Iran

**Acknowledgements**

cooperation in this research work.

, Sadjad Gharehbaghi2

\*Address all correspondence to: abbas.moustafa@yahoo.com

1 Department of Civil Engineering, University of Birjand, Birjand, Iran

3 Department of Civil Engineering, Minia University, Minia, Egypt

design more quickly than the BGSA model.


[27] Gharehbaghi, S. & Khatibinia, M. (2015). Optimal seismic design of reinforced con‐ crete structures subjected to time-history earthquake loads using an intelligent hy‐ brid algorithm*. Earthquake Engineering and Engineering Vibration*, Vol. 14(1), pp. 97-109.

[14] Khatibinia, K.; Salajegheh, J.; Fadaee, M.J. & Salajegheh, E. (2013). Seismic reliability assessment of RC structures including soil-structure interaction using wavelet weighted least squares support vector machine, *Reliability Engineering and System*

[15] Salajegheh, E.; Gholizadeh, S. & Khatibinia, M. (2008). Optimal design of structures for earthquake loads by a hybrid RBF-BPSO method, *Earthquake Engineering and Vi‐*

[16] Gholizadeh, S. & Salajegheh, E. (2009). Optimal design of structures subjected to time history loading by swarm intelligence and an advanced meta-model*, Computer Meth‐*

[17] Chen, K.Y. (2007). Forecasting systems reliability based on support vector regression with genetic algorithms, *Reliability Engineering and System Safety*, Vol. 92, pp. 423-32.

[18] Zhiwei, G. & Guangchen B. (2009). Application of least squares support vector ma‐ chine for regression to reliability analysis, *Chinese Journal of Aeronautics*, Vol. 22, pp.

[19] Tan, X.H.; Bi, W.H.; Hou, X.L. & Wang, W. (2011). Reliability analysis using radial basis function networks and support vector machines, *Computers and Geotechnics*, Vol.

[20] Moura, M.C.; Zio, E.; Lins, I.D. & Droguett, E. (2011). Failure and reliability predic‐ tion by support vector machine regression of time series data, *Reliability Engineering*

[21] Dai, H.; Zhang, H. & Wang, W. (2012). A support vector density-based importance sampling for reliability assessment, *Reliability Engineering and System Safety*, Vol. 106,

[22] Khatibinia, M.; Salajegheh, E.; Salajegheh, J. & Fadaee M.J. (2013). Reliability-based design optimization of RC structures including soil-structure interaction using a dis‐ crete gravitational search algorithm and a proposed meta-model, *Engineering Optimi‐*

[23] Rashedi, E.; Nezamabadi-pour, H. & Saryazdi, S. (2009). GSA: a gravitational search

[24] He, S.; Wu, Q.H.; Wen, J.Y.; Saunders, J.R. & Paton, R.C. (2004). A particle swarm op‐

[25] Rashedi, E.; Nezamabadi-pour, H. & Saryazdi, S. (2010). BGSA: binary gravitational

[26] Suykens, J.A.K.; Brabanter, J.D.; Lukas, L. & Vandewalle, J. (2002). Weighted least squares support vector machines: robustness and sparse approximation, *Neurocom‐*

algorithm, *Information Science*, Vol. 179(13), pp. 2232-2248.

search algorithm, *Natural Computing*, Vol. 9(3), pp. 727-745.

timizer with passive congregation, *BioSystems*, Vol. 78, pp. 135-147.

*ods in Applied Mechanics and Engineering*, Vol. 198, pp. 2936-2949.

300 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

*Safety*, Vol. 110, pp. 22-33.

160-166.

38, pp. 178-86.

pp. 86-93.

*bration Engineering*, Vol. 7, pp. 14-24.

*and System Safety*, Vol. 96, pp. 1527-34.

*zation*, Vol. 45 (10), pp. 1147-1165.

*puting*, Vol. 48, pp. 85-105.


[39] Khatibinia, M. & Naseralavi, S. S. (2014). Truss optimization on shape and sizing with frequency constraints based on orthogonal multi-gravitational search algorithm,

302 Earthquake Engineering - From Engineering Seismology to Optimal Seismic Design of Engineering Structures

[40] Shinozuka, M. & Sato, Y. (1967). Simulation of nonstationary random processes. *Jour‐*

[41] Mckay, M.D.; Beckman, R.J. & Conover, W.J. (1979). A comparison of three methods for selecting values on input variables in the analysis of output from a computer

[42] ACI318 Committee. (2005). *Building code requirements for structural concrete and com‐*

[43] Kunnath, S.K.; Reinhorn, A.M. & Lobo, R.F. (1992*). IDARC Version 3: A program for the inelastic damage analysis of RC structures*. Technical Report NCEER-92-0022, National Center for Earthquake Engineering Research, State University of New York, Buffalo.

[44] Park, Y.J. & Ang, A.H. (1985). Mechanistic seismic damage model for reinforced con‐

[45] Moustafa, A. (2011). Damage-based design earthquake loads for SDOF inelastic

[46] Moustafa, A. & Mahmoud, S. (2014). Damage assessment of adjacent buildings under

[47] Heidari, A. & Gharehbaghi, S. (2015). Seismic performance improvement of special truss moment frames using damage and energy concepts, *Earthquake Engineering and*

[48] Gharehbaghi, S.; Salajegheh, E.; Shojaee, S. & Khatibinia, M. (2009). Optimum seismic design of RC structures considering damage index, *proceeding of the 1st international*

[49] Mazzoni, S.; McKenna, F.; Scott, M.H. & Fenves G.L., (2010). *OpenSEES: Open system for earthquake engineering simulation*. Pacific Earthquake Engineering Research Centre

[50] Kent, D.C. & Park, R. (1997). Flexural members with confined concrete*, Structural Di‐*

[51] Saatcioglu, M. & Razvi, S.R. (1992). Strength and ductility of confined concrete, *Jour‐*

[52] Paulay, T. & Priestley, M.J.N. (1992). *Seismic design of reinforced concrete and masonry*

*conference of lightweight structures and earthquake*, Kerman, Iran (In Persian).

structures*, Journal of Structural Engineering*, ASCE, Vol. 137, pp. 456–467.

*Sound and Vibration*, Vol. 333, pp. 6349–6369.

code, *Technometrics,* Vol. 21, pp. 439-445.

*nal of Engineering Mechanics,* ASCE, Vol. 93, pp. 11-40.

*mentary*, American Concrete Institute, Michigan, USA.

crete*, Journal of Structural Engineering,* Vol. 111, pp. 722-39.

earthquake loads, *Engineering Structures*, Vol. 61, pp. 153–165.

*Structural Dynamics*, DOI: 10.1002/eqe.2499, (in press).

*nal of Structural Engineering,* ASCE, Vol. 118, pp. 1590-1607.

*buildings*, 1st Ed. New York: Wiley-Interscience.

(PEER): http://opensees.berkeley.edu/.

*vision*, ASCE, Vol. 97, pp. 1969-1990.

