Generating Fragility Curves for Seismic Risk Assessment

In the field of performance-based earthquake engineering, fragility curves are indispensable tools. These probabilistic functions describe the likelihood that a structure or component will exceed a certain damage state when subjected to varying levels of seismic intensity—typically measured using Spectral Acceleration (Sa) or Peak Ground Acceleration (PGA).

At QuakeLogic, we help researchers, engineers, and emergency planners develop and automate fragility analysis as part of seismic resilience studies. This blog post introduces a fully open-source Python script developed by our team that reads structural response data and generates fragility curves for multiple damage states.

What Are Fragility Curves?

A fragility curve is a cumulative distribution function (CDF) that defines the conditional probability:

P(Damage ≥ State | Intensity Measure)

This means that, for a given level of shaking (e.g., 0.3g Sa), the curve tells you the probability that a structure will exceed a limit-state such as Moderate or Complete damage.

We typically use lognormal distributions to model fragility because they capture the exponential rise in damage probability with increasing shaking intensity.

Engineering Demand Parameters (EDPs) and Intensity Measures (IMs)

In seismic fragility analysis, the relationship is between:

• EDP (Engineering Demand Parameter): A measurable response of the structure.
• IM (Intensity Measure): A scalar value representing the strength of ground shaking.

Common EDPs:

• Maximum interstory drift ratio (MIDR) – most widely used for buildings
• Peak floor acceleration
• Residual drift
• Relative displacement
• Roof displacement
• Member strains or rotations (for local component fragility)

Common IMs:

• PGA (Peak Ground Acceleration)
• Sa(T1) – Spectral acceleration at fundamental period (T1)
• PGV (Peak Ground Velocity)
• Arias Intensity
• CAV (Cumulative Absolute Velocity)

Choosing the right IM and EDP depends on the structure type, the available data, and the analysis objective.

Understanding the Fragility Curve Plot

The plot generated by the code shows the fragility curves for multiple damage states. Here’s how to read it:

• X-axis: Intensity Measure (IM)
  In our example, this is Sa(T1) (spectral acceleration), but it can be PGA, PGV, or any other meaningful measure of shaking intensity.
• Y-axis: Probability of Exceedance
  This is the probability that the structure’s response (EDP) will exceed the defined threshold for each damage state.

Interpretation:

A steeper curve implies a rapid transition from low to high damage probability — common for brittle systems or clearly defined performance thresholds.Methodology Overview

As the intensity of shaking (IM) increases, the probability of damage naturally increases.

Fragility curves slope upwards — from 0 to 1 — showing this increasing probability.

Robust structures have curves that shift rightward, indicating a lower probability of damage at a given IM.

1. Input Data: Simulated or observed structural responses (e.g., drift ratios) vs. seismic intensities.
2. Damage States: Thresholds for damage classification (e.g., 1%, 3%, 4%, 5% drift).
3. Binning: The Sa values are grouped and used to calculate the fraction of samples exceeding each damage threshold.
4. Curve Fitting: We fit a lognormal CDF to the exceedance probabilities using non-linear regression.
5. Visualization: Fragility curves are plotted for each damage state.

Example Output

Input Format (input.txt)

A simple tab-separated file with the following columns:

Sa_T1_g    MaxDrift_ratio
0.1        0.0043
0.15       0.0098
...        ...

You can download an example file here.

Python Code: Fragility Curve Generator

################################################################################
# Fragility Curve Generator for Seismic Risk Analysis
#
# This script calculates and visualizes fragility curves for structural
# components or systems based on nonlinear response data from a suite of
# ground motions.
#
# Fragility curves represent the conditional probability that a structure
# will exceed a predefined damage state for a given ground motion intensity
# measure (IM), typically spectral acceleration (Sa). These curves are
# essential tools for probabilistic seismic risk assessment and
# performance-based earthquake engineering.
#
# The script uses a lognormal cumulative distribution function (CDF) to fit
# observed probabilities of exceedance derived from drift thresholds and
# structural response data. Drift ratios are compared against damage state
# thresholds to calculate these probabilities at various IM bins.
#
# Generated by: QuakeLogic Inc.
# Contact: support@quakelogic.net
################################################################################

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import lognorm
from scipy.optimize import curve_fit

def run():
    try:
        df = pd.read_csv("input.txt", sep='t')
    except FileNotFoundError:
        print("❌ input.txt not found. Please ensure it exists in the working directory.")
        return

    IMs = df['Sa_T1_g'].values
    EDPs = df['MaxDrift_ratio'].values

    damage_states = {
        'Slight': 0.01,
        'Moderate': 0.03,
        'Extensive': 0.04,
        'Complete': 0.05
    }

    fragility_results = {}
    unique_IMs = np.linspace(min(IMs), max(IMs), 20)
    tolerance = 0.05

    for ds, threshold in damage_states.items():
        prob_exceed = []
        for im in unique_IMs:
            edp_subset = EDPs[np.abs(IMs - im) < tolerance]
            if len(edp_subset) > 0:
                prob = np.mean(edp_subset >= threshold)
            else:
                prob = 0
            prob_exceed.append(prob)

        def lognorm_cdf(x, mu, sigma):
            return lognorm.cdf(x, s=sigma, scale=np.exp(mu))

        prob_array = np.array(prob_exceed)
        if np.all(prob_array == 0) or np.all(prob_array == 1):
            print(f"⚠️ Skipping '{ds}' — flat probability distribution (all 0s or all 1s)")
            fragility_results[ds] = (np.nan, np.nan)
            continue

        try:
            popt, _ = curve_fit(lognorm_cdf, unique_IMs, prob_exceed, p0=[np.log(np.mean(unique_IMs)), 0.5])
            mu, sigma = popt
            fragility_results[ds] = (mu, sigma)
            print(f"✔️ Fitted '{ds}': mu = {mu:.3f}, sigma = {sigma:.3f}")
        except Exception as e:
            print(f"❌ Failed to fit for '{ds}': {e}")
            fragility_results[ds] = (np.nan, np.nan)

    plt.figure(figsize=(10, 6))
    x_vals = np.linspace(min(IMs), max(IMs), 300)
    for ds, (mu, sigma) in fragility_results.items():
        if np.isnan(mu): continue
        plt.plot(x_vals, lognorm.cdf(x_vals, s=sigma, scale=np.exp(mu)), label=ds)

    plt.ylim(0, 1)
    plt.xlabel("Spectral Acceleration Sa(T1) [g]", fontsize=12)
    plt.ylabel("Probability of Exceedance", fontsize=12)
    plt.title("Seismic Fragility Curves by Damage State", fontsize=14)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.legend(title="Damage State")
    plt.tight_layout()
    plt.savefig("fragility_curves_combined.png", dpi=300)
    plt.show()
    print("📈 Fragility plot saved as 'fragility_curves_combined.png'")

if __name__ == '__main__':
    run()

Summary

This fragility analysis tool allows engineers and researchers to:

• Automate fragility curve generation from drift or response data
• Visualize performance thresholds for multiple damage states
• Incorporate results into risk models or decision-making tools

Need Help?

QuakeLogic Inc. provides advanced tools and consulting for:

• Earthquake early warning systems
• Structural health monitoring (SHM)
• Fragility model development
• Shake table testing and seismic instrumentation

📧 Reach out to us anytime: support@quakelogic.net