SOBOL Sensitivity Analysis

The following example converts Parameter Study to perform a global variance-based sensitivity analysis following the method of Saltelli (2002).

Input File

To perform a the sensitivity analysis only two minor changes to the input file from Parameter Study need to occur: the sampler needs be changed and another statistics calculation object needs to be added, the remainder of the file remains the same. The exception being that the sampler name in a few of the objects was updated to use the "sobol" Sampler name, as defined in the following sub-section.

Sobol Sampler

The Sobol sampling scheme consists of using a sample and re-sample matrix to create a series of matrices that can be used to compute first-order, second-order, and total-effect sensitivity indices. The first-order indices are the portion of the variance in the quantity of interest (e.g., TavgT_{avg}) due to the variance of the uncertain parameters (e.g., γ\gamma). The second-order indices are the portion of the variance in the quantity of interest due to the variance two uncertain parameters interacting (e.g., γ\gamma as ss). Finally, the total-effect indices are the portion of the variance of the quantity of interest due to the variance of an uncertain parameter and all interactions with the other parameters. For complete details of the method please refer to Saltelli (2002).

The Sobol object requires two input samplers to form the sample and re-sample matrix. Thus, the Samplers block contains three sample objects. The first two are used by the third, the "sobol" sampler, which is the Sampler used by the other objects in the simulation.

[Samplers]
  [hypercube_a]
    type = LatinHypercube
    num_rows = 10000
    distributions = 'gamma q_0 T_0 s'
    seed = 2011
  []
  [hypercube_b]
    type = LatinHypercube
    num_rows = 10000
    distributions = 'gamma q_0 T_0 s'
    seed = 2013
  []
  [sobol]
    type = Sobol
    sampler_a = hypercube_a
    sampler_b = hypercube_b
  []
[]
(contrib/moose/modules/stochastic_tools/examples/sobol/main.i)

The sobol method implemented here requires n(2k+2)n(2k+2) model evaluations, where kk is the number of uncertain parameters (4) and nn is the number of replicates (10,000). Therefore, for this example 100,000 model evaluations were performed to compute the indices.

Sobol Statistics

The values of the first-order, second-order, and total-effect are computed by the SobolReporter object. This object is a Reporter, as such it is added to the Reporters block. The output of the object includes all the indices for each of the vectors supplied.

[Reporters]
  [results]
    type = StochasticReporter
    outputs = none
  []
  [stats]
    type = StatisticsReporter
    reporters = 'results/results:T_avg:value results/results:q_left:value'
    compute = 'mean'
    ci_method = 'percentile'
    ci_levels = '0.05 0.95'
  []
  [sobol]
    type = SobolReporter
    sampler = sobol
    reporters = 'results/results:T_avg:value results/results:q_left:value'
    ci_levels = '0.05 0.95'
  []
[]
(contrib/moose/modules/stochastic_tools/examples/sobol/main.i)

Results

If JSON output is enabled, the SobolReporter object will write a file that contains columns of data. Each column comprises of the computed indices for the quantities of interest. For example, Listing 1 is and example output from the SobolReporter object for this example problem.

Listing 1: Computed Sobol indices for the example heat conduction problem.

