Stochastic Emulation using Generalized Stratified Sampling for Performance-Based Risk Optimization of Structures

Stochastic Emulation using Generalized Stratified Sampling for Performance-Based Risk Optimization of Structures Isabela D. Rodriguesa , Seymour M.J. Spenceb,∗, Henrique M. Kroetzc , André T. Becka a

Department of Structural Engineering, São Carlos School of Engineering, University of São Paulo, São Carlos, 13566-590, São Paulo, Brazil b Department of Civil and Environmental Engineering, University of Michigan, Ann Arbor, 48109, Michigan, USA c Center for Marine Studies, Federal University of Paraná, Pontal do Paraná, 13566-590, Paraná, Brazil

arXiv:2608.05006v1 [cs.LG] 5 Aug 2026

Abstract Metamodels are instrumental in reducing the computational burden associated with nested reliability analyses and optimization loops in Performance-Based Risk Optimization (PBRO) of structures under stochastic loads. In this context, stochastic emulators are particularly useful because they approximate response distributions while accounting for the intrinsic stochasticity of the simulator. Among these methods, Stochastic Polynomial Chaos Expansion (SPCE) is especially attractive because it does not require replications of nonlinear analyses at fixed input conditions. However, SPCE may present limitations in accurately representing extreme responses in the tails of structural response distributions. To address this limitation, this study proposes a framework that combines Generalized Stratified Sampling (GSS) with SPCE. The GSS scheme partitions the input space into strata according to the intensity of the hazard, improving the representation of extreme responses, while independent SPCE emulators are trained within each stratum. The conditional exceedance probabilities estimated in each stratum are then recombined using the total probability theorem to evaluate the probabilistic constraints in the optimization problem. The proposed GSS-SPCE framework is applied to the optimal design of the cross-sectional areas of buckling-restrained braces in a two-story steel building. The objective is to minimize the initial construction cost while satisfying prescribed probabilistic performance constraints. In the case study, the strata are defined based on the spectral acceleration associated with the seismic excitation, ∗

Corresponding author Email address: smjs@umich.edu (Seymour M.J. Spence)

Preprint submitted to Elsevier

August 6, 2026

while the design variables are included as additional inputs to the SPCE emulators. This allows the response distributions associated with each design candidate to be efficiently and accurately estimated. Results show that the proposed framework accurately estimates structural response distributions, including their tail regions, while substantially reducing the number of nonlinear model evaluations required for PBRO. Keywords: Stochastic Emulators, Generalized Stratified Sampling, Performance-Based Risk Optimization, Uncertainty quantification, Structural reliability 1. Introduction The quantification of structural performance in terms of risk and reliability has become increasingly common in the assessment and design of structures exposed to natural hazards, including earthquakes [1] and extreme winds [2]. By explicitly accounting for uncertainties in hazard intensity and structural parameters, risk- and reliability-based performance assessments provide information that helps designers and stakeholders make more informed decisions regarding structural design, assessment, and management, with the aim of reducing economic losses and other adverse consequences. In this context, structural performance assessment can be integrated naturally into optimization frameworks, in which design alternatives are systematically explored to achieve an appropriate balance among structural safety, functionality, and economy. This class of problems has been referred to in the literature as Performance-Based Design Optimization (PBDO) [3]; in the present study, following Rodrigues et al. [4], the term Performance-Based Risk Optimization (PBRO) is adopted. The PBRO framework can be interpreted within the broader context of Reliability-Based Design Optimization (RBDO) and Risk Optimization (RO) [5]. In both formulations, reliability analysis is nested within the optimization and must be repeated throughout the search for the optimum. This nested structure introduces a significant computational burden, which may become prohibitive for high-dimensional problems or complex structural systems [6, 7, 8, 9, 10]. Among the various alternatives proposed to address this issue, metamodels (or surrogate models) have shown significant promise in approximating the input–output relationships of high-fidelity numerical models at substantially reduced computational cost [11, 12, 13, 14, 15, 16, 17]. 2

Although metamodels provide an efficient data-driven approach, their development for PBRO problems is further complicated by the stochastic nature of seismic and wind excitations. The estimation of performance metrics typically relies on repeated nonlinear simulations to obtain structural response quantities, commonly referred to as Engineering Demand Parameters (EDPs), such as displacements and accelerations. Because the response depends on the random process representing the excitation, its evaluation becomes particularly demanding within an optimization framework, where each iteration requires response estimates for the current design candidate. Representative PBRO studies using metamodels include Gidaris and Taflanidis [18], who incorporated Kriging into an optimization framework to approximate repair costs and reference velocities from stochastic simulations as functions of the design variables. Other studies have approximated the conditional distributions of EDPs. Gidaris et al. [19] developed a Kriging metamodel to approximate the median and standard deviation of the logarithmic structural response for use in fragility and risk calculations. An alternative formulation was later proposed by Kyprioti and Taflanidis [20], in which the EDP distribution was approximated directly using Kriging, thereby reducing the number of replications of the nonlinear simulation. Kim et al. [21] approximated the moments of the logarithmic structural response using a heteroscedastic Gaussian process metamodel within an active-learning-based RBDO framework. More recently, Rodrigues et al. [4] proposed a noisy Kriging framework trained to predict the median, dispersion, and correlation of the EDP distribution, while incorporating statistical uncertainty arising from a limited number of ground motions through a noise term. Further methodologies for RBDO problems under stochastic excitation were reviewed by Jerez et al. [22]. Despite the efficiency of the aforementioned metamodeling frameworks, their application still faces challenges related to the direct propagation of high-dimensional uncertainties through nonlinear simulations. These challenges have motivated the recent development of stochastic emulators, which aim to approximate entire response distributions rather than deterministic model outputs, with applications involving seismic [23, 24, 25, 26] and wind excitations [27, 28, 29]. In this context, a nonlinear structural model subjected to stochastic excitation can be interpreted as a stochastic simulator: for the same set of input parameters, different realizations of the stochastic excitation may lead to different structural response 3

values. As a result, the structural response is a random output for fixed inputs, thereby accounting for the variability induced by the stochastic excitation. Recent formulations of stochastic emulators include, among others, Generalized Lambda Models (GLMs) [30, 31] and Stochastic Polynomial Chaos Expansions (SPCEs) [32]. In this paper, SPCE is employed as the stochastic emulator, following its recent application to reliability analysis problems [33, 34]. SPCE introduces an artificial latent variable to represent the intrinsic stochasticity of the simulator, together with a noise term that improves numerical stability. This formulation allows the emulator to approximate general, nonparametric conditional response distributions without requiring replications of the original simulator at the same input point [32]. More recently, Moustapha and Sudret [35] incorporated SPCE into an RBDO framework to approximate conditional response distributions and compared its performance with that of Kriging metamodels. However, as pointed out by Zhu et al. [23], SPCE still faces limitations in seismic fragility analysis, particularly when failure data are scarce at large damage thresholds, because it is not specifically designed to accurately represent the tail of the EDP distribution. To address this limitation, Zhu et al. [23] suggested using advanced sampling schemes to better represent extreme events, and this direction is explored in the present work. Advanced sampling schemes, such as importance sampling [36, 37] and subset simulation [38, 39], have been widely studied and applied in reliability analysis to reduce the large number of model evaluations required to estimate small failure probabilities using Monte Carlo methods [40, 41, 42, 43, 44]. These schemes generate samples that more effectively explore the failure region, thereby increasing the proportion of realizations that contribute to the failure probability estimate. Several studies have also combined advanced sampling schemes with metamodels to further reduce the computational burden associated with reliability evaluations involving very small failure probabilities [45, 46, 47, 48]. With this background, this work combines the Generalized Stratified Sampling (GSS) scheme, formally introduced by Arunachalam and Spence [49], with SPCE to improve the approximation of the tails of conditional EDP distributions within a Performance-Based Engineering (PBE) framework. The proposed framework, referred to as GSS-SPCE, is subsequently used to solve a PBRO problem. In this framework, GSS partitions the input space 4

according to hazard intensity, improving the representation of rare events, while separate SPCE emulators are trained within each stratum with the design variables included as inputs. The main contribution of this work is the integration of an advanced sampling scheme with stochastic emulation to improve tail representation and enable the efficient estimation of probabilistic constraints within PBRO. The remainder of the paper is organized as follows. Section 2 presents the problem definition, including the PBE and PBRO concepts adopted in this work. Section 3 provides background on the stochastic simulation, SPCE, and GSS methods underlying the proposed framework. Section 4 presents the proposed GSS-SPCE framework and its implementation within PBRO. Section 5 demonstrates the framework through a case study involving the optimal design of the cross-sectional areas of Buckling-Restrained Braces (BRBs) in a twostory steel building, with the objective of minimizing the initial construction cost while satisfying prescribed probabilistic performance constraints. 2. Problem definition 2.1. Assessment of structural performance The assessment of building performance generally requires the characterization of a structural response quantity, represented by the random variable Y , under uncertainties in the system parameters and load excitation. For a given response level y, the exceedance probability is defined as: p (y) = P (Y > y) = 1 − FY (y)

