SENSITIVITY: MULTIPLE OUTPUTS¶

This example showcases an application of sensitivity analysis to a simply supported beam model. The model computes the deflections at several points along the length of the beam subjected to a uniform random load.

Package imports¶

In [1]:
from uqpylab import sessions
import pandas as pd

Start a remote UQCloud session¶

In [2]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (d87125ba3d3f4ad69039616ee8348730) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [3]:
uq.rng(50,'twister');

Computational model¶

The simply supported beam model is shown in the following figure:

No description has been provided for this image

The (negative) deflection of the beam at any longitudinal coordinate $s$ is given by:

$$V(s) = -\frac{p \,s (L^3 - 2\, s^2 L + s^3) }{2E b h^3}$$

This computation is carried out by the function simply_supported_beam_9points. The function evaluates the inputs gathered in the $N \times M$ matrix X, where $N$ and $M$ are the numbers of realizations and inputs, respectively. The inputs are given in the following order:

  • $b$: beam width $(m)$
  • $h$: beam height $(m)$
  • $L$: beam length $(m)$
  • $E$: Young's modulus $(Pa)$
  • $p$: uniform load $(N/m)$

The function returns the beam deflection $V(s_i)$ at nine equally-spaced points along the length $s_i = i \cdot L/10, \; i=1,\ldots,9.$

Create a MODEL object from the simply_supported_beam_9points function:

In [4]:
ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'simply_supported_beam_9points.model',
    'isVectorized': 'true'
}

myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The simply supported beam model has five independent input parameters modeled by lognormal random variables. The parameters of the distributions are given in the following table:

Variable Description Distribution Mean Std. deviation
b Beam width Lognormal 0.15 m 7.5 mm
h Beam height Lognormal 0.3 m 15 mm
L Length Lognormal 5 m 50 mm
E Young modulus Lognormal 30000 MPa 4500 MPa
p Uniform load Lognormal 10 kN/m 2 kN/m

Define an INPUT object with the following marginals:

In [5]:
InputOpts = {
    'Marginals': [
        {
        'Name': 'b', # beam width
        'Type': 'Lognormal',
        'Moments': [0.15, 0.0075] # (m)
        },
        {
        'Name': 'h', # beam height
        'Type': 'Lognormal',
        'Moments': [0.3, 0.015] # (m)
        },
        {
        'Name': 'L', # beam length
        'Type': 'Lognormal',
        'Moments': [5, 0.05] # (m)
        },
        {
        'Name': 'E', # Young's modulus
        'Type': 'Lognormal',
        'Moments': [3e10, 4.5e9] # (Pa)
        },
        {
        'Name': 'p', # uniform load
        'Type': 'Lognormal',
        'Moments': [1e4, 2e3] # (N/m)
        }]
}

myInput = uq.createInput(InputOpts)

SENSITIVITY ANALYSIS¶

Sensitivity analysis is performed by calculating the Sobol' indices for each of the output components separately.

Select the Sobol' sensitivity tool:

In [6]:
SobolOpts = {
    "Type": "Sensitivity",
    "Method": "Sobol"
}

Specify the sample size of each variable:

In [7]:
SobolOpts["Sobol"] = {
    "SampleSize": 1e4
}

Note that the total cost of computation is $(M+2) \times N$, where $M$ is the input dimension and $N$ is the sample size. Therefore, the total cost for the current setup is $(5+2) \times 10^4 = 7 \times 10^4$ evaluations of the full computational model.

Run the sensitivity analysis:

In [8]:
mySobolAnalysis = uq.createAnalysis(SobolOpts)

Processing .
 done!

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: simply_supported_beam_9points.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((70000, 9))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

RESULTS VISUALIZATION¶

Print out a report of the results:

In [9]:
uq.print(mySobolAnalysis)

--------------------------------------------------
     Total Sobol' indices for output component 1
--------------------------------------------------
b           h           L           E           p           
0.030236    0.265587    0.019782    0.263724    0.456656    
--------------------------------------------------
--------------------------------------------------
    First Order Sobol' indices for output component 1
--------------------------------------------------
b           h           L           E           p           
0.018281    0.251242    0.013435    0.240050    0.421727    
--------------------------------------------------
Total cost (model evaluations): 70000

Retrieve the analysis results (the total and first-order indices):

In [10]:
SobolResults = mySobolAnalysis['Results']
TotalSobolIndices = SobolResults['Total']
FirstOrderIndices = SobolResults['FirstOrder']

Plot the total Sobol' indices for all the output components:

In [11]:
df = pd.DataFrame(TotalSobolIndices)
df.columns = [f'Y<sub>{i+1}' for i in range(len(TotalSobolIndices[0]))]
df.index = SobolResults['VariableNames']

uq.display_bar(
    df, 
    xaxis_title='Input variable', 
    yaxis_title='S<sub>i</sub><sup>T', 
    title='Total Sobol\' indices', 
    showlegend=True, 
    width=900
)

Plot the first-order Sobol' indices for all the output components:

In [12]:
df = pd.DataFrame(FirstOrderIndices)
df.columns = [f'Y<sub>{i+1}' for i in range(len(FirstOrderIndices[0]))]
df.index = SobolResults['VariableNames']

uq.display_bar(
    df, 
    xaxis_title='Input variable',
    yaxis_title='S<sub>i', 
    title='First-order Sobol\' indices', 
    showlegend=True, 
    width=900
)

Terminate the remote UQCloud session¶

In [13]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session d87125ba3d3f4ad69039616ee8348730 terminated.
Out[13]:
True