SENSITIVITY: COMPARISON OF MC- VS. PCE- AND LRA-BASED SOBOL' INDICES¶

In this example, Sobol' sensitivity indices for the borehole function are calculated with three different methods: Monte Carlo (MC) simulation, polynomial chaos expansion (PCE), and canonical low-rank approximation (LRA).

Package imports¶

In [1]:
from uqpylab import sessions
import numpy as np

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 (263581be9baf480fa4c43fe124dad25c) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

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

Computational model¶

The computational model is an $8$-dimensional analytical formula that is used to model the water flow through a borehole. The borehole function uq_borehole is supplied with UQLab.

Create a MODEL object from the function file:

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

myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The probabilistic input model consists of eight independent random variables.

Specify the marginals as follows:

In [5]:
InputOpts = {
    'Marginals': [
        {
            'Name': 'rw', # Radius of the borehole
            'Type': 'Gaussian',
            'Parameters': [0.10, 0.0161812] # % (m)
        },
        {
            'Name': 'r', # Radius of influence
            'Type': 'Lognormal',
            'Parameters': [7.71, 1.0056] # % (m)
        },
        {
            'Name': 'Tu', # Transmissivity, upper aquifer
            'Type': 'Uniform',
            'Parameters': [63070, 115600] # % (m^2/yr)
        },
        {
            'Name': 'Hu', # Potentiometric head, upper aquifer
            'Type': 'Uniform',
            'Parameters': [990, 1110] # % (m)
        },
        {
            'Name': 'Tl', # Transmissivity, lower aquifer
            'Type': 'Uniform',
            'Parameters': [63.1, 116] # % (m^2/yr)
        },
        {
            'Name': 'Hl', # Potentiometric head , lower aquifer
            'Type': 'Uniform',
            'Parameters': [700, 820] # % (m)
        },
        {
            'Name': 'L', # Length of the borehole
            'Type': 'Uniform',
            'Parameters': [1120, 1680] # % (m)
        },
        {
            'Name': 'Kw', # Borehole hydraulic conductivity
            'Type': 'Uniform',
            'Parameters': [9855, 12045] # % (m/yr)
        },
    ]
}

Create an INPUT object based on the specified marginals:

In [6]:
myInput = uq.createInput(InputOpts)

Sensitivity analysis¶

Sobol' indices are calculated first with a direct MC simulation of the model and subsequently through post-processing of the coefficients of its PCE and LRA approximation.

MC-based Sobol' indices¶

Select the sensitivity analysis module in UQLab and the Sobol' analysis method:

In [7]:
SobolOpts = {
    'Type': 'Sensitivity',
    'Method': 'Sobol'
}

Specify the maximum order of the Sobol' indices to be calculated:

In [8]:
SobolOpts['Sobol'] = {
    'Order': 1
}

Specify the sample size for the MC simulation:

In [9]:
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 $(8+2) \times 10^4 = 10^5$ evaluations of the full computational model.

Run the sensitivity analysis:

In [10]:
mySobolAnalysisMC = uq.createAnalysis(SobolOpts)

Processing .
 done!

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

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

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

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

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

Retrieve the analysis results for comparison:

In [11]:
mySobolResultsMC = mySobolAnalysisMC['Results']

PCE-based Sobol' indices¶

Select the metamodeling tool in UQLab and the polynomial chaos expansion (PCE) type:

In [12]:
PCEOpts = {
    'Type': 'Metamodel',
    'MetaType': 'PCE',
    'Method': 'LARS'
}

Assign the full computational model:

In [13]:
PCEOpts['FullModel'] = myModel['Name']

The full model is used to generate an experimental design for the PCE.

Specify the maximum polynomial degree (default: sparse PCE):

In [14]:
PCEOpts['Degree'] = 5

Specify the size of the experimental design (i.e., the total computational cost of constructing the metamodel):

In [15]:
PCEOpts['ExpDesign'] = {
    'NSamples': 200
}

Calculate the PCE:

In [16]:
myPCE = uq.createModel(PCEOpts)

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

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

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

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

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

The same options structure SobolOpts can be re-used to create a new analysis on the PCE model:

In [17]:
mySobolAnalysisPCE = uq.createAnalysis(SobolOpts)

Retrieve the results for comparison:

In [18]:
mySobolResultsPCE = mySobolAnalysisPCE['Results']

Note: The PCE metamodel does not need to be explicitly specified as the last created model will be used in the analysis. The sensitivity analysis module in UQLab recognizes PCE models and calculates the Sobol' indices accordingly by post-processing the coefficients of a PCE model without sampling.

