.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_04_sensitivity_analysis.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_04_sensitivity_analysis.py: ======================= 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. .. GENERATED FROM PYTHON SOURCE LINES 12-14 Step 1: Import Required Libraries ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 14-24 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 25-31 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 :func:`~pycfast.parsers.parse_cfast_file` to load the `USN_Hawaii_Test_03.in `_ input file. .. GENERATED FROM PYTHON SOURCE LINES 31-34 .. code-block:: Python model = parse_cfast_file("data/USN_Hawaii_Test_03.in") .. GENERATED FROM PYTHON SOURCE LINES 35-36 The parsed model is displayed below. .. GENERATED FROM PYTHON SOURCE LINES 36-39 .. code-block:: Python model .. raw:: html

🔥

CFAST Model

USN_Hawaii_Test_03_parsed.in

Total Components: 5

Simulation: HawaiiTest 3

Duration: 10 min

Compartments (1)

Bay 1: 97.6×74.0×14.8 m

Fires (1)

Hawaii_03 in Bay 1: Hawaii_03_Fire

Devices (1)

Targ 1 in Bay 1: PLATE

Materials (2)

CONCRETE: Concrete Normal Weight (6 in)
STEELSHT: Steel Plain Carbon (1/16 in)

.. GENERATED FROM PYTHON SOURCE LINES 40-51 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. .. GENERATED FROM PYTHON SOURCE LINES 51-73 .. code-block:: Python 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]}]") .. rst-class:: sphx-glr-script-out .. code-block:: none 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] .. GENERATED FROM PYTHON SOURCE LINES 74-84 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. .. GENERATED FROM PYTHON SOURCE LINES 84-95 .. code-block:: Python 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()) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 96-106 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 :meth:`~pycfast.CFASTModel.update_fire_params` and :meth:`~pycfast.CFASTModel.update_material_params` - Run the CFAST simulation with :meth:`~pycfast.CFASTModel.run` - Extract and store the output values Progress indicators help track completion of all evaluations. .. GENERATED FROM PYTHON SOURCE LINES 106-137 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none 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,) .. GENERATED FROM PYTHON SOURCE LINES 138-147 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. .. GENERATED FROM PYTHON SOURCE LINES 147-150 .. code-block:: Python Si = sobol_analyze.analyze(problem, np.array(Y)) .. GENERATED FROM PYTHON SOURCE LINES 151-161 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 162-172 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. .. GENERATED FROM PYTHON SOURCE LINES 172-266 .. code-block:: Python 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}") .. image-sg:: /auto_examples/images/sphx_glr_plot_04_sensitivity_analysis_001.png :alt: Sobol Sensitivity Indices with 95% Confidence Intervals, First-order Sensitivity Indices (S1), Total-order Sensitivity Indices (ST) :srcset: /auto_examples/images/sphx_glr_plot_04_sensitivity_analysis_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 267-271 Cleanup ------- CFAST generates temporary output files during simulation runs. We clean these up to keep the workspace tidy and avoid confusion with future runs. .. GENERATED FROM PYTHON SOURCE LINES 271-286 .. code-block:: Python 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.") .. rst-class:: sphx-glr-script-out .. code-block:: none 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! .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 44.422 seconds) .. _sphx_glr_download_auto_examples_plot_04_sensitivity_analysis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_04_sensitivity_analysis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_04_sensitivity_analysis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_04_sensitivity_analysis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_