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

1. Introduction

Past destructive earthquakes (e.g. the 1994 Northridge earthquake and the 1995 Kobe earthquake) have left a clear signature on the engineering community worldwide, changing thinking of structural engineers [1-2]. As such, after holding several workshops and conferences, an innovative approach namely Performance-Based Design (PBD) was presented by modern guidelines [3-5]. In principle, a structure designed using PBD approach should meet performance objectives in accordance with a set of specified reliabilities over its service life. This is aimed to reach structural design candidates associated with more predictable seismic behavior, quantifying and controlling the risk at an engineered acceptable level.

Both seismic demands and capacity parameters, that are inherently uncertain, are highly influential on the acceptable performance level of a structure. Furthermore, due to the fact that a structure on underlying soil is not rigid, soil-structure interaction (SSI) affects the responses of structures during an earthquake. Obviously, ignoring the SSI effects could lead to unrealistic structural responses and seismic demands. Hence, the effects of SSI should be considered in the seismic responses of structures [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 manufacturing and operating processes may lead to unreliable design resulting in inappropriate balance between cost and safety. A proper design procedure must reasonably account for the inherent uncertain nature of a structural system and associated external load [7]. In structural optimization, non-deterministic performance of structures can be taken into account using 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 problems 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 computational effort. This drawback may accentuate when the nonlinear dynamic structural responses are required in RBDO of the structure using the Monte-Carlo Simulation (MCS) method, the importance sampling technique and the response surface method. In order to obtain an acceptable confidence within probabilities of the order close to 10^-4 - 10^-6, the MCS method requires a large number of structural analyses. Based on Lagaros et al. [11], for such results, the large number of analyses ranges from $$1.6\times {10}^6$$ to $$1.6\times {10}^9$$ 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

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 probability of non-performance with target values for each performance level. The objective of the RBDO problem is to minimize the total cost whereas the deterministic and probabilistic constraints should not exceed a specified target.

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

$$\begin{array}{l}Minimize{C}_{TOT}\\ Subjectto{g}_i\left(X,{\overline{X}}_{rand}\right)\le 0\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}i=1,2,...,N_g\\ {\beta }_{annual}\left[G_j\left(X,{\overline{X}}_{rand}\right)\ge 0\right]\ge {\beta }_{\begin{array}{l}j\\ \end{array}}^{Target}j=1,2,3\\ X=\{x_1,x_2,...,x_k,...,x_n\}\in R^d\end{array}$$E1

where C_TOT is the total cost; $${\beta }_{annual}$$ is the system reliability index corresponding to the jth performance level (performance function), $$G_j(X,{\overline{X}}_{rand})$$ ; $${\beta }_j^{Target}$$ is the prescribed target value of the reliability index; N_g is the number of deterministic constraints, $$g_i(X)$$ ; A given set of discrete values, X, is expressed by R^d and design variables, x_k, can take values only from this set. The vector $${\overline{X}}_{rand}$$ represents the random variables.

2.2. Life- cycle cost assessment of RC structure

The total cost, C_TOT, 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:

$$C_{TOT}\left(X,{\overline{X}}_{rand},t\right)=C_{IC}\left(X,{\overline{X}}_{rand}\right)+C_{RC}\left(X,{\overline{X}}_{rand},t\right)$$E2

where C_IC is the initial cost of structure; C_RC 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, C_C, and the total cost of reinforcing steel, C_S, is given [14]:

$$C_{IC}=C_C+C_S={\displaystyle \sum _{i=1}^{Ne}{w}_C{b}_i{h}_i{L}_i}+{\displaystyle \sum _{i=1}^{Ne}{w}_S{A}_{Si}{L}_i}$$E3

where $$w_C$$ and $$w_S$$ are the unit cost coefficients of each material; b_i, h_i, and L_i are the section dimensions of ith element and its length; and A_Si is the total reinforcement in the section of element. N_e 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, DI_overall, 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]:

$$C_{RC}={\displaystyle {\int }_0^1{C}_{RC}|{}_{D{I}_{overall}}\cdot f_{D{I}_{overall}}\left(D{I}_{overall}\right)}\cdot d\left(D{I}_{overall}\right)$$E4

where $$f_{D{I}_{overall}}$$ is the probability density function for the index, DI_overall. $$C_{RC}|{}_{D{I}_{overall}}$$ is the expected present cost which is defined as [34]:

$$C_{RC}|{}_{D{I}_{overall}}={\displaystyle {\int }_0^{\infty }{C}_{f0}\left(D{I}_{overrall}\right)f\left(t\right))dt=}{C}_f\left(D{I}_{overrall}\right)\frac{\upsilon }{r+\upsilon }$$E5

