【MATLAB第64期】【保姆级教程】基于MATLAB的SOBOL全局敏感性分析模型运用（含无目标函数，考虑代理模型）

版本更新：

2023/8/5：
1.因BP作为代理模型不稳定，经过测试，libsvm比rf /bp 效果稳定且精度较高。故用libsvm替换原来的bp，并增加选择libsvm的原因。
2.增加用libsvm作为代理模型的sobol敏感结果对比分析及验证内容。
3.增加遍历来筛选sobol样本数量，进行结果比对。
4.单独以sobol作为一章。因为内容比较多，为了便于观看，后期会更新其他的全局敏感性分析方法。(PAWN,GSA等)

引言

在前面几期，介绍了局部敏感性分析法，本期来介绍sobol全局敏感性分析模型，因还在摸索中，其他全局敏感性模型敬请期待。
【MATLAB第31期】基于MATLAB的降维/全局敏感性分析/特征排序/数据处理回归问题MATLAB代码实现(持续更新)
【MATLAB第32期】【更新中】基于MATLAB的降维/全局敏感性分析/特征排序/数据处理分类问题MATLAB代码实现
【MATLAB第63期】基于MATLAB的改进敏感性分析方法IPCC，拥挤距离与皮尔逊系数法结合实现回归与分类预测

一、SOBOL（有目标函数）

（1）评价指标

评价指标包括：一阶影响指数S，总效应指数ST**

*一阶影响指数S：*显示由各个输入变量的方差产生的因变量的方差，根据一阶影响指数可以量化单个变量对模型的敏感程度

总效应指数ST：显示由每个输入变量的方差及其与其他输入变量的相互作用而产生的因变量的方差。

每个因变量和所有变量的 Sobol 指数都显示在专用的 Sobol 图中
，其中直方图按总效应指数ST排序。因变量对具有最高总效应指数ST的输入变量最敏感。

输入变量的总效应指数ST和一阶影响指数S之间的差异可以衡量该输入与其他输入变量之间相互作用的效果。

（2）运行思路

A、设定目标函数（3个变量，即维度D=3）
Y=X1^2+2*X2+X3-1

y=x(1)^2+2*x(2)+x(3)-1;

B、设定变量上下限

VarMin=[0 0 0];%各个参数下限
VarMax=[10  10 10];%各个参数上限

C、设定sobol其他参数

M=D*2 %2倍D数量的矩阵，提高样本丰富度
%则此时相当于共6个输入变量
VarMin=[0 0 0 0 0 0];%各个参数下限
VarMax=[10  10 10 10 10 10];%各个参数上限
nPop=4;%采样数量，样本数量（数量设置越大，准确率越高。为了方便展示数值，选取nPop=4）

增加对nPop样本数量筛选功能。

对nPop采用遍历形式，即4:50:1000，50为间隔数。
在这里插入图片描述
其中，第四条线代表所有变量S /ST之和。

通过nPop=4结果可以看出，S /ST值相对稳定性较差。基本在nPop=200左右保持稳定。为了方便展示数值，选取nPop=4。

D、生成sobol序列样本数据

这里要说，除了sobol序列函数可以生成样本数据，其他也可以，比如正交设计、超立方抽样等等。很关键！！！举一反三即可

1、生成多组N*M（即N行6列）的样本矩阵p。用自带函数sobolset生成。

p= sobolset(M)

**p  矩阵形式： 9007199254740992x6 sobolset**

2、筛选nPop*M（即4行6列）的样本矩阵R。
两种思路，第一种直接选取前nPop行的p采样数据，优点是方便快捷，但是缺点是样本不随机，并没有考虑上下限对样本的影响。

% 第一种
R=p(1:nPop,:);%选取前nPop行

在这里插入图片描述

% 第二种
for i=1:nPop% 选取前nPop行被上下限空间处理后的样本
    r=p(i,:);
    r=VarMin+r.*(VarMax-VarMin);
   R=[R; r]; 
end

本例因为最小和最大值一样，如果最小值和最大值均为0/1，则两种方法结果一致。
第二种方法R矩阵为：
在这里插入图片描述

明显第二种方法更符合逻辑。

E、R样本拆分变换（提高样本丰富度）

1.将矩阵的前D列设置为矩阵A，后D列设置为B列，在我们的例子中就是矩阵m的前3列设置为矩阵A，后3列设置为矩阵B。

A=R(:,1:D);% 每行代表一组参数，其中每列代表每组参数的一个参数；行数就代表共有几组参数
B=R(:,D+1:end);

在这里插入图片描述

