Quick Start: Basic Simulation Example#

This example demonstrates how to calibrate a stochastic model to a recorded ground motion. We utilize the ModelInverter class to obtain a fitted StochasticModel, which is then used to generate synthetic ground motion simulations.

Step 1: Import Libraries and Load the Target Ground Motion#

Import necessary sgsim classes along with necessary libraries. Load an existing accelerogram record to serve as the target for the simulation.

# %% import libraries
import numpy as np
import matplotlib.pyplot as plt
import time
from sgsim import GroundMotion, ModelInverter, Functions

# %% Prepare the target ground motion
# Change the path to the file
gm = GroundMotion.load_from(source='nga', file='RSN123_Example.AT2')
# If necessary, preprocess the ground motion (e.g., trimming, baseline correction, tapering and filtering), otherwise skip this step
gm = gm.trim_by_energy((0.001, 0.999)).taper(0.05).butterworth_filter((0.05, 100.0))
# An alterantive to triming gm for just model inversion is to use fit_range in the inverter.fit() method to specify the time range for fitting (e.g., fit_range=(0.01, 0.99) for 1% to 99% cumulative energy)

Step 2: Perform Model Inversion#

Use the ModelInverter class to fit a stochastic model to the target ground motion. Use Functions to define the functional forms for the model parameters (amplitude, frequency, and damping evolution with time).

# %% Specify the stochastic model functional forms
q = Functions.BetaCentroidSpread()  # Modulating function
wu = Functions.Constant()           # Upper frequency
zu = Functions.Constant()           # Upper damping
wl = Functions.Constant()           # Lower frequency
zl = Functions.Constant()           # Lower damping

inverter = ModelInverter(ground_motion=gm, modulating=q,
                         upper_frequency=wu, upper_damping=zu,
                         lower_frequency=wl, lower_damping=zl)

# Optional: Provide initial guesses and fix parameters using tight bounds
# Here we fix upper_damping to 0.707 and lower_damping to 1.0 by making min/max bounds equal
ig = {
    "upper_frequency": [10.0],
    "upper_damping": [0.707],
    "lower_frequency": [0.3],
    "lower_damping": [1.0],
}

bs = {
    "upper_frequency": [(0.1, 40.0)],
    "upper_damping": [(0.707, 0.707)], # Fixed parameter
    "lower_frequency": [(0.01, 0.99)],
    "lower_damping": [(1.0, 1.0)],     # Fixed parameter
}

# Fit to the target ground motion and measure execution time
start = time.time()
model = inverter.fit(initial_guess=ig, bounds=bs, fit_range=(0.01, 0.99), mode='stepwise')  # stepwide mode is faster and lead to unique solution, while joint mode is more robust but slower and may lead to multiple local minima. Adjust fit_range as needed to focus on specific time windows of the ground motion.
end = time.time()

print(f"Execution time: {end - start:.2f} seconds")

# Print a summary of the fitted parameters
model.summary()

Step 3: Generate Simulations#

Generate synthetic ground motions using the fitted stochastic model parameters. The resulting simulations are returned as a GroundMotion object.

Step 4: Visualize the Results#

You can use standard libraries like matplotlib to visualize the model, simulations, and the target.

Tip

See Quick Start: Basic Ground Motion Analysis Example for plotting response spectra and other intensity measures.

# 1. Compare Acceleration Time Series
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(10, 6), gridspec_kw={'hspace': 0.15})

# Target Ground Motion
axes[0].plot(gm.t, gm.ac, color='black', lw=1, label='Target')
axes[0].plot(model.t, model.q, color='red', ls='--', lw=1.2, label='Model Envelope')
axes[0].plot(model.t, -model.q, color='red', ls='--', lw=1.2)
axes[0].set_ylabel('Acceleration (g)')
axes[0].legend(frameon=False)
axes[0].set_title('Target vs Simulation - Acceleration')

# First Simulation
axes[1].plot(sm.t, sm.ac[0], color='tab:blue', lw=1, label='Simulation #1')
axes[1].set_ylabel('Acceleration (g)')
axes[1].set_xlabel('Time (s)')
axes[1].legend(frameon=False)

plt.show()

# 2. Compare Displacement Time Series
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(10, 6), gridspec_kw={'hspace': 0.15})

# Target Ground Motion
axes[0].plot(gm.t, gm.disp, color='black', lw=1, label='Target')
axes[0].set_ylabel('Displacement (cm)')
axes[0].legend(frameon=False)
axes[0].set_title('Target vs Simulation - Displacement')

# First Simulation
axes[1].plot(sm.t, sm.disp[0], color='tab:blue', lw=1, label='Simulation #1')
axes[1].set_ylabel('Displacement (cm)')
axes[1].set_xlabel('Time (s)')
axes[1].legend(frameon=False)

plt.show()

# 3. Compare Fourier Amplitude Spectra (FAS)
plt.figure(figsize=(7, 4))
plt.loglog(sm.freq, sm.fas.T, color='gray', alpha=0.3, lw=0.7)  # 20 FAS simulations
plt.loglog(sm.freq, np.mean(sm.fas, axis=0), color='red', ls='--', lw=2, label='Simulation Mean')
plt.loglog(gm.freq, gm.fas, color='black', lw=1.5, label='Target')
plt.loglog(model.freq / (2 * np.pi), model.fas, color='Blue', lw=1.5, label='Model')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Fourier Amplitude')
plt.legend(frameon=False)
plt.grid(True, which="both", ls=":", alpha=0.3)
plt.title('Fourier Amplitude Spectra')
plt.ylim(1e-4, 1e1)  # if necessary, adjust y-axis limits for better visualization
plt.show()

# 4. Compare Response Spectra (5% Damping)
tp = np.arange(0.0, 10.0, 0.1)
_, _, sa = gm.response_spectra(tp)
_, _, sm_sa = sm.response_spectra(tp)

plt.figure(figsize=(7, 4))
plt.loglog(tp, sm_sa.T, color='gray', alpha=0.3, lw=0.7)
plt.loglog(tp, np.mean(sm_sa, axis=0), color='red', ls='--', lw=2, label='Simulation Mean')
plt.loglog(tp, sa, color='black', lw=1.5, label='Target')
plt.xlabel('Period (s)')
plt.ylabel('Spectral Acceleration (g)')
plt.legend(frameon=False)
plt.grid(True, which="both", ls=":", alpha=0.3)
plt.title('Response Spectra (5% Damping)')
plt.show()

Step 6: Export Simulations#

Tip

See Quick Start: Basic Ground Motion Analysis Example for saving and plotting results. The simulated ground motions are GroundMotion objects and can be saved or processed similarly.