(1)

where p (y) is the probability that the structural response Y exceeds the response level y, and FY is the cumulative distribution function (CDF) of Y . Structural performance can then be assessed by evaluating this probability at a prescribed critical value ycr , which gives rise to the design constraint: p (ycr ) = P (Y > ycr ) ≤ pt

(2)

where pt is the acceptable exceedance probability specified by the designer. Therefore, the following probabilistic integral must be evaluated: p (ycr ) = 1 − FY (ycr ) = 1 − 5

Z

Ωycr

fX (x) dx

(3)

where X is a random vector defined over the input space DX ⊆ RnRV and represents the uncertainties associated with the system parameters and load excitation, nRV is the number of random variables, fX is the joint probability density function (PDF) of X, and Ωycr = {x : Y (x) ≤ ycr } is the domain in which the structural response does not exceed the critical value. When the stochastic nature of the load excitation is explicitly modeled, as in this work, the integral in Eq. (3) becomes high-dimensional and generally must be estimated using simulation methods [11]. For instance, if Monte Carlo simulation is considered, the following steps are required to obtain FY (ycr ), considering the structural response as a nonlinear function y = gNL (x): 1. Generate a realization x(i) of the uncertain vector X; 2. Generate a realization of the stochastic load using the appropriate components of x(i) ; 3. Evaluate the corresponding realization of the structural response, y (i) , using gNL ; 4. If the target total number of samples has been generated, proceed to the next step. Otherwise, return to Step 1; 5. Estimate FY (ycr ) from the samples of Y . In Step 5, FY (ycr ) can be estimated directly from the samples using a nonparametric approach or by fitting a parametric probability distribution. If the acceptable exceedance probability is small, a large number of system evaluations, nsys , must be performed. For example, under crude Monte Carlo simulation, if pt = 10−k , approximately nsys = 10k+2 model evaluations are required to estimate p (ycr ) with a coefficient of variation of 10%. For computationally intensive nonlinear models, this requirement is generally prohibitive [49]. 2.2. Optimization problem Considering the structural performance assessment defined in Section 2.1, the design of a system can be formulated as a PBRO problem. The optimization consists of selecting values for a set of parameters, denoted by d, that satisfy a set of performance constraints while minimizing the system cost C. The vector d may represent structural element sizes, while the performance constraints may be associated with the requirements defined in Eq. (2). The

6

optimization problem can be written as: Find d = {d1 , . . . , dNd }T to minimize C(d) subject to P (Yj (d, X) > ycr,j ) ≤ pt,j

for j = 1, . . . , Nc

(4)

d∈D X ∈ DX where d is an Nd -dimensional vector of design variables belonging to the design space D ⊆

RNd , Yj (d, X) is the jth structural response of interest, ycr,j is the corresponding critical value, pt,j is the associated acceptable exceedance probability, and Nc is the number of constraints. Solving the optimization problem is not straightforward due to the reliability analyses nested within the optimization loop. Since the constraints are implicit functions of d, their evaluation requires a full probabilistic assessment, as described in Section 2.1, for each design candidate. This computational burden limits the use of evolutionary optimizers, while the implicit dependence of the constraints on d also makes gradient-based optimization challenging due to the difficulty of evaluating constraint gradients. Although metamodeling schemes have been proposed to approximate the structural response, limitations remain when highdimensional uncertainties must be propagated directly through nonlinear simulations and when the relevant tail probabilities, as in Eq. (1), are very small. To address these limitations, this work proposes combining an advanced sampling scheme with a stochastic emulation approach to accelerate the evaluation of the probabilistic constraints for any given design d. The GSS-SPCE method uses SPCE to emulate the structural response distributions associated with the constraints, while GSS partitions the input space according to hazard intensity to improve the approximation of the response distribution tails. The proposed approach facilitates the efficient solution of the PBRO problem defined in Eq. (4) by allowing the probabilistic constraints to be evaluated rapidly without additional calls to the nonlinear structural model.

7

3. Background methods The implementation of the GSS-SPCE framework requires the structural performance assessment to be formulated as a stochastic simulation problem. This formulation enables the structural response distribution to be emulated using SPCE, while the input space is partitioned through the GSS scheme. The background formulations of these methods are presented in the following sections. 3.1. Stochastic simulators To formulate the nonlinear structural response model y = gNL (x) as a stochastic simulator, Yx = gNL (x; ω) ≡ Ms (x, ω), the uncertainty associated with stochastic excitation can be separated into two main components. The first component is described by parametric explanatory variables that characterize the hazard intensity, such as seismic spectral acceleration or wind velocity. These variables are commonly referred to as intensity measures (IMs) and are explicitly included in the input vector x. The second component is associated with the complex physical mechanisms involved in the generation and propagation of the stochastic excitation. In seismic applications, this component is commonly referred to as record-to-record (RTR) variability. It represents the different possible load time histories that may occur for the same IM value, limiting the extent to which the excitation can be completely described by the selected parametric explanatory variables [25, 50]. This separation naturally leads to the formulation of the structural response as a stochastic simulator. While a deterministic simulator produces a unique scalar- or vector-valued output for a fixed input vector, a stochastic simulator produces random outputs even when the explicit input parameters are fixed. In this context, the output randomness represents the intrinsic variability of the simulator, which may account, for instance, for the effects of RTR variability. A stochastic simulator can be written as: Ms : DX × Ω → R

(5)

(x, ω) 7→ Ms (x, ω) where Ω is the sample space associated with the internal stochasticity of the simulator and ω ∈ Ω. Therefore, for a fixed input vector x, different realizations of ω may lead to different 8

values of the structural response. For a fixed x ∈ DX , the response Yx = Ms (x, ·) : Ω → R is therefore a random variable. Accordingly, when the input vector is fixed at x(0) , each simulator run corresponds to a
different event ω (j) and produces a realization Ms x(0) , ω (j) .

This behavior is common in structural engineering applications involving stochastic loads.

For ground motions, the response may remain random even when event parameters such as magnitude, source-to-site distance, and site conditions are fixed, due to RTR variability [23, 25]. Similarly, under wind loads, the structural response may vary due to the turbulent component of the excitation, even for fixed IM values [34, 28]. 3.2. Stochastic Polynomial Chaos Expansion SPCE is a stochastic emulator designed to reproduce the conditional distribution of the response Y of a stochastic simulator given the input vector X, that is, Yx ≡ Y | X = x. Originally introduced by Zhu and Sudret [32], SPCE extends the framework of Polynomial Chaos Expansions (PCEs) to represent the intrinsic stochasticity of the original simulator Ms . PCEs constitute a well-established and widely used metamodeling technique for deterministic simulators. Considering the random input vector X ∈ RnRV described by the joint PDF fX , a deterministic model output y = Md (x) can be approximated by an orthogonal polynomial expansion, expressed as: Md (x) ≈

X

cα ψα (x)

(6)

α∈A

where ψα denotes a multivariate polynomial basis function that is orthonormal with respect to the joint probability density function fX , α = [α1 , α2 , . . . , αnRV ] is the corresponding multiindex, and cα is the associated expansion coefficient to be estimated. The set A ⊂ NnRV is a finite set of multi-indices obtained by applying a truncation rule. The standard truncation scheme consists of selecting all polynomials whose total degree is less than or equal to a n o P RV prescribed value pdg , such that Apdg ,nRV = α ∈ NnRV : nj=1 αj ≤ pdg .

To construct the SPCE, a latent random variable Z is introduced to represent the intrinsic

stochasticity of the stochastic simulator. Let FY |X (y | x) denote the conditional CDF of the 9

response for a fixed input vector x. Using the probability integral transform and assuming a prescribed distribution for the latent variable Z, the conditional response can be expressed as Yx = FY−1|X (FZ (Z) | x), where FZ is the CDF of Z. Assuming that Y has finite variance,

FY−1|X (FZ (Z) | x) can be represented by a PCE in the augmented input space (X, Z). By adding a noise term ϵ ∼ N (0, σ 2 ) to improve numerical stability, the SPCE can be written as: Yx ≈ Yex =

X

cα ψα (x, Z) + ϵ

(7)

α∈A

The SPCE defined in Eq. (7) is constructed by estimating the coefficients cα and the standard deviation of the noise term, σ. In this work, the calibration follows the procedure described in [32, 51]. This procedure also includes the adaptive selection of the latent-variable distribution and the truncation scheme A ⊂ NnRV +1 . For this purpose, a Design of Experiments (DoE) is defined as:

  XS = x(1) , . . . , x(S) 

 YS = Ms x(1) , ω (1) , . . . , Ms x(S) , ω (S)

(8)

where the matrix XS of support points contains S realizations of X sampled in DX , and YS contains the corresponding stochastic simulator evaluations. Each support point is evaluated only once. As discussed in the Introduction, SPCE is particularly suitable for the present problem because it emulates the conditional response distribution without imposing a predefined parametric form and without requiring replications of the original simulator at the same input point. Once the SPCE is built, its realizations can be used to estimate response distributions and exceedance probabilities, as defined in Eq. (3). Nevertheless, the estimation of small exceedance probabilities requires an accurate representation of the distribution tails, which may be challenging when rare response levels are poorly represented in the training set [23]. For this reason, SPCE is combined in this work with the GSS scheme. The GSS procedure is used to generate stratum-specific DoE samples and the corresponding model evaluations, which are then used to train independent SPCEs within each stratum.

