Gaussian Process (GP) Surrogate Modeling

Note

Surrogate modeling functionality of R2D is built upon GPy library (available under BSD 3-clause license), an opensource Python framework for Gaussian process modeling developed in the Sheffield machine learning group.

The Train GP Surrogate Model module is used to construct a Gaussian process (GP) based surrogate model that substitutes expensive computational simulation models. The challenge of surrogate modeling in earthquake engineering arrives from the stochasticity in the ground motion time history and corresponding stochastic output.

Surrogate Modeling in EE-UQ

• Output of surrogate model: Following the performance-based earthquake engineering practice, the response estimation module in SimCenter workflow concerns only the peak responses (e.g. peak acceleration, peak drift ratio at each floor) rather than the whole time history. Therefore, the surrogate modeling algorithms in this tool are designed to predict the non-time history responses.

• Input of surrogate model: To keep the input dimension of the surrogate model to a reasonable scale, the ground motion time history is first parameterized in the algorithm into several key features. This can be, for example, the source properties (e.g. M, R) or time/frequency domain excitation property at the site of interest (e.g. intensity measures) ([Kyprioti2021], [Zhong2023]).

The surrogate prediction of EDP given the ground motion parameters typically inherits significant randomness, as the amount of information provided by these parameters is not sufficient to determine the exact model peak response. Therefore, the key task in surrogate modeling is to capture the remaining amount of uncertainty in the EDP predictions, as shown in the above figure.

The trained surrogate model can be saved in a portable file (.json format) and later imported in place of the original simulation model for UQ analysis or optimization purposes. See user guide section 2.2.5 for how to import the surrogate model in R2D.

Kyprioti2021

Kyprioti, A. P., & Taflanidis, A. A. (2021). “Kriging metamodeling for seismic response distribution estimation”. Earthquake Engineering & Structural Dynamics, 50(13), 3550-3576.

Zhong2023

Zhong, K., Navarro, J.G., Govindjee, S., and Deierlein, G.G. (2023). “Surrogate Modeling of Structural Seismic Response Using Probabilistic Learning on Manifolds,” Earthquake Engineering & Structural Dynamics, 1-22, https://doi.org/10.1002/eqe.3839

Input description

Case 1: Sampling and Simulation

When the Training Dataset option is set to the Sampling and Simulation, a simulation model should be presented in the later tabs. The training points are sampled adaptively by the design of experiments.

Input panel for surrogate modeling

• Maximum Number of Model Runs: When the number of simulation runs reaches the limit, the analysis will be terminated.

• Maximum Computation Time (in minutes): When the tolerance limit of the computation time is reached, the analysis will be terminated. There will be a few minutes of error.

• Target Accuracy (Normalized Error): The target accuracy is defined in terms of normalized root-mean-squared error (NRMSE) estimated by leave-out-one cross-validation (LOOCV).

  ()\[\begin{align*} \rm{NRMSE} ~ &= \frac{\sqrt{\frac{1}{N} \sum^{N}_{k=1} (y_k-\hat{y}_k)^2}}{\max_{k=1,...,N}(y_k)-\min_{k=1,...,N}(y_k)} \end{align*}\]

  where

  \(y_k\): exact response from the model simulation

  \(\hat{y}_k\): estimated response by LOOCV surrogate model prediction

  \(N\): number of samples used to train the surrogate model

• Random Seed: Seed of the random number generator

• Parallel execution: This engine implemented multiprocessing (local) or mpi4py (remote) python packages to run parallel execution.

Note that the results from the parallel and serial run may not be exactly the same because parallel execution sets the number of batch design of experiments (DoE) in order to maximize the use of resources (Default DoE interval: 5). User can also activate the following advanced options:

Advanced Options for Gaussian Process Model

Sampling and Simulation - Case 1

Start with Existing Dataset

RV (Random Variables) Tab

The bounds of RVs need to cover the domain of interest in future applications, while it should not be unnecessarily stretched. Input type and Distribution should be set to Parameters and Uniform. When a dataset is provided, make sure to match the order of RVs in the RV tab to the order of data columns. Any correlation values will be ignored.

Example of RV tab

Output description

Goodness-of-Fit

Once the training is completed, the following three verification measures are presented based on leave-one-out cross-validation (LOOCV) error estimation.

• Leave-one-out cross-validation (LOOCV)

  LOOCV prediction \(\hat{\boldsymbol{y}}_k\) at each sample location \(\boldsymbol{x}_k\) is obtained by the following procedure: A temporary surrogate model \(\hat{\boldsymbol{y}}=f^{sur}_{loo,k}(\boldsymbol{\boldsymbol{x}})\) is constructed using the samples \(\{\boldsymbol{x}_1,\boldsymbol{x}_2,...,\boldsymbol{x}_{k-1},\boldsymbol{x}_{k+1},...,\boldsymbol{x}_N\}\) and the calibrated parameters, and the prediction \(\hat{\boldsymbol{y}}_k=f^{sur}_{loo,k}(\boldsymbol{x}_k)\) is compared with the exact outcome.

We provide different verification measures for two different cases.

