SENSITIVITY: BOREHOLE MODEL¶

This example showcases the application of different sensitivity analysis techniques available in UQpyLab to the borehole function.

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 (02f75d466023435bbd3626b777413327) 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. This computation is carried out by the function model supplied within borehole.py. The input parameters of this function are gathered into the vector X.

Create a MODEL object from the function file:

In [4]:
ModelOpts = {
    'Type': 'Model', 
    'ModelFun':'borehole.model'
}
In [5]:
myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The probabilistic input model consists of eight independent random variables.

Specify the marginals as follows:

In [6]:
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 [7]:
myInput = uq.createInput(InputOpts)

Sensitivity Analysis¶

Sensitivity analysis on the borehole model is performed with the following methods:

  • Input/output correlation
  • Standard Regression Coefficients
  • Perturbation method
  • Cotter sensitivity indices
  • Morris elementary effects
  • Sobol' sensitivity indices
  • Borgonovo sensitivity indices

Input/output correlation analysis¶

Select the sensitivity tool and the correlation method:

In [8]:
CorrSensOpts = {
    "Type": "Sensitivity",
    "Method": "Correlation"
}

Specify the sample size used to calculate the correlation-based indices:

In [9]:
CorrSensOpts["Correlation"] = {
    "SampleSize": 1e4
}

Run the sensitivity analysis:

In [10]:
CorrAnalysis = uq.createAnalysis(CorrSensOpts)

 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 ((10000,))...

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

Print the results of the analysis:

In [11]:
uq.print(CorrAnalysis)

-------------------------------------------
Correlation-based sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.807255    -0.019360   -0.006952   0.311201    0.003291    -0.312184   -0.303220   0.157166    

-------------------------------------------
Rank-Correlation-based sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.815478    -0.010427   -0.005394   0.311967    0.005590    -0.311895   -0.292013   0.150018    
Total cost (model evaluations): 10000

Display a graphical representation of the results:

In [12]:
uq.display(CorrAnalysis);

Standard Regression Coefficients (SRC)¶

Select the sensitivity tool and the SRC method:

In [13]:
SRCSensOpts = {
    "Type": "Sensitivity",
    "Method": "SRC"
}

Specify the sample size used to calculate the regression-based indices:

In [14]:
SRCSensOpts["SRC"] = {
    "SampleSize": 1e4
}

Run the sensitivity analysis:

In [15]:
SRCAnalysis = uq.createAnalysis(SRCSensOpts)

 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 ((10000,))...

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

Print the results of the analysis:

In [16]:
uq.print(SRCAnalysis)

-------------------------------------------
Standard Regression sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.817350    0.001432    0.002735    0.306399    0.004482    -0.309804   -0.301170   0.146809    

-------------------------------------------
Standard Rank-Regression sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.823424    -0.005044   0.004579    0.306992    0.001742    -0.305730   -0.290278   0.143663    
-------------------------------------------

Display a graphical representation of the results:

In [17]:
uq.display(SRCAnalysis);

Perturbation-based indices¶

Select the sensitivity tool and the perturbation method:

In [18]:
PerturbationSensOpts = {
    "Type": "Sensitivity",
    "Method": "Perturbation"
}

Run the sensitivity analysis:

In [19]:
PerturbationAnalysis = uq.createAnalysis(PerturbationSensOpts)

 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 ((9,))...

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

Print the results of the analysis:

In [20]:
uq.print(PerturbationAnalysis)

-------------------------------------------
Perturbation-based sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.697424    0.000002    0.000000    0.095836    0.000004    0.095836    0.088715    0.022184    
-------------------------------------------

Display a graphical representation of the results:

In [21]:
uq.display(PerturbationAnalysis);

Cotter sensitivity indices¶

Select the sensitivity tool and the Cotter method:

In [22]:
CotterSensOpts = {
    "Type": "Sensitivity",
    "Method": "Cotter"
}

Specify the boundaries for the factorial design:

In [23]:
CotterSensOpts["Factors"] = {
    "Boundaries": 0.5
}

Run the sensitivity analysis:

In [24]:
CotterAnalysis = uq.createAnalysis(CotterSensOpts)

 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 ((18,))...

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

Print the results of the analysis:

In [25]:
uq.print(CotterAnalysis)

-------------------------------------------
Cotter sensitivity indices:
-------------------------------------------