10

3.3. Generalized Stratified Sampling scheme The GSS scheme adopted in this work follows the methodologies proposed by Arunachalam and Spence [42] and Arunachalam and Spence [49] for reliability assessment problems involving natural hazards and the estimation of multiple small failure probabilities. In the present study, GSS is used to generate the DoE samples within each stratum required for the calibration of SPCE. Accordingly, only Phase I of the generalized stratified sampling framework proposed by Arunachalam and Spence [49] is considered, since the objective is to construct stratum-specific samples rather than to estimate stratum-specific failure probabilities directly. The basic idea is to partition the input space DX into a collection of nS disjoint strata: S {Si }ni=1

such that:

nS [

i=1

Si = DX

(9)

and Si ∩ Sj = ∅ for i ̸= j

(10)

The partition is defined in terms of a scalar stratification variable V, which may be either an input random variable or the output of an auxiliary computational model H that depends on a subset of the input uncertainties. In the latter case, V is evaluated as V = H(τ ), where τ denotes the subset of input components required to compute the stratification variable. The use of H is particularly convenient when its computational cost is much lower than that of the nonlinear structural model gNL , allowing a large number n̂ of realizations of V to be generated. After selecting the number of strata nS and the probability constant pGSS , the stratum boundaries are defined adaptively from the empirical distribution of V. Following Arunachalam and Spence [49], values of pGSS commonly lie in the range [0.1, 0.3]. Equivalently, if the n̂ realizations of V are sorted in increasing order, Vi is taken as the (1 − piGSS ) n̂th ordered value for i = 1, . . . , nS − 1. Setting V0 = −∞ and VnS = ∞, the strata are induced in the input space by: Si = {τ : V ∈ (Vi−1 , Vi ]} ,

i = 1, . . . , nS

(11)

Once the strata are defined, samples are generated within each stratum to form the DoE. 11

The components of τ are sampled conditionally on the corresponding stratum, while the remaining input variables are sampled according to their prescribed distributions. When the stratification variable and the stratum boundaries are appropriately defined, GSS provides a more efficient representation of rare-event regions than standard Monte Carlo sampling using the same number of model evaluations. This motivates its use in the present work to improve the coverage of the input space used to calibrate the SPCE, particularly in regions associated with the tails of the structural response distribution. 4. Proposed framework 4.1. Overview The GSS-SPCE framework proposed in this work is used to emulate the distribution of the structural response generated by the stochastic simulator associated with the nonlinear structural model and the stochastic load excitation. Instead of approximating only a deterministic response value for a given input vector x, the proposed framework generates realizations of the conditional emulated response Yex following the SPCE formulation defined

in Eq. (7). By evaluating the trained stratum-specific SPCEs over samples of X, the resulting realizations are combined to estimate the unconditional CDF of the structural response, denoted by FYe , and the corresponding exceedance probability: 

pe (y) = P Ye > y



(12)

as an approximation of the performance metric defined in Eq. (1). Once the stochastic emulators are trained and their accuracy is assessed in representing the unconditional response distributions, they are used within the optimization framework defined in Eq. (4) to estimate the probabilistic constraints for each design candidate. The following subsections describe the construction of the GSS-SPCE in the augmented input space and its subsequent use in the optimization problem. 4.2. Construction of the SPCE in the augmented input space The complete set of uncertainties involved in the structural performance assessment can be conceptually divided into model-related uncertainties, Xm , load-intensity uncertainties, Xh , 12

and the stochastic variability of the load excitation, Xw . These sources may be represented  T by Ξ = XTm , XTh , XTw .

As defined in Eq. (7), the SPCE is formulated in the augmented input space (X, Z), T  where X = XTm , XTh contains the variables explicitly sampled and provided as inputs to the emulator. The latent random variable Z represents the remaining stochastic variability

associated with Xw , which is not explicitly included in X. To avoid repeated nonlinear-model evaluations when estimating the probabilistic constraints within the optimization, the design variables d are also included as emulator inputs. Therefore, the GSS-SPCE is constructed in  T the fully augmented input space (U, Z), where U = dT , XTm , XTh . For the purpose of DoE

generation and SPCE training, the design variables are sampled from uniform distributions over the design domain D. As a result, the explicit input vector of the stochastic emulator  T is U, while Z remains a latent variable. For a realization u = dT , xT , the structural response estimated by the SPCE is conditioned on the design variables and the explicitly sampled random inputs and is denoted by Yeu = Ye (d, x). For a fixed design vector d(ℓ) ∈ D

within the optimization process, random realizations of X are propagated through the trained
SPCEs to generate samples of Ye d(ℓ) , X . Alternative formulations are also possible. For instance, part of the uncertainty associated

with the explicitly sampled random inputs X could instead be represented in the latent space through Z. In the limiting case, Z would represent all stochastic variability associated with Ξ, leaving only the design variables as explicit emulator inputs. The resulting mapping from the design variables to the structural response would still constitute a stochastic simulator, with the latent variable accounting for all output randomness. 4.3. Stratified scheme for constructing the GSS-SPCE Following the GSS formulation described in Section 3.3, the stratification variable is nS evaluated from the subset τ of the hazard-related variables Xh . Once the strata {Si }i=1 are

defined, samples of Xh are selected within each stratum by conditioning the components τ on the corresponding stratum. If P (SnS ) ≈ 10−k , approximately 10k+2 evaluations of H are required to obtain an expected 102 samples in the last stratum, corresponding approximately to a coefficient of variation of 10% in the estimated stratum probability [49].

13

The remaining components of Xh , if any, are sampled according to their prescribed distributions. These samples are then combined with uniformly distributed realizations of the design variables d and realizations of the model-related uncertainties Xm to form the stratumspecific input vectors U(i) used to train the SPCEs. This leads to a stratum-specific emulated response, denoted by Ye (i) . The exceedance probability in Eq. (12) can then be expressed within each stratum as:


(i) (i) e > y | Si pe (y) = P Y U 

(i)

(13)

Following the law of total probability, the unconditional exceedance probability is obtained by recombining the stratum-specific conditional probabilities as: pe (y) =

nS X i=1

pe(i) (y) P (Si )

(14)

where P (Si ) denotes the probability of occurrence of the ith stratum. This recombination S is valid because the strata {Si }ni=1 form a mutually exclusive and collectively exhaustive

partition of the input space, as defined in Eq. (10). In the general SPCE formulation, the DoE was denoted by XS in Eq. (8) and defined in terms of the input vector X. In the proposed GSS-SPCE framework, the corresponding DoE (i)

is defined in the augmented input space and denoted by US for each stratum Si as:   (i) US = u(1,i) , . . . , u(S,i)  

 (j,i) T T (j,i) (j,i) T (j,i) T u = d , xm , xh , (j,i)

where j indexes the support points and xh

(15) j = 1, . . . , S

denotes a realization of the hazard-related

variables associated with stratum Si . 4.4. Algorithm The overall procedure of the proposed framework is summarized in the flowchart shown in Fig. 1. The following steps are implemented to partition the input space, select samples within each stratum, train the SPCE emulators, and estimate the exceedance probabilities within the optimization loop:

 T 1. Define the augmented input vector U = dT , XTm , XTh . 14

GSS-SPCE Define input space U = [d,Xm,Xh]T Define t, nS, pGSS and 𝒱=ℋ(t)

Sample X and assign samples to strata Si

Perform (S x nS) analysis to obtain 𝒴S(i)=gNL(𝒰S(i))

Define Si={t: 𝒱 ∈ (𝒱i-1,𝒱i]}

~

Yes

i =nS-1

ℓ=1

Combine Xh(i) with random samples of d and Xm to form 𝒰S(i)

i=1

No

Design option d(ℓ)

Select S samples Xh(i) of t from each stratum Si

^ MC Generate n samples of t

i = i+1

Optimization

Evaluate trained SPCEs ~ Y (i) ~

Train nS SPCEs Y (i) with 𝒰S(i) and 𝒴S(i)

Estimate p(i) (ycr|d(ℓ) ) ~ and p (ycr|d(ℓ) )

ℓ = ℓ+1 No

Minimum cost?

Yes

Done

Figure 1: Flowchart of the proposed methodology.

2. Define the subset τ of the hazard-related variables Xh used to evaluate the stratification variable V = H(τ ), together with the number of strata nS and the probability constant pGSS ∈ [0.1, 0.3]. Generate a sufficiently large number n̂ of realizations of τ to define the stratum boundaries adaptively, as described in Section 3.3. This procedure partitions S the input space into nS disjoint strata {Si }ni=1 .

3. For each stratum Si , select S realizations of the hazard-related variables, denoted by (j,i)

{xh }Sj=1 . Combine these samples with realizations of the design variables {d(j,i) }Sj=1 (i)

(j,i)

and the model-related uncertainties {xm }Sj=1 to form the stratum-specific DoE US in the augmented input space. 4. Perform nS × S evaluations of the stochastic simulator Ms to obtain the corresponding (i)

response vectors YS for i = 1, . . . , nS .

5. Train one SPCE for each stratum. Each emulator Ye (i) is calibrated using the stratum(i)

(i)

specific DoE US and the corresponding response vector YS .

6. Within the optimization loop, for each design candidate d(ℓ) , generate a sufficiently large number of realizations of X. For each realization, identify the corresponding
stratum Si from the value of τ and evaluate the associated emulator Ye (i) d(ℓ) , xm , xh .
Estimate the exceedance probability conditioned on each stratum, pe(i) y | d(ℓ) , using 15

Eq. (13), and recombine the stratum-specific probabilities using Eq. (14) to obtain the
unconditional exceedance probability pe y | d(ℓ) .

7. Repeat Step 6 for the successive design candidates generated by the optimization algorithm until its termination criteria are satisfied and the optimal design d⋆ is obtained by minimizing the cost function while satisfying the probabilistic constraints. 5. Case study 5.1. Overview

The proposed framework is demonstrated through a case study involving a two-story buckling-restrained braced-frame (BRBF) archetype building designed in accordance with [52, 53]. Buckling-restrained braces (BRBs) are concentric structural elements designed to prevent global brace buckling and concentrate inelastic deformations within a yielding steel core [54]. This mechanism enhances lateral stiffness and energy dissipation capacity, thereby improving the seismic performance of the building by limiting damage to the surrounding structural members. The selected archetype is documented in the National Institute of Standards and Technology (NIST) GCR 10-917-8 report [55] and is assumed to be located in urban California. Plan and elevation views of the building are presented in Figs. 2a and 2b, respectively.

4.00 m

38.10 m

5 bays @ 7.62 m

4.00 m

38.10 m

7.62 m

5 bays @ 7.62 m

(a)

(b)

Figure 2: Case-study structure: (a) typical plan view of the building; (b) elevation of the steel BRBF.

16

Table 1: Member sizes for the BRBF.

Story

Beam size

Column size

1

W21×73

W14×82

2

W21×50

W14×82

The building performance in this case study is assessed with respect to its response under seismic excitation. The seismic inputs are generated using the stochastic ground-motion model described in Appendix A. The optimization problem considers the cross-sectional areas of the BRBs on the first and second stories, denoted by d1 and d2 , respectively, as design variables. The objective is to determine the design vector d that minimizes the construction cost of the frame while satisfying prescribed probabilistic performance constraints. The PBRO problem is formulated as: Find d = [d1 , d2 ]T to minimize C(d)

  subject to pe1 (ycr | d) = P Ye1 (d, X) > ycr ≤ pt   e pe2 (ycr | d) = P Y2 (d, X) > ycr ≤ pt

(16)

20 ≤ d1 ≤ 70 cm2 20 ≤ d2 ≤ 70 cm2

where Ye1 (d, X) and Ye2 (d, X) denote the emulated first- and second-story interstory drifts,

respectively, while pe1 (ycr | d) and pe2 (ycr | d) are their corresponding exceedance probabilities.

The optimal design obtained by solving Eq. (16) is denoted by d⋆ .

For demonstration purposes, the critical interstory drift threshold and acceptable ex-

ceedance probability are set to ycr = 15% and pt = 4 × 10−3 , respectively. The building is located on Site Class C and has a typical story height of 4.0 m. The beams and columns are made of ASTM A992 Grade 50 steel, while the BRBs are made of ASTM A36 steel. Additional design details can be found in [55]. Table 1 summarizes the geometric specifications of the structural elements. The construction cost of the steel frame, denoted by C(d) in Eq. (16), is estimated as: C(d) = Cs (Wb + Wc ) + CBRB WBRB 17

(17)

where Cs is the cost per unit weight of structural steel, CBRB is the cost per unit weight of the BRBs, and Wb , Wc , and WBRB are the total weights of the beams, columns, and BRBs, respectively. The unit costs adopted in this study are Cs = $4.72/kg [56] and CBRB = $8.99/kg [57]. The nominal weights of the beams and columns are obtained from the American Institute of Steel Construction (AISC) database, version 16.0 [58]. The total weight of the BRBs P BRB is calculated as WBRB = ρ N b=1 Ab Lb , where NBRB is the total number of BRBs, and Ab

and Lb are the cross-sectional area and length of the bth BRB, respectively. The material density is taken as ρ = 7850 kg/m3 for ASTM A36 steel. 5.2. Finite element model of the structural system A two-dimensional nonlinear dynamic finite element (FE) model of the two-story BRBF archetype is developed in OpenSees [59]. The model represents one exterior braced frame in the direction of analysis. Only the bare steel components of the BRBF are explicitly

modeled; the composite floor slabs and gravity framing system are not directly represented. The floor slabs are assumed to act as rigid diaphragms, such that the seismic mass tributary to the analyzed frame is lumped at the floor levels. Beams and columns are modeled using elastic beam-column elements with concentrated flexural springs at their ends to represent the plastic hinge regions. The nonlinear response of these springs is defined through the moment–rotation relationships for each member section. The cyclic behavior of the steel springs is represented using the Steel02 material model available in OpenSees, originally based on the Giuffrè–Menegotto–Pinto formulation with isotropic hardening parameters [60]. The rotational springs are implemented using zeroLength elements. A Young’s modulus of 200 GPa and a yield strength of 345 MPa are adopted for the beams and columns, together with a post-yield strain-hardening ratio of 1%. The elastic portions of the beams and columns between the end springs are represented by elasticBeamColumn elements using the corresponding cross-sectional properties. The BRBs are modeled using twoNodeLink elements with axial response only. The cyclic behavior of the BRB steel core is also represented using the Steel02 material model. A yield strength of 290 MPa and a post-yield strain-hardening ratio of 2% are adopted for the BRB core.

18

Rayleigh damping is assigned to reproduce a damping ratio of 2% in the first and third vibration modes. Geometric nonlinearities are considered by explicitly including P –∆ effects in the beam and column elements. Gravity loads are applied at the floor levels consistently with the tributary vertical loads of the exterior braced bay, while the seismic masses are assigned according to the tributary floor area associated with the analyzed frame under the rigid-diaphragm assumption. The gravity-load and seismic-mass assumptions are defined based on the archetype design information provided in [55]. 5.3. GSS-SPCE implementation 5.3.1. Input space and probability models To emulate the structural response Ye and estimate the exceedance probability defined

in Eq. (12), the input space must first be defined. In this example, the input vector of the GSS-SPCE is given by:

 T U = dT , XTh

(18)

where d = [d1 , d2 ]T is the vector of design variables defined in Eq. (16), and Xh = [Mw , r]T is the vector of hazard-related random variables. Here, Mw is the earthquake moment magnitude and r is the epicentral distance. The remaining stochastic variability associated with the ground-motion excitation is represented implicitly through the latent variable Z of the SPCE. For the construction of the DoE, the design variables are sampled uniformly within their admissible ranges, while the hazard-related variables are sampled according to their prescribed probability distributions. The earthquake magnitude Mw is modeled using a truncated Gutenberg–Richter distribution [49], with PDF given by: fMw (Mw ) =

β exp [−β (Mw − Mw,min )] 1 − exp [−β (Mw,max − Mw,min )]

(19)

where Mw,min = 6, Mw,max = 8, and β = 0.9 ln (10). The epicentral distance r is modeled as a lognormal random variable with a median of 5 km and a coefficient of variation of 0.4. For the construction of the DoE, the BRB areas d1 and d2 are sampled independently from uniform distributions over their admissible ranges, such that d1 , d2 ∈ [20, 70] cm2 . The number of support points within each stratum is set 19

to S = 2000, with each support point corresponding to one realization of the explicit input (i)

vector U. These support points form the stratum-specific DoE: US . 5.3.2. Stratum-specific samples To construct the stratum-specific DoEs, the number of strata is set to nS = 5, and the probability constant is set to pGSS = 0.2. The stratification variable V is chosen as the spectral acceleration at the first-mode period, Sa (T̄ ), considering a damping ratio of 5%. This choice is based on Arunachalam and Spence [49], who observed an increasing trend in drift ratios with increasing spectral acceleration, thereby supporting the use of Sa (T̄ ) as an effective stratification variable. Accordingly, the auxiliary computational model H used to evaluate V consists of the stochastic ground-motion model, parameterized by the hazard-related variables, followed by the evaluation of the response of a linear oscillator. Thus, V = Sa (T̄ ) = H(τ ), where τ = Xh = [Mw , r]T in this case study. The reference period T̄ = 0.624 s is computed for the

design d1 = d2 = 45 cm2 , corresponding to the midpoint of each admissible design-variable range, and is used throughout the stratification process. As discussed in Section 4.3, a region with a probability of the order of 10−3 requires approximately 105 evaluations of H to obtain an expected 102 samples within that region. In this study, n̂ = 107 realizations of τ are generated using Monte Carlo simulation to provide a sufficiently large candidate pool for selecting S = 2000 support points within the final stratum. Because spectral acceleration is nonnegative, the lower boundary is set to V0 = 0. The resulting stratum boundaries are {V0 , V1 , V2 , V3 , V4 , V5 } = {0, 1.207, 2.223, 3.404, 4.710, ∞} in units of g. The estimated empirical spectral-acceleration hazard curve and the corresponding stratum boundaries are shown in Fig. 3a, assuming λMw 6 = 0.6 as the occurrence rate associated with the lower magnitude threshold Mw,min = 6. Figure 3b presents the stratum-specific scatter of the seismic hazard parameters Mw and r.

20

10-2

10-4

10-6

10-1

100

101

(a)

(b)

Figure 3: Stratum-specific samples: (a) estimated empirical spectral-acceleration hazard curve λ(V) and the corresponding stratum boundaries; (b) stratum-specific scatter of Mw and r.

5.4. Results The results are presented in the following subsections: (1) construction and validation of the GSS-SPCE over the augmented input space U; and (2) optimization results obtained using the trained GSS-SPCEs. 5.4.1. GSS-SPCE results From the candidate samples generated in Section 5.3.2, S = 2000 realizations of Mw and r are randomly selected within each stratum to serve as hazard-related support points for the construction of the SPCEs. The selected samples within each stratum are shown in Fig. 4. For each stratum, S support points of the design vector d are generated using Latin hypercube sampling and combined with the hazard-related support points to form the stratum(i)

specific DoE: US . The accuracy of the GSS-SPCE is first assessed within each stratum. For this purpose, an independent reference dataset containing 5000 realizations of the seismic hazard variables Mw and r is generated for each stratum. The design variables d1 and d2 are also sampled within their admissible ranges to form the corresponding input vectors U(i) . For each input vector, a stochastic ground-motion realization is generated and a nonlinear time-history analysis (NLTHA) is performed using the stochastic simulator Ms to obtain the (i)

reference response samples YR . Using the same input vectors U(i) , the trained GSS-SPCEs are evaluated within each 21

101

100 6

6.5

7

7.5

8

Figure 4: Selected samples of Mw and r for the S = 2000 support points within each stratum.

stratum to generate the corresponding emulated response samples, denoted by Ye (i) . The resulting conditional complementary cumulative distribution functions (CCDFs) are presented (i)

in Fig. 5 for Y1 and Fig. 6 for Y2 , together with the reference results YR and their associated confidence bounds. The 2000 structural responses used to train each stratum-specific (i)

GSS-SPCE are also shown and denoted by YT . The conditional results within each stratum are plotted according to their contribution to the recombined exceedance probability pe(y) defined in Eq. (14). Because the contribution

of the ith stratum is weighted by P (Si ), the conditional exceedance probability pe(i) (y) con-

tributes to the unconditional estimate through the product pe(i) (y)P (Si ). Therefore, to assess

the accuracy of the GSS-SPCE down to a prescribed unconditional probability level plim , the conditional CCDF within each stratum must be evaluated approximately down to: (i)

pelim =

plim P (Si )

(20)

The plotted range of each conditional CCDF is adjusted according to the probability of occurrence of the corresponding stratum. For the stratification adopted in this study, with pGSS = 0.2 and nS = 5, the stratum probabilities are {P (Si )}5i=1 = {0.80, 0.16, 0.032, 0.0064, 0.0016}. Thus, the conditional CCDFs in the tail strata do not need to be represented down to the same probability level as those in the more frequent strata because their smaller stratum weights already reduce their contributions to the unconditional exceedance probability. For instance, adopting plim = 10−3 for the 22

P Ỹ (2) > y | S2

P Ỹ (1) > y | S1

(1)

YR Ỹ

10−2

10

5

(2)

YR

10−1

(1)

10−3 0.0

(2)

YT

P Ỹ (3) > y | S3

(1)

YT

10−1

100

100

100

Ỹ (2)

(3)

YT

10−1

−2

10

5

y1 (%)

(3)

YR

Ỹ (3)

10

5

y1 (%)

15

P Ỹ (5) > y | S5

1

P Ỹ (4) > y | S4

1 0.8 0.6

10

y1 (%)

0.9

0.4

0.8

(4)

YT

(4) YR (4)

0.2

(5)

YR

0.7

Ỹ

5

(5)

YT

10

Ỹ (5)

15

5

y1 (%)

10

15

y1 (%)

Figure 5: Conditional exceedance probabilities for the first-story interstory drift Y1 within each stratum.

(1) YR (1)

Ỹ

10−2 10−3

10

(2) YR (2)

−1

Ỹ

10

10−2

2

4

0.0

2

y2 (%)

(3)

YT

(3)

YR

Ỹ (3) −1

4

0.0

2

y2 (%) (4)

P Ỹ (4) > y | S4

YT

1

1 0.8 0.6

5

y2 (%)

0.9

(4)

YR Ỹ

0.4

(5)

YT

P Ỹ (5) > y | S5

0.0

100

(2)

YT

P Ỹ (3) > y | S3

−1

P Ỹ (2) > y | S2

P Ỹ (1) > y | S1

10

(1)

YT

100

100

(4)

(5)

YR

Ỹ (5)

0.8 0.7

0.2 0.0

2

4

0.2

y2 (%)

0.5

0.8

y2 (%)

Figure 6: Conditional exceedance probabilities for the second-story interstory drift Y2 within each stratum.

23

100

100

10-1

10-1

10-2

10-2

10-3

10-3

10-4

0

10

20

30

40

10-4

50

(a) First-story interstory drift Y1

0

5

10

15

20

(b) Second-story interstory drift Y2

Figure 7: Unconditional exceedance probabilities of the first- and second-story interstory drifts. (i)

graphical assessment leads to conditional probability limits pelim of approximately 1.25×10−3 , 6.25 × 10−3 , 3.13 × 10−2 , 1.56 × 10−1 , and 6.25 × 10−1 for strata S1 through S5 , respectively.

The stratum-specific results show good agreement between the reference responses and

the samples generated by the GSS-SPCEs over the plotted probability ranges, including the tail regions represented within each stratum. These results support the subsequent use of the weighted recombination to estimate the unconditional exceedance probabilities. To assess the accuracy of the GSS-SPCE after recombination, an independent reference dataset containing 50,000 NLTHAs is generated using random realizations of the input vector  T U = dT , XTh . The corresponding reference response samples are denoted by YR . The same

explicit input realizations are then evaluated using the trained stratum-specific SPCEs. For

each input realization, the corresponding stratum is identified from the values of Mw and r, and the associated SPCE is used to generate one emulated response sample. The resulting emulated response samples are denoted by Ye . The empirical unconditional exceedance prob-

ability curves are compared in Fig. 7. The figure also includes the weighted recombination of the stratum-specific empirical distributions obtained from the 10,000 NLTHAs used to generate the training responses, denoted by YT . In addition, a lognormal distribution is fitted

to the reference response samples and plotted for each component of Y . The results indicate that the GSS-SPCEs generate structural response samples that accurately represent the reference distributions down to exceedance probabilities of approximately

24

10−3 . This probability level is consistent with the stratification process because the lower boundary of the final stratum shown in Fig. 3a corresponds to an exceedance rate close to 10−3 . If greater accuracy at lower probabilities is required, additional strata could be introduced to extend the range covered by the stratification. The GSS-SPCE result requires the 10,000 NLTHAs used for training, corresponding to one-fifth of the 50,000 NLTHAs used to construct the reference distribution. The results also indicate that the lognormal approximation, which is commonly adopted in performance-based engineering applications, begins to deviate from the reference response distributions at exceedance probabilities below approximately 10−1 . This limitation may become more pronounced when the distribution is fitted using fewer response samples because of the computational cost of the nonlinear structural analyses. To assess the accuracy of GSS-SPCEs trained using fewer support points within each stratum, an error metric focused on the upper tail of the response distribution is defined. Because the reliability assessment is governed by small exceedance probabilities, the discrepancy between the reference and emulated distributions is evaluated only above the 90% quantile of the reference distribution. The dimensionless tail error metric is defined as: h  i 2 E W2,q YR , Ỹ εtail = , q = 0.90 (21) Var [YR ] 2 denotes the squared Wasserstein distance restricted to the upper tail of the where W2,q

distributions, defined as: 2 W2,q (Ya , Yb ) =

1 1−q

Z 1 q

[Qa (u) − Qb (u)]2 du

(22)

where Qa and Qb are the quantile functions of the random variables Ya and Yb , respectively. The factor 1/(1 − q) normalizes the metric by the length of the tail interval, so that the error represents the average squared discrepancy between the two quantile functions over the upper tail. For this assessment, each trained GSS-SPCE is evaluated once for each input vector, consistent with the standard use of the emulator after training. The number of support points S within each stratum is varied, and a new set of stratum-specific SPCEs is trained for each value of S. To account for variability arising from the selection of the DoE and the 25

100

10-1

10-2

10-3

200

500

1000

1500

2000

Figure 8: Tail error εtail for Y1 as a function of the number of support points S used to train each stratumspecific emulator Ye (i) .

intrinsic stochasticity of the simulator, each scenario is repeated 20 times using independent DoEs. The resulting values of εtail are summarized using box plots in Fig. 8 for the firststory interstory drift Y1 . The dash-dotted line connects the mean errors obtained from the 20 repetitions for each value of S, while the horizontal line represents the error obtained by approximating the reference response distribution as lognormal. The results indicate that using S = 1500, or even S = 1000, does not materially degrade the tail accuracy relative to that obtained using S = 2000. This suggests that the number of support points used to train each stratum-specific emulator could be reduced by as much as one-half while maintaining an adequate level of accuracy. Such a reduction would further improve the computational efficiency of the proposed framework. 5.4.2. Optimization results The results obtained by solving the optimization problem defined in Eq. (16) are presented in this section. For each design candidate d(ℓ) , the probabilistic constraints associated with

pe1 ycr | d(ℓ) and pe2 ycr | d(ℓ) are evaluated using the trained stratum-specific SPCEs. This

process is repeated until the termination criteria of the optimization algorithm are satisfied and the optimal design d⋆ is obtained.

The optimization problem is solved using a genetic algorithm (GA) [61], which explores the design space through selection, crossover, and mutation operations. Although a gradientbased algorithm could also be considered, the probabilistic constraints are implicit functions 26

100

100

10-1

10-1

10-2

10-2

10-3

10-3

10-4

0

10

20

30

40

10-4

50

(a) First-story interstory drift Y1

0

10

20

30

40

50

(b) Second-story interstory drift Y2

Figure 9: Unconditional exceedance probabilities of the first- and second-story interstory drifts for d⋆ = T

[39.14, 20] cm2 .

of the design variables, as discussed in Section 2.2, and their gradients are not readily available. The GA is therefore adopted because it does not require constraint gradients and is suitable for optimization problems with potentially nonconvex constraint functions. A population size of 40 candidates is used. For each generation, 50,000 random realizations of the hazard-related variables Xh are generated and combined with each design h iT
(ℓ) (ℓ) T T candidate d to form the corresponding input vectors U = d , Xh . For each re-

alization, the corresponding stratum Si is identified from the values of Mw and r, and the

associated stratum-specific SPCEs are evaluated to generate realizations of the emulated re
(i) (i) sponses Yej , for j = 1, 2. The conditional exceedance probabilities pej ycr | d(ℓ) are then

estimated within each stratum and recombined according to Eq. (14) to obtain the uncondi
tional exceedance probabilities pej ycr | d(ℓ) .

The optimal design obtained is d⋆ = [39.14, 20]T cm2 . For this design, an independent

reference dataset containing 50,000 response realizations is generated through NLTHA using random realizations of the hazard-related variables Mw and r. The resulting reference responses are denoted by YR . The same 50,000 explicit input realizations are then evaluated using the trained stratum-specific SPCEs to generate the corresponding emulated responses, denoted by Ye . The resulting unconditional exceedance probability curves are presented in Fig. 9 for the first- and second-story interstory drifts, shown in Figs. 9a and 9b, respec-

tively. The results show good agreement between the GSS-SPCE estimates and the reference 27

distributions for the optimal design d⋆ down to exceedance probabilities of approximately 10−3 . This probability level is below the acceptable exceedance probability adopted in the optimization problem, pt = 4 × 10−3 , indicating that the probabilistic constraints are accurately evaluated over the probability range of interest. These results support the use of the GSS-SPCE approach for optimization problems involving acceptable exceedance probabilities within the range validated in this study, which would otherwise require a large number of evaluations of the nonlinear structural model. Regarding computational effort, the GA terminated after 870 evaluations of the objective function and probabilistic constraints. The termination criterion was satisfied when the average change in the fitness value fell below the prescribed function tolerance while the maximum constraint violation remained below the constraint tolerance of 10−5 . Within the proposed framework, the probabilistic constraints are evaluated using the trained GSS-SPCEs. Therefore, no additional evaluations of the nonlinear structural model gNL are required during the optimization loop. To provide an estimate of the computational savings achieved by the proposed framework, two hypothetical direct optimization strategies are considered. For this comparison, the same 870 design-candidate evaluations are assumed for all strategies. This assumption enables a consistent comparison of the required numbers of nonlinear model evaluations, although the convergence behavior of the optimization algorithm could differ when the probabilistic constraints are estimated directly through simulation. Without the emulators, each design candidate would require the direct propagation of samples through gNL to estimate the exceedance probabilities associated with the probabilistic constraints. If a direct MC-based optimization were performed using 50,000 samples per design candidate, the optimization would require approximately 870 × 50,000 = 4.35 × 107 evaluations of gNL . A direct GSSbased optimization would reduce this number by estimating the exceedance probabilities through stratified sampling. Ideally, the number of samples assigned to each stratum would be determined using an optimal allocation scheme, as proposed by Arunachalam and Spence [49]. Such an allocation depends on the conditional exceedance probabilities within the individual strata, which are not known a priori and may vary with the design candidate. Therefore, for the purpose of this illustrative comparison, 100 nonlinear model evaluations 28

Table 2: Estimated numbers of nonlinear model evaluations required by the different optimization strategies. Optimization approach

Evaluations per design candidate

Design-candidate evaluations

One-time training evaluations

Total evaluations of gNL

Total relative to GSS-SPCE

Direct MC

50,000

870

0

4.35 × 107

4350

5 × 100 = 500

870

0

4.35 × 105

43.5

870

1.00 × 104

1.00 × 104

1.0

GSS GSS-SPCE

0

Table 3: Comparison of the original NIST design and the optimized design obtained from the GSS-SPCE framework. Here, yej,pt is the interstory drift corresponding to an exceedance probability of pt for story j. Design

d1 (cm2 )

d2 (cm2 )

NIST design

58.06

38.71

Optimized design

39.14

20.00

pe1 (ycr )

2.67 × 10

−3

3.84 × 10−3

pe2 (ycr )

5.40 × 10

−5

8.76 × 10−4

ye1,pt (%) 13.22

14.69

ye2,pt (%) 3.14

7.56

C (d) ($)

2.16 × 104

1.94 × 104

are assumed within each of the five strata. This results in 5 × 100 = 500 evaluations of gNL per design candidate and approximately 870 × 500 = 4.35 × 105 evaluations over the optimization process. In contrast, constructing the GSS-SPCEs requires a one-time computational cost of S = 2000 support points within each of the nS = 5 strata, resulting in 5 × 2000 = 10,000 evaluations of gNL . No additional nonlinear model evaluations are required during the subsequent optimization. Under the assumptions adopted for this comparison, the GSS-SPCE framework reduces the number of nonlinear model evaluations by an estimated factor of 4.35 × 107 /10,000 = 4350 relative to direct MC-based optimization and by an estimated factor of 4.35 × 105 /10,000 = 43.5 relative to direct GSS-based optimization. The corresponding computational requirements are summarized in Table 2. It is also informative to compare the optimization results obtained in this study with the original design reported in NIST GCR 10-917-8 [55]. Because the design objectives and performance criteria adopted in the two procedures differ, this comparison is intended only to provide a reference point for interpreting the optimized design. As shown in Table 3, the optimized design has smaller BRB areas than the original NIST design while satisfying the probabilistic constraints within the tolerance adopted by the optimization algorithm. The optimized design also has a lower estimated initial construction cost under the cost model adopted in this study. These results indicate that the proposed PBRO formulation can identify a more economical design while maintaining the prescribed probabilistic seismic performance requirements considered in this study. 29

6. Conclusions This paper proposes a framework that combines an advanced stratified sampling scheme with stochastic emulation to estimate exceedance probabilities of structural responses under stochastic excitation, with particular emphasis on accurately representing distribution tails. The proposed GSS-SPCE framework is intended for use within PBRO problems in which the constraints are defined in terms of acceptable exceedance probabilities. GSS partitions the input space according to hazard intensity, while independent SPCEs are trained within the resulting strata to emulate the corresponding conditional response distributions. The stratum-specific conditional exceedance probabilities are then recombined using the law of total probability, providing an efficient means of evaluating the probabilistic constraints throughout the optimization process. The proposed framework is evaluated through a case study involving the optimization of the BRB areas in a two-story BRBF. The objective is to minimize the estimated initial construction cost while satisfying prescribed probabilistic constraints on the first- and second-story interstory drifts under seismic excitation. The results show that the GSS-SPCE accurately represents the response distributions over the probability range relevant to the optimization, including their upper tails. The optimized design satisfies the prescribed probabilistic constraints while requiring smaller BRB areas and having a lower estimated initial construction cost than the original reference design. Under the assumptions adopted for the computational comparison, the proposed framework requires approximately 4000 times fewer nonlinear model evaluations than direct Monte Carlo-based optimization and 40 times fewer than direct GSS-based optimization. The proposed framework provides a computationally efficient approach for optimizing structures under stochastic excitation when response distributions and exceedance probabilities must be evaluated repeatedly for different design candidates. Its use may also be extended to broader performance-based engineering applications involving probabilistic performance assessment under uncertainty. Future developments may include applications to higher-dimensional design problems, the consideration of multiple performance objectives and constraints, the use of gradient-based optimization methods, and the assessment of alternative stratification variables for different natural hazards. 30

Nomenclature Y

Structural response

Ye

Emulated structural response

YR

Reference structural response

YT

Training structural response

y

Response level

ycr

Critical response threshold

p(y)

Exceedance probability associated with response level y

pe(y)

Emulated exceedance probability

pt

Acceptable exceedance probability

plim

Prescribed unconditional probability level

F

Cumulative distribution function

f

Probability density function

d

Vector of design variables

D

Design space

C(d)

Construction cost

Wb

Total weight of the beams

Wc

Total weight of the columns

WBRB

Total weight of the BRBs

X

Vector of explicit random input variables

DX

Input space of X

Xm

Vector of model-related random variables

Xh

Vector of hazard-related random variables

Xw

Vector representing the stochastic variability of the load excitation

Ξ

Complete vector of uncertainties involved in the structural performance assessment

nRV

Number of random variables

Mw

Earthquake moment magnitude

r

Epicentral distance

gNL

Nonlinear structural response model 31

Ms

Stochastic simulator

Md

Deterministic simulator

H

Auxiliary model used to evaluate the stratification variable

Ω

Sample space associated with the internal stochasticity of the simulator

ω

Event in Ω

U

Augmented explicit input vector composed of design and random variables

Z

Latent random variable representing intrinsic simulator variability

ϵ

Additive SPCE noise term

σ

Standard deviation of the SPCE noise term

α

Polynomial multi-index

cα

SPCE expansion coefficient

ψα

Multivariate polynomial basis function

A

Finite set of polynomial multi-indices

S

Number of support points within each stratum

XS

DoE matrix of explicit random inputs

YS

Simulator evaluations associated with XS

US

Stratum-specific DoE in the augmented input space

τ

Subset of hazard-related variables used for stratification

V

Stratification variable

Vi

Boundary of the ith stratum

(i)

Si

ith stratum

P (Si )

Probability of occurrence of the ith stratum

nS

Number of strata

pGSS

Probability constant used in the GSS scheme

n̂

Number of samples used to determine the stratum boundaries

εtail

Dimensionless upper-tail error metric

2 W2,q

Squared Wasserstein distance restricted to the upper tail

q

Quantile defining the lower boundary of the upper tail

Qa , Qb

Quantile functions of the random variables Ya and Yb

32

Declaration of competing interest The authors have no conflict of interest in submitting this manuscript. Data availability Data will be made available on request. Acknowledgment This research was supported in part by the Brazilian agencies FAPESP (São Paulo Research Foundation, grants 2020/14072-7, 2023/02504-8, and 2025/13211-7), CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior, Finance Code 001), and CNPq (National Council for Scientific and Technological Development, grant No. 309107/2020-2). The contribution of the second author is based upon work supported in part by the U.S. National Science Foundation under Award No. CMMI-2118488. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the U.S. National Science Foundation. Appendix A. Seismic hazard and stochastic ground-motion model The stochastic ground-motion model used to generate synthetic records that are statistically consistent with the target spectrum employs a two-corner point-source formulation that captures the physics of fault rupture and wave propagation [62, 49, 63]. Nonstationarity in time is incorporated through a temporal envelope function, while the frequency content is controlled by a parametric radiation spectrum A(f ; Mw , r), where f denotes frequency. Both depend on the earthquake magnitude Mw and epicentral distance r. The frequency-domain radiation spectrum is expressed as: A(f ; Mw , r) = (2πf )2 E(f ; Mw )P (f ; r)G(f )