LRA-based Sobol' indices¶

Select the metamodeling tool in UQLab and the canonical low-rank approximation (LRA) type:

In [19]:
LRAOpts = {
    'Type': 'Metamodel',
    'MetaType': 'LRA'
}

Assign the full computational model:

In [20]:
LRAOpts['FullModel'] = myModel['Name']

Specify the rank range:

In [21]:
LRAOpts['Rank'] = np.arange(1,20).tolist()

Specify the degree range:

In [22]:
LRAOpts['Degree'] = np.arange(1,20).tolist()

Specify the size of the experimental design (i.e., the total computation cost of the metamodel):

In [23]:
LRAOpts['ExpDesign'] = {
    'NSamples': 200
}

Calculate the LRA:

In [24]:
myLRA = uq.createModel(LRAOpts)

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

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

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

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

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

As before, the same options structure SobolOpts can be re-used to create a new analysis on the LRA model:

In [25]:
mySobolAnalysisLRA = uq.createAnalysis(SobolOpts)

Retrieve the results for comparison:

In [26]:
mySobolResultsLRA = mySobolAnalysisLRA['Results']

Note: The LRA metamodel does not need to be explicitly specified as the last created model will be used in the analysis. The sensitivity analysis module in UQLab recognizes LRA models and calculates the Sobol' indices accordingly by post-processing the coefficients of a LRA model without sampling.

Comparison of the results¶

Print the results of the Sobol' indices calculation based on the MC simulation:

In [27]:
uq.print(mySobolAnalysisMC)

--------------------------------------------------
     Total Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.681016    0.000003    0.000000    0.106299    0.000008    0.105510    0.100767    0.025361    
--------------------------------------------------
--------------------------------------------------
    First Order Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.668886    0.021934    0.021918    0.119483    0.021888    0.114352    0.110502    0.039414    
--------------------------------------------------
Total cost (model evaluations): 100000

Print the results of the Sobol' indices calculation based on PCE:

In [28]:
uq.print(mySobolAnalysisPCE)

--------------------------------------------------
     Total Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.695782    0.000103    0.000202    0.107328    0.000056    0.103624    0.100941    0.023204    
--------------------------------------------------
--------------------------------------------------
    First Order Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.668395    0.000000    0.000000    0.096756    0.000000    0.093267    0.090301    0.021007    
--------------------------------------------------
Total cost (model evaluations): 0

Print the results of the Sobol' indices calculation based on LRA:

In [29]:
uq.print(mySobolAnalysisLRA)

--------------------------------------------------
     Total Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.696226    0.000095    0.000028    0.106453    0.000016    0.102549    0.104212    0.026409    
--------------------------------------------------
--------------------------------------------------
    First Order Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.665599    0.000082    0.000025    0.093762    0.000014    0.090277    0.091761    0.023015    
--------------------------------------------------
Total cost (model evaluations): 0

Create a bar plot to compare the total Sobol' indices:

In [30]:
data = {
    'MC-based ({:.2E} simulations)'.format(mySobolResultsMC['Cost']): mySobolResultsMC['Total'],
    'PCE-based ({} simulations)'.format(myPCE['ExpDesign']['NSamples']): mySobolResultsPCE['Total'],
    'LRA-based ({} simulations)'.format(myLRA['ExpDesign']['NSamples']): mySobolResultsLRA['Total']
}

uq.display_bar(
    data, 
    VarNames=mySobolAnalysisMC['Results']['VariableNames'], 
    xaxis_title='Variable name', 
    yaxis_title='Total Sobol\' indices',
    showlegend=True, 
    width=900
)

Create a bar plot to compare the first-order Sobol' indices:

In [31]:
data = {
    'MC-based ({:.2E} simulations)'.format(mySobolResultsMC['Cost']): mySobolResultsMC['FirstOrder'],
    'PCE-based ({} simulations)'.format(myPCE['ExpDesign']['NSamples']): mySobolResultsPCE['FirstOrder'],
    'LRA-based ({} simulations)'.format(myLRA['ExpDesign']['NSamples']): mySobolResultsLRA['FirstOrder']
}

uq.display_bar(
    data, 
    VarNames=mySobolAnalysisMC['Results']['VariableNames'], 
    xaxis_title='Variable name', 
    yaxis_title='First-order Sobol\' indices',
    showlegend=True, 
    width=900
)

Terminate the remote UQCloud session¶

In [32]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 263581be9baf480fa4c43fe124dad25c terminated.
Out[32]:
True