计算一阶灵敏度与总灵敏度
参考链接
代码如下：
因为随机抽样，每次计算结果都不同
主函数
% 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_pPS(:,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_baset=t+1;kp(t,:)=PS(i,:);output(i,:)=Sobol_obj(kp(t,:)); %代入目标函数/模型计算for j=1:n_pt=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_basef0 = f0+output(i,:)/n_base; % 模型积分D = D+output(i,:).^2/n_base; % 计算模型输出方差
endD=D-f0.^2; %模型输出方差
Dj=ones(n_p,1)*D;
Dtotj=zeros(n_p,1);
for i = 1:n_basefor j = 1:n_pDj(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）
一阶敏感指数表示单参数对模型输出方差的贡献率。
高阶敏感指数表示多个参数耦合效应对模型输出方差的贡献率。
全局敏感度指数：全局敏感度度量在修改一个元组时查询结果的最大变化，只与查询函数相关，并且独立于数据集本身。
)
