【MATLAB第119期】基于MATLAB的KRR多输入多输出全局敏感性分析模型运用（无目标函数，考虑代理模型）

下一期研究SHAP的多输入多输出敏感性分析方法

一、SOBOL（无目标函数）

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

1.运行思路

A、设定KRR代理模型和变量上下限（6个变量，2个因变量，维度D=6）
(1).选用KRR模型作为代理模型*
代理模型讲究运行效率快、精度高、模型简单，适用于多输入多输出，与BP训练效果进行对比
(2).数据设置：常用的案例数据，103*8 ，前6列代表输入变量，最后2列代表因变量。

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

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

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

=== 输出1评价指标 ===
训练集: R2=0.9975, MAE=2.0306, MBE=-0.0019
测试集: R2=0.9887, MAE=4.5770, MBE=-0.7719

=== 输出2评价指标 ===
训练集: R2=0.9969, MAE=0.3082, MBE=-0.0001
测试集: R2=0.9937, MAE=0.4756, MBE=0.0310

BP:
在这里插入图片描述

=== 输出1评价指标 ===
训练集: R2=0.9901, MAE=5.0527, MBE=1.4881
测试集: R2=0.9813, MAE=5.4096, MBE=-0.8009

=== 输出2评价指标 ===
训练集: R2=0.9728, MAE=1.0644, MBE=0.1991
测试集: R2=0.9416, MAE=1.3128, MBE=0.1372

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

B、生成sobol序列样本数据

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

2、筛选nPop*M（即400行12列）的样本矩阵R。

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

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

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

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

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

AB = zeros(nPop, D, D_out)  % 400*600*2

AB=
在这里插入图片描述

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

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

[YA YB YAB1 YAB2 ]组合起来
依次各列数据代表YA YB YAB1 YAB2 值

E、一阶影响指数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:Dfor j=1:nPopVarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));endVarX(i)=1/nPop*VarX(i);S(i)=VarX(i)/VarY;  %一阶影响指数
end

3.总效应指数ST值

for i = 1:Dfor k = 1:D_outsum_term = 0;for j = 1:nPopsum_term = sum_term + (YA(j, k) - AB(j, i, k))^2;endST0(i, k) = (1/(2*nPop)) * sum_term / VarY(k);end
end