Global Sensitivity Analysis#

This section provides an implementation for global sensitivity analysis using the Sobol’ method. The Sobol’ method is implemented using SALib, a Python library built specifically for performing sensitivity analysis. The method is demonstrated on the Welded Beam Function and Water Flow Function that are available in the benchmarking problems section of smt. The first-order and total Sobol’ indices are found using the true function and a Kriging model of the true function.

NOTE: SALib will have to be installed using pip install SALib in your Conda environment before running the code in this section.

1. Welded Beam Function#

The three-dimensional welded beam function can be represented as:

\[\begin{split} \begin{gathered} f(\textbf{x}) = \sqrt{\frac{\tau'+\tau''^2+l\tau'\tau''}{\sqrt{0.25(l^2+(h+t)^2)}}} \\ \textrm{where} \quad \tau' = \frac{6000}{\sqrt{2}hl}, \tau'' = \frac{6000(14+0.5l)\sqrt{0.25(l^2+(h+t)^2)}}{2[0.707hl(\frac{l^2}{12}+0.25(h+t)^2)]}, \quad \textrm{for} \quad h \in [0.125,1.0], \quad l,t \in [5,10] \\ \end{gathered} \end{split}\]

The block of code below imports the required packages and defines the Welded Beam Function given in smt.

# Imports
import numpy as np
import matplotlib.pyplot as plt
from smt.sampling_methods import LHS
from smt.surrogate_models import KRG
from smt.problems import WaterFlow, WeldedBeam
from SALib.sample import sobol, latin
from SALib.analyze import sobol as SOBOL

# Defining welded beam function problem
ndim = 3
weldedbeam = WeldedBeam(ndim=ndim)

The next block of code defines the sensitivity analysis problem for SALib. The problem is defined as a dictionary with the number of variables, names and bounds of the variables.

# Defining sensitivity analysis problem
xlimits = np.array([[0.125, 1.0], [5.0, 10.0], [5.0, 10.0]])
problem = {
    'num_vars': ndim,
    'names': ['$h$', '$l$', '$t$'],
    'bounds': xlimits
}

In the next block of code, Sobol’ sequences are used to create a set of samples for performing Monte Carlo simulations to approximate the relevant integrals for calculating the Sobol’ indices. Essentially, the block of code uses functionality in SALib to generate the samples using Sobol’ sequences and calculate the corresponding function values using the true function. Refer to the lecture notes for further details.

param_values = sobol.sample(problem, 1024)
Y = weldedbeam(param_values)
Y = Y.reshape(-1)

Now, the Sobol’ indices can be calculated using the samples and function values generated in the previous code block. Both the first order (\(S_1\)) and total indices (\(S_T\)) are calculated.

Si = SOBOL.analyze(problem, Y)

print("First-order indices:", Si['S1'])
print("\nTotal indices:", Si['ST'])

First-order indices: [0.00992783 0.95928776 0.025763  ]

Total indices: [0.01476005 0.96516655 0.02891984]

Next, the code blocks will plot graphs to visualize the Sobol’ indices for the different variables of the welded beam function.

# Plotting first order indices
vars = problem['names']
s1_values = Si['S1']
st_values = Si['ST']
r = np.arange(len(s1_values))
width = 0.3
fig, ax = plt.subplots(figsize=(7,5))
ax.bar(r, s1_values, color ='blue', width = width, label = "$S_1$")
ax.bar(r+width, st_values, color ='green', width = width, label = "$S_T$")
ax.set_xlabel("Design Variables", fontsize = 14)
ax.set_ylabel("Sobol' Indices", fontsize = 14)
ax.set_xticks(r + width/2, vars, fontsize = 14) 
ax.legend(fontsize = 14)
ax.grid()
ax.set_ylim([0.0, 1.0])

(0.0, 1.0)
../_images/2713dd9cdba9aa7bf0e43edbdf09a991e467301fd31b6dc8c190d2d70f93bfe9.png

From the values of the Sobol’ indices, it is clear that the welded beam function is highly sensitive to the \(l\) variable and far less sensitive to the other variables. Additionally, the values of the total indices are similar to the first-order indices which means that there are few interaction effects between the design variables of this function.

In the next few code blocks, the use of surrogates to perform sensitivity analysis is demonstrated. A Kriging model is created of the Welded Beam function using various number of samples and the Monte Carlo simulations needed for calculating the Sobol’ indices are performed using the Kriging model instead of the true function. The idea is to reduce the computational cost of calculating the Sobol’ indices. The Monte Carlo simulations require 1000s of evaluations of the true function which can be a high-fidelity physics simulation that is expensive to evaluate. Using a surrogate in place of the high-fidelity analysis will greatly reduce the computational cost.

Plots are created to visualize the values of the first-order and total Sobol’ indices and a comparison is drawn between the indices calculated using the true function and the Kriging model.

# Defining sample sizes
samples = np.arange(5,35,5)

for size in samples:

    # Generate training samples using LHS
    sampling = LHS(xlimits=xlimits, criterion="ese")
    xtrain = sampling(size)
    ytrain = weldedbeam(xtrain)

    # Create kriging model
    corr = 'squar_exp'
    sm = KRG(theta0=[1e-2], corr=corr, theta_bounds=[1e-6, 1e2], print_global=False)
    sm.set_training_values(xtrain, ytrain)
    sm.train()

    # Generate values for Sobol sequence samples using the surrogate
    surrogate_Y = sm.predict_values(param_values)
    surrogate_Y = surrogate_Y.reshape(-1)
    
    Si_surrogate = SOBOL.analyze(problem, surrogate_Y)
    print("\nNumber of samples:", size)
    print("First-order indices:", Si_surrogate['S1'])
    print("\nTotal indices:", Si_surrogate['ST'])
    
    # Plotting first order indices 
    vars = problem['names']
    s1_surrogate = Si_surrogate['S1']
    r = np.arange(len(s1_surrogate))
    width = 0.4
    fig, ax = plt.subplots(1, 2, figsize=(12,5))
    ax[0].bar(r, s1_values, color ='blue', width = width, label = "True function")
    ax[0].bar(r+width, s1_surrogate, color ='red', width = width, label = "Kriging model")
    ax[0].set_xlabel("Design Variables", fontsize = 12)
    ax[0].set_ylabel("First-order Sobol' Indices", fontsize = 12)
    ax[0].set_xticks(r + width/2, vars, fontsize = 12) 
    ax[0].legend()
    ax[0].grid()
    ax[0].set_ylim([0.0, 1.0])

    # Plotting total indices
    st_surrogate = Si_surrogate['ST']
    ax[1].bar(r, st_values, color ='green', width = width, label = "True function")
    ax[1].bar(r+width, st_surrogate, color ='red', width = width, label = "Kriging model")
    ax[1].set_xlabel("Design Variables", fontsize = 12)
    ax[1].set_ylabel("Total Sobol' Indices", fontsize = 12)
    ax[1].set_xticks(r + width/2, vars, fontsize = 12) 
    ax[1].grid()
    ax[1].legend(fontsize = 10)
    ax[1].set_ylim([0.0, 1.0])
    fig.suptitle("Number of samples: {}".format(size))

Number of samples: 5
First-order indices: [1.45187146e-07 9.35368622e-01 6.28880660e-02]

Total indices: [1.14404722e-07 9.39716948e-01 6.44601575e-02]

Number of samples: 10
First-order indices: [0.00828477 0.97146836 0.01237414]

Total indices: [0.01633838 0.98070572 0.01331145]

Number of samples: 15
First-order indices: [0.00883056 0.96285217 0.02542413]

Total indices: [0.01216097 0.96688958 0.02705144]

Number of samples: 20
First-order indices: [0.00984579 0.96188797 0.02412994]

Total indices: [0.01342875 0.96764912 0.02681691]

Number of samples: 25
First-order indices: [0.00950821 0.95981969 0.02641781]

Total indices: [0.01343733 0.96506613 0.0289949 ]

Number of samples: 30
First-order indices: [0.00959373 0.96031721 0.02538151]

Total indices: [0.01401013 0.96591804 0.02815098]
../_images/6879a7ee5fd41b42f588107e8cabdeaacd717ea1a7a8a609ad8581429f5e9894.png ../_images/fc93489cadfddf107b73e0f54b4a4cb8eab398360e1a86bc144ad67e8d9528b8.png ../_images/6167f857517feb1fdeb0a0690ff1603766615c1a35e508c5412efaa6900188ae.png ../_images/bb9c2c53443776ace6cdbf29fc0a1084896656255a730c7e347f5d5c0d23c92b.png ../_images/3a886fc024258d17c75208346c996938d63eed5cf960b526a4ec968ee4735d12.png ../_images/1716f6574f9dc91c95bae82a589d77b754bc39ebe42e596ecd31f7b3491cca54.png

With a few number of samples, the Kriging model cannot model the Welded Beam function well enough to accurately calculate the Sobol’ indices. As the number of samples increases, the Kriging model becomes more accurate and the Sobol’ indices calculated using the Kriging model match closely with those calculated using the true function. However, calculating the Sobol’ indices using the Kriging model only required 20 - 30 evaluations of the true function which makes it effective in reducing the computational cost of calculating global sensitivities.

2. Water Flow Function#

The eight-dimensional water flow function can be represented as:

\[\begin{split} f(\textbf{x}) = \frac{2\pi T_u (H_u - H_l)}{\textrm{ln}(\frac{r}{r_w})[1 + \frac{2LT_u}{\textrm{ln}(\frac{r}{r_w})r_w^2K_w}+\frac{T_u}{T_l}]} \\ \end{split}\]
\[\begin{split} \textrm{where} \quad r_w \in [0.05,0.15], \quad r \in [100, 50000], \quad T_u \in [63070, 115600], \quad H_u \in [990, 1110], \quad T_l \in [63.10,116], \quad H_l \in [700, 820], \quad L \in [1120, 1680], \quad K_w \in [9855, 12045] \\ \end{split}\]

The blocks of code below defines the Water Flow Function given in smt and calculates the Sobol’ indices for the variables in a similar manner as the Welded Beam Function.

# Defining the water flow function
ndim = 8
waterflow = WaterFlow(ndim=ndim)

# Defining sensitivity analysis problem
xlimits = np.array([[0.05, 0.15], [100.0, 50000.0], [63070, 115600], [990, 1110], [63.10, 116], [700, 820], [1120, 1680], [9855, 12045]])
problem = {
    'num_vars': ndim,
    'names': ['$r_w$', '$r$', '$T_u$', '$H_u$', '$T_l$', '$H_l$', '$L$', '$K_w$'],
    'bounds': xlimits
}

Now, the Sobol’ indices are calculated using the true function and plots are created to visualize the values of the Sobol’ indices.

param_values = sobol.sample(problem, 1024)
Y = waterflow(param_values)
Y = Y.reshape(-1)

Si = SOBOL.analyze(problem, Y)

print("First-order indices:", Si['S1'])
print("\nTotal indices:", Si['ST'])

First-order indices: [ 8.18402669e-01 -6.08736000e-05  2.78690418e-08  4.15112337e-02
  2.77024069e-05  4.91288814e-02  4.49470052e-02  8.38858193e-03]

Total indices: [8.70999389e-01 2.56800170e-06 1.07823405e-11 5.41311167e-02
 1.05217942e-05 5.44311283e-02 5.17587036e-02 1.27692866e-02]

# Plotting first order indices
vars = problem['names']
s1_values = Si['S1']
st_values = Si['ST']
r = np.arange(len(s1_values))
width = 0.3
fig, ax = plt.subplots(figsize=(7,5))
ax.bar(r, s1_values, color ='blue', width = width, label = "$S_1$")
ax.bar(r+width, st_values, color ='green', width = width, label = "$S_T$")
ax.set_xlabel("Design Variables", fontsize = 14)
ax.set_ylabel("Sobol' Indices", fontsize = 14)
ax.set_xticks(r + width/2, vars, fontsize = 14) 
ax.legend(fontsize = 14)
ax.grid()
ax.set_ylim([1e-6, 1.0])

(1e-06, 1.0)
../_images/d788134651e384cdbc6c434a3d3946db076e4e929ba7e5173d5bed92418528c8.png

From the values of the Sobol’ indices, it is clear that the water flow function is highly sensitive to the \(r_w\) variable and slightly sensitive to the \(H_u\), \(H_l\), \(L\), and \(K_w\) variables. It can also be seen that the total indices are larger than the first-order indices indicating the presence of interaction effects between the variables of this function.

The next few code block utilize a Kriging model to conduct the sensitivity analysis in a similar manner to the analysis conducted for the welded beam function.

# Defining sample sizes
samples = np.arange(5,45,5)

for size in samples:

    # Generate training samples using LHS
    sampling = LHS(xlimits=xlimits, criterion="ese")
    xtrain = sampling(size)
    ytrain = waterflow(xtrain)

    # Create kriging model
    corr = 'squar_exp'
    sm = KRG(theta0=[1e-2], corr=corr, theta_bounds=[1e-6, 1e2], print_global=False)
    sm.set_training_values(xtrain, ytrain)
    sm.train()

    # Generate values for Sobol sequence samples using the surrogate
    surrogate_Y = sm.predict_values(param_values)
    surrogate_Y = surrogate_Y.reshape(-1)
    
    Si_surrogate = SOBOL.analyze(problem, surrogate_Y)
    print("\nNumber of samples:", size)
    print("First-order indices:", Si_surrogate['S1'])
    print("\nTotal indices:", Si_surrogate['ST'])
    
    # Plotting first order indices 
    vars = problem['names']
    s1_surrogate = Si_surrogate['S1']
    r = np.arange(len(s1_surrogate))
    width = 0.4
    fig, ax = plt.subplots(1, 2, figsize=(12,5))
    ax[0].bar(r, s1_values, color ='blue', width = width, label = "True function")
    ax[0].bar(r+width, s1_surrogate, color ='red', width = width, label = "Kriging model")
    ax[0].set_xlabel("Design Variables", fontsize = 12)
    ax[0].set_ylabel("First-order Sobol' Indices", fontsize = 12)
    ax[0].set_xticks(r + width/2, vars, fontsize = 12) 
    ax[0].legend(fontsize = 10)
    ax[0].grid()
    ax[0].set_ylim([0.0, 1.0])

    # Plotting total indices
    st_surrogate = Si_surrogate['ST']
    ax[1].bar(r, st_values, color ='green', width = width, label = "True function")
    ax[1].bar(r+width, st_surrogate, color ='red', width = width, label = "Kriging model")
    ax[1].set_xlabel("Design Variables", fontsize = 12)
    ax[1].set_ylabel("Total Sobol' Indices", fontsize = 12)
    ax[1].set_xticks(r + width/2, vars, fontsize = 12) 
    ax[1].grid()
    ax[1].legend(fontsize = 10)
    ax[1].set_ylim([0.0, 1.0])
    fig.suptitle("Number of samples: {}".format(size))

Number of samples: 5
First-order indices: [ 8.44423805e-01  5.37484118e-06  2.05471941e-05  1.24192699e-03
  3.38859581e-03  1.35813097e-01  1.48746358e-03 -4.89705185e-06]

Total indices: [8.62082397e-01 1.98924544e-06 3.75125274e-05 1.56336097e-03
 3.69293343e-03 1.42601013e-01 1.43981070e-03 1.69405005e-06]

Number of samples: 10
First-order indices: [8.35561003e-01 4.96902726e-04 1.78652937e-03 7.69974030e-02
 4.30161649e-06 1.26526642e-03 8.15101801e-02 7.21941070e-08]

Total indices: [8.38883777e-01 5.02172947e-04 1.71870783e-03 7.71203149e-02
 2.75969474e-06 1.24206144e-03 8.14787489e-02 3.44898333e-08]

Number of samples: 15
First-order indices: [8.44709980e-01 1.06700504e-06 3.24133283e-04 2.58552242e-02
 2.23341376e-03 3.86293202e-02 6.01782145e-02 5.34915388e-03]

Total indices: [8.74604841e-01 3.71130584e-07 2.55879659e-04 3.01456625e-02
 5.29969067e-03 3.77321321e-02 6.85605400e-02 5.67403775e-03]

Number of samples: 20
First-order indices: [8.19892797e-01 6.31571231e-06 7.03944889e-04 4.51452468e-02
 1.90837343e-05 4.36003556e-02 4.87573143e-02 7.73569634e-03]

Total indices: [8.67075125e-01 9.80602097e-06 8.35516378e-04 5.52205345e-02
 8.60621761e-06 4.94937591e-02 5.47902672e-02 8.77088441e-03]

Number of samples: 25
First-order indices: [ 8.02887662e-01  1.74078006e-04 -1.43141177e-06  4.93418341e-02
 -1.96451403e-05  4.97651557e-02  4.16374778e-02  9.06475717e-03]

Total indices: [8.63180007e-01 1.96386623e-04 6.31011552e-06 6.45979479e-02
 4.67704986e-04 5.76288295e-02 4.96236828e-02 1.20119411e-02]

Number of samples: 30
First-order indices: [ 8.02221590e-01  6.72364141e-06  2.34418924e-05  4.24600731e-02
 -5.82003505e-06  5.15332400e-02  5.04175093e-02  8.07873797e-03]

Total indices: [8.60553632e-01 5.28448495e-06 6.52269168e-05 5.58684586e-02
 6.88091676e-07 5.91501359e-02 6.02273886e-02 1.19472721e-02]

Number of samples: 35
First-order indices: [ 8.22678633e-01  2.20211882e-05 -2.94498646e-05  3.88931737e-02
  7.17476494e-06  4.81593111e-02  4.37957354e-02  7.66566672e-03]

Total indices: [8.72433374e-01 7.30928459e-06 7.36740106e-06 5.00135221e-02
 1.15173542e-05 5.34646444e-02 5.22718464e-02 1.16243393e-02]

Number of samples: 40
First-order indices: [ 8.25248289e-01  2.69045540e-05 -3.23954277e-05  4.22595624e-02
  4.24453179e-05  4.65729453e-02  4.22731887e-02  8.24864990e-03]

Total indices: [8.72544574e-01 2.07191606e-05 2.01344175e-05 5.36659637e-02
 5.36527662e-05 5.15890532e-02 4.84399706e-02 1.24735957e-02]
../_images/800c2b96a233524c21910a299bd6eb544cc5071f3e977ba6d380768fd240e0df.png ../_images/6076bb51b76140f35dd788a9d692f0dde92f4506f5fa40d50d0cc2f87e367d12.png ../_images/be3fe2fa7675de015284f19c321035ff7e8bdb1ad8a7211a6cf6febf7084671e.png ../_images/d0e36f81e7e73ab2b0887b6dbc16077aa62e27b13fc7b23053851adf709c8509.png ../_images/cd500e5543ef3ddf223058a27f90732874d355b6a3ef31b5d2db4615767d4131.png ../_images/4323dc6504f918afc1059004978d9a06f689f9ed57200bfe1b11aecf8c05e906.png ../_images/e7a64e77b9a186898cdd296c231bc3b22faff3c34739b04d7d2236db511eac51.png ../_images/3a50e2686eebbfa0960a520ff865d3874d2168aad077dc0c6260ba464e821bd1.png

The Kriging model cannot model the Water Flow function well enough to accurately calculate the Sobol’ indices with a few number of samples. As the number of samples increases, the Kriging model becomes more accurate and the Sobol’ indices calculated using the Kriging model match closely with those calculated using the true function. In this case, almost 40 samples are required to accurately calculate the Sobol’ indices as there are more number of variables in this function and there are more interaction effects to capture. However, the computational cost is still much cheaper than calculating the Sobol’ indices using the true function.