$$C_f\left(D{I}_{overall}\right)=\{\begin{array}{l}k{C}_0\left(\frac{D{I}_{overall}}{0.60}\right)\text{if}D{I}_{overall}\le 0.60\\ k{C}_0\text{if}D{I}_{overall}>0.60\end{array}$$E6

where $$\upsilon$$ and r are the mean occurrence of earthquakes and the discount rate, respectively. C₀ 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

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:

$$fit\left(X,{\overline{X}}_{rand}\right)=C{}_{TOT}\left(1+r_{p}PF\right)$$E7

where $$fit(X)$$, PF and r_p 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={\displaystyle \sum _i{\left[\max \left(g_i,0\right)\right]}^2}+{\displaystyle \sum _j{\left[\max \left(1-\frac{{\beta }_{annual,j}}{{\beta }_j^{Target}},0\right)\right]}^2}$$E8

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 earthquakes, 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/sec²). 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 longitudinal bars of structural elements’ cross-sections, $$\rho$$, 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, u_max, maximum inter-storey drift, DR_max, maximum local damage index, DIL_max, and overall damage index, DI_overall, are considered as the structural seismic responses. After that, the mean, $${\overline{R}}_i$$, and the standard deviation, $${\sigma }_{Ri}$$, of ith seismic response, R_i, 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 DIL_max and DI_overall.

3.2. Limit state functions

The operational, life safety and collapse prevention levels have been defined as the performance levels. A performance level depends on some limit state functions. A limit state function, $$G(\overline{X})$$, is determined by capacity, $$R_{{}_{LIM}}$$, and demand, $$R(\overline{X})$$, as follows [13, 22]:

$$G\left(\overline{X}\right)=R_{LIM}\left(\mu ,\delta \right)-R\left(\overline{X}\right)$$E9

where $$R_{{}_{LIM}}$$ is the limiting value for a seismic response $$R(\overline{X})$$ at a given performance level, with mean value of $$\mu$$ and coefficient of variation $$\delta$$.

The limit state functions and their probability distribution function (PDF) for the performance levels, according to SEAOC guidelines (2000), are shown in Table 1.

Performance level	Limit sate function	G $$(\overline{\text{X}})$$	PDF
Operational
	Elastic roof displacement	$$({\overline{u}}_y,0.1)-u_{\max }(\overline{X})$$	Lognormal
	Inter-storey drift	$$(0.005,0.1)-D{R}_{max}(\overline{X})$$	Lognormal
Life safety
	Inter-storey drift	$$(0.015,0.1)-DR{}_{max}(\overline{X})$$	Lognormal
	Max. local damage index	$$(0.6,0.1)-DI{L}_{max}(\overline{X})$$	Beta
	Global damage index	$$(0.4,0.1)-D{I}_{overall}(\overline{X})$$	Beta
Collapse
	Inter-storey drift	$$(0.025,0.1)-D{R}_{max}(\overline{X})$$	Lognormal
	Max. local damage index	$$(1,0.1)-DI{L}_{max}(\overline{X})$$	Beta
	Global damage index	$$(0.4,0.1)-D{I}_{overall}(\overline{X})$$	Beta

Table 1.

The limit sate functions of the performance levels [14, 22]

In Table 1, $${\overline{u}}_y$$ is the yielding horizontal displacement at top storey of the frame which is determined by a pushover analysis. Furthermore, the demand, $$R(\overline{X})$$, corresponding to each of seismic responses are defined using the mean, $$\overline{R}$$, and the standard deviation, $${\sigma }_R.$$

3.3. Annual probability of non-performance

The non-performance probability, $$P_f$$, 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(\overline{X})\le 0$$, $$P_f$$ is calculated as follows:

$$P_f={\displaystyle \underset{G(\overline{X})\le 0\text{}}{\int }{\displaystyle \int ...{\displaystyle \int f_{\overline{X}}\left(\overline{X}\right))d\text{}\overline{X}}}}$$E10

where $$f_{\overline{X}}(\overline{X})$$ is the joint probability density function of $$\overline{X}$$.

For each performance level, the total exceeding probability, $$P{f}_E$$, is considered as a series system. Determination of the total exceeding probability, $$P{f}_E$$, 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]:

$$P{f}_{annual}=1-\exp \left(-\upsilon P{f}_E\right)\to {\beta }_{annual}=-{\Phi }^{-1}\left(P{f}_{annual}\right)$$E11

where $${\beta }_{annual}$$ can be expressed as an reliability index as shown in Eq. (11), using the standard normal cumulative distribution function, $$\Phi (\cdot )$$.

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:

$$M\Delta \ddot{u}+C\Delta \dot{u}+K_{\text{T}}\Delta u=-m_x{\ddot{u}}_{g,x}\left(t+\Delta t\right)-F\left(t\right)$$E12

where $$M,C$$ and $$K_T$$ are mass, damping, and tangent stiffness matrices of SSI model, respectively; $$\Delta u\text{}\text{}$$ is the incremental vector of the relative displacements for SSI system between times $$t$$ and $$t+\Delta t$$ ; and $$F(t)$$ is the vector of internal forces at time t. The term $${\ddot{u}}_{g,x}(t+\Delta t)$$ is the free-field component of acceleration in x direction. The column matrix, $$m_x$$, 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 }_m-{\theta }_y}{{\theta }_u-{\theta }_y}+\frac{\beta }{{\theta }_u{M}_y}{\displaystyle \int d{E}_H}$$E13