(A.1)

where E(f ; Mw ) represents the source term, P (f ; r) models path effects, and G(f ) accounts for site amplification.

33

The temporal envelope function, which governs the evolution of the ground motion over time, is given by [64]: e(t; Mw , r) = at



t tn

bt

  t exp −ct tn

(A.2)

where t denotes time, tn is the ground-motion duration, and the constants at , bt , and ct are selected such that e(t; Mw , r) peaks at t = λt tn with a normalized maximum of 1 and satisfies e(tn ; Mw , r) = ηt . Synthetic ground motions are generated by modulating the white-noise sequence Xw = [Xw (i∆t)], where ∆t is the time increment, as follows: 1. Multiply Xw by the envelope function e(t; Mw , r). 2. Transform the resulting sequence to the frequency domain. 3. Normalize by the square root of the mean-square amplitude spectrum. 4. Multiply by the radiation spectrum A(f ; Mw , r). 5. Transform back to the time domain to obtain the acceleration time history. The model parameters used for ground-motion generation follow [62, 49]: radiation pattern RΦ = 0.55, shear-wave velocity V̄s = 3.5 km/s, rock density ρs = 2.8 g/cm3 , seismic velocity cQ = 3.5 km/s, and anelastic attenuation function Q(f ) = 180f 0.45 . Geometric spreading is modeled as Z(R) = 1/R for R < 70 km and Z(R) = 1/70 for R ≥ 70 km, where R is the source-to-site distance. Site amplification follows the generic rock-site model of Boore and Joyner [65]. The envelope-function parameters are set to λt = 0.2 and ηt = 0.05 [64]. The high-dimensional vector Xw represents RTR variability, while Mw and r constitute the hazard-related random vector Xh and represent the dominant seismic hazard parameters [62, 49]. References [1] H. Krawinkler, Challenges and progress in performance-based earthquake engineering, in: International Seminar on Seismic Engineering for Tomorrow–In Honor of Professor Hiroshi Akiyama, volume 26, 1999. [2] S. M. Spence, S. Arunachalam, Performance-based wind engineering: Background and state of the art, Frontiers in Built Environment 8 (2022) 830207. 34

