【MATLAB第106期】#源码分享 | 基于MATLAB的有限差分算法的全局敏感性分析模型

【MATLAB第106期】#源码分享 | 基于MATLAB的有限差分法的全局敏感性分析模型

一、原理

有限差分法是一种数值方法，用于估计函数对输入参数的敏感性。在全局敏感性分析中，这种方法特别有用，因为它可以评估模型输出对所有输入参数变化的整体响应。以下是根据提供的MATLAB代码分析有限差分法实现全局敏感性分析的原理：

1、初始化和参数设置：清除所有变量（clear all），确保没有旧数据干扰。
初始化随机数生成器（rng(‘shuffle’)），确保随机抽样的可重复性。

2、定义样本数量和步长：
N：定义每个参数的样本数量，即要评估的模型运行次数。
h：定义有限差分的步长，即在每个参数上进行微小变化的量。

3、预分配内存：
y：存储模型输出（如增长率）的数组。
y_plus：存储扰动后模型输出的数组。
grad_y：存储模型输出梯度的数组。
Xs：存储归一化参数的数组。

4、设置参数的初始值和变化范围：
setvals：定义参数的初始值。
var：定义参数变化的百分比。
xl和xu：分别计算参数的下限和上限，用于随机抽样。

5、计算基准参数的模型输出：
使用for循环，对每个样本进行迭代。
在每次迭代中，随机抽样参数（Xs(jj,:)），并将其转换为实际参数值（params）。
计算模型输出（y(jj)），通常通过调用一个函数（如func_test1）实现。

6、移除不合理输出（进行约束）：
如果模型输出超出预期范围（如大于15），则重新抽样参数并重新计算，直到获得合理的输出。

7、计算扰动参数的模型输出：
对每个参数进行扰动，即在每个参数上增加或减少步长h。
计算扰动后的模型输出（y_plus），并同样移除不物理的输出。

8、计算梯度：
使用有限差分法计算每个参数的梯度（grad_y），即模型输出对每个参数变化的敏感度。梯度计算公式为：
grad

9、奇异值分解（SVD）：
对梯度矩阵进行奇异值分解，得到权重（U）、奇异值（S）和特征向量（V）。
第一列权重（w）和第二列权重（w2）分别表示模型输出对每个参数的相对重要性。

10、计算特征值和权重：
计算特征值（evalues），这些值表示模型输出对参数变化的敏感性。
计算权重（W和W2），这些值表示每个参数对模型输出的贡献度。

通过这种方法，可以全面了解模型输出对各个输入参数的敏感性，从而识别出对模型结果影响最大的参数。这对于优化模型参数、理解模型行为以及进行不确定性分析都非常有用。以下为目标函数：

function gamma = func_test1(x)
gamma=x(1)^2 + x(2)^2 + x(3)^2  + x(4)^2 +x(5)^2;                    
end

二、代码展示

%使用有限差分法计算参数变化的敏感性。通过计算梯度并使用SVD分析结果，可以了解模型输出对每个参数的敏感程度。
clear all;%清除所有变量
rng('shuffle');%初始化随机数生成器：将随机数生成器设置为“shuffle”，以确保每次运行时生成的随机数序列不同。
N = 1000;   %每个参数的样本数量
h = 1e-5;  %有限差分的步长y = zeros(N,1);                        %存储输出数组
Nparams = 5;
y_plus = zeros(N,Nparams);             %存储扰动后输出的数组
grad_y = zeros(Nparams,N);             %存储输出的梯度数组 
Xs = zeros(N,Nparams);                    %存储输出的梯度数组setvals = [1; 0.5; 0.4; 0.3; 1.5]; % 设置参数的初始值
var = 0.05; % 参数变化的百分比  5
xl = (1-var)*setvals;%参数的下限
xu = (1+var)*setvals;%参数的上限tic
%使用for循环计算基准参数和扰动参数的y。
for jj = 1:Njjrng(sum(100*clock)+pi*jj); % 在可接受的范围内随机抽样参数Xs(jj,:) = 2*rand(1,Nparams) - 1;params = 1/2*(diag(xu - xl)*Xs(jj,:)' + (xu + xl));% 计算函数的根，func_test函数用于计算给定参数的函数y。y(jj) = func_test1(params);% 设置阈值移除必要的根，%如果计算得到的y大于15，则重新随机采样参数直到获得合理的y。while y(jj) > 15                                                        Xs(jj,:) = 2*rand(1,Nparams) - 1;params = 1/2*(diag(xu - xl)*Xs(jj,:)' + (xu + xl));y(jj) = func_test1(params);end 
end%扰动参数
%略for jj = 1:N%使用有限差分计算appx梯度，%使用有限差分法计算增长率对每个参数的梯度。grad_y(:,jj) = (y_plus(jj, :) - y(jj))/h;
end
toc% 计算权重和特征值
%计算矩阵C的奇异值分解
[U,S,V] = svd(1/sqrt(N)*grad_y);
w = abs(U(:,1));
W=w/sum(w);
w2 = abs(U(:,2));W2=w2/sum(w2);   
%计算C的特征值
evalues = diag(S.^2);
%使用奇异值分解（SVD）计算grad_y矩阵的权重和特征值。
%找出最大和最小grad_y的差值以检查错误
diff_y = max(max(grad_y)) - min(min(grad_y));
%计算grad_y矩阵中最大值和最小值的差异，以检查可能的错误。
figure()
bar([W   W2]);xlabel('输入变量')ylabel('权重')legend('初始结果','扰动结果')title('有限差分敏感性分析')