where ${\theta }_m$ is the maximum rotation during loading history; ${\theta }_u$ and ${\theta }_y$ are the ultimate and yield rotation, respectively; $M_y$ is the yield moment; and $${\displaystyle \int d{E}_H}$$ 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:

$$D{I}_{story}={\displaystyle \sum _{i=1}^{Ne}D{I}_i.\cdot {\lambda }_{i,component};}{\lambda }_{i,component}={\left[E_i/{\displaystyle \sum E_i}\right]}_{component}$$E14

$$D{I}_{overall}={\displaystyle \sum _{i=1}^{Ns}D{I}_{i,story}.{\lambda }_{i,\text{}story};}{\lambda }_{i,story}={\left[E_i/{\displaystyle \sum E_i}\right]}_{story}$$E15

where ${\lambda }_i$ is energy weighting factor; and $E_i$ 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

OpenSEES [49], as an open-source computational software framework, is used for by simulation of SSI system, and performing nonlinear dynamic analyses of SSI system depicted in Fig. 1. Assuming materials of constant properties over its depth, soil encompasses different layers, and the foundation is considered as rigid strip footing. Beams and columns of structure are modeled using force-based nonlinear beam-column element with considering the spread plasticity along the element’s length. The integration along each element is based on Gauss-Lobatto quadrature rule. Also, the infinite boundaries of soil are modeled using the artificial boundaries (Fig. 1). The model of soil-structure system shown in Fig. 1 was successfully used by [14, 22, 31].

The Kent-Scott-Park model [50] is utilized for modeling the confined and unconfined concrete of cross-sections of structural elements. The constitutive parameters of this model are: f_c=concrete peak strength in compression, f_u=residual strength, $${\epsilon }_0$$ =strain at peak strength, and $${\epsilon }_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:

$${\epsilon }_{uc}={\epsilon }_{uo}+1.4\frac{\rho {}_v{f}_{yh}{\epsilon }_{us}}{f_{cc}}$$E16

where $${\epsilon }_{uc}$$, $${\epsilon }_{uo}$$ and $${\epsilon }_{us}$$ are the ultimate compressive strain of confined and unconfined concrete, and the ultimate strain of longitudinal steel reinforcement in tensile stress, respectively; $$\rho {}_v,f_{yh}$$ and $$f_{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₂ plasticity model with linear hardening is utilized for modeling the constitutive behavior of the steel reinforcement. The material parameters defining J₂ plasticity model are: f_y=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]

Figure 2.

Material constitutive models; (a) Concrete, (b) Steel [22]

Soil layers are modeled using isoperimetric four-node quadrilateral finite elements and assuming bilinear displacement interpolation. The plane strain condition is assumed for the soil domain with considering a constant soil thickness corresponding to the inter-frame distance. The material of the soil is modeled using a modified pressure-independent multi-yield-surface J₂ 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_n=a\rho V_p;V_p=\sqrt{\frac{2G\left(1-v\right)}{\rho \left(1-2v\right)}}$$E17

$$C_s=b\rho V_s;V_s=\sqrt{\frac{G}{\rho }}$$E18

where $$C_n$$ and $$C_s$$ are the normal and shear damping, respectively; $$\rho$$ and $$v$$ are the mass density and Poisson ratio of soil, respectively; a and b are dimensionless parameters to be determined, and $$V_p$$ and $$V_s$$ 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 $$G_i$$ as low strain shear modulus, $$B_i$$ as bulk modulus, and $${\tau }_i$$ that is shear strength, with i=1,2,...,n representing the soil layers’ number. Other parameters of soil material depend on $$B_i,G_i$$, and $$V_{p,i}$$.

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 characteristics with a desired site. It is often difficult or impossible in some cases to choose a proper record for a site, since historically recorded accelerograms for a given site could be limited or scare. Hence, artificial earthquakes that are statistically influenced by desired properties of the given site are very useful for seismic design of structures. In this work, spectral representation method based on time domain procedure is used for the generation of synthetic ground motion records. The non-stationary ground motion is simulated using this method as [13]:

$$a\left(t\right)=I_m\left(t\right){{\displaystyle \sum _{n=1}^{NFR}\left\{4S_{KT}\left(n\Delta \text{}f\right)\left[1+{\delta }_S{R}_N\right]\right\}\Delta \text{}f}}^{1/2}\sin \left(2\pi n\Delta \text{}ft+{\theta }_n\right)$$E19