[3] A. Hassanzadeh, S. Moradi, H. V. Burton, Performance-based design optimization of structures: State-of-the-art review, Journal of Structural Engineering 150 (2024a) 03124001. [4] I. D. Rodrigues, S. M. Spence, A. T. Beck, Incorporating statistical uncertainty into metamodel-driven performance-based risk optimization for seismic design, Engineering Structures 348 (2026) 121744. [5] A. T. Beck, W. J. de Santana Gomes, A comparison of deterministic, reliability-based and risk-based structural optimization under uncertainty, Probabilistic Engineering Mechanics 28 (2012) 18–29. [6] S. M. Spence, A. Kareem, Performance-based design and optimization of uncertain wind-excited dynamic building systems, Engineering Structures 78 (2014) 133–144. [7] M. Rastegaran, S. B. Aval, E. Sangalaki, Multi-objective reliability-based seismic performance design optimization of SMRFs considering various sources of uncertainty, Engineering Structures 261 (2022) 114219. [8] P. Ni, Z. Yuan, J. Fu, Y. Bai, L. Liu, Stochastic design optimization of nonlinear structures under random seismic excitations using incremental dynamic analysis, Probabilistic Engineering Mechanics 78 (2024) 103707. [9] R. H. Lopez, K. R. dos Santos, L. F. F. Miguel, An efficient approach for taking into account uncertainties in structural parameters in performance-based design optimization, Engineering Structures 336 (2025) 120382. [10] P. T. Movaghar, A. A. Taflanidis, A. Giaralis, D. Vamvatsikos, Sustainability-driven risk-based design of inerter vibration absorbers for multistory hysteretic buildings under seismic hazard, Journal of Engineering Mechanics 151 (2025) 04025058. [11] S. M. J. Spence, Optimization of uncertain and dynamic high-rise structures for occupant comfort: An adaptive kriging approach, Structural Safety 75 (2018) 57–66. [12] J. Kim, J. Song, Reliability-based design optimization using quantile surrogates by adaptive gaussian process, Journal of Engineering Mechanics 147 (2021) 04021020. 35

