【Matlab】Sobol灵敏度分析

计算一阶灵敏度与总灵敏度

参考链接
代码如下：
因为随机抽样，每次计算结果都不同

主函数

% function main
clc
clear
close all
%% 参数设置
n_p = 4;  % 待分析参数数目
n_base= 2000; % 参数样本数目
N = n_base*(n_p*2+1); % 模型/函数运行总次数
% 设置参数上下限，此处1<x1<1000,10<x2<100,1<x3<10,0<x4<2*pi
datamin=[1,10,1,0];%所有参数的下线
datamax=[1000,100,10,2*pi];%所有参数的上限
PS=PS_generate(n_base,n_p,datamin,datamax);% 参数样本空间（一）
comp_PS=PS_generate(n_base,n_p,datamin,datamax);% 参数样本空间（二），用于Monte Carlo采样
%% Sobol主程序
[Sob_1,Sob_t]=Sobol_cal(n_p,n_base,PS,comp_PS);
disp(Sob_1); %一阶灵敏度
disp(Sob_t); %总灵敏度
%负值或总和大于1为截断误差、非线性模型等因素造成

Sobol_obj.m

%% 样例函数
function out=Sobol_obj(kp)
% a function that works with the trial Sobol Method program
% the trial function is fx: x1+5*x2*x3^2;
x1=kp(1);
x2=kp(2);
x3=kp(3);
x4=kp(4);
out= x1+ x2*x3^2+sin(x4);
end

PS_generate.m

function PS=PS_generate(n_base,n_p,datamin,datamax)
PS=zeros(n_base,n_p);
for i=1:n_p
    PS(:,i)=ceil((datamin(i) + datamax(i).*rand(n_base,1)));
end

Sobol_cal.m

function [Sob_1,Sob_t]=Sobol_cal(n_p,n_base,PS,comp_PS)
% 计算模型输出
t=0;
for i=1:n_base
   t=t+1;
   kp(t,:)=PS(i,:);
   output(i,:)=Sobol_obj(kp(t,:)); %代入目标函数/模型计算
   for j=1:n_p
      t=t+1;
      kp(t,:)=[comp_PS(i,1:j-1),PS(i,j),comp_PS(i,j+1:n_p)]; % 构造Sobol抽样
      c_out_1(i,:,j)=Sobol_obj(kp(t,:)); %代入目标函数/模型计算
    
      t=t+1;
      kp(t,:) = [PS(i,1:j-1),comp_PS(i,j),PS(i,j+1:n_p)]; % 构造Sobol抽样
      c_out_t(i,:,j)=Sobol_obj(kp(t,:)); %代入目标函数/模型计算
   end
end
% t=N here;
%蒙特卡洛积分
n_out = size(output,2); 
f0 = zeros(1,n_out); 
D=zeros(1,n_out);  

for i = 1:n_base
 f0 = f0+output(i,:)/n_base; % 模型积分
 D = D+output(i,:).^2/n_base; % 计算模型输出方差
end
 
D=D-f0.^2; %模型输出方差
Dj=ones(n_p,1)*D;  
Dtotj=zeros(n_p,1); 
for i = 1:n_base
 for j = 1:n_p
 Dj(j,:)=Dj(j,:)-(output(i,:)-c_out_1(i,:,j)).^2/(2*n_base); %计算偏方差
 Dtotj(j,:)=Dtotj(j,:)+(output(i,:)-c_out_t(i,:,j)).^2/(2*n_base); %计算参数j的总方差
 end
end
%计算敏感度
Sob_1 = Dj./(ones(n_p,1)*D); %first order effect 一阶敏感度
Sob_t = Dtotj./(ones(n_p,1)*D); % total effect 总敏感度

注：一阶灵敏度与总灵敏度
（1）
一阶影响指数：显示由各个输入变量的方差产生的因变量的方差，根据一阶影响指数可以量化单个变量对模型的敏感程度
总效应指数：显示由每个输入变量的方差及其与其他输入变量的相互作用而产生的因变量的方差。
（2）
一阶影响指数：度量单个模型输入对输出方差的贡献
总阶影响指数：度量模型输入对输出方差的贡献，包括一阶与更高阶
（3）
一阶敏感指数表示单参数对模型输出方差的贡献率。
高阶敏感指数表示多个参数耦合效应对模型输出方差的贡献率。
全局敏感度指数：全局敏感度度量在修改一个元组时查询结果的最大变化，只与查询函数相关，并且独立于数据集本身。