2.构造nPop*D的矩阵ABi(i = 1,2,…,D)，即用矩阵B中的第i列替换矩阵A的第i列，以本体为例：

for i=1:D
    tempA=A;
    tempA(:,i)=B(:,i);
    AB(1:nPop,1:D,i)=tempA;
end

AB=
在这里插入图片描述

经过这三步我们构造了A、B、AB1、AB2、AB3这五个矩阵，这样我们就有(D + 2) * nPop (即20)组x1、x2、x3输入数据，因此我们将有20组Y值。将上述的数据带入函数，这里详细的计算过程就不描述了。根据输入我们得出对应的Y值矩阵。
F、计算所有样本对应的Y值

for i=1:nPop
    YA(i)=myfun(A(i,:));  %A矩阵对应的YA值
    YB(i)=myfun(B(i,:));%B矩阵对应的YB值
    for j=1:D
        YAB(i,j)=myfun(AB(i,:,j));%YAB矩阵对应的YAB值
    end
end

[YA YB YAB1 YAB2 YAB3 ]组合起来
依次各列数据代表YA YB YAB1 YAB2 YAB3值
在这里插入图片描述

G、一阶影响指数S值、总效应指数ST值计算

1.计算公式：
在这里插入图片描述
var方差函数为matlab自带

2.一阶影响指数S值

VarX=zeros(D,1);% S的分子
S=zeros(D,1);
VarY=var([YA;YB],1);% S的分母。 计算基于给定的样本总体的方差(EXCEL var.p())
for i=1:D
    for j=1:nPop
         VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));
    end
     VarX(i)=1/nPop*VarX(i);
     S(i)=VarX(i)/VarY;  %一阶影响指数
end

在这里插入图片描述

3.总效应指数ST值
在这里插入图片描述绘图：

H、分析

1.回到目标函数：y=x(1)^2+2*x(2)+x(3)-1;
可根据数学所学知识，得X1项为幂函数，X2项为系数=2的一次函数，X3项为系数=1的一次函数
根据常识即理论可知，敏感度排序X1>X2>X3
通过SOBOL的总效应指数ST柱状图结果也可以证实以上结论。

2.其次，一阶影响指数S中，第二个变量对应的S为负值，表示单个变量对因变量的敏感度，即所谓的局部敏感性分析法。
|X1|>|X2|>|X3|
而全局要考虑不同变量对因变量的影响，即ST定义——每个输入变量的方差及其与其他输入变量的相互作用而产生的因变量的方差。

3.输入变量的总效应指数ST和一阶影响指数S之间的差异可以衡量该输入与其他输入变量之间相互作用的效果。

二、SOBOL（无目标函数）

1.解决思路

（1）针对简单线性数据及非线性数据，用函数拟合得到公式，随后思路与上面一致。
（2）无法拟合得到公式，即复杂非线性函数，需要通过借用机器学习模型，作为训练学习模型（黑箱子模型）
本文具体研究攻克第二种情况
有个前提（模型拟合性较好，对应数据较好）
即训练学习模型，训练集和测试集拟合效果很棒。
如果拟合效果差，SOBOL分析结果一定存在较大误差。

2.模型选择

1、选用libsvm模型作为代理模型
原因：
（1）代理模型讲究运行效率快、精度高、模型简单。libsvm符合以上情况，仅有的两超参数c、g经验值结果普遍较好，基本不用调参。
（2）进行对比以bp为代表的神经网络模型，因其机理中涉及随机初始的权值阈值等参数，会让模型不够稳定。
（3）进行对比的rf随机森林模型，训练效果远差于bp /libsvm ，且参数调整较为复杂。
（4）深度学习模型更适合大数据模型，对于平时用的小数据，传统模型不见得效果比深度模型差，其次深度学习运行时间、模型复杂程度，调参难度等问题明显无法与传统方法相比。
故综合以上原因，选择libsvm作为代理模型。

libsvm运行插件在往期文章分享过，可直接下载。

【MATLAB第8期】#源码分享//基于MATLAB的最简易且不用安装的支持向量机LIBSVM函数及SVM分类回归模型参数设置

2.数据设置：常用的案例数据，103*8 ，前7列代表输入变量，最后1列代表因变量。

3.选用模型后，几个点需要注意：
（1）数据固定，即训练样本/测试样本固定，所代表的模型评价才够稳定。

（2）使用固定算子函数代码（神经网络代理模型是必要的），即开头代码为： rng default 或者rng（M）等，M根据实际测试效果确定。可固定输出结果，保证运行结果一致。此一致代表此刻你打开的matlab，在不关闭情况下每次运行结果一致。跟matlab版本有关，系统版本，以及电脑有关。

