【MATLAB第111期】基于MATLAB的sobol全局敏感性分析方法二阶指数计算

一、简介

在MATLAB中计算Sobol二阶效应指数通常涉及到全局敏感性分析（Global Sensitivity Analysis, GSA），其中Sobol方法是一种流行的技术，用于评估模型输入参数的敏感性。Sobol二阶效应指数衡量的是两个参数之间的交互作用对模型输出的影响。

Sobol二阶效应指数的计算涉及到以下步骤：

1、生成Sobol序列，并在一阶的基础上，生成N=(2D+2)*npop个样本集，其中D为变量数， npop为采样的数量， N为总样本数。
在一阶基础上，N=(D+2)*nPop样本，包含A、AB和B矩阵。
二阶需要生成BA矩阵，用来评估二阶指数。

2、模型计算
可参考64期文章，利用sobol函数进行抽样，得到的X值，通过bp组成的代理模型进行计算。

3、计算Sobol指数：使用Sobol序列和模型输出，计算每个参数的一阶、二阶效应指数和总效应指数，其中，一阶和总效应指数较为好计算，二阶效应指数可以参考python的Salib库进行研究，

Sobol二阶效应指数的计算公式如下：
在这里插入图片描述
其中i为1：D，j为i+1：D，
4、计算效果

在这里插入图片描述

二、部分源码

S2计算代码如下：

 Vjk = mean(BAj .* ABk - A .* B) / var(y);
    Sj = first_order(A, ABj, B);
    Sk = first_order(A, ABk, B);
    S2 = Vjk - Sj - Sk;

核心参考代码如下（需要自行二次编译）：



    % Normalize the model output
    Y = (Y - mean(Y)) / std(Y);
    % Separate output values
    A = Y(1:2*D+2:end);
    B = Y((end-1):-(2*D+1):1);
    AB = zeros(length(Y)/ (2*D+2), D);
        BA = zeros(length(Y)/ (2*D+2), D);
        for j = 1:D
            AB(:, j) = Y((j+1):2*D+2:end);
            BA(:, j) = Y((j+1+D):2*D+2:end);
        end
    end

    % Calculate second order indices if required

        for j = 1:D
            for k = j+1:D
                Si.S2(j, k) = second_order(A, AB(:, j), AB(:, k), BA(:, j), B);
end
end
function S = first_order(A, AB, B)
    % First order estimator following Saltelli et al. 2010 CPC, normalized by
    % sample variance
    y = [A, B];
    if range(y) == 0
        warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
        S = 0;
        return;
    end
    S = mean(B .* (AB - A)) / var(y);
end

function S = total_order(A, AB, B)
    % Total order estimator following Saltelli et al. 2010 CPC, normalized by
    % sample variance
    y = [A, B];
    if range(y) == 0
        warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
        S = 0;
        return;
    end
    S = 0.5 * mean((A - AB) .^ 2) / var(y);
end

function S = second_order(A, ABj, ABk, BAj, B)
    % Second order estimator f ollowing Saltelli 2002
    y = [A, B];
    if range(y) == 0
        warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
        S = 0;
        return;
    end
    Vjk = mean(BAj .* ABk - A .* B) / var(y);
    Sj = first_order(A, ABj, B);
    Sk = first_order(A, ABk, B);
    S = Vjk - Sj - Sk;
end