1. When nugget variance is low: The LOOCV prediction \(\hat{\boldsymbol{y}}_k\) is expected to match the exact outcome \(\boldsymbol{y_k}=f(\boldsymbol{x}_k)\) when the surrogate model is well-trained. To quantify the goodness, R2 error, normalized root-mean-squared-error (NRMSE), and correlation coefficient are provided:

   • R2 error

     R2 error is defined in terms of the total sum of squares over the residual sum of squares

     ()\[\begin{align*} R^2 &= 1 - \frac{\sum^N_{k=1} (\hat{y}_k-\mu_\hat{y})^2}{\sum^N_{k=1} (\hat{y}_k-y_k)^2} \end{align*}\]

     The surrogate model is considered well-trained when the R2 (<1) approaches 1

   • Normalized root-mean-squared-error (NRMSE)

     ()\[\begin{align*} \rm{NRMSE} ~ &= \frac{\sqrt{\frac{1}{N_t} \sum^{N_t}_{k=1} (y_k-\hat{y}_k)^2}}{\max_{k=1,...,N_t}(y_k)-\min_{k=1,...,N_t}(y_k)} \end{align*}\]

     The surrogate model is considered well-trained when the NRMSE (>0) approaches 0

   • Correlation coefficient

     Correlation coefficient is a statistic that measures the linear correlation between two variables
     ()\[ \rho_{y,\hat{y}} = \frac{\sum^N_{k=1}(y_k-\mu_{y})(\hat{y}_k-\mu_{\hat{y}})} {\sigma_y \sigma_\hat{y}}\]
     where

     \(\mu_{y}\) : mean of \(\{y_k\}\)

     \(\mu_{\hat{y}}\): mean of \(\{\hat{y}_k\}\)

     \(\sigma_{y}\): standard deviation of \(\{y_k\}\)

     \(\sigma_{\hat{y}}\): standard deviation of \(\{\hat{y}_k\}\)
     The surrogate model is considered well-trained when the correlation coefficient ( \(-1<\rho<1\) ) approaches 1

2. When nugget variance is high: The distance between LOOCV prediction \(\hat{\boldsymbol{y}}_k\) and the exact outcome \(\boldsymbol{y_k}=f(\boldsymbol{x}_k)\) is expected to follow a normal distribution when the surrogate model is well-trained. To quantify the goodness, inter-quartile ratio (IQR) and Cramer-Von Mises statistics can be evaluated:

• Inter-quartile ratio (IQR): IQR provides the ratio of the sample QoIs that lies in 25-75% LOOCV prediction bounds (interquartile range). The IQR values should theoretically approach 0.5 if the prediction is accurate.

• Cramer-Von Mises statistics: Cramer-Von Mises statistics calculate the normality score. A GP assumes that the observations follow a normal distribution conditional on the input parameters. To assess the normality of the model predictions, the difference between the mean prediction \(\hat{y}_k\) and the sample observation \(y_k\) value is divided by the standard deviation prediction from surrogate \(\hat{\sigma}_{y,k}\):

  ()\[ u_k = \frac{y_k-\hat{y}_k} {\hat{\sigma}_{y,k}}\]

If the values of \({u_k}\) follow the standard normal distribution, the resulting surrogate model may be considered well-constructed. The Cramer-Von Mises test is calculated using the scipy.stats.cramervonmises function in the Python package Scipy, and the resulting p-value is displayed. Conventionally, if the p-value exceeds a significance threshold, e.g. 0.05, the null hypothesis that the samples are from a normal distribution is not rejected, meaning the samples may be considered to follow a Gaussian distribution.

Additionally, scatter plot between the predicted and exact responses is presented: Well-trained model will form a clear diagonal line in case (i), or be distributed normally around the diagonal line in case (ii). A poorly trained model will have points that are more scattered around.

Well-trained surrogate (left) and poorly trained surrogate (right) models

Note

Since these validation measures are calculated from the cross-validation predictions, they can be biased, particularly when a highly localized nonlinear range exists in the actual response surface, and those regions are not covered by the training samples. The introduction of adaptive design of experiments helps to suppress the bias by enabling the targeted selection of simulation points around potentially faulty regions.

Warning

Note that GP-based surrogate models can be used to fit only smooth, continuous functions. When the surrogate model is poorly trained, a parametric study is highly recommended to check for any possible discontinuity presented in the simulation model.

Saving Options

• Save GP Model: The constructed surrogate model is saved. One main file and one auxiliary folder will be saved.

  • SurroateGP Info File (default name: SimGpModel.json): This file contains information required to quickly reconstruct the surrogate model and predict the response for different input realizations. This can be later imported into EEUQ.

  • Simulation template directory (default name: tmplatedir_SIM): This folder contains all the scripts and commands to run the original dynamic time history analysis. This folder can later be imported into EEUQ along with the surrogate model to alternate between original simulations and surrogate predictions or compare the surrogate predictions to the response of the original model.

  • Only for the multi-fidelity case in quoFEM, additional SurroateGP model file (default name: SimGpModel.pkl) will be written along with the above two.

Warning

Do not place the above surrogate model files in your root, downloads, or desktop folder as when the application runs it will copy the contents on the directories and subdirectories containing these files multiple times. If you are like us, your root, Downloads, or Documents folders contain an awful lot of files and when the backend workflow runs you will slowly find you will run out of disk space!

• Save GP Info: This is a report generated for user reference. It contains the GP model parameter and other information. The default file name is GPresults.out.

• RV Data, QoI Data:It saves the samples of RV and QoI. The default file names are X.txt and Y.txt, respectively. IMPORTANT: To continue surrogate modeling with additional simulations, save the RV and QoI sample files using this button and import them as initial points. Refer to the ‘Start with Existing Dataset’ option in Case 1.

Saving options

Example outputs from saving options

Patsialis2021

Patsialis, D., and A. A. Taflanidis. (2021). Multi-fidelity Monte Carlo for seismic risk assessment applications. Structural Safety 93: 102129.