[13] L. F. F. Miguel, R. H. Lopez, H. Carvalho, A. T. Beck, Performance-based optimization of nonlinear friction-folded ptmds of structures subjected to stochastic excitation, Mechanical Systems and Signal Processing 195 (2023) 110291. [14] L. F. F. Miguel, S. Elias, A. T. Beck, Reliability-based optimization of supported pendulum tmds’ nonlinear track shape using padé approximants, Engineering Structures 306 (2024) 117861. [15] K. R. dos Santos, A. T. Beck, R. H. Lopez, Sequential simulated annealing for life-cycle optimization of nonlinear stochastic systems via arbitrary polynomial chaos expansion, Engineering Structures 304 (2024) 117675. [16] D. Gu, W. Han, J. Guo, H. Guo, S. Gao, X. Liu, A kriging-based adaptive adding point strategy for structural reliability analysis, Probabilistic Engineering Mechanics 74 (2023) 103514. [17] S. Misra, P. Bocchini, Metaimnet: A physics-informed neural network architecture for surrogate response and fragility modeling of structures subjected to time-varying hazard loads, Structural Safety (2025) 102650. [18] I. Gidaris, A. A. Taflanidis, Performance assessment and optimization of fluid viscous dampers through life-cycle cost criteria and comparison to alternative design approaches, Bulletin of Earthquake Engineering 13 (2015) 1003–1028. [19] I. Gidaris, A. A. Taflanidis, G. P. Mavroeidis, Kriging metamodeling in seismic risk assessment based on stochastic ground motion models, Earthquake Engineering & Structural Dynamics 44 (2015) 2377–2399. [20] A. P. Kyprioti, A. A. Taflanidis, Kriging metamodeling for seismic response distribution estimation, Earthquake Engineering & Structural Dynamics 50 (2021) 3550–3576. [21] J. Kim, S.-r. Yi, J. Song, Active learning-based optimization of structures under stochastic excitations with first-passage probability constraints, Engineering Structures 307 (2024) 117873.