{
    "reporters": {
        "sobol": {
            "confidence_intervals": {
                "levels": [
                    0.05,
                    0.95
                ],
                "method": "percentile",
                "replicates": 10000,
                "seed": 1
            },
            "indices": [
                "FIRST_ORDER",
                "TOTAL",
                "SECOND_ORDER"
            ],
            "num_params": 4,
            "type": "SobolReporter",
            "values": {
                "results_results:T_avg:value": {
                    "type": "SobolIndices<double>"
                },
                "results_results:q_left:value": {
                    "type": "SobolIndices<double>"
                }
            }
        },
        "stats": {
            "confidence_intervals": {
                "levels": [
                    0.05,
                    0.95
                ],
                "method": "percentile",
                "replicates": 10000,
                "seed": 1
            },
            "type": "StatisticsReporter",
            "values": {
                "results_results:T_avg:value_MEAN": {
                    "stat": "MEAN",
                    "type": "std::pair<double, std::vector<double> >"
                },
                "results_results:q_left:value_MEAN": {
                    "stat": "MEAN",
                    "type": "std::pair<double, std::vector<double> >"
                }
            }
        }
    },
    "time_steps": [
        {
            "sobol": {
                "results_results:T_avg:value": {
                    "FIRST_ORDER": [
                        [
                            0.7079682315460529,
                            0.022003810752935962,
                            0.06629141722528406,
                            0.12306349437144783
                        ],
                        [
                            [
                                -2.778576837390401,
                                -0.0923840362252813,
                                -0.33953641462917444,
                                -0.6272646085746221
                            ],
                            [
                                3.1575055384065687,
                                0.10810117542142149,
                                0.44730763561354525,
                                0.7136139871038889
                            ]
                        ]
                    ],
                    "SECOND_ORDER": [
                        [
                            [
                                0.7079682315460529,
                                5.3391885740933,
                                2.213176168087708,
                                3.2444306898982767
                            ],
                            [
                                5.3391885740933,
                                0.022003810752935962,
                                0.10021434300595027,
                                0.03146757769484784
                            ],
                            [
                                2.213176168087708,
                                0.10021434300595027,
                                0.06629141722528406,
                                0.011329509079781405
                            ],
                            [
                                3.2444306898982767,
                                0.03146757769484784,
                                0.011329509079781405,
                                0.12306349437144783
                            ]
                        ],
                        [
                            [
                                [
                                    -2.778576837390401,
                                    -6.454549229380251,
                                    -8.013830000296329,
                                    -8.121033794023335
                                ],
                                [
                                    -6.454549229380251,
                                    -0.0923840362252813,
                                    -0.8743948178853919,
                                    -1.1316087615874393
                                ],
                                [
                                    -8.013830000296329,
                                    -0.8743948178853919,
                                    -0.33953641462917444,
                                    -1.6300124540179968
                                ],
                                [
                                    -8.121033794023335,
                                    -1.1316087615874393,
                                    -1.6300124540179968,
                                    -0.6272646085746221
                                ]
                            ],
                            [
                                [
                                    3.1575055384065687,
                                    6.594936922311948,
                                    8.192699613090781,
                                    8.316284425832205
                                ],
                                [
                                    6.594936922311948,
                                    0.10810117542142149,
                                    0.9183264421363888,
                                    1.2152704067270645
                                ],
                                [
                                    8.192699613090781,
                                    0.9183264421363888,
                                    0.44730763561354525,
                                    1.8568944936372327
                                ],
                                [
                                    8.316284425832205,
                                    1.2152704067270645,
                                    1.8568944936372327,
                                    0.7136139871038889
                                ]
                            ]
                        ]
                    ],
                    "TOTAL": [
                        [
                            0.7170011239748011,
                            -0.23410520806105195,
                            0.014050245820542373,
                            -0.06683002867472965
                        ],
                        [
                            [
                                0.049115893640706565,
                                -2.240148385235004,
                                -1.5694929791243908,
                                -1.7111043447276533
                            ],
                            [
                                1.8641209326473396,
                                3.9534370745538534,
                                3.4019193072456204,
                                3.585576242089134
                            ]
                        ]
                    ]
                },
                "results_results:q_left:value": {
                    "FIRST_ORDER": [
                        [
                            0.7389200088104518,
                            0.020349594867369668,
                            0.22549127987918816,
                            0.025694489604817107
                        ],
                        [
                            [
                                -4.091739873702765,
                                -0.16637293825878882,
                                -1.3830155445581565,
                                -0.21054326591717482
                            ],
                            [
                                3.947523492747353,
                                0.1665429587899574,
                                1.3267423877123623,
                                0.19340501736478963
                            ]
                        ]
                    ],
                    "SECOND_ORDER": [
                        [
                            [
                                0.7389200088104518,
                                1.1482851633655828,
                                1.3355470568848986,
                                1.0958597073025564
                            ],
                            [
                                1.1482851633655828,
                                0.020349594867369668,
                                -0.08344803820392069,
                                -0.04418713548380046
                            ],
                            [
                                1.3355470568848986,
                                -0.08344803820392069,
                                0.22549127987918816,
                                -0.1269011054809108
                            ],
                            [
                                1.0958597073025564,
                                -0.04418713548380046,
                                -0.1269011054809108,
                                0.025694489604817107
                            ]
                        ],
                        [
                            [
                                [
                                    -4.091739873702765,
                                    -9.139881163051552,
                                    -11.458210899796187,
                                    -8.682853541799657
                                ],
                                [
                                    -9.139881163051552,
                                    -0.16637293825878882,
                                    -1.656948876741212,
                                    -0.3918542256440673
                                ],
                                [
                                    -11.458210899796187,
                                    -1.656948876741212,
                                    -1.3830155445581565,
                                    -1.6848825449903648
                                ],
                                [
                                    -8.682853541799657,
                                    -0.3918542256440673,
                                    -1.6848825449903648,
                                    -0.21054326591717482
                                ]
                            ],
                            [
                                [
                                    3.947523492747353,
                                    9.818740534575209,
                                    11.68784380386844,
                                    9.752682930726047
                                ],
                                [
                                    9.818740534575209,
                                    0.1665429587899574,
                                    2.1155879208584008,
                                    0.5543252646398215
                                ],
                                [
                                    11.68784380386844,
                                    2.1155879208584008,
                                    1.3267423877123623,
                                    2.119622204584661
                                ],
                                [
                                    9.752682930726047,
                                    0.5543252646398215,
                                    2.119622204584661,
                                    0.19340501736478963
                                ]
                            ]
                        ]
                    ],
                    "TOTAL": [
                        [
                            0.7167820377142026,
                            -0.06820553511007565,
                            0.28209547467004326,
                            -0.027625629805629437
                        ],
                        [
                            [
                                -0.8931482758355531,
                                -5.2654851300237615,
                                -3.2439900542445317,
                                -5.111217986149409
                            ],
                            [
                                2.5115384112506294,
                                7.2346296356909665,
                                4.793000313371506,
                                6.889820155566966
                            ]
                        ]
                    ]
                }
            },
            "stats": {
                "results_results:T_avg:value_MEAN": [
                    197.08718601660595,
                    [
                        188.41345732966371,
                        205.97591428995048
                    ]
                ],
                "results_results:q_left:value_MEAN": [
                    174.12450852747952,
                    [
                        161.95826742103895,
                        186.5833699460805
                    ]
                ]
            },
            "time": 2.0,
            "time_step": 2
        }
    ]
}
(contrib/moose/modules/stochastic_tools/examples/sobol/gold/main_out.json)

