Surrogate probabilistic seismic demand modelling of inelastic single‐degree‐of‐freedom systems for efficient earthquake risk applications

This paper proposes surrogate models (or metamodels) mapping the parameters controlling the dynamic behaviour of inelastic single‐degree‐of‐freedom (SDoF) systems (i.e., force‐displacement capacity curve, hysteretic behaviour) and the parameters of their probabilistic seismic demand model (PSDM, i.e., conditional distribution of an engineering demand parameter [EDP] given a ground‐motion intensity measure [IM]). These metamodels allow the rapid derivation of fragility curves of SDoF representation of structures. Gaussian Process (GP) regression is selected as the metamodelling approach because of their flexibility in implementation, the resulting accuracy and computational efficiency. The metamodel training dataset includes 10,000 SDoF systems analysed via cloud‐based non‐linear time‐history analysis (NLTHA) using recorded ground motions. The proposed GP regressions are tested in predicting the PSDM of both the SDoF database (through ten‐fold cross validation) and eight realistic reinforced concrete (RC) frames, benchmarking the results against NLTHA. An application is conducted to propagate such modelling uncertainty to both fragility and vulnerability/loss estimations. Error levels are deemed satisfactory for practical applications, especially considering the low required modelling effort and analysis time. Regarding single‐building applications enabled by the proposed metamodel, this paper presents a first attempt at a direct loss‐based design procedure, which allows setting a target loss level for the designed structure (shown for a realistic RC frame). An earthquake risk model involving dynamic exposure and vulnerability modules is illustrated as an example of building portfolio applications. Specifically, the proposed application considers a retrofit‐based seismic risk‐reduction policy for a synthetic building portfolio, for which it is possible estimating the loss evolution over time.