36

[22] D. Jerez, H. Jensen, M. Beer, Reliability-based design optimization of structural systems under stochastic excitation: An overview, Mechanical Systems and Signal Processing 166 (2022) 108397. [23] X. Zhu, M. Broccardo, B. Sudret, Seismic fragility analysis using stochastic polynomial chaos expansions, Probabilistic Engineering Mechanics 72 (2023) 103413. [24] S.-r. Yi, A. A. Taflanidis, Stochastic emulation with enhanced partial-and no-replication strategies for seismic response distribution estimation, Earthquake Engineering & Structural Dynamics 53 (2024) 2354–2381. [25] S.-r. Yi, A. A. Taflanidis, Multi-output stochastic emulation with applications to seismic response correlation estimation, Structural Safety (2025) 102578. [26] J. Kim, Z. Wang, Uncertainty quantification for seismic response using dimensionality reduction-based stochastic simulator, Earthquake Engineering & Structural Dynamics 54 (2025) 471–490. [27] X. Zhu, B. Sudret, Replication-based emulation of the response distribution of stochastic simulators using generalized lambda distributions, International Journal for Uncertainty Quantification 10 (2020). [28] H. Kroetz, A. Beck, L. Costa, F. Macedo, S. Marelli, B. Sudret, Fragility analysis of structures under non-synoptic winds using stochastic emulators, Engineering Structures 357 (2026) 122515. [29] F. C. Macedo, A. T. Beck, K. Giannoukou, L. G. Costa, H. M. Kroetz, L. F. Fadel Miguel, S. Marelli, B. Sudret, Fragility assessment of transmission line supports under non-synoptic wind loads using multi-fidelity stochastic emulators, Probabilistic Engineering Mechanics 85 (2026) 103971. [30] X. Zhu, B. Sudret, Emulation of stochastic simulators using generalized lambda models, SIAM/ASA Journal on Uncertainty Quantification 9 (2021) 1345–1380.

