sobol_indices_salib.py

"""
This script uses SALib to estimate Sobol' indices from the Beam model.
The beam model is implemented as a Python function.
Moreover, we measure the time necessary to evaluate the model and estimate 
the Sobol' indices.

Output
------
N = 100000
Elapsed = 218.30088353157043 (s)

References
----------
http://openturns.github.io/openturns/latest/user_manual/_generated/openturns.usecases.cantilever_beam.CantileverBeam.html
"""

# %%
import openturns as ot
import numpy as np
import matplotlib.pyplot as plt
from SALib.sample import sobol
from SALib.analyze import sobol as sobol_analyze
import time

# %%
t1 = time.time()

# %%
def function_of_interest(U):
    # Problem definition for the cantilever beam model
    distribution_E = ot.Beta(0.9, 3.1, 2.8e7, 4.8e7)
    F_parameters = ot.LogNormalMuSigma(3.0e4, 9.0e3, 15.0e3)
    distribution_F = ot.ParametrizedDistribution(F_parameters)
    distribution_L = ot.Uniform(250.0, 260.0)
    distribution_II = ot.Beta(2.5, 4, 310.0, 450.0)
    distribution_X = ot.JointDistribution([distribution_E, distribution_F, distribution_L, distribution_II])
    # Transform uniform input to actual distribution
    distribution_U = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * 4)
    mapping = ot.DistributionTransformation(distribution_U, distribution_X)
    E, F, L, II = mapping(U)
    #
    Y = F * L**3 / (3 * E * II)
    return [Y]


# %%
# Generate bounds for SALib
problem = {
    "num_vars": 4,
    "names": ["E", "F", "L", "I"],
    "bounds": [[0.0, 1.0]] * 4,
}
problem

# %%
# Generate samples using Sobol' sequence
N = 100000
param_values = sobol.sample(problem, N, calc_second_order=False)

# %%
# Evaluate the model
Y = np.array([function_of_interest(U)[0] for U in param_values])
Y

# %%
# Perform Sobol analysis
Si = sobol_analyze.analyze(problem, Y, calc_second_order=False, print_to_console=False)
Si

# %%
t2 = time.time()
print(f"Elapsed = {t2 - t1} (s)")

# %%
# Plot
ind = np.arange(problem["num_vars"])
fig = plt.figure()
plt.bar(ind-0.125, Si["S1"], width=0.25, yerr=Si["S1_conf"], capsize=5)
plt.bar(ind + 0.125, Si["ST"], width=0.25, yerr=Si["ST_conf"], capsize=5)
plt.xticks(ind, problem["names"], rotation=45)
plt.title(f"Sobol indices, n = {N}")
plt.ylabel("Sensitivity Index")
plt.xlabel("Input Variables")
plt.legend(["First-order", "Total order"])
plt.grid()
plt.ylim(-0.1, 1.1)
fig.savefig("beam_sobol_salib.png")


# %%