04 Sensitivity Analysis

This example demonstrates how to perform sensitivity analysis on CFAST fire simulation models using the SALib library. We’ll use Sobol indices to quantify parameter importance and visualize the results.

Step 1: Import Required Libraries

We’ll import:

• NumPy: For generating parameter ranges and arrays

• Pandas: For organizing and analyzing simulation results

• Matplotlib: For visualizing the generated data

• SALib: For performing sensitivity analysis (see sample and analyze)

• PyCFAST: For parsing and running CFAST models (see parse_cfast_file and run)

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.ticker import FormatStrFormatter
from SALib.analyze import sobol as sobol_analyze
from SALib.sample import sobol as sobol_sample

from pycfast.parsers import parse_cfast_file

Step 2: Load Base Model

For this example, we’ll use a predefined CFAST input file as our base model. You can replace this with your own model file. We’ll use parse_cfast_file to load the USN_Hawaii_Test_03.in input file.

model = parse_cfast_file("data/USN_Hawaii_Test_03.in")

The parsed model is displayed below.

print(model.summary())

Model: USN_Hawaii_Test_03_parsed.in
Simulation: 'HawaiiTest 3' (570s)

Components:
  Material Properties (2):
    Material 'CONCRETE' (Concrete Normal Weight (6 in)): k=1.6, ρ=2400.0, c=0.75, t=0.15, ε=0.94
    Material 'STEELSHT' (Steel Plain Carbon (1/16 in)): k=60.0, ρ=7850.0, c=0.48, t=0.0015, ε=0.9
  Compartment (1):
    Compartment 'Bay 1': 97.6x74.0x14.8 m, volume: 106891.52 m³ (ceiling: CONCRETE, wall: CONCRETE, floor: CONCRETE)
  Fire (1):
    Fire 'Hawaii_03' (Hawaii_03_Fire) in 'Bay 1' at (36.7, 39.9) (peak: 1 kW, duration: 11min, χr: 0.4)
  Device (1):
    Target 'Targ 1' (PLATE) in 'Bay 1' at (36.7, 39.9, 14.7) (material: STEELSHT, depth: 0.00075m)

Step 3: Define the Problem for Sensitivity Analysis

We specify which parameters to analyze and their realistic ranges. For this example, we focus on four key parameters:

1. Heat of combustion: Energy released per unit mass of fuel (affects fire intensity)

2. Radiative fraction: Fraction of fire energy released as radiation (affects heat transfer)

3. Target thickness: Material thickness of temperature measurement targets

4. Target emissivity: Surface radiation properties of targets

Each parameter gets a range of realistic values to explore during the analysis.

problem = {
    "num_vars": 4,
    "names": [
        "heat_of_combustion",
        "radiative_fraction",
        "target_thickness",
        "target_emissivity",
    ],
    "bounds": [
        [100, 50000],
        [0.1, 0.8],
        [0.10, 0.60],
        [0.8, 0.95],
    ],
}

Below you can see the defined parameters and their ranges for the sensitivity analysis.

for i, name in enumerate(problem["names"]):
    bounds = problem["bounds"][i]
    print(f"{name}: [{bounds[0]}, {bounds[1]}]")

heat_of_combustion: [100, 50000]
radiative_fraction: [0.1, 0.8]
target_thickness: [0.1, 0.6]
target_emissivity: [0.8, 0.95]

Step 4: Generate Parameter Samples

We use Sobol sequences to create well-distributed parameter combinations. The number of samples affects accuracy but also computational time.

N = 64  # Number of samples (will generate N*(2*num_vars+2) total samples)
param_values = sobol_sample.sample(problem, N)

print(f"Generated {len(param_values)} parameter combinations")
print(f"Sample shape: {param_values.shape}")

df_samples = pd.DataFrame(param_values, columns=problem["names"])

Generated 640 parameter combinations
Sample shape: (640, 4)

First 5 parameter combinations

df_samples.head()

	heat_of_combustion	radiative_fraction	target_thickness	target_emissivity
0	28417.103877	0.301596	0.151390	0.822629
1	30359.977065	0.301596	0.151390	0.822629
2	28417.103877	0.215610	0.151390	0.822629
3	28417.103877	0.301596	0.563827	0.822629
4	28417.103877	0.301596	0.151390	0.873541

Step 5: Run Model Evaluations

Now we run CFAST simulations for each parameter combination. This is the most time-consuming step as we need to:

• Modify the model with each parameter set using update_fire_params and update_material_params

• Run the CFAST simulation with run

• Extract and store the output values

outputs = []