For each set of indices (FIRST_ORDER, SECOND_ORDER, and TOTAL) contains a pair entries: the first is the values computed and the second corresponds to the 5% and 95% confidence intervals. This problem examines four uncertain parameters, so each element of the set of indices corresponds to the parameter indicated in the input file (γ\gamma, q0q_0, T0T_0, and ss). See SobolReporter for further information regarding the output.

For the problem at hand, the first-order and second-order indices for the two quantities of interest are presented in Table 1 and Table 2. The diagonal entries are the first-order indices and the off-diagonal terms are the second-order indices. For example for TavgT_{avg} the first order-index for γ\gamma is S1=0.763S_1 = 0.763 and the second-order index S1,2=0.014S_{1,2} = 0.014 for γ\gamma interacting with q0q_0. The negative values are essentially zero, if more replicates were executed these numbers would become closer to zero.

python ../../python/visualize_sobol.py main_out.json  --markdown-table \
--values results_results:T_avg:value --stat second_order \
--param-names '$\gamma$' '$q_0$' '$T_0$' '$s$' \
--number-format .3g

Table 1: First-order and second-order Sobol indices for TavgT_{avg}.

Si,jS_{i,j} (5.0%, 95.0%) CIγ\gammaq0q_0T0T_0ss
γ\gamma0.69 (0.626, 0.77)---
q0q_00.00108 (-0.122, 0.122)0.00737 (0.00598, 0.00892)--
T0T_00.0114 (-0.126, 0.148)0.00412 (-0.003, 0.0111)0.0902 (0.0807, 0.102)-
ss0.00538 (-0.153, 0.16)0.0039 (-0.0101, 0.018)0.00353 (-0.0207, 0.0283)0.201 (0.181, 0.225)
python ../../python/visualize_sobol.py main_out.json  --markdown-table \
--values results_results:q_left:value --stat second_order \
--param-names '$\gamma$' '$q_0$' '$T_0$' '$s$' \
--number-format .3g