INTRODUCTION AND MOTIVATION
Building-level seismic fragility is defined as the probability of reaching or exceeding various damage states (DSs) as a function of a hazard intensity measure (IM). DSs are usually expressed in terms of engineering demand parameter (EDP) thresholds, such as inter-storey drift limits. In the realm of analytical (or numerical) fragility analysis methods, nonlinear time-history analysis (NLTHA) of refined structural models is the best practice when it comes to building-specific applications. 1,2 On the other hand, NLTHA of equivalent single degree of freedom (SDoF) systems is commonly used to characterise a building vulnerability class, 2 with the results being used for building-portfolio earthquake risk analysis. 3 This reflects the trade-off between simplicity and accuracy that different stakeholders generally tolerate on the matter: 4 private owners likely need a detailed assessment of individual buildings or small building portfolios, while governmental agencies or (re)insurance companies generally look at large portfolios accepting larger uncertainties and a lower analysis refinement level. The above-mentioned analysis methods, respectively targeted at single-building and building-portfolio applications, allow an analyst to easily overcome computational issues, especially considering the excellent performance of ordinary modern computers. However, at least arguably, some applications may still be computationally unaffordable if using the above approaches. An example of such applications related to single buildings may be loss-based conceptual seismic design. This would require calculating a loss metric (e.g., the expected annual loss [EAL]) for (several) tentative structural configurations, considering a specific seismic hazard profile: only the configurations complying with the target loss level will be used to continue the detailed design process. Even considering SDoF-based NLTHA, the required analysis time may not be compatible with such a preliminary/conceptual design process, which requires pseudo-instantaneous analysis results. Lower-refinement results may suffice, and more refined analysis methods can be applied in the subsequent phases of the design. Considering building portfolios, the design of regional policies for seismic risk reduction (e.g., based on structural retrofit and/or seismic insurance) may be computationally unaffordable using SDoF NLTHA. In fact, this may require an earthquake risk model with many exposure scenarios (which may also be time-dependent), allowing selecting the optimal policy using one or more loss metrics.
Metamodelling approaches may be suitable for the applications mentioned above. A metamodel is "a model of a model": it defines a (statistical) relationship between a given set of inputs and outputs (obtained through numerical simulations). In this particular case, metamodels mapping the parameters controlling the dynamic behaviour of inelastic SDoF systems (e.g., force-displacement capacity curve, hysteretic behaviour) and their probabilistic seismic demand model (PSDM; i.e., conditional probability distribution of an EDP given an IM) seems suitable. In fact, those can easily lead to building-level fragility curves for a set of structure-specific DSs, enabling the computation of various risk metrics/decision variables (e.g., economic losses, casualties, downtime). Different metamodels were proposed in the literature for the seismic response of structures, often referring to refined multi-degree-of-freedom (MDoF) structures (e.g., ref. 5,6). For the scope of this paper, the most relevant example of such metamodels is SPO2IDA, 7 surrogating the different percentiles of incremental dynamic analysis (IDA) curves of SDoF systems from their piecewise static pushover (SPO) curve, which include elastic, hardening and degrading branches. This parametric model is based on regressions, and it is reasonably accurate within the range of its training dataset. However, the authors highlight its limitations related to the low number of ground-motion records (i.e., 30) and the single hysteresis model (moderately pinching) adopted for the training of the model. The authors provide a methodology to extend the training dataset to other periods by using their proposed functional forms. Although a similar methodology may be applied to extend the range of applicability of SPO2IDA to different parameters (e.g., considering other hysteresis models), a fixed functional form may fail to capture the complex interaction between the input parameters and the parameters of the resulting PSDMs. An example of this process is presented in ref. 8. In this paper, a metamodelling approach is proposed to surrogate SDoF PSDMs while addressing the above issues, enabling efficient seismic risk assessment applications. First, it is recognised that most widely adopted DS descriptions (e.g., ref. [9][10][11] often cover DS levels up to the so-called "near collapse." This generally involves a structural deformation level corresponding to the peak of the pushover curve or to a 20% post-peak strength drop (which can be reasonably close to the peak, depending on the softening branch of the pushover). In this context, characterising PSDMs above such deformation levels may be avoided by accepting an approximation of the results. Restricting the PSDM characterisation to pre-peak deformation levels allows drastically reduce the number of SDoF systems in the training set. Moreover, a cloud-based analysis approach 12 is herein preferred to IDA, to use a higher number of seed ground-motion records at a lower computational cost. Most importantly, a Gaussian Process (GP) regression is adopted since it does not require any prior definition of the output functional form (a GP regression is a non-parametric model). Not only this approach would F I G U R E 1 Steps for the calibration of the proposed surrogate model ensure high accuracy and flexibility in implementation for earthquake engineering applications (e.g., ref. [13][14][15], but it will also result in an infinitely scalable surrogate model. In fact, it will be possible to add new case studies to the training dataset (e.g., adding more hysteresis models or damping ratios) and re-fit the GP regressions without the need to change the metamodel structure.
After describing the metamodelling strategy (Section 2), its accuracy is shown in Section 3, considering both SDoF systems and eight reinforced concrete (RC) frame case studies. Finally, Section 4 shows two use cases for the provided metamodel: direct loss-based seismic design (DLBD) and scenario-based dynamic earthquake loss modelling for riskinformed decision making.

METAMODELLING STRATEGY
The mapping between the parameters controlling the dynamic behaviour of the considered inelastic SDoF systems and their PSDM consists of three phases ( Figure 1). First, a number of input parameters directly affecting the PSDM are selected, together with their desired range (described in detail in Section 2.1). Plain Monte Carlo sampling is performed to generate a database of case-study SDoF systems. Second, each SDoF system is analysed via cloud analysis adopting a suite of 100 ground-motion records. A PSDM is constructed for each SDoF sample in the training dataset (described in detail in Section 2.2). A GP regression is finally fit to surrogate the complex relationship between the input SDoF parameters and the output PSDM ones (Section 2.3). It is worth mentioning that the GP regressions are tested with a ten-fold cross validation (Section 3). The code implementing the proposed metamodel is freely available (https://github.com/robgen/surrogatedPSDM). In addition, the same repository includes the adopted training datasets and the code to train the GP regressions so that users can filter, extend or modify the training dataset and update the fitting.

Generation of the SDoF database
The selected input parameters defining each SDoF system ( Figure 1A) are the considered hysteresis model "ℎ ," governing stiffness degradation under unloading-reloading conditions; the fundamental period, , controlling the elastic stiffness; the yield shear strength (normalised to the total weight), ; the hardening ratio. ℎ. It is worth mentioning that no within-cycle strength degradation is considered since the surrogated PSDMs will not be used for predictions exceeding the ductility at the peak pushover strength. For the same reason, no explicit limit for the ductility (or displacement) at peak capacity is considered.
Five different hysteresis models are adopted for the model training ( Figure 1A): • Kinematic hardening, KIN, showing no stiffness degradation, mimicking the response of steel structures. Such an interpretation implicitly neglects the Bauschinger effect, which is likely not relevant for deformation levels smaller than the peak pushover response; • Modified Takeda "fat," 16 MTf, appropriate for structures with a non-linear response dominated by members with low axial load (e.g., newly-designed RC frames dominated by beam behaviour). The unloading ( ) and reloading ( ) stiffness coefficients are respectively equal to 0.3 and 0.6; 17 • Modified Takeda "thin," MTt, appropriate for structures with a non-linear response dominated by members with large axial load (e.g., soft storey-like RC frames; RC bridge piers). The and stiffness coefficients are respectively equal to 0.5 and 0; 17 • Modified Sina, MS, appropriate for existing structures exhibiting pinching response (e.g., RC frames dominated by joint failure; timber structures dominated by the failure of their connections). The unloading and reloading stiffness coefficients are respectively equal to 0.3 and 0.6, while pinching is defined by a closing force equal to 25% of the yield force; 17 • Flag shape, FS, appropriate to simulate the behaviour of hybrid pre-stressed structures with re-centring behaviour. The energy dissipation coefficient ( ), controlling the force level where the re-centring action starts, is equal to 0.5.
In fitting the GP regression, the hysteresis type is considered a categorical (i.e., non-numerical) variable, as opposed to explicitly mapping the response of the SDoF systems to the numerical parameters of the hysteresis models (e.g., and for the MT). This is done for three reasons: (1) because the hysteresis-parameter types are different for each hysteresis model; (2) to avoid excessive interpolation error (e.g., between the two sets of MT hysteresis parameters); (3) to avoid a high increase in the number of required SDoF case studies (to consider more hysteresis-parameters sets).
For each hysteresis model, 2,000 realisation of the parameters ( , , ℎ) are sampled via a plain Monte Carlo approach, adopting uniform distributions. As discussed in Rasmussen and Williams,18 this is an appropriate assumption when adopting GP regressions as the focus is the input-output mapping rather than on modelling the specific distribution of the input parameters. The specific parameters of the assumed distributions, clearly defining the scope of the final metamodels, are ∼ (0.2 , 1.5 ); ∼ (0.05, 0.6); ℎ ∼ (0, 0.3).

Seismic response analysis and PSDM
Each SDoF system in the dataset is analysed to derive a cloud of points in the EDP vs IM space. The selected EDP in this study is the ductility demand . The selected IM is the pseudo-spectral acceleration at the SDoF period, normalised to the yield base shear coefficient, = ∕ . Although more advanced IMs are available (e.g., ref. 19), such a simple choice may simplify the hazard analysis in practical applications or even exploit existing hazard models.
A set of 100 natural (i.e., recorded) ground motions is selected from the SIMBAD database, "Selected Input Motions for displacement-Based Assessment and Design." 20 As per ref. 13, the three-component 467 records in the database are ranked according to their peak ground acceleration, PGA, values (by using the geometric mean of the two horizontal components) and then keeping the horizontal component with the largest PGA value. The first 100 records are arbitrarily selected, characterised by moment magnitudes in the range 5-7.3, station-to-source distance smaller than 35 km; and PGA ranging between 0.29 and 1.77 g. Such an ground-motion selection approach is consistent with the adopted response analysis method (i.e., cloud analysis). Hazard-consistent, site-specific record selection is outside the scope of this work since the aim is to calibrate a flexible metamodel of the PSDM for a large number of hazard/site conditions. NLTHAs are conducted scaling the records to ensure non-linear response ( > 1) for every case study (which show a wide range of period of vibration and yield strength), considering 100 guesses target values of ductility demand, equally spaced between 1 and 6. The equal displacement rule ( ≈ ) is used to derive reasonable guesses of the scale factors ( ): (1) a record and a guess ductility demand are randomly selected (without repetition) from the relevant sets; (2) the scale factor is calculated assuming that the scaled spectral acceleration is equal to , and using the spectral acceleration for the unscaled record at the SDoF period. The adopted scale factors range between 0.20 and 5.6, keeping the bias introduced in the response analyses to acceptable levels. 21 It is worth mentioning that a recent study 22 suggests that this level of scaling might create bias in estimating collpase probability. Therefore, for some specific applications, this aspect should be further checked and possibly corrected by using procedures available in the literature (e.g., ref. 23).
The ductility demand of each analysis is checked against the guessed value. Analysis results leading to a ductility demand outside of the target range are excluded from the PSDM fit described below, ensuring that at least 70 records are adopted for each cloud.
The adopted bi-linear PSDM ( Figure 1, step 2) is defined according to Equation 1. Bi-linear PSDMs have been demonstrated effective in previous work (e.g., ref. 24). Clearly, the PSDM in Equation 1 does not imply the equal displacement rule. Characterising the behaviour in the elastic range is trivial and follows from the definition of an elastic SDoF ( = ), not requiring any analysis for its calibration. The inelastic range is obtained performing a linear regression in the logarithmic space, where ln( −1)| −1 (henceforth simply called ) is the logarithmic standard deviation of the pairs − 1 versus − 1, and is a standard normal variable. Therefore, the median relationship is the line = ( − 1) + 1, where is its slope. Such a model choice, compatible with other literature studies (e.g., ref. 7), implies a lognormal distribution of the residuals, which is desirable in calculating lognormal fragility curves (see Section 2.4). Such model also implies homoscedasticity for > 1. As demonstrated previously, 12 this assumption is deemed satisfactory for a variety of engineering systems (including SDoF systems).
Depending on the considered practical applications, the SDoF parameters may be affected by a degree of variability. Although this is out of the scope of this paper, both the provided GP regressions and/or the training data can be used to appropriately account for the effect of such variability on the predicted seismic response (e.g., ref. 25).

2.3
Gaussian process regressions for the ( , , , ) -( , ) map For each SDoF realisation depending on a vector of input parameters = {ℎ , , , ℎ} , the methodology in Section 2.2 allows defining two PSDM output parameters ( , ). Two independent training datasets, ( , ), are composed by the matrix collecting the input vectors for all the SDoF realisations, and the vector which includes the related or outputs. Based on such training datasets, a GP regression is fitted for each PSDM parameter. This takes a vector of unique inputs ( ) and produces an output/target = ( ) using a statistical model, while being computationally cheaper than the methodology shown in Section 2.2. Based on a training dataset made of inputs (covariates) and known outputs, a GP regression is fitted so that it is possible to make predictions for any input vector outside the training dataset. Although a detailed mathematical description of GP regressions (and their fitting) is outside the scope of this paper, it is worth providing a general, high-level perspective on the matter. Rasmussen and Williams 18 present an exhaustive mathematical description/derivation of GP regression. A step-by-step mathematical description for an earthquake engineering application is provided, for example, ref. 13 In a GP regression, ( ) is regarded as a realised value of a GP. A GP is a generalisation of the Gaussian probability distribution model, describing the distribution of functions ( ) rather than scalars or vectors. Its mean and covariance functions fully specify a GP. According to a Bayesian framework, the first step in a GP regression is to set a prior distribution for all the possible functions ( ), reflecting the starting knowledge about the output before having any data. This is done by assigning some properties of the mean and the covariance functions (e.g., smoothness). Then, the prior distribution is converted into a posterior distribution (over functions) based on the observed data, and such a posterior distribution is used for predictions. The properties of the output function ( )-with particular reference to its smoothness-are governed by the covariance function, which captures the correlation among different input vectors and reflects it in the output. The training of a generic, scalar-input function is shown in Figure 2.
The structure of the covariance function is selected by the GP user and should reflect the expected behaviour of the output. A popular choice of covariance function is the squared exponential covariance-also adopted in this studysince it reflects the "stability" of the involved physical quantities (i.e., a small perturbation of the input SDoF parameters produces a slight change in the output PSDM parameters). The covariance function parameters are called hyperparameters since they are not assigned (a GP regression is a non-parametric model) but inferred from the training dataset (usually with a maximum likelihood estimator). Given the hyperparameters, the posterior predictive distribution of the GP is obtained by conditioning the prior distribution to the training dataset. The expected value of the posterior distribution is adopted to make predictions.
With regard to the GP regression fitting for this particular study, a constant basis function is assumed for the posterior distribution of the mean function. In addition, a squared exponential covariance function with separate length scales is adopted. This is done to have a reasonable amount of hyperparameters while reflecting the mechanics of the phenomenon under investigation. A quasi-Newton method 26 is adopted to optimise the hyperparameters, including the noise variance.

Seismic fragility and loss analysis
Once the PSDM parameters ( , ) are evaluated using the trained GP regressions, it is possible to derive fragility functions to perform seismic risk and/or loss analysis. Building-level fragility curves are calculated for a set of structure-specific DSs, identified by the thresholds . One possibility involves choosing four DSs, slight, moderate, extensive and complete damage, as defined according to HAZUS. 9 Other definitions of the DSs are possible, 11,27 and the proposed framework is independent of their particular choice.
According to the properties of the adopted PSDM (Equation 1), lognormal fragility curves for each DS, representing the DSs' exceeding probability, = ( ≥ ), are completely specified by their median and logarithmic standard deviation (simply called dispersion), which are specified in Equation 2 both for the elastic and inelastic ranges. It is worth mentioning that collapse cases (corresponding to ground motions leading to dynamic instability of the analysis or exceedance of a conventional 10% drift threshold) are not expected since the SDoF systems are only subjected to pre-peak ductility demand levels.
Loss analysis is based on vulnerability curves, which can be derived using a building-level consequence model relating the repair-to-replacement cost to structural and non-structural DSs. This involves defining the expected building-level damage-to-loss ratios ( ) for each DS. The (mean) loss ratio ( ) for a given IM value is defined according to Equation 3, allowing deriving a mean vulnerability curve (involving the difference between the exceeding probability of the ( + 1) ℎ and ℎ DSs). More advanced, component-based fragility/loss methodologies are available (e.g., ref. 28), which are appropriate for more refined applications, such as the detailed seismic loss estimation of individual buildings, for instance.
One of the relevant loss metrics commonly adopted in practice is the EAL. Assuming a hazard curve, this is calculated according to Equation 3, where is the mean annual frequency of exceeding a given value of the IM.

Surrogated-vs-modelled error and cross-validation of the GP regressions
The trained GP regressions are subjected to various tests to investigate their predictive power. First, the PSDM parameters are predicted for the entire SDoF database, measuring the surrogate-vs-modelled normalised root mean squared error, ). Figure 3A shows that the surrogated-vs-modelled points are particularly close to the = line, with the NRMSE being equal to 2.5% for the slope and 6.2% for the logarithmic standard deviation of the PSDMs. Considering the various sources of uncertainty commonly involved in seismic performance assessment or risk/loss models, such error levels are deemed satisfactory for practical applications.
To test the predictive power of the GP regressions outside the training set, ten-fold cross validation is carried out. This first involves randomly dividing the training dataset into ten subsets made of 1,000 SDoF samples each. Then, ten new GP regressions are fitted alternatively, excluding one subset from the new training dataset. For each subset, the fitted GP regressions are used to make (in-fold) predictions for the subset kept out of the training dataset. The in-fold prediction errors for each of the ten sub-sets are used to calculate the in-fold NRMSE, which is 0.1% higher than the previously calculated one for both PSDM parameters, thus confirming the validity of the error estimates. It is worth mentioning that such error estimates are consistent with the results in ref. 7, which, for the case studies most similar to those in this paper, report a 3%-6% error on the PSDM estimation (although those are based on 30 ground-motion records).
As expected for GP regressions, Figure 3B shows that the prediction errors (herein quantified as surrogated-vs-modelled ratios) follow a Normal distribution with a mean equal to 1.0001 and 1.0002 (for and ) and standard deviation essentially equivalent to the NRMSE. This result is particularly important to infer the propagation of such error in fragility or risk/loss estimations, as discussed below in this Section.
The surrogated PSDMs for the entire SDoF dataset are used to estimate seismic fragility curves. Those are derived for two DSs related to ductility thresholds, respectively, equal to 3 and 4, representing DS3 (extensive damage according to HAZUS) and DS4 (complete damage according to HAZUS). Elastic DSs are not considered since, according to Equation 2, no error is expected for those. Figure 4A shows the surrogated vs modelled plot for the median of the fragility curves, also indicating a satisfactory NRMSE, respectively equal to 2.0% and 2.2% for DS3 and DS4. Consistently with the definition of the PSDM (Figure 1 and Equation 1), a specified error in the slope parameter propagates to higher error levels for the fragility median as the ductility threshold increases. It is worth mentioning that no further discussion is needed for the fragility dispersion since this is exactly equal to the PSDM logarithmic standard deviation. Also in this case, as expected, the surrogated-vs-modelled ratios for the median of the fragility curves ( Figure 4B) follow a Normal distribution (with mean equal to 0.9993 and standard deviations equal to 0.018 and 0.02 for DS3 and DS4, respectively). Quantifying how the uncertainty on the SDoF parameters affects the prediction of the PSDM and fragility ones is considered out of scope herein. However, after characterising the probability distribution of the SDoF parameters (and their correlation), a user can adopt the proposed GP regressions to carry out the above uncertainty quantification. The last step of this error-propagation exercise involves calculating both vulnerability curve and EAL for two illustrative (arbitrary) case-study SDoF systems with an elastic period equal to 0.5s. Once again, both systems are characterised according to the four structure-specific DSs qualitatively defined according to HAZUS and quantitatively measured on the SDoFs backbones. The median of the fragilities for the elastic-range DSs is the same for both case studies and is equal to 0.04 g and 0.1 g for DS1 and DS2, respectively. The first case study represents a structure with poor seismic performance; thus, the DS3 and DS4 medians are equal to 0.12 g and 0.16 g, respectively. For the second case study, showing considerably higher performance, the DS3 and DS4 medians are equal to 0.36 g and 0.44 g, respectively. A fragility dispersion equal to 0.24 is adopted for all DSs (both elastic and inelastic) for both case studies. A detailed motivation to assume the same dispersion for both the elastic and inelastic ranges is given in Section 3.2.
DLRs equal to 2%, 10%, 43.5%, and 100% of the total reconstruction cost (DS1 to DS4) are adopted for this example, consistently with the assumptions in ref. 13. Vulnerability curves are calculated according to Equation 3. Finally, the EAL is calculated according to Equation 4, assuming the hazard curve appropriate for L'Aquila, a high-seismicity town in central Italy. Such a curve is consistent with the official Italian hazard model (Stucchi et al. 29 ), implemented in the current Italian building code. 30 Although the absolute values of the hazard curve are deemed less relevant for this error-propagation exercise, the adopted hazard model is consistent with those used in Section 4.1.
Starting from the "baseline" fragility curves described above, a plain Monte Carlo approach is used to simulate 30,000 new realisations of the fragility parameters such that = * , = * . The variables , are consistent with the empirical distributions of the surrogated-vs-modelled ratios shown in Figure 4B (for the medians) and Figure 3B (for the dispersion). By re-calculating both the vulnerability for each realisation, the fragility error is propagated to the EAL. For both case studies, Figure 5 shows the empirical distribution of the EAL as a ratio of the baseline EAL (the insertions in each panel show the baseline fragility curves for the case studies). The overall distribution F I G U R E 6 Selected incremental retrofit solutions: (A) plastic mechanism at the life-safety damage state, DS3; (B) pushover curves (modified after Gentile et al. 32 ) of the error is a particularly narrow Normal distribution with a mean respectively equal to 1.0007 and 1.0006 and standard deviation equal to 0.015 and 0.01 for the lower-and higher-performance case studies, respectively. In the Authors' opinion, such low error levels confirm the suitability of the proposed surrogate model in loss assessments involving SDoF-based fragility curves and building-level damage-to-loss ratios (e.g., ref. 2).

Application to realistic RC frame case studies
Most practical applications involve more complex MDoF structures rather than SDoF systems. Therefore, it is crucial to test the predictive power of the trained GP regression against realistic structures. The selected case study is the central longitudinal frame of a three-storey RC building with rectangular plan geometry and structural details of beams, columns and joints consistent with pre-1976 design according to an older Italian building code. 31 Apart from this as-built configuration, consistent with gravity-only design and neglecting any capacity-design provision, seven retrofitted configurations providing incremental seismic performance are designed implementing the RC jacketing technique ( Figure 6). Detailed descriptions and illustrations of geometry, load analysis and structural details of the as-built and retrofitted members are provided in ref. 32. They are not repeated here for brevity. A 2D lumped-plasticity model (bare frame) is developed using the finite element software Ruaumoko 33 for each configuration. The adopted numerical modelling strategy was extensively validated against experimental results. 34 Floor diaphragms are modelled as rigid in their plane, and fully fixed boundary conditions are considered at the base. P-Delta effects are not modelled since they are deemed negligible for three-storey frames. A 5% tangent stiffness-proportional elastic damping is assigned to all frequencies. The flexural capacity of the RC members is derived using moment-curvature analysis. The flexural response is checked against other failure mechanisms that may significantly modify the lateral member response. The capacity model of the members included in the analysis is modified (if appropriate) to include those failure mechanisms accordingly, considering slab-related flange effect for the negative beam moment capacity, lap splice failure, shear capacity, bar buckling. 13 The modified Takeda hysteresis model 16 is used for beams and columns, with the columns having a thinner loop. The hysteretic behaviour of the beam-column joints is modelled using the Modified Sina model, 16 which can capture their pinching behaviour. Each frame configuration is first analysed via displacement-control pushover with a linear force profile ( Figure 6). This allows quantifying four displacement thresholds (Δ ) compatible with the DS definitions in HAZUS (described in Section 2.4), to be adopted for fragility derivation. Those refer to the first member in the frame reaching first cracking, yielding, ¾ of the near collapse drift and the near collapse drift respectively for DS1-DS4. It is worth mentioning that the member "causing" a given DS may change for different DSs.
The as-built configuration shows a storey-level failure mode developing at the first storey ( Figure 6A), characterised by plastic hinging of the columns and shear failure of the exterior beam-column joints. Although the shear failure of the exterior joints "avoids" the soft-storey behaviour, this failure mode results in a low strength ( Figure 6A) approximately equal to 0.2 g spectral acceleration capacity for the equivalent SDoF system, whose elastic period is 0.77s. The accelerationdisplacement capacity ( Figure 6B) of the equivalent system is derived from the pushover force-displacement curve assuming a (first-mode) effective mass equal to 90% of the total building mass (obtained averaging the participating mass of the case studies, ranging between 88% and 93%). Moreover, the effective height displacement is adopted as an EDP, using the effective height formulation by Priestley et al. 17 The adopted retrofit strategy has the objective of inverting the local hierarchy of strength at the sub-assembly level and transforming such a localised failure mode into a more reliable Beam-Sway global mechanism (with flexural plastic hinges forming in all beam ends and at the ground section of the columns). Such a goal is obtained "incrementally," adopting concrete column jacketing as the selected retrofit technique. This results in seven retrofit solutions, for which Figure 6A shows the plastic mechanism at DS3 (extensive damage according to HAZUS; generally associated with life safety). For the solutions I1, I2 and I3, only the interior columns are jacketed (respectively up to the first, second or third storey). Similarly, the solutions IE1, IE2, and IE3 include column jacketing for both interior and exterior columns. Finally, the IE3+ retrofit solution improves IE3 by involving enhanced jacketing for the first-storey columns to provide higher strength for the frame. It is worth mentioning that the beam-column joints located between two jacketed columns are reinforced with horizontal stirrups having the same layout as the jacketed columns, thus significantly enhancing their shear capacity (and avoiding shear hinging).
The response of each frame configuration is independently analysed using two analysis methods with increased refinement. The cloud capacity spectrum method (Cloud CSM, 35 ) is first adopted: this involves applying the CSM, 36 assuming a suite of as-recorded ground-motion spectra. In this case, the set of 100 records described in Section 2.2 is adopted without using any scale factors. The response (effective height displacement) for each ground motion is calculated by using the pushover curves in Figure 6B. Moreover, by adopting the same record suite, cloud-based NLTHA 12 is carried out for each frame configuration, registering the peak effective height displacement. The spectral acceleration at the fundamental period, ( 1 ), is finally calculated for each record and adopted as an IM. 1 ranges between 0.77s for the as built configuration and 0.55s for the IE3+ one.
The two alternative (EDP vs IM) sets are adopted to derive the conditional mean and standard deviation of EDP given IM and derive PSDMs for each frame configuration. Consistently with the common practice, the power-law model EDP = 1 2 is obtained via the least square method, where 1 and 2 are the parameters of the regression. This allows deriving fragility curves compatible with those derived using the GPs, where = ( ln( It is worth mentioning that some ground motions led to collapse, which is herein defined as a global dynamic instability (i.e., non convergence) of the numerical analysis, likely corresponding to a plastic mechanism (i.e., the structure is underdetermined) or exceeding the nominal threshold of 10% maximum inter-storey drift. The information carried out in these analyses with non-numerical EDP values are included in the definition of the fragility curves according to the procedure in ref. 37, which is described in detail in ref. 13. The parameters of the fragility curves are finally computed using the trained GP regressions, after bi-linearising the pushover curves, 38 adopting the MTf hysteresis type (deemed appropriate for RC frames not developing a soft-storey mechanism), and using DS ductility thresholds consistent with the above-mentioned Δ ( Table 1). Based on the GP results, literature formulations can be used to calculate other EDPs such as displacement profiles 17 or peak floor accelerations. 28 Before comparing fragility estimates of a surrogated SDoF model and an explicit MDoF model, it is worth highlighting the various sources of error affecting it. First, the two models involve a different functional form of the PSDM, respectively suited for SDoF and MDoF models. Moreover, the comparison adopts the spectral acceleration at the fundamental period as IM, although more sufficient and efficient IMs are available (e.g., ref. 19). Moreover, although the same set of records are used, those are scaled in amplitude in the metamodel training, possibly leading to further bias. 21

F I G U R E 7 Surrogated versus modelled fragility median for the MDoF case study. Modelling is alternatively based on (A) capacity spectrum method and (B) non-linear time history analysis
As shown in Figure 7A, comparing the surrogated medians with those calculated via the Cloud CSM highlights a very low prediction error: the NRMSE is equal to 1.0%, 0.6%, 10.4%, and 11.2% for DS1, DS2, DS3, and DS4, respectively. Such error level increase to 24.0%, 26.0%, 17.7%, and 19.0% if the benchmark analysis method is the NLTHA of the MDoF numerical models ( Figure 7B). Although such an increase is substantial, it is still deemed consistent with the simplicity of the GP metamodel as opposed to the required modelling effort and time involved in the MDoF NLTHA.
From a qualitative point of view, the GP regressions provide a conservative underestimation of the medians (with respect to the NLTHA), except for the DS4 estimation for three case studies. Consequently, the GP regressions would provide conservative overestimations in risk/loss analyses. Confirming the discussion in Section 3.1, the error levels are generally higher for severe DSs and higher-performing structures (higher fragility medians). However, in risk/loss estimations, such error levels will likely have a low impact given the lower hazard frequency related to high intensity (e.g., in Equation 4 they will be weighted down by the lower values as IM increases). This effect, valid for any combination of building models and sites, is further discussed in Section 4.1, also showing an example.
On the other hand, the GP regressions show a lower predictive power in estimating the fragility dispersion. In fact, the fragility dispersion predicted with the GP regressions for the non-linear range lies within 0.35 and 0.4 for all the considered frame configurations. Despite an NRMSE equal to 23.6% and 29.0%, respectively, assuming the Cloud CSM and the NLTHA as benchmark models, in both cases, the error for the single datapoints approximately ranges between 9% and 59%, and it is generally an underestimation.
It is worth highlighting that the fragility dispersion, in turn depending on the EDP|IM standard deviation, is strongly affected by the refinement of the model (e.g., a higher dispersion is expectable for MDoFs than for SDoF systems, e.g., due to higher-mode effects and lower sufficiency of the IMs). Accordingly, the adopted PSDM model for MDoFs allows accounting for the EDP|IM uncertainty also in the elastic range, while this is not present in the SDoF response. To partially compensate for this deficiency, it is suggested to use the GP-based fragility dispersion also for the elastic DSs (which would have = 0 according to Equation 2). On the other hand, considering collapse cases in the MDoF model usually reduces the fragility dispersion, especially if the structure develops a plastic mechanism for low seismic intensity levels. The combination of such effects, not adequately captured in the GP regression approach, is case-dependent, and it is arguably challenging to predict a general trend.
For this reason, although the fragility dispersion errors are numerically higher than those of the medians, those are still comparable to those obtained explicitly modelling an SDoF system (arguably the standard for portfolio analyses) subjected to cloud-based NLTHA.

Direct loss-based seismic design: a first attempt
The first application relying on the proposed surrogate model is the DLBD. A first attempt DLBD procedure, among other possible alternatives (e.g., refs. 39-41), is proposed in this Section, along with an application to a realistic RC frame structure analysed along one horizontal direction. This procedure is not intended as definitive, and future research work is suggested, and underway, to refine and validate it, possibly modifying some of its steps. By analogy with the words of Priestley, 17 the fundamental philosophy of the proposed procedure is to design a structure that would achieve, rather than be bounded by, a given loss-related metric (e.g., EAL) under the relevant site-specific seismic hazard. Moreover, the adjective "direct" refers to the ability to specify the desired level of loss as a fundamental input parameter before performing all the relevant steps to achieve such loss level with a reasonable tolerance. The loss estimation module in this procedure is based on building-level vulnerability curves defined in terms of IM (Section 2.4). Although more refined approaches with vulnerability defined in terms of EDP are available, possibly at a component-based level, the selected approach is arguably deemed more appropriate for a preliminary design phase. For example, an EDP-based loss estimation would allow explicitly considering acceleration-sensitive non-structural components. However, since their exact typology and number is most likely not known in the preliminary design phase, it is preferred to make a more generic (and flexible) assumption to embed them in the definition of the DLRs. Moreover, assumptions to include indirect losses can be embedded in the definition of DLRs, rather than needing to add specific vulnerability functions.
The proposed DLBD can be outlined in the following steps, selecting the EAL as the fundamental loss metric for design: 1. Select (e.g., 0.5% of the total reconstruction cost). This is the level of EAL that the designed structure will achieve; 2. Retrieve an appropriate set of site-specific hazard curves in terms of SA in a wide range of periods. Also, select a set of DSs (e.g., DS1-DS4 as defined in Section 2.4) relevant to the considered structural typology (e.g., RC frames).
Relatively to the ductility capacity at peak strength , which is an intermediate design parameter calculated in step 5, select reasonable guesses for the DS thresholds consistent with the qualitative definition of the DSs (e.g., = [0. 5 1 0.75 ] may be consistent with the DS definition above). The main idea is to provide values as close as possible to those obtainable from the numerical pushover analysis of the final design case. Select a relevant damageto-loss model consistent with the involved structural typology and the adopted DSs (e.g., = [7 15 50 100] % of the total reconstruction cost); 3. Select the basic geometric and material properties of the considered structure, also according to gravity-load preliminary design. For a geometrically-regular RC frame building, the parameters to be (tentatively) set are the yield stress and elastic modulus of steel ( , ), the number of storeys and bays ( , ), the centre-to-centre inter-storey height , the centre-to-centre length of the beams , the depth of beams and columns (ℎ , ℎ ), the seismic storey mass directly affecting the considered frame; 4. Select a number of combinations of yield strength and ductility capacity to define seed SDoF systems. For example, those can be 1 equally-spaced points within , and , and 2 equally-spaced points between , and , , leading to 1 2 seed SDoF systems. The range for such parameters should be carefully selected based on engineering judgement, to ensure including possible design cases. According to the relevant structural type, complete the definition of the seed SDoF systems consistent with the GP regression inputs (i.e., bilinear backbone and expected hysteresis type). For RC frames, for example, the yield displacement (Δ ) only depends on geometry, 17 and it can be calculated according to Equation 5, where is the height of storey from the ground and is the displacement shape ( = 1) assumed according to Priestley et al. 17 Therefore, yield strength ( , ) is linearly related to stiffness, and the elastic period ( ) of the seed SDoF systems is calculated via Eq. 6. The hardening ratio of the seed SDoF systems can be reasonably guessed (e.g., ℎ = 0.05 for RC frames) or alternatively included in the initial grid-based definition. Therefore, the peak normalised based shear of the SDoF system is equal to = (1 + ℎ). Finally select the most appropriate hysteresis model for the analysed structure (e.g., MTf for newly-designed RC frames); Using the herein-proposed GP regressions, calculate fragility curves for each seed SDoF system consistent with the selected DS thresholds . According to the selected and the hazard curve (interpolated based on ), use Equations 3 and 4 to calculate the building-level vulnerability curve and the EAL for each seed (see Figure 8A); 6. Select all the SDoF seeds that meet the target EAL level within a set tolerance (e.g., = 0.01 ). Clearly, the finer the ( , ) grid is, the smaller such tolerance can be. For each selected seed, run the CSM using spectra for each DS demand (code-based, as per Figure 8B, or site-specific ones) and disregard the cases not meeting the seismic demand for any DS (i.e.,

, >
). In addition, calculate the frequency of exceeding the complete damage DS, by integrating the complete damage fragility with , analogously to Equation 4, and disregard the cases above a conventionally-established threshold (e.g., between 10 -5 and 10 -4 ) 42 . The seed SDoF systems meeting the target EAL, the complete damage exceeding frequency bound and complying with the code-based displacement check are equally valid candidate design SDoF systems. One of those can be arbitrarily selected as the final design SDoF system, possibly according to design requirements not related to seismic actions; 7. Among many other alternatives, the principles of displacement-based design are herein suggested to detail each member of the structure to comply with the design SDoF's backbone and achieve a favourable plastic mechanism (e.g., plastic hinges for the base column sections and all beams' end sections, for an RC frame). For a geometrically-regular RC frame, the beams' plastic drift demand ( , ) compatible with the displacement capacity at peak force (Δ ) of the design SDoF is calculated according to Equations 7 and 8, where is the inter-storey drift profile compatible with the displacement profile Δ . On the other hand, the drift demand on the first-storey columns ( ) is simply equal to the inter-storey drift at that level. Such an approach involves a degree of error with respect to more refined, NLTHAbased approaches. However, pros and cons of displacement-based design, including the calibration of all the involved equations and assumptions are provided in ref. 17; 8. The strength demand for the members developing inelastic behaviour is computed via an equilibrium approach. For RC frames, this involves the overturning moment (OTM) equilibrium, 17 according to Equation 9, where (Δ ) = is the SDoF base shear, = ∑ Δ Δ is the effective mass, , is the base moment of column , is the frame length, , is the shear of the beam at storey and bay . Any allocation of strength to the members developing inelastic behaviour that complies with Equation 11 would satisfy the design objective. One possible design choice may involve assuming (Equation 10) that the contra-flexure point of the first-storey columns locates at 60% of their height so that capacity design of the first-storey beam-column joints is ensured. Moreover, if the beams have the same detailing in both end sections and within the frame, Equation 11 holds, where , ( , ) is the beam moment demand corresponding to its plastic drift demand. With such assumptions, it is possible to invert Equation 9 and calculate the beams' moment demand , corresponding to the target base shear (Δ ). The moment demand of the base columns, at (Δ ), can be reasonably proportioned assuming that the interior columns will be approximately twice as strong as the exterior ones, thus , = 0.6 1 (Δ ) 2( −1)+2 and , = 2 , .
9. Once the theoretical deformation and strength demand of the members developing inelastic behaviour are obtained, the structural details of such members are designed. For an RC frame, this involves the demand for beams ( , ; , ) and first-storey columns ( , ; , ). The structural details of such members can be designed via a moment-curvature approach. Generally, the provided detailing will lead to a degree of deviation from the theoretical values above. Thus, Equation 9 can be used to re-calculate the peak base shear (Δ ) provided by the design. Moreover, with a reasonable assumption for the hardening (e.g., ℎ = 0.05 for RC frames), and by knowing Δ , the entire backbone can be calculated and compared to that of the design SDoF system ( Figure 9A). Such match is particularly important since it ensures the match between the EAL of the designed frame and the target EAL, at least under the assumptions adopted for the fragility and vulnerability models. Possibly, a pushover analysis can be conducted to further verify the simplified equilibrium/compatibility-based calculations; 10. Any undesired mismatch between the SDoF and the frame backbone curves may be corrected via: (a) iterations in the member design, thus repeating step 8; (b) revising Δ or ℎ, thus restarting from step 3; (c) based on the results of a pushover analysis, revising the relative definition of with respect to the ductility capacity at peak strength , thus restarting from step 1; 11. As per any seismic design procedure involving non-linear behaviour, members intended to remain elastic must be capacity protected. Moreover, some member characteristics may be modified to comply with code-based requirements.
This procedure is demonstrated for a four-storey RC case-study building in the city of L'Aquila and designed to achieve a 0.5% EAL. According to the Italian seismic risk classification, this would correspond to a building in the A+ risk class. 43 Both and exemplified in step 1 are adopted for the design (consistently with those proposed in ref. 43 for Italian concrete buildings). The code-based Italian seismic hazard model is adopted 29 ; thus, no ad hoc site-specific probabilistic seismic hazard analysis is explicitly performed. For each point of a 5 km-spaced grid, the above model provides both PGA and spectral acceleration values for nine probabilities of exceedance in 50 years. The model provides the 16th, 50th, and 84th percentiles of the spectral accelerations for 10 period values between 0.1s and 2s, assuming rock conditions. The Italian code 30 also provides analytical approximations of the uniform hazard spectra consistent with the above model, including correction factors to obtain, among others, different soil and topography conditions. For the present illustration, a C-type soil category according to the Italian code is assumed (shear wave velocity in the first 30 m within the range 180-360 m/s).
The rectangular building is composed of four 9 m-spaced parallel frames with = 3, = 6 , = 3.6 , ℎ = 0.6 , ℎ = 0.6 , = 300 , = 200 . A storey mass approximately equal to 400 is calculated by adopting, apart from the self-weight of the structure, a superimposed dead load equal to 2 and a factored live load equal to 0.9 . Thus, the storey mass affecting one frame is = 100 . Figure 8A shows the EAL of the seed SDoF systems, defined using a vector of 100 equally-spaced points between 0.1 and 0.4 for and 100 equally-spaced points within 1.5 and 6 for , thus leading to 10000 seed SDoF systems. A linear displacement shape is assumed for this four-storey frame, 17 and the calculated yield displacement is Δ = 0.065 , with H = 10.5 . The MTf hysteresis model is used for this case study. Among the SDoF seeds meeting the target EAL (with a tolerance equal to 0.01 ), those complying with the code-based demand as per step 5 (design candidates) are represented in Figure 8B, together with the seismic demand corresponding to respectively, 30, 50, 475, and 975 year mean return periods. Those are related to the operability, damage limitation, life safety and near collapse DSs, defined consistently with HAZUS. The final design SDoF system is arbitrarily selected among those, and it is characterised by =0.97s (secant-to-yielding), = 0.27, = 2.65, ℎ = 0.05. The frequency of exceeding the complete damage DS for the design SDoF system is equal to 1.4 × 10 -6 , which is deemed appropriate according to the above-mentioned limits.
According to the displacement capacity (at DS4) of the design SDoF system (Δ = 0.17 ), steps 6 and 7 of the above procedure returned a demand equal to 325kNm(@0.8% plastic drift) for the beams, 347kNm(@0.8% plastic drift) for the exterior columns and 694kNm(@0.8% plastic drift) for the interior ones. Based on such theoretical values, cross-section design is performed characterising the flexural capacity of beams and columns via moment-curvature analysis (including gravity axial load). The model by Mander 44 is used for concrete, including calculating its ultimate strain (accounting for confinement). King et al. 45 is used for steel. Plastic hinge length is calculated according to ref. 46 Such a cross-section design required four iterations to reasonably match the SDoF design curve by appropriately balancing the (combined) strength of the members to achieve the desired OTM.
The RC frame is equipped with 0.6 × 0.3 m beams reinforced with 4 22 mm bars both in the tension and the compression sides of the cross section (340kNm yield moment). The interior columns have a 0.6 × 0.4 m cross section reinforced with 4 26mm bars on both sides of the section and 2 26 mm bars in the section centreline (729kNm yield moment). Finally, the exterior columns have a 0.5 × 0.3 m cross section equipped with 8 22mm bars equally spaced along the perimeter (336kNm yield moment). The transverse reinforcement of the listed members is composed of 12mm stirrups with 100 mm spacing, and the clear cover is equal to 30 mm. The unconfined compressive strength of the adopted concrete is equal to 30 MPa.
A pushover analysis (as per Section 3.2) is first carried out ( Figure 9A). The comparison of the numerical pushover curve (both in the refined or bi-linearised forms) to the design SDoF curve is particularly satisfactory since the curves are essentially superimposed. The eigenvalue analysis of the frame indicates a fundamental period equal to 0.93s (secant-toyielding), which is satisfactorily close to the SDoF one (0.98s). As discussed in Section 3.2, the member-level pushover results are used to identify refined DS thresholds. Figure 9A shows that the comparison of such thresholds (circular marks) with those assumed in step 1 of the procedure (triangular marks) is satisfactory since the maximum registered (displacement-wise) error is 11% (for DS4).
To thoroughly test the DLBD of this case study, NLTHA is conducted adopting the set of 100 records described in Section 2.2 (without using any scale factors). As per the frames in Section 3.2, fragility and vulnerability curves are herein derived based on the NLTHA results. In Figure 9A, such fragility curves (solid lines) are compared with those obtained via the GP regressions (dashed lines), respectively adopting the refined and simplified DS thresholds. Similarly, to the discussion in Section 3.2, the GP regressions generally under-estimate the median of the fragilities, with higher errors for more severe DSs (15.3%, 11.2%, 21.2%, and 19.3% for DS1 to DS4). The fragility dispersion obtained with the GP regressions underestimates the NLTHA-based one by 29.2% for this particular case. The fragility estimation error propagates to the vulnerability curves, with higher errors as IM increased. However, as anticipated in Section 3.2, such error levels are less relevant in the loss estimation since they correspond to hazard levels with low MAF of exceedance. This is reflected in the final calculation of EAL based on NLTHA, which is equal to 0.46% (conservatively predicted within 8% from the target value).
The result of this practical illustration, although encouraging, cannot be regarded as general. The tentative DLBD procedure proposed here is, at least arguably, worth a more comprehensive investigation/refinement/modification involving F I G U R E 9 Design SDoF case versus design MDoF case: (A) pushover and fragility curves; (B) vulnerability and hazard curves different case studies (in terms of material, lateral resisting system, geometry, etc.), different hazard levels, different loss metrics (e.g., the tail value at risk) and different target values of such metrics. Detailed comparison of design cases against refined NLTHA is required (also including 3D effects). If the procedure is proved successful, a final verification against more refined, component-by-component methodologies for loss analysis (e.g., FEMA P58, as opposed to the adopted building-level approach) is needed to estimate the confidence with which DLBD allows to set a target loss level.

Dynamic earthquake risk modelling
A flexible and reliable surrogate model for the PSDM of structures (and hence for the derivation of fragility and vulnerability curves) can enable the fast development of a high number of earthquake risk models with different scenario exposure. Among all the possible examples falling in this category, this Section demonstrates a dynamic earthquake risk model with exposure changing over time, representing the implementation of a retrofit-based seismic risk reduction policy for a region. 47,48 This illustrative example involves a synthetic portfolio of 100 RC buildings. In the as-built condition (time = 0), the buildings SDoF parameters are simulated based on the uniform distributions of ∼ (0.2 , 1.5 ), ∼ (0.1, 0.25), ℎ ∼ (0.01, 0.1), and ∼ (1.2, 3). All the buildings are characterised by the MTf hysteresis model. An earthquake risk model can be used to simulate the effects, in terms of portfolio economic losses, of a risk-mitigation policy of this kind: each householder in the area must retrofit at a time ∼ (0, 15 ), and they must increase the yield strength such that Δ ∼ (0, 0.15) and ductility capacity such that Δ ∼ (0.2, 1.5). Given the inherent freedom of the householders, the implementation of this policy may be regarded as a random retrofit process. Using the provided GP regressions, simulating such a process to obtain the distribution of the evolving portfolio losses becomes computationally feasible. This can be done according to Algorithm 1.

Algorithm 1. Dynamic earthquake risk model to simulate a portfolio retrofit process.
: number of simulations of the retrofit process; : number of buildings in the portfolio; : time;  Figure 10 shows an application of Algorithm 1, considering one single simulation of the retrofit process described above ( = 1). The details of the earthquake risk modelling methodology and related assumptions are deemed less relevant for this illustrative application. They are only briefly introduced here, while they are described in detail in ref. 49. The buildings in the portfolio are equally spaced in a rectangular grid ( Figure 10A) in the vicinity of a case-study strike-slip line fault. A 10,000-year stochastic catalogue is considered; for each event, 500 realisations of the ground motion fields are F I G U R E 1 0 (A) Simulated portfolio loss curves at t = 0 (as-built portfolio) and t = 15years (completed retrofit process); (B) Evolution of the median portfolio loss curves during the retrofit process considered. Those are expressed in terms of SA in the period range (0.2s, 1.5s) with 0.1s spacing, appropriately considering the correlation coefficients at different vibration periods. The considered fault follows a log-linear recurrence relationship (i.e. frequency of occurrence vs moment magnitude ) for moderate (non-characteristic) events (5 ≤ < 6.5), and a constant probability branch of occurrence for events with 6.5 ≤ ≤ 7. Any pulse-like feature of the ground motions is neglected.
The fragility curves for each building in the portfolio are derived by adopting the fitted GP regressions. Those are related to DSs qualitatively defined according to HAZUS and quantified by the ductility thresholds = [0.5, 1, 0.75 , ]. Conservatively, the dispersion related to the non-linear range (Eq. 2) has been assigned to every DS. DLRs equal to 2%, 10%, 43.5%, and 100% of the total reconstruction cost (DS1 to DS4) are adopted. Figure 10A shows 500 simulations of the portfolio loss curves (mean annual frequency, MAF, of loss exceedance vs loss ratio) for the as-built portfolio ( = 0) and after the full simulation of the retrofit process ( = 15 ), showing the significant reduction of losses. Figure 10B shows the time evolution of the median portfolio loss curve, together with the evolution of the portfolio EAL.

CONCLUSIONS
This paper proposed a metamodelling approach mapping the parameters controlling the dynamic behaviour of SDoF systems (i.e., force-displacement capacity curve, hysteretic behaviour) and their probabilistic seismic demand model (i.e., EDP vs IM distribution). This metamodel allows rapidly calculating fragility curves of SDoF representation of structures to be used in seismic risk and/or loss analyses. The selected approach is the GP regression since it does not require any a priori definition of the output functional form (a GP is a non-parametric model), and therefore, it results in an infinitely-scalable surrogate model. The dataset used to train the GP regressions is based on a Monte Carlo sampling of 10,000 SDoF systems, each analysed via a cloud-based NLTHA using 100 ground-motion records. The code to make predictions using the proposed set of GP regressions is freely available (https://github.com/robgen/surrogatedPSDM), together with the adopted training datasets. With the same code, the users are also allowed to filter, extend or modify the training dataset and customise the fitting.
The results of this study can be highlighted as follows: • Within the training dataset, the GP predictions for the slope of the PSDMs show an NRMSE equal to 2.5%, while this is equal to 6.2% for the logarithmic standard deviation. Such error estimates increase by 0.1% if a ten-fold cross validation test is conducted. Estimates of the median of the fragility curves consistent with DS3 and DS4 show an NRMSE respectively equal to 2.0% and 2.2%, while the error on the fragility dispersion remains equal to 6.2%. • By selecting two baseline case-study SDoF systems, the distribution of the above fragility errors is propagated to the EAL. Simulating 30,000 realisations of the fragility parameters affected by such errors, the resulting simulated-over-baseline EAL error follows a Normal distribution with mean equal to 1 and standard deviation equal to 0.015 and 0.01, respectively for the two case studies. • The GP-based fragility predictions are tested against eight realistic RC frames with different levels of seismic performance. The NRMSE of such predictions is equal to 24.0%, 26.0%, 17.7%, and 19.0% for the DS1, DS2, DS3, and DS4 medians, respectively, if compared to NLTHA-based results. On the other hand, the inherent simplicity of the metamodels leads to a higher estimation error for the fragility dispersion, despite the NRMSE is equal to 29%. • Considering the various sources of uncertainty typically involved in the seismic performance or risk/loss models, often not captured due to the simplified nature of the models themselves, the error levels introduced by using the proposed GP regressions are deemed satisfactory for practical applications, especially considering the low modelling effort and time required for the GP-based predictions. • The proposed metamodels enable the development of a direct loss based design, for which a tentative procedure is proposed. This allows designing structures complying with a given target EAL, which is a fundamental input selected by the designer. The procedure consists of two phases: (1) the GP regressions are used to quickly select a design SDoF system complying with both the selected EAL target and the code-based demand at different DSs; (2) the members composing the lateral resisting system of the structure are designed (via equilibrium and compatibility principles) by matching the force-displacement curve of the designed structure with a target SDoF backbone. A realistic 2D RC frame is designed to achieve 0.5% EAL. Adopting an NLTHA-based approach and a building-level damage-to-loss model, the same parameter is measured to be 0.46%, thus showing an 8% error. Although this result is encouraging, the procedure requires further validations, refinement and/or modifications before it can be regarded generality. • The proposed GP-regression approach also enables the fast development of a high number of earthquake risk models to be used in scenario-based decision making. This paper demonstrates a dynamic-exposure earthquake risk model representing the implementation of a retrofit-based seismic risk reduction policy for a region, in which householders are required to improve the performance of their buildings within a given time frame. By using the provided GP regressions, simulating such a random retrofit process becomes computationally feasible, thus allowing to obtain the distribution of the evolving portfolio losses. A seismic risk-mitigation policy implementation for a synthetic portfolio of 100 buildings is analysed for illustration purposes.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available at https://github.com/robgen/surrogatedPSDM.