rw          r           Tu          Hu          Tl          Hl          L           Kw          
11.776240   0.029201    0.000037    4.812646    0.037382    4.812646    4.909657    2.250203    
-------------------------------------------

Display a graphical representation of the results:

In [26]:
uq.display(CotterAnalysis);

Morris' elementary effects¶

Select the sensitivity tool and the Morris method:

In [27]:
MorrisSensOpts = {
    "Type": "Sensitivity",
    "Method": "Morris"
}

Specify the boundaries for the Morris method:

In [28]:
MorrisSensOpts["Factors"] = {
    "Boundaries": 0.5
}

Make sure there are no unphysical values (e.g., with the positive-only lognormal variable #2).

Specify the maximum cost (in terms of model evaluations) to calculate the Morris elementary effects:

In [29]:
MorrisSensOpts["Morris"] = {
    "Cost": 1e4
}

Run the sensitivity analysis:

In [30]:
MorrisAnalysis = uq.createAnalysis(MorrisSensOpts)

 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 ((9999,))...

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

Print the results of the analysis:

In [31]:
uq.print(MorrisAnalysis)

--------------------------------------------------
     Morris sensitivity indices:
--------------------------------------------------
       rw          r           Tu          Hu          Tl          Hl          L           Kw          
mu:    22.857933   -0.046465   0.000056    8.484848    0.055494    -8.495304   -8.168279   4.076042    
mu*:   22.857933   0.046465    0.000056    8.484848    0.055494    8.495304    8.168279    4.076042    
sigma: 1.679861    0.015583    0.000013    0.954075    0.012958    0.924936    1.016367    0.499855    
--------------------------------------------------
Total cost (model evaluations): 9999

Display a graphical representation of the results:

In [32]:
uq.display(MorrisAnalysis);

Sobol' indices¶

Select the sensitivity tool and the Sobol' method:

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

Specify the maximum order of the Sobol' indices calculation:

In [34]:
SobolOpts["Sobol"] = {
    "Order": 1
}

Specify the sample size for the indices estimation of each variable

In [35]:
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 [36]:
SobolAnalysis = 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.

Print the results of the analysis:

In [37]:
uq.print(SobolAnalysis)

--------------------------------------------------
     Total Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.689004    0.000003    0.000000    0.102545    0.000008    0.098916    0.099439    0.024472    
--------------------------------------------------
--------------------------------------------------
    First Order Sobol' indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.669092    0.004567    0.004554    0.092456    0.004577    0.090096    0.096124    0.028173    
--------------------------------------------------
Total cost (model evaluations): 100000

Create a graphical representation of the results:

In [38]:
uq.display(SobolAnalysis);

Borgonovo indices¶

Select the sensitivity tool and the Borgonovo method:

In [39]:
BorgonovoOpts = {
    "Type": "Sensitivity",
    "Method": "Borgonovo"
}

Specify the sample size:

In [40]:
BorgonovoOpts["Borgonovo"] = {
    "SampleSize": 1e4
}

A relatively large sample size is recommended for Borgonovo indices estimation, especially for complex functions.

Specify the amount of classes in Xi direction:

In [41]:
BorgonovoOpts["Borgonovo"]["NClasses"] = 20

By default, UQLab will then create classes that contain the same amount of sample points.

Run the sensitivity analysis:

In [42]:
BorgonovoAnalysis = uq.createAnalysis(BorgonovoOpts)

 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 ((10000,))...

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

Print the results of the analysis:

In [43]:
uq.print(BorgonovoAnalysis)

--------------------------------------------------
     Borgonovo indices for output component 1
--------------------------------------------------
rw          r           Tu          Hu          Tl          Hl          L           Kw          
0.410648    0.050435    0.052653    0.121902    0.056921    0.124577    0.120598    0.077573    
--------------------------------------------------
Total cost (model evaluations): 10000

Create a graphical representation of the results:

In [44]:
uq.display(BorgonovoAnalysis);

In order to assess the accuracy of the results, it is possible to inspect the 2D histogram estimation of the joint distribution used in the calculation of an index:

In [45]:
uq.display(BorgonovoAnalysis, outidx=1, Joint_PDF=True, inidx=1);

Terminate the remote UQCloud session¶

In [46]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 02f75d466023435bbd3626b777413327 terminated.
Out[46]:
True