Table 2: First-order and second-order Sobol indices for qleftq_{left}.

Si,jS_{i,j} (5.0%, 95.0%) CIγ\gammaq0q_0T0T_0ss
γ\gamma0.815 (0.765, 0.874)---
q0q_00.00787 (-0.094, 0.108)0.0267 (0.0245, 0.0292)--
T0T_00.0188 (-0.0959, 0.132)-0.000773 (-0.00685, 0.00532)0.135 (0.126, 0.145)-
ss0.011 (-0.0887, 0.109)0.000268 (-0.0021, 0.00265)-0.000539 (-0.00844, 0.00736)0.00423 (0.00312, 0.00538)

The data in these two tables clearly indicates that a majority of the variance of both quantities of interest are due to the variance of γ\gamma, ss also affects TavgT_{avg} and T0T_0 affects qleftq_{left}. Additionally, a small contribution of the variance is from a second-order interaction, S1,2S_{1,2}, between γ\gamma and T0T_0. The importance of γ\gamma, T0T_0, and ss is further evident by the total-effect indices, as shown in Table 3.

python ../../python/visualize_sobol.py main_out.json --markdown-table --stat total \
--names '{"results_results:T_avg:value":"$T_{avg}$", "results_results:q_left:value":"$q_{left}$"}' \
--param-names '$\gamma$' '$q_0$' '$T_0$' '$s$' \
--number-format .3g

Table 3: Total-effect Sobol indices for TavgT_{avg} and qleftq_{left}.

STS_T (5.0%, 95.0%) CIγ\gammaq0q_0T0T_0ss
TavgT_{avg}0.707 (0.673, 0.736)0.0217 (-0.0861, 0.111)0.115 (0.016, 0.196)0.208 (0.12, 0.28)
qleftq_{left}0.837 (0.825, 0.848)0.0415 (-0.0242, 0.0999)0.157 (0.0988, 0.209)0.0137 (-0.0537, 0.0743)

To help visualize the sensitivities, visualize_sobol.py can also represent it as a bar pot or heat map:

python ../../python/visualize_sobol.py main_out.json --bar-plot --log-scale --stat total \
--names '{"results_results:T_avg:value":"$T_{avg}$", "results_results:q_left:value":"$q_{left}$"}' \
--param-names '$\gamma$' '$q_0$' '$T_0$' '$s$'
python ../../python/visualize_sobol.py main_out.json --heatmap --log-scale --stat second_order \
--names '{"results_results:T_avg:value":"$T_{avg}$", "results_results:q_left:value":"$q_{left}$"}' \
--param-names '$\gamma$' '$q_0$' '$T_0$' '$s$'

References

  1. Andrea Saltelli. Making best use of model evaluations to compute sensitivity indices. Computer Physics Communications, 145(2):280–297, 2002. URL: https://doi.org/10.1016/S0010-4655(02)00280-1.[BibTeX]