37

[31] X. Zhu, B. Sudret, Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models, Reliability Engineering & System Safety 214 (2021) 107815. [32] X. Zhu, B. Sudret, Stochastic polynomial chaos expansions to emulate stochastic simulators, International Journal for Uncertainty Quantification 13 (2023). [33] A. V. Pires, M. Moustapha, S. Marelli, B. Sudret, Reliability analysis for data-driven noisy models using active learning, Structural Safety 112 (2025) 102543. [34] A. V. Pires, M. Moustapha, S. Marelli, B. Sudret, Reliability analysis for nondeterministic models using stochastic polynomial chaos expansion and active learning, Structural Safety (submitted) (2025). [35] M. Moustapha, B. Sudret, High-dimensional reliability-based design optimization using stochastic emulators, arXiv preprint arXiv:2604.05759 (2026). [36] R. Melchers, Importance sampling in structural systems, Structural safety 6 (1989) 3–10. [37] S.-K. Au, J. L. Beck, Important sampling in high dimensions, Structural Safety 25 (2003) 139–163. [38] S.-K. Au, J. L. Beck, Estimation of small failure probabilities in high dimensions by subset simulation, Probabilistic Engineering Mechanics 16 (2001) 263–277. [39] S.-K. Au, J. L. Beck, Subset simulation and its application to seismic risk based on dynamic analysis, Journal of Engineering Mechanics 129 (2003) 901–917. [40] Z. Wang, M. Broccardo, J. Song, Hamiltonian monte carlo methods for subset simulation in reliability analysis, Structural Safety 76 (2019) 51–67. [41] W. Chen, Z. Wang, M. Broccardo, J. Song, Riemannian manifold hamiltonian monte carlo based subset simulation for reliability analysis in non-gaussian space, Structural Safety 94 (2022) 102134.

38

[42] S. Arunachalam, S. M. J. Spence, An efficient stratified sampling scheme for the simultaneous estimation of small failure probabilities in wind engineering applications, Structural Safety 101 (2023) 102310. [43] J. Xian, Z. Wang, Relaxation-based importance sampling for structural reliability analysis, Structural Safety 106 (2024) 102393. [44] D. Lee, Z. Wang, J. Song, Efficient seismic reliability and fragility analysis of lifeline networks using subset simulation, Reliability Engineering & System Safety 260 (2025) 110947. [45] V. Dubourg, B. Sudret, F. Deheeger, Metamodel-based importance sampling for structural reliability analysis, Probabilistic Engineering Mechanics 33 (2013) 47–57. [46] X. Zhu, Z. Lu, W. Yun, An efficient method for estimating failure probability of the structure with multiple implicit failure domains by combining meta-is with is-ak, Reliability Engineering & System Safety 193 (2020) 106644. [47] N.-C. Xiao, H. Zhan, K. Yuan, A new reliability method for small failure probability problems by combining the adaptive importance sampling and surrogate models, Computer Methods in Applied Mechanics and Engineering 372 (2020) 113336. [48] T. Yu, Z. Lu, A novel single-loop kriging importance sampling method for estimating failure probability upper bound under random-interval mixed uncertainties, Aerospace Science and Technology 143 (2023) 108722. [49] S. Arunachalam, S. M. J. Spence, Generalized stratified sampling for efficient reliability assessment of structures against natural hazards, Journal of Engineering Mechanics 149 (2023) 04023042. [50] J. W. Baker, C. A. Cornell, Uncertainty propagation in probabilistic seismic loss estimation, Structural Safety 30 (2008) 236–252. [51] N. Lüthen, X. Zhu, S. Marelli, B. Sudret, UQLab user manual – Stochastic polynomial chaos expansions (SPCE), Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, Zurich, 2026. Report UQLab-V2.2-121. 39

[52] AISC, ANSI/AISC 341-10: Seismic provisions for structural steel buildings, American Institute of Steel Construction (2010). [53] ASCE, ASCE/SEI 7-10: Minimum design loads for buildings and other structures, American Society of Civil Engineers (2013). [54] C.-M. Uang, M. Nakashima, K.-C. Tsai, Research and application of buckling-restrained braced frames, International Journal of Steel Structures 4 (2004) 301–313. [55] NIST, NIST GCR 10-917-8: Evaluation of the FEMA P695 methodology for quantification of building seismic performance factors, National Institute of Standards and Technology (2010). [56] A. Ghasemof, M. Mirtaheri, R. K. Mohammadi, Multi-objective optimization for probabilistic performance-based design of buildings using fema p-58 methodology, Engineering Structures 254 (2022) 113856. [57] H. Guerrero, A. Terán-Gilmore, T. Ji, J. A. Escobar, Evaluation of the economic benefits of using buckling-restrained braces in hospital structures located in very soft soils, Engineering Structures 136 (2017) 406–419. [58] AISC, Specification for structural steel buildings, ANSI/AISC (2022). [59] F. McKenna, G. Fenves, M. Scott, Opensees: Open system for earthquake engineering simulation. pacific earthquake engineering research center, uc berkeley, ca, 2006. [60] F. C. FILIPPOU, E. P. POPOV, V. V. BERTERO, Effects of bond deterioration on hysteretic behavior of reinforced concrete joints, Earthquake Engineering Research Center, University of California, College of Engineering, Berkeley, 1983. (Report no UCB/EERC/19, August 1983). [61] J. H. Holland, Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence, MIT Press, Cambridge, 1992.

40

[62] C. Vetter, A. A. Taflanidis, Global sensitivity analysis for stochastic ground motion modeling in seismic-risk assessment, Soil Dynamics and Earthquake Engineering 38 (2012) 128–143. [63] G. M. Atkinson, W. Silva, Stochastic modeling of california ground motions, Bulletin of the Seismological Society of America 90 (2000) 255–274. [64] D. M. Boore, Simulation of ground motion using the stochastic method, Pure and applied geophysics 160 (2003) 635–676. [65] D. M. Boore, W. B. Joyner, Site amplifications for generic rock sites, Bulletin of the seismological society of America 87 (1997) 327–341.

41

Related documents

• Optimal experimental design: Formulations and computations#183155
• Exact Algorithms for Network Fortification and Design Problems under Stochastic Interdiction#244652
• A Sensitivity Analysis Methodology for Rule-Based Stochastic Chemical Systems#627583
• Economy of scale in algorithm design : approximation algorithms for submodular optimization and network design#195727
• A Comprehensive Review of Archimedes Optimization Algorithm with its Theory, Variants, Hybridization, and Applications#250495