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	47410.363561	0.766250	0.144942	0.805988
1	14522.579276	0.766250	0.144942	0.805988
2	47410.363561	0.261726	0.144942	0.805988
3	47410.363561	0.766250	0.443230	0.805988
4	47410.363561	0.766250	0.144942	0.833259

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=47410.36356110126, rf=0.766249826643616, thickness=0.14494219878688455, emissivity=0.8059879719745369
Sample 50/640: hoc=11987.648245692253, rf=0.28936069132760167, thickness=0.19125614119693637, emissivity=0.8465721073560417
Sample 100/640: hoc=8030.70919001475, rf=0.7464817758649588, thickness=0.35063072564080355, emissivity=0.8359183150343598
Sample 150/640: hoc=45106.93948259577, rf=0.3475165076553822, thickness=0.4631066582165658, emissivity=0.8703478946816177
Sample 200/640: hoc=34593.38555680588, rf=0.7257574131712318, thickness=0.5088625754229724, emissivity=0.8613909552339465
Sample 250/640: hoc=14369.037607312202, rf=0.24015680784359575, thickness=0.15233415449038148, emissivity=0.8677431277465075
Sample 300/640: hoc=20942.112887464464, rf=0.7866323976777494, thickness=0.23676777994260192, emissivity=0.8560732133686543
Sample 350/640: hoc=25510.251733846962, rf=0.12212109677493573, thickness=0.26938829887658355, emissivity=0.8647939141839742
Sample 400/640: hoc=40124.19003220275, rf=0.7351479401811958, thickness=0.32676704730838535, emissivity=0.9069193908944726
Sample 450/640: hoc=708.9922304265201, rf=0.35874115619808433, thickness=0.2516225766390562, emissivity=0.9425162628758699
Sample 500/640: hoc=2893.7357547692955, rf=0.5636289807967843, thickness=0.4307295041158795, emissivity=0.8428680005948991
Sample 550/640: hoc=37548.86142062023, rf=0.18585570296272635, thickness=0.38307219222187994, emissivity=0.8072757772170007
Sample 600/640: hoc=27991.00147727877, rf=0.7978781790472568, thickness=0.4482516389340162, emissivity=0.9245299776084721
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.4604
radiative_fraction: 0.5368
target_thickness: 0.0573
target_emissivity: 0.0077

Total-order indices (ST):

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

heat_of_combustion: 0.4605
radiative_fraction: 0.4574
target_thickness: 0.1960
target_emissivity: 0.0281

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.000
radiative_fraction: -0.079
target_thickness: 0.139(interactions dominant)
target_emissivity: 0.020

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.status
Removed USN_Hawaii_Test_03_parsed_masses.csv
Removed USN_Hawaii_Test_03_parsed_compartments.csv
Removed USN_Hawaii_Test_03_parsed.log
Removed USN_Hawaii_Test_03_parsed_devices.csv
Removed USN_Hawaii_Test_03_parsed_zone.csv
Removed USN_Hawaii_Test_03_parsed.in
Removed USN_Hawaii_Test_03_parsed_walls.csv
Removed USN_Hawaii_Test_03_parsed.smv
Removed USN_Hawaii_Test_03_parsed_vents.csv
Removed USN_Hawaii_Test_03_parsed.plt
Removed USN_Hawaii_Test_03_parsed.out
Cleanup completed!