Note
Go to the end to download the full example code.
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#
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.
model
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:
Heat of combustion: Energy released per unit mass of fuel (affects fire intensity)
Radiative fraction: Fraction of fire energy released as radiation (affects heat transfer)
Target thickness: Material thickness of temperature measurement targets
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", # Heat of combustion (MJ/kg)
"radiative_fraction", # Fire radiative fraction (adimensional)
"target_thickness", # Target material thickness (m)
"target_emissivity", # Target material emissivity (adimensional)
],
"bounds": [
[100, 50000], # Heat of combustion: 100-50000 MJ/kg
[0.1, 1.0], # Radiative fraction: 0.1-1.0
[0.15, 0.60], # Target thickness: 0.15-0.60 m
[0.8, 0.95], # Target emissivity: 0.8-0.95
],
}
print("Sensitivity analysis problem defined:")
for i, name in enumerate(problem["names"]):
bounds = problem["bounds"][i]
print(f" {name}: [{bounds[0]}, {bounds[1]}]")
Sensitivity analysis problem defined:
heat_of_combustion: [100, 50000]
radiative_fraction: [0.1, 1.0]
target_thickness: [0.15, 0.6]
target_emissivity: [0.8, 0.95]
Step 4: Generate Parameter Samples#
We use Sobol sequences to create well-distributed parameter combinations. Sobol sampling ensures:
Good coverage of the parameter space
Efficient sampling for sensitivity analysis
Proper computation of interaction effects
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"])
print("\nFirst 5 parameter combinations:")
print(df_samples.head())
Generated 640 parameter combinations
Sample shape: (640, 4)
First 5 parameter combinations:
heat_of_combustion radiative_fraction target_thickness target_emissivity
0 31492.521572 0.571776 0.513563 0.814845
1 6615.572754 0.571776 0.513563 0.814845
2 31492.521572 0.120317 0.513563 0.814845
3 31492.521572 0.571776 0.223460 0.814845
4 31492.521572 0.571776 0.513563 0.874361
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()andupdate_material_params()Run the CFAST simulation with
run()Extract and store the output values
Progress indicators help track completion of all evaluations.
print("Running model evaluations...")
outputs = []
for i, params in enumerate(param_values):
if i % 50 == 0:
print(f"Processing sample {i}/{len(param_values)}")
print(
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], # heat of combustion in MJ/kg
radiative_fraction=params[1], # radiative fraction (0-1)
)
temp_model = temp_model.update_material_params(
material="STEELSHT",
thickness=params[2], # thickness in meters
emissivity=params[3], # emissivity (0-1)
)
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}")
Running model evaluations...
Processing sample 0/640
hoc=31492.521572485566, rf=0.5717756270430983, thickness=0.5135625824332237, emissivity=0.8148452483583242
Processing sample 50/640
hoc=9666.265978943557, rf=0.2575571526773274, thickness=0.5482137045823037, emissivity=0.8459060095250607
Processing sample 100/640
hoc=9145.329202152789, rf=0.6222266593016684, thickness=0.1785783180501312, emissivity=0.8263050523120911
Processing sample 150/640
hoc=37227.449157927185, rf=0.31942721037194133, thickness=0.2556590145919472, emissivity=0.8726047772914172
Processing sample 200/640
hoc=38890.82333482802, rf=0.6520382239483297, thickness=0.3485228758770972, emissivity=0.8603055193088949
Processing sample 250/640
hoc=19871.727732289582, rf=0.4109266906045378, thickness=0.4929217284545302, emissivity=0.8670865799766034
Processing sample 300/640
hoc=14215.965873934329, rf=0.6053294532932341, thickness=0.45672696651890865, emissivity=0.8488405632786453
Processing sample 350/640
hoc=47400.74555622414, rf=0.4638157720677555, thickness=0.4750896418467163, emissivity=0.8562633258290588
Processing sample 400/640
hoc=30484.69004649669, rf=0.6103650408796966, thickness=0.3960473406594246, emissivity=0.8879477511625737
Processing sample 450/640
hoc=2438.445478025824, rf=0.3102090531028807, thickness=0.4446570959407836, emissivity=0.9131355569697917
Processing sample 500/640
hoc=1020.0531443580985, rf=0.8546855271793902, thickness=0.2248543360736221, emissivity=0.8426053887233138
Processing sample 550/640
hoc=29075.427332986146, rf=0.536978883575648, thickness=0.16195285082794725, emissivity=0.8185939020011574
Processing sample 600/640
hoc=49595.82871813327, rf=0.5855433824472129, thickness=0.2445194221101701, emissivity=0.9369054885581135
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.
print("Sobol Sensitivity Analysis Results:")
print("\nFirst-order indices (S1):")
for i, param_name in enumerate(problem["names"]):
print(f" {param_name}: {Si['S1'][i]:.4f}")
print("\nTotal-order indices (ST):")
for i, param_name in enumerate(problem["names"]):
print(f" {param_name}: {Si['ST'][i]:.4f}")
Sobol Sensitivity Analysis Results:
First-order indices (S1):
heat_of_combustion: 0.0279
radiative_fraction: 0.7355
target_thickness: 0.0661
target_emissivity: 0.0355
Total-order indices (ST):
heat_of_combustion: 0.5501
radiative_fraction: 0.8324
target_thickness: 0.0406
target_emissivity: 0.0574
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.522(interactions dominant)
radiative_fraction: 0.097(interactions dominant)
target_thickness: -0.025
target_emissivity: 0.022
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_masses.csv
Removed USN_Hawaii_Test_03_parsed.out
Removed USN_Hawaii_Test_03_parsed.log
Removed USN_Hawaii_Test_03_parsed.status
Removed USN_Hawaii_Test_03_parsed.smv
Removed USN_Hawaii_Test_03_parsed.plt
Removed USN_Hawaii_Test_03_parsed.in
Removed USN_Hawaii_Test_03_parsed_compartments.csv
Removed USN_Hawaii_Test_03_parsed_vents.csv
Removed USN_Hawaii_Test_03_parsed_walls.csv
Removed USN_Hawaii_Test_03_parsed_devices.csv
Removed USN_Hawaii_Test_03_parsed_zone.csv
Cleanup completed!
Total running time of the script: (1 minutes 44.422 seconds)