for i, params in enumerate(param_values):
    if i % 50 == 0:
        print(
            f"Sample {i}/{len(param_values)}: "
            f"hoc={params[0]}, rf={params[1]}, thickness={params[2]}, emissivity={params[3]}"
        )

    temp_model = model.update_fire_params(
        fire="Hawaii_03_Fire",
        heat_of_combustion=params[0],
        radiative_fraction=params[1],
    )
    temp_model = temp_model.update_material_params(
        material="STEELSHT",
        thickness=params[2],
        emissivity=params[3],
    )

    results = temp_model.run()
    max_target_temp = results["devices"]["TRGSURT_1"].max()

    outputs.append(max_target_temp)

Y = np.array(outputs)
print(f"Completed {len(outputs)} model evaluations")
print(f"Output shape: {Y.shape}")

Sample 0/640: hoc=28417.103877384216, rf=0.3015960895456374, thickness=0.15139032937586308, emissivity=0.8226289795245976
Sample 50/640: hoc=16678.96753353998, rf=0.6153631545603275, thickness=0.19618909191340209, emissivity=0.8606753397732974
Sample 100/640: hoc=15285.871427785605, rf=0.3359390974976123, thickness=0.5748151544481516, emissivity=0.8039458022918553
Sample 150/640: hoc=26928.842971939594, rf=0.5471831217408181, thickness=0.5226884795352816, emissivity=0.8419867265969515
Sample 200/640: hoc=46281.35332502425, rf=0.35854666652157907, thickness=0.4575213006697595, emissivity=0.8481106677092612
Sample 250/640: hoc=3360.743639897555, rf=0.6461753271520138, thickness=0.13369938489049674, emissivity=0.8453330488875509
Sample 300/640: hoc=9506.77813719958, rf=0.2790872951038182, thickness=0.2842662594281137, emissivity=0.8669409015215933
Sample 350/640: hoc=40886.01847086102, rf=0.7488805322907865, thickness=0.24182693092152477, emissivity=0.8511400885414332
Sample 400/640: hoc=32423.039924446493, rf=0.32799937920644884, thickness=0.34130285158753393, emissivity=0.8993061103392392
Sample 450/640: hoc=19123.613521363586, rf=0.5556386437267066, thickness=0.25817580092698333, emissivity=0.9361760112456977
Sample 500/640: hoc=20542.33846226707, rf=0.14300404693931343, thickness=0.5036803366616368, emissivity=0.8579021774698049
Sample 550/640: hoc=33739.35492550954, rf=0.784896035026759, thickness=0.5489695839583874, emissivity=0.820952822547406
Sample 600/640: hoc=42295.51629126072, rf=0.2868569494225085, thickness=0.5195466446690261, emissivity=0.9300749773159622
Completed 640 model evaluations
Output shape: (640,)

Step 6: Perform Sobol Sensitivity Analysis

We calculate Sobol indices to quantify parameter importance:

• First-order indices (S1): Direct effect of each parameter alone

• Total-order indices (ST): Total effect including interactions with other parameters

• Interaction effects (ST - S1): How much parameters interact with each other

Higher values indicate greater influence on the simulation output.

Si = sobol_analyze.analyze(problem, np.array(Y))

First-order indices (S1):

for i, param_name in enumerate(problem["names"]):
    print(f"{param_name}: {Si['S1'][i]:.4f}")

heat_of_combustion: 0.3237
radiative_fraction: 0.2578
target_thickness: 0.1758
target_emissivity: 0.0008

Total-order indices (ST):

for i, param_name in enumerate(problem["names"]):
    print(f"{param_name}: {Si['ST'][i]:.4f}")

heat_of_combustion: 0.5178
radiative_fraction: 0.2998
target_thickness: 0.1940
target_emissivity: 0.0261

Step 7: Visualize Sensitivity Results

The bar charts show:

• Left panel: First-order effects (direct parameter influence)

• Right panel: Total effects (including parameter interactions)

• Error bars: Confidence intervals for the estimates

Parameters with higher bars have greater influence on the simulation output. Large differences between S1 and ST indicate significant parameter interactions.

plt.style.use("seaborn-v0_8-whitegrid")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)

names = problem["names"]
x = np.arange(len(names))

S1 = np.asarray(Si["S1"])
S1_conf = np.asarray(Si["S1_conf"])
ST = np.asarray(Si["ST"])
ST_conf = np.asarray(Si["ST_conf"])

cmap = plt.get_cmap("tab10")
bar_width = 0.6

