sobol_indices_openturns.py

"""
This script uses OpenTURNS to estimate Sobol' indices from the Beam model implemented as an OpenTURNS's function.
Moreover, we measure the time necessary to evaluate the model and estimate the Sobol' indices.

Output
------
Elapsed = 1.0558149814605713 (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
import openturns.viewer as otv
import time

# %%
t1 = time.time()


# %%
def function_of_interest(X):
    E, F, L, I = X
    Y = F * L**3 / (3 * E * I)
    return [Y]


# %%
# Problem definition for the cantilever beam model
E = ot.Beta(0.9, 3.1, 2.8e7, 4.8e7)
F_parameters = ot.LogNormalMuSigma(3.0e4, 9.0e3, 15.0e3)
F = ot.ParametrizedDistribution(F_parameters)
L = ot.Uniform(250.0, 260.0)
II = ot.Beta(2.5, 4, 310.0, 450.0)
distribution = ot.JointDistribution([E, F, L, II])

# %%
dimension = distribution.getDimension()
model = ot.PythonFunction(dimension, 1, function_of_interest)

# %%
sizeSobol = 100000
computeSecondOrder = False
sie = ot.SobolIndicesExperiment(distribution, sizeSobol, computeSecondOrder)
inputDesignSobol = sie.generate()
print("Sample size = ", inputDesignSobol.getSize())
inputNames = distribution.getDescription()
inputDesignSobol.setDescription(inputNames)
inputDesignSobol.getSize()
outputDesignSobol = model(inputDesignSobol)

# %%
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(
    inputDesignSobol, outputDesignSobol, sizeSobol
)
S1 = sensitivityAnalysis.getFirstOrderIndices()
ST = sensitivityAnalysis.getTotalOrderIndices()
S1 = np.array(S1)
ST = np.array(ST)
print("First order indices:", S1)
print("Total order indices:", ST)

# %%
if computeSecondOrder:
    S2 = sensitivityAnalysis.getSecondOrderIndices()
    S2 = np.array(S2)
    print("S2 =")
    print(S2)


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

# %%
graph = sensitivityAnalysis.draw()
graph.setTitle(f"Sobol indices, n = {sizeSobol}")
view = otv.View(graph)
view.save("beam_sobol_openturns.png")