where $$a(t),I_m(t)$$ and $$S_{KT}(.)$$ 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 }$$, $${\delta }_S$$ and $$R_N$$ are the coefficient of variation and a standard normal variable that used in ordinates of PSDF, $$\Delta f$$ is frequency step, and $${\theta }_n$$ are random phase angles with a uniform distribution between 0 and $$2\pi$$. In this work, the modulation function expressed in [13] is used:

$$I_m\left(t\right)=\{\begin{array}{l}{\left(t/T_1\right)}^{d}0\le t\le T_1\\ 1T_1\le t\le T_2\\ e^{-c\left(t-T_2\right)}{T}_2\le t\le T_{}\end{array}$$E20

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:

$$S_{KT}\left(f\right)=S_0\left[\frac{1+4{\xi }_g^2{\left(f/f_g\right)}^2}{{\left(1-{\left(f/f_g\right)}^2\right)}^2+4{\xi }_g^2{\left(f/f_g\right)}^2}\right]\left[\frac{{\left(f/f_f\right)}^4}{{\left(1-{\left(f/f_f\right)}^2\right)}^2+4{\xi }_f^2{\left(f/f_f\right)}^2}\right]$$E21

where $$S_0$$ is the constant PSDF of input white-noise random process; $$f_g$$ and $${\xi }_g^{}$$ are the characteristic ground frequency and the ground damping ratio; $$f_f$$ and $${\xi }_f^{}$$ are parameters for a high-pass filter to attenuate low frequency components. As listed in Table 2, the parameters for the generation of simulated ground motion are selected according to the values proposed by Möller et al. [13].

Parameter	$$a_G\le 350cm/s^2$$ *	$$350\le a_G\le 700cm/s^2$$	$$a_G\ge 700cm/s^2$$
T (sec)	5.12	10.24	20.48
T1 (sec)	0.50	1.50	2.00
T2 (sec)	4.00	8.00	16.00
c	2.0	1.0	0.7
d	2.0	2.0	2.0
NFR	100	200	300
fmax (Hz)	12	15	15
$${\delta }_S$$	0.40	0.40	0.40

Table 2.

Parameters for the generation of simulated ground motion (* $$a_G$$ = $$\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.

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

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 ith agent in tth iteration is associated with a position vector, $$X_i^{}(t)=\{x_i^1,...,x_i^d,...,x_i^D\}$$, and a velocity vector, $$V_i^{}(t)=\{v_i^1,...,v_i^d,...,v_i^D\}$$. 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 minimization problem as follows:

$$M_i\left(t\right)=\frac{fi{t}_i\left(t\right)-worst\left(t\right)}{{\displaystyle \sum _{j=1}^N\left(fi{t}_j\left(t\right)-worst\left(t\right)\right)}}$$E22

where N, M_i(t) and $$fi{t}_i(t)$$ represent the population size, the mass and the fitness value of agent i at tth 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\left(t\right)={\displaystyle \sum _{j\in kbest,j\ne i}ran{d}_{j}G\left(t\right)}\frac{M_j\left(t\right)}{R_{ij}\left(t\right)+\epsilon }\left(x_j^d\left(t\right)-x_i^d\left(t\right)\right)$$E23

$$v_i^d\left(t+1\right)=ran{d}_i{v}_i^d\left(t\right)+a_i^d\left(t\right)$$E24

$$x_i^d\left(t+1\right)=x_i^d\left(t\right)+v_i^d\left(t+1\right)$$E25

where $$a_i^d$$, $$v_i^d$$ and $$x_i^d$$ present the acceleration, velocity and position of ith agent in dimension d, respectively. $$ran{d}_i$$ and $$ran{d}_j$$ are two uniformly distributed random numbers in the interval [0, 1]; $$\epsilon$$ is a small value; and $$R_{i,j}(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_0$$ at the beginning and decreased with time. Here, $$k_0$$ is set to N and is decreased linearly to 1.0. G is a decreasing function of time.

6.2. The proposed discrete GSA

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(\text{}{v}_i^d(t))$$ [25]. For updating the position of the agents, Eq. (25) is re-defined as follows [25]:

$$x_i^d\left(t+1\right)=\{\begin{array}{l}complement\left(x_i^d\left(t\right)\right)\text{for}rand<S\left(v_i^d\left(t\right)\right)\\ x_i^d\left(t\right)\text{for}rand\ge S\left(v_i^d\left(t\right)\right)\end{array}$$E26

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{array}{l}{v}_i^d\left(t+1\right)=r_{i,1}{v}_i^d\left(t\right)+a_i^d\left(t\right)+r_{i,2}\left(p_{best,i}^d\left(t\right)-x_i^d\left(t\right)\right)+\\ r_{i,3}\left(p_g^d\left(t\right)-x_i^d\left(t\right)\right)+r_{i,4}\left(p_i^d(t)-x_i^d\left(t\right)\right)\end{array}$$E27

where $$p_{best}^d$$, $$p_g^d$$ and $$p_i^d$$ are the best previous position of the ith agent, the best previous position among all the agents and a randomly selected agent from the population, respectively. r_i,1 and r_i,2 are the random number in the interval [0, 1], respectively; r_i,3 is the uniform random number in the interval (0, 1).

In DGSA, the scalar $$x_i^d\in \{1,2,\mathrm{...},n\}$$ corresponds to discrete values of the set A={A₁, A₂,...,A_n}. Hence, the position of agents is updated by the following equation instead of Eq. (26):

$$x_i^d\left(t+1\right)=INT\left(x_i^d\left(t\right)+v_i^d\left(t\right)\right)$$E28

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:

$$x_{i\text{}}^d\left(t+1\right)=\{\begin{array}{l}{P}_g^{d}iffit\left(p_g^d\left(t\right)\right)=fit\left(p_g^d\left(t-1\right)\right)\\ P_{i,best\text{}}^d(t)iffit\left(p_g^d\left(t\right)\right)\ne fit\left(p_g^d\left(t-1\right)\right)\\ INT\left(x_{i\text{}}^d\left(t\right)+v_{i\text{}}^d\left(t\right)\right)otherwise\end{array}$$E29

where $$INT(\cdot )$$ 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 meta-model 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

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\left(w,e\right)=\frac{1}{2}{\Vert w\Vert }^2+\frac{1}{2}\gamma {\displaystyle \sum _{i=1}^n{\overline{v}}_i{e}_i^2}$$E30

Subject to the following equality constraints:

$$y_i=w^T\phi \left(x_i\right)+b+e_i,i=1,2,...,n$$E31

where $${\{x_i,\text{}{y}_i\}}_{i=1}^n$$ is a training data set; $$x_i\in R^n$$ and $$y_i\in R$$ represent the input and output data. Operator $$\phi \text{}(.):R^n\to R^d$$ is a function which maps the input space into a higher dimensional space. The vector, $$w\in R^d$$, represents weight vector in primal weight space. The symbols, $$e_i\in R$$ and $$b\in 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:

$$y\left(x\right)=w^T\phi \left(x\right)+b$$E32

It is impossible to indirectly compute $$w$$ from Eq. (30), for the structure of the function $$\phi (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\text{};x\right)=J\left(w,e\right)-{\displaystyle \sum _{i=1}^n{\alpha }_i\left(w^T\phi \left(x_i\right)+b+e_i-y_i\right)}$$E33

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

$$\left[\begin{array}{cc}\Omega +V_{\gamma }& {\text{1}}_n^T\\ {\text{1}}_n& \text{0}\end{array}\right]\left[\begin{array}{c}\text{α}\\ b\end{array}\right]=\left[\begin{array}{c}y\\ 0\end{array}\right]$$E34

where

$$\begin{array}{l}{V}_{\gamma }=diag\left\{1/\gamma \text{}{\overline{v}}_1,...,1/\gamma \text{}{\overline{v}}_n\right\};{\Omega }_{i,j}={⟨\phi \left(x_i\right),\phi \left(x_j\right)⟩}_{H}i,j=1,...,n\\ y={\left[y_1,...,y_n\right]}^T;{\text{1}}_n^T=\left[1,...,1\right];\alpha =\left[{\alpha }_1,...,{\alpha }_n\right]\end{array}$$E35

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

$$K\left(x,\overline{x}\right)={⟨\phi \left(x\right),\phi \left(\overline{x}\right)⟩}_H$$E36

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

$$y\left(x\right)={\displaystyle \sum _{i=1}^n{\alpha }_{i}K\left(x_i,x\right)+b}$$E37

Weight $${\overline{v}}_k$$ is estimated as follows [57-58]:

$${\overline{v}}_k=\{\begin{array}{c}1if\left|e_i/\widehat{s}\right|\le c_1\\ \frac{c_2-\left|e_i/\widehat{s}\right|}{c_2-c_1}if{c}_1<\left|e_i/\widehat{s}\right|\le c_2\\ {10}^{-4}otherwise\end{array}$$E38

where $$\widehat{s}$$ is a robust estimation of the standard deviation for the error variables ($$e_i=a_i/D_{ii}^{-1}$$); constants $$c_1$$ and $$c_2$$ are typically chosen as $$c_1=2.5$$ and $$c_2=3$$. Here $$D_{ii}^{-1}$$ denotes the ith 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 $${\overline{v}}_k$$ 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\left(x,\overline{x}\right)=\exp \left(-\frac{{\Vert x-\overline{x}\Vert }^2}{{\sigma }^2}\right)$$E39

where $${\sigma }^2$$ 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_k,\text{}{y}_k\}}_{k=1}^{Ntot}$$, set N=Ntot.

Step 2. Find an optimum $$(\gamma ,\sigma )$$ 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 $$\left|\alpha \right|$$.

Step 4. Remove a small number of M points (typically 5% of the N points) which has the smallest values in the sorted $$\left|\alpha \right|$$.

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 \left(t\right)=\frac{1}{\sqrt{a}}\cos \left({\omega }_0\left(\frac{t-b}{a}\right)\right)\exp \left(-0.5{\left(\frac{t-b}{a}\right)}^2\right)$$E40

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(x,\overline{x}\right)={\displaystyle \prod _{i=1}^n\Psi \left(\frac{x_i-{\overline{x}}_i}{a}\right)}$$E41

where n is the number of samples; $$x$$ and $$\overline{x}\in R^n$$

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

$$K\left(x,\overline{x}\right)={\displaystyle \prod _{i=1}^n\frac{1}{\sqrt{a}}\cos \left({\omega }_0\left(\frac{x_i-{\overline{x}}_i}{a}\right)\right)\exp \left(-0.5\frac{{\Vert x_i-{\overline{x}}_i\Vert }^2}{a^2}\right)}$$E42

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, $$\gamma$$, and the wavelet kernel parameters, a and $${\omega }_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}{\displaystyle \sum _{i=1}^{n}100\times \left|\frac{y_i-{\overline{y}}_i}{y_i}\right|}$$E43

where y and $$\overline{y}$$ are the actual value and the predicted value, respectively.

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

$$\begin{array}{l}MinimizeMAP{E}_{k-foldCV}\\ Subjectto{\gamma }_{\min }\le \gamma \le {\gamma }_{\max }\\ a_{\min }\le a\le a_{\max }\\ {\omega }_{0\min }\le {\omega }_0\le {\omega }_{0\max }\end{array}$$E44

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

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 $$X_{FCM}$$. This optimum solution is copied N_FCM times to create the some part of the initial swarm FPSO. Other particles of the initial swarm, i.e. $${\text{X}}_{rnd,j}(j=1,2,\mathrm{...},N_{PSO}-{\text{N}}_{\text{FCM}})$$, 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 parameters 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 least-squares 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

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/m³ for all soil layers. The values of the dead and live loads are considered to be DL=5.884 N/mm² (600 kg/m²) and LL=1.961 N/mm² (200 kg/m²), 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 $$w_C=60$$ and $$w_S=700$$ Euros. To calculate the total expected cost of repair, first, the cumulative distribution is obtained using MCS and the proposed meta-model for the response DI_overall 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, $${\epsilon }_{0o}$$, and the residual strength, f_uo, for unconfined concrete are selected as 0.002 and 0.0, respectively. Shear-wave velocity, V_s, and friction angle, $$\phi$$, 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, a_g, and the central frequency, f_g, for soil filter are considered as random variables in the RBDO process. The PGA value is also taken into account as follows [13, 22]:

$$a_g={\overline{a}}_g\left(1+{\sigma }_{\overline{a}g}\right)$$E45

where $${\overline{a}}_g$$ and $${\sigma }_{\overline{a}g}$$ are the mean and the standard deviation of PGA, respectively. The values of these parameters are shown in Table 6.

Number	Column (h $$\times$$ h)^a			Beam (h $$\times$$ b)^a
Type 1	Type 2	Type 3	Type 1	Type 2	Type 3
1	65 $$\times$$ 65	55 $$\times$$ 55	45 $$\times$$ 45	55 $$\times$$ 45	50 $$\times$$ 45	45 $$\times$$ 40
2	60 $$\times$$ 60	50 $$\times$$ 50	40 $$\times$$ 40	55 $$\times$$ 40	50 $$\times$$ 40	45 $$\times$$ 35
3	55 $$\times$$ 55	50 $$\times$$ 50	40 $$\times$$ 40	55 $$\times$$ 35	50 $$\times$$ 35	45 $$\times$$ 30
4	55 $$\times$$ 55	45 $$\times$$ 45	40 $$\times$$ 40	55 $$\times$$ 30	50 $$\times$$ 30	45 $$\times$$ 25
5	55 $$\times$$ 55	45 $$\times$$ 45	35 $$\times$$ 35	50 $$\times$$ 40	45 $$\times$$ 40	40 $$\times$$ 35
6	50 $$\times$$ 50	45 $$\times$$ 45	40 $$\times$$ 40	50 $$\times$$ 35	45 $$\times$$ 35	40 $$\times$$ 30
7	55 $$\times$$ 55	45 $$\times$$ 45	35 $$\times$$ 35	50 $$\times$$ 30	45 $$\times$$ 30	40 $$\times$$ 30
8	-	-	-	50 $$\times$$ 30	45 $$\times$$ 30	35 $$\times$$ 30

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 G₀ = 50 is utilized in this study. In order to consider the stochastic nature of the optimization process, ten independent optimization runs are performed and the best solution is considered as the final results.

9.1. Example 1: Six-storey RC frame

Six-storey RC frame is shown in Fig. 8(a). In the frame, the length of each bay and the height of stories are 5 m and 3 m, respectively. The members of the structure are divided into four groups for the columns C₁, C₂, C₃, C₄ and four groups B₁, B₂, B₃ and B₄ for the beams. The groups of structural elements are presented in Fig. 8(a).

Material	Variable	PDF	Mean	SD
Concrete
	fco (Mpa)	Lognormal	28	2
	$${\epsilon }_{uo}$$	Normal	0.0035	0.00035
Steel
	fy (Mpa)	Lognormal	340	25
	E (Mpa)	Lognormal	210000	8000
	H	Normal	0.015	0.0015
Soil
	Vs,1 (m/s)	Normal	375	10
$${\phi }_1$$	Normal	37.5	1.0
Vs,2 (m/s)	Normal	300	20
$${\phi }_2$$	Normal	37.5	1.0
Vs,3 (m/s)	Normal	200	10
$${\phi }_3$$	Normal	32.5	1.0

Table 5.

The marginal probability distribution, mean and standard deviation of materials

Variable	PDF	Mean	SD
$${\overline{a}}_g$$ (cm/s^2)	Lognormal	300	100
$${\sigma }_{\overline{a}g}$$	Normal	-	0.15
$$f_g$$ (Hz)	Normal	2.5	0.375

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, $${\overline{R}}_i$$, and the standard deviation, $${\sigma }_{Ri}$$, of ith 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 $$\gamma \in [1.0,500],a\in [0.5,5.0]$$ and $${\omega }_0\in [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 root-mean-squared error, i.e. RRMSE, and R² as the absolute fraction of variance, during testing the meta-model and WLS-SVM, are defined using the following equations:

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

$$R^2=1-\left(\frac{{\displaystyle \sum _{i=1}^n{(y_i-{\overline{y}}_i)}^{\text{2}}}}{{\displaystyle \sum _{i=1}^n{\overline{y}}_i^2}}\right)$$E47

Parameter	Mean			Standard deviation
$$\gamma$$	$$a$$	$${\omega }_0$$	$$\gamma$$	$$a$$	$${\omega }_0$$
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

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 R² are indicative of better performance generality. The comparison of WWLS-SVM and WLS-SVM, with respect to MAPE, RRMSE, and R² is shown in Table 8.

Method	Statistical parameters	Mean				Standard deviation
umax	DRmax	DILmax	DIoverall	umax	DRmax	DILmax	DIoverall
WWLS-SVM
	MAPE	2.008	2.417	2.337	2.172	2.163	2.075	2.532	3.406
	RRMSE	0.018	0.022	0.020	0.021	0.023	0.024	0.045	0.028
	R^2	0.9999	0.9998	0.9998	0.9999	0.9999	0.9998	0.9998	0.9997
WLS-SVM
	MAPE	5.328	6.302	5.392	5.862	5.386	6.017	6.007	5.737
	RRMSE	0.057	0.074	0.059	0.061	0.079	0.094	0.091	0.064
	R^2	0.9988	0.9979	0.9987	0.9987	0.9985	0.9984	0.9988	0.9988

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

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, $${\beta }_{annual}$$, are estimated using WWLS-SVM-based MCS with 10⁶ 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, $${\beta }_{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.

Element groups no.	DGSA		BGSA
Cross-section	$$\rho (\%)$$	Cross-section	$$\rho (\%)$$
C1	55 $$\times$$ 55	2.67	55 $$\times$$ 55	3.23
C2	55 $$\times$$ 55	2.64	55 $$\times$$ 55	3.20
C3	45 $$\times$$ 45	2.56	45 $$\times$$ 45	2.60
C4	45 $$\times$$ 45	2.33	45 $$\times$$ 45	2.48
B1	50 $$\times$$ 45	1.98	50 $$\times$$ 45	2.38
B2	50 $$\times$$ 40	2.13	50 $$\times$$ 40	2.30
B3	45 $$\times$$ 40	1.81	45 $$\times$$ 40	1.88
B4	45 $$\times$$ 35	1.69	45 $$\times$$ 35	1.78
CIC (Euro)	3448		3542
CRC (Euro)	1080		1045
CTOT (Euro)	4528		4587
Iterations	89		118
$${\beta }_{annual}^{Operational}$$	1.4435		1.6901
$${\beta }_{annual}^{Lifesafety}$$	2.4885		2.7342
$${\beta }_{annual}^{Collapse}$$	2.7986		3.1384

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

Nine-storey RC frame is shown in Fig. 8(b). In the frame, the length of each bay and the height of stories are 5 m and 3 m, respectively. The members of the structure are divided into six groups for the columns and six groups for the beams. The groups of structural elements are presented in Fig. 8(b).

9.2.1. Training and testing the meta-model

After training database using the presented WWLS-SVM optimal parameters of the meta-model associated with the mean and the standard deviation of seismic responses are shown in Table 10. Furthermore, the performance generality of the proposed meta-model and WLS-SVM is given in Table 10 in terms of MAPE, RRMSE and R².

Method	Statistical parameters	Mean				Standard deviation
umax	DRmax	DILmax	DIoverall	umax	DRmax	DILmax	DIoverall
WWLS-SVM
	MAPE	2.837	3.028	3.127	2.689	3.024	2.682	3.008	3.105
	RRMSE	0.026	0.038	0.035	0.023	0.032	0.024	0.031	0.036
	R^2	0.9999	0.9998	0.9998	0.9999	0.9998	0.9999	0.9998	0.9998
WLS-SVM
	MAPE	5.538	6.346	6.483	6.006	6.305	5.396	5.843	5.579
	RRMSE	0.094	0.138	0.162	0.105	0.1057	0.0987	0.0998	0.0924
	R^2	0.9987	0.9985	0.9981	0.9988	0.9987	0.9988	0.9988	0.9988

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, $${\beta }_{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 approximate analyses by the meta-model, respectively.

Element groups No.	DGSA		BGSA
Cross-section	$$\rho (\%)$$	Cross-section	$$\rho (\%)$$
C1	65 $$\times$$ 65	2.49	65 $$\times$$ 65	2.90
C2	60 $$\times$$ 60	2.50	60 $$\times$$ 60	3.08
C3	55 $$\times$$ 55	2.53	55 $$\times$$ 55	2.87
C4	50 $$\times$$ 50	2.49	50 $$\times$$ 50	2.79
C5	45 $$\times$$ 45	2.27	45 $$\times$$ 45	2.68
C6	40 $$\times$$ 40	2.29	40 $$\times$$ 40	2.89
B1	55 $$\times$$ 45	1.88	55 $$\times$$ 45	2.18
B2	55 $$\times$$ 40	1.80	55 $$\times$$ 40	1.98
B3	50 $$\times$$ 45	1.82	50 $$\times$$ 45	2.01
B4	50 $$\times$$ 40	1.87	50 $$\times$$ 40	1.98
B5	45 $$\times$$ 40	1.68	45 $$\times$$ 40	1.79
B6	45 $$\times$$ 35	1.66	45 $$\times$$ 35	1.85
CIC (Euro)	5571		5736
CRC (Euro)	1004		987
CTOT (Euro)	6575		6723
Iterations	83		123
$${\beta }_{annual}^{Operational}$$	1.4985		1.7101
$${\beta }_{annual}^{Lifesafety}$$	2.4953		2.7996
$${\beta }_{annual}^{Collapse}$$	2.8514		3.1837

Table 11.

Optimum designs obtained by DGSA and BGSA

Figure 10.

Convergence histories of the best solution of DGSA and BGSA

10. Conclusions

In general, the optimum design of structures depends on a number of parameters that are inherently uncertain. Reliability-based design optimization (RBDO) has been employed as the only method that assesses the influence of uncertain parameters and balance both cost and safety of structures. To account for all necessary uncertain and random parameters in RBDO of RC structures and to achieve the realistic optimum design of RC structures, the uncertain material properties of soil and structure, as well as the characteristics of ground motions should be considered as random parameters. Furthermore, the realistic seismic responses of RC structures can be account by consideration of soil-structure interaction (SSI) effects. In this work, a new discrete gravitational search algorithm (DGSA) and a new meta-modeling framework were incorporated for RBDO of RC structures with Performance-Based Design (PBD) under seismic loading. The objective of the RBDO problem was to minimize the total cost whereas the deterministic constraints and the system reliability index corresponding to each of the performance levels should not exceed a specified target. Based on this study, the following conclusions can be drawn:

• To reduce the computational effort and computational cost of the Monte-Carlo Simulation (MCS) method, a new meta-model based on a wavelet weighted least squares support vector machine (WWLS-SVM) and gravitational search algorithm (GSA) was utilized in the RBDO procedure. Therefore, the proposed meta-model, as a substitute for the nonlinear dynamic analysis of SSI system, can estimate the reliability index through MCS with a small computational cost.

• The WWLS-SVM and kernel parameters were simultaneously optimized in the proposed meta-model in order to improve performance generality of WWLS-SVM. Numerical results of training and testing the meta-model indicated that performance generality of the meta-model was higher in comparison to WLS-SVM. Hence, the proposed meta-model can predict the nonlinear dynamic analysis of SSI system in terms of accuracy and flexibility.

• The proposed DGSA was presented based on the standard GSA with passive congregation. The passive congregation strategy can be considered as perturbations operator in the optimization procedure. Therefore, the presented DGSA using the passive congregation can transfer information among agents avoiding local minima. Furthermore, the coding and encoding of the position of agents as a time consuming process is omitted in DGSA. To eliminate this drawback, the position of agents was calculated as the integer value. The optimum designs obtained by DGSA were compared with those produced by BGSA model. Numerical examples showed that the proposed DGSA can converge and reach the optimum design more quickly than the BGSA model.

Future extension of current research could include reducing the computations involved in the PBD by replacing MCS with the response surface method or the importance sampling technique. The constraints imposed on the objective function could be also treated as random quantities (see [68]).

Acknowledgments

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