ymax_candidate = np.nanmax([S1 + S1_conf, ST + ST_conf])
ymax = (
    float(np.clip(ymax_candidate + 0.06, 0, 1.0))
    if not np.isnan(ymax_candidate)
    else 1.0
)
ymin = 0.0

ax1.bar(x, S1, width=bar_width, color=cmap(0), edgecolor="black", alpha=0.9)
ax1.errorbar(
    x, S1, yerr=S1_conf, fmt="none", ecolor="black", capsize=5, elinewidth=1.25
)
ax1.set_title("First-order Sensitivity Indices (S1)", fontsize=14)
ax1.set_ylabel("Sensitivity Index", fontsize=12)
ax1.set_xlabel("Parameter", fontsize=12)
ax1.set_xticks(x)
ax1.set_xticklabels(names, rotation=45, ha="right", fontsize=10)
ax1.set_ylim(ymin, ymax)
ax1.grid(which="major", axis="y", linestyle="--", alpha=0.4)

for i, v in enumerate(S1):
    if np.isnan(v):
        txt = "nan"
    else:
        txt = f"{v:.2f} ± {S1_conf[i]:.2f}"
    ax1.text(
        i,
        (v if not np.isnan(v) else 0) + (0.02 if not np.isnan(v) else 0.01),
        txt,
        ha="center",
        fontsize=9,
    )

ax2.bar(x, ST, width=bar_width, color=cmap(1), edgecolor="black", alpha=0.9)
ax2.errorbar(
    x, ST, yerr=ST_conf, fmt="none", ecolor="black", capsize=5, elinewidth=1.25
)
ax2.set_title("Total-order Sensitivity Indices (ST)", fontsize=14)
ax2.set_ylabel("Sensitivity Index", fontsize=12)
ax2.set_xlabel("Parameter", fontsize=12)
ax2.set_xticks(x)
ax2.set_xticklabels(names, rotation=45, ha="right", fontsize=10)
ax2.set_ylim(ymin, ymax)
ax2.grid(which="major", axis="y", linestyle="--", alpha=0.4)

for i, v in enumerate(ST):
    if np.isnan(v):
        txt = "nan"
    else:
        txt = f"{v:.2f} ± {ST_conf[i]:.2f}"
    ax2.text(
        i,
        (v if not np.isnan(v) else 0) + (0.02 if not np.isnan(v) else 0.01),
        txt,
        ha="center",
        fontsize=9,
    )

ax1.yaxis.set_major_formatter(FormatStrFormatter("%.2f"))
ax2.yaxis.set_major_formatter(FormatStrFormatter("%.2f"))

for ax in (ax1, ax2):
    ax.axhline(0, color="gray", linewidth=0.8)
    ax.set_ylim(bottom=ymin)

plt.suptitle("Sobol Sensitivity Indices with 95% Confidence Intervals", fontsize=16)
plt.show()

interaction_effects = ST - S1
print("Interaction Effects (ST - S1):")
for i, param_name in enumerate(names):
    val = interaction_effects[i]
    if np.isnan(val):
        print(f"{param_name}: nan")
    else:
        tag = "(interactions dominant)" if val > 0.05 else ""
        print(f"{param_name}: {val:.3f}{tag}")

Interaction Effects (ST - S1):
heat_of_combustion: 0.194(interactions dominant)
radiative_fraction: 0.042
target_thickness: 0.018
target_emissivity: 0.025

Cleanup

CFAST generates temporary output files during simulation runs. We clean these up to keep the workspace tidy and avoid confusion with future runs.

import glob
import os

files_to_remove = glob.glob("USN_Hawaii_Test_03*")

for file in files_to_remove:
    if os.path.exists(file):
        os.remove(file)
        print(f"Removed {file}")

if files_to_remove:
    print("Cleanup completed!")
else:
    print("No files to clean up.")

Removed USN_Hawaii_Test_03_parsed_zone.csv
Removed USN_Hawaii_Test_03_parsed.in
Removed USN_Hawaii_Test_03_parsed_vents.csv
Removed USN_Hawaii_Test_03_parsed.smv
Removed USN_Hawaii_Test_03_parsed.out
Removed USN_Hawaii_Test_03_parsed.status
Removed USN_Hawaii_Test_03_parsed_compartments.csv
Removed USN_Hawaii_Test_03_parsed_devices.csv
Removed USN_Hawaii_Test_03_parsed.plt
Removed USN_Hawaii_Test_03_parsed_masses.csv
Removed USN_Hawaii_Test_03_parsed.log
Removed USN_Hawaii_Test_03_parsed_walls.csv
Cleanup completed!