Sobol’ sensitivity indices from chaos

In this example we are going to compute global sensitivity indices from a functional chaos decomposition.

We study the Borehole function that models water flow through a borehole:

\frac{2 \pi T_u (H_u - H_l)}{\ln{r/r_w}(1+\frac{2 L T_u}{\ln(r/r_w) r^2_w K_w}\frac{T_u}{T_l})}

With parameters:

  • r_w: radius of borehole (m)

  • r: radius of influence (m)

  • T_u: transmissivity of upper aquifer (m^2/yr)

  • H_u: potentiometric head of upper aquifer (m)

  • T_l: transmissivity of lower aquifer (m^2/yr)

  • H_l: potentiometric head of lower aquifer (m)

  • L: length of borehole (m)

  • K_w: hydraulic conductivity of borehole (m/yr)

import openturns as ot
from operator import itemgetter
import openturns.viewer as viewer
from matplotlib import pylab as plt

ot.Log.Show(ot.Log.NONE)

borehole model

dimension = 8
input_names = ["rw", "r", "Tu", "Hu", "Tl", "Hl", "L", "Kw"]
model = ot.SymbolicFunction(
    input_names, ["(2*pi_*Tu*(Hu-Hl))/(ln(r/rw)*(1+(2*L*Tu)/(ln(r/rw)*rw^2*Kw)+Tu/Tl))"]
)
coll = [
    ot.Normal(0.1, 0.0161812),
    ot.LogNormal(7.71, 1.0056),
    ot.Uniform(63070.0, 115600.0),
    ot.Uniform(990.0, 1110.0),
    ot.Uniform(63.1, 116.0),
    ot.Uniform(700.0, 820.0),
    ot.Uniform(1120.0, 1680.0),
    ot.Uniform(9855.0, 12045.0),
]
distribution = ot.JointDistribution(coll)
distribution.setDescription(input_names)

Freeze r, Tu, Tl from model to go faster

selection = [1, 2, 4]
complement = ot.Indices(selection).complement(dimension)
distribution = distribution.getMarginal(complement)
model = ot.ParametricFunction(
    model, selection, distribution.getMarginal(selection).getMean()
)
input_names_copy = list(input_names)
input_names = itemgetter(*complement)(input_names)
dimension = len(complement)

design of experiment

size = 1000
X = distribution.getSample(size)
Y = model(X)

create a functional chaos model

algo = ot.FunctionalChaosAlgorithm(X, Y)
algo.run()
result = algo.getResult()
print(result.getResiduals())
print(result.getRelativeErrors())
[0.00175769]
[3.70801e-06]

Quick summary of sensitivity analysis

sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
sensitivityAnalysis
FunctionalChaosSobolIndices
  • input dimension: 5
  • output dimension: 1
  • basis size: 181
  • mean: [75.5283]
  • std-dev: [30.9634]
Input Variable Sobol' index Total index
0 rw 0.698383 0.726572
1 Hu 0.084703 0.095073
2 Hl 0.084747 0.095132
3 L 0.081092 0.092178
4 Kw 0.019748 0.022672
Index Multi-index Part of variance
1 [1,0,0,0,0] 0.623730
3 [0,0,1,0,0] 0.084747
2 [0,1,0,0,0] 0.084703
4 [0,0,0,1,0] 0.080217
5 [0,0,0,0,1] 0.019748
56 [6,0,0,0,0] 0.019191
31 [5,0,0,0,0] 0.018039
61 [7,0,0,0,0] 0.011068
6 [2,0,0,0,0] 0.010331


draw Sobol’ indices

first_order = [sensitivityAnalysis.getSobolIndex(i) for i in range(dimension)]
total_order = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(dimension)]
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, first_order, total_order)
view = viewer.View(graph)
Sobol' indices

We saw that total order indices are close to first order, so the higher order indices must be all quite close to 0

for i in range(dimension):
    for j in range(i):
        sij = sensitivityAnalysis.getSobolIndex([i, j])
        print(f"{input_names[i]} & {input_names[j]}: {sij:.4f}")

plt.show()
Hu & rw: 0.0088
Hl & rw: 0.0088
Hl & Hu: 0.0000
L & rw: 0.0083
L & Hu: 0.0012
L & Hl: 0.0011
Kw & rw: 0.0020
Kw & Hu: 0.0003
Kw & Hl: 0.0003
Kw & L: 0.0003