（3）最为关键的一点，变量的上下限不能超过案例数据的上下限，为了保证模型的普适性和有效性！！！
比如案例数据的训练样本中， X1-X7的最小值为：
[137 0 0 160 4.4 708 650]
X1-X7的最大值为：
[374 193 260 240 19 1049.90 902]
那么你的sobol序列生成的数据也只能在这个范围，才能保证代理模型的有效性。
（4）生成样本的数量当然以多为好，但不能跟案例数据样本数量差距太大，减少偶然性。
（5）代理模型效果（以libsvm为例）

在这里插入图片描述
训练集数据的R2为：0.99787
测试集数据的R2为：0.96186
训练集数据的MAE为：0.32795
测试集数据的MAE为：1.2748
训练集数据的MBE为：0.019637
测试集数据的MBE为：-1.1294

个人认为，训练集和测试集R2如果均大于0.9还是可以的，评价指标好坏全凭主观意思。包括评价指标的选择，不一定是R2，R2更适合这样的波动的曲线。

(6)保存模型所需要的变量

save svmnet model ps_output ps_input
通过sobol生成样本进行仿真预测。

3.SOBOL模型分析

（1）sobol参数设置

%% 设定：给定参数个数和各个参数的范围
D=7;% 7个参数
M=D*2;%
nPop=80;% 采样点个数,跟训练样本数量大概一致
VarMin=[137	0	0	160	4.4	708	650];%各个参数下限
VarMax=[374	193	260	240	19 1049.90	902];%各个参数上限

nPop=10:50：500
nPop数量遍历结果：
在这里插入图片描述
nPop=10 / 60时，S/ST值结果不稳定，样本数=200时偏于稳定
本文便于分析，选取nPop=80。

（2）运行结果

敏感程度（libsbm作为代理模型）：X4>X3≈X1>X7>X6>X5>X2
在这里插入图片描述
跟原先使用BP模型，分析结果进行对比：
敏感程度（BP作为代理模型）： X3≈X1>X4>X2>X7>X6>X5

在这里插入图片描述
最显著区别是，关于X4变量的敏感程度的区别，其次是X6-X7变量的敏感程度的区别。
两者结果不同，需要通过控制变量法，剔除部分变量，看代理模型的训练效果是否能够印证sobol分析的结果。

4.SOBOL结果验证

验证方法有很多，其中极差分析法是相对比较理想的方法。当然极差分析法也可以直接取代sobol进行分析，原理就是通过控制变量改变X1-X7参数固定比例的值，然后看对Y结果的影响。比如对于X1来说，每个样本的X1增加-10% -5% 5% 10% 四种情况，来看对Y结果的影响。当然计算量比较大，不过结果是非常可观的，可以直接通过Y变化百分比来显示，不像sobol 的S /ST结果那么抽象。

本文为了分析简便，以结果为导向，通过筛选变量，对代理模型重新预测，看预测效果。当然结果也只能验证本文选用的代理模型不错，但无法证明最优。

（1）第一招：主打的就是和平（求同存异）
X1 / X3 /X4 三个变量相对都比较重要。用这3个变量分别测试libsvm /bp的预测效果。**
A、libsvm结果
在这里插入图片描述

B、BP结果
在这里插入图片描述
分析可得，libsvm结果符合逻辑，比所有参数作为输入结果差，但结果也基本满足要求，R2均在0.9左右。而BP结果不合理，测试集R2甚至为负数。

（2）第二招：是来找茬的
找到结果区别最大的X4变量作为研究对象，在筛除X4作为输入后，看两者结果。
A、libsvm结果
在这里插入图片描述
B、BP结果

可以看出，在筛除X4之后， libsvm明显测试结果下来了，不过R2=0.75还是交了个及格答卷。而BP结果仍然比较差。

（3）第三招：兵戎相见
通过分别把各自模型中比较满意的敏感度较高的变量作为输入，看效果。
libsvm：X1 X3 X4 X6 X7
BP：X1 X2 X3 X4
A、libsvm结果
在这里插入图片描述
B、BP结果

以上结果有两个关键结论：
a、 libsvm筛选的这5个变量结果的确好， R2均在0.98以上。而从BP这里可以看出，结果真不行。这里不是针对BP，而是指在场的神经网络模型。
b、根据 libsvm剔除X2、X5的测试集R2为0.98，大于X1-X7作为输入时测试的结果0.96，这里可以真正体现使用sobol的意义。对于S值和ST值均较小的变量，剔除后结果有所改善。