% 参数敏感性分析
clear;
%进行样本采样 M为参数样本
N=3000;%采样次数
Np=4;%参数数目
Par={'a','percent',[0 1];
'b','percent',[0 100];
'c','percent',[0 1];
'd','percent',[0,20]};
Rnom = zeros(2,Np);
for k=1:Np
Rnom(:,k) = [Par{k,3}(1) ; Par{k,3}(2)];
end
SampleMethod= 'Uniform';
% SampleMethod= 'LatinHypercube';
switch SampleMethod
case 'Uniform'
M = zeros(N,Np);
for k = 1:Np
pdfun = makedist('Uniform','lower',Rnom(1,k),'upper',Rnom(2,k));
M(:,k) = random(pdfun,1,N);
end
case 'LatinHypercube'
M = lhsdesign(N,Np).*( ones(N,1)*(Rnom(2,:)-Rnom(1,:)) ) + ones(N,1)*Rnom(1,:);
end
% A,B:每一列为1类参数，每一行为一个样本
% A = [p1(1) p2(1) ... pNp(1);
% p1(2) p2(2) ... pNp(2);
% ...;
% p1(N/2) p2(N/2) ... pNp(N/2)]
%%
%以文档中的M进行举例
clear;
M=[0.5 0.5 0.5;...
0.75 0.25 0.25;...
0.25 0.75 0.75;...
0.375 0.375 0.625;...
0.5 0.5 0.5;...
0.25 0.75 0.75;...
0.75 0.25 0.25;...
0.875 0.375 0.125];
N=size(M,1);
Np=size(M,2);
Ohandle=@test1;
A = M(1:N/2,:);
B = M(N/2+1:N,:);
% BAi B中的第i列来自于A
BAi = cell(1,Np);
for k = 1:Np
BAi{k} = B;
BAi{k}(:,k) = A(:,k);
end
% ABi A中的第i列来自于B
ABi = cell(1,Np);
for k = 1:Np
ABi{k} = A;
ABi{k}(:,k) = B(:,k);
end
YA = sens_monte(Ohandle,A,1);
YB = sens_monte(Ohandle,B,1);
Y=[YA;YB];
VY=var(Y,1);
SY = zeros(Np,1);
STY = zeros(Np,1);
YABi = cell(1,Np);
YBAi = cell(1,Np);
sensmethod='Saltelli';
switch sensmethod
case 'sobol'
f02 = mean(Y)^2;
for k = 1:Np
YBAi{k} = sens_monte(Ohandle,BAi{k});
YABi{k}= sens_monte(Ohandle,ABi{k});
SY(k) = (mean( YA.*YBAi{k}) - f02 )/VY;%一阶敏感度
STY(k) = mean( YA.*(YA - YABi{k}) )/VY;%全局敏感度
end
case 'Saltelli'
for k = 1:Np
YABi{k} = sens_monte(Ohandle,ABi{k},1);
SY(k) = mean(YB.*(YABi{k} - YA),1)./VY;
STY(k) = mean((YA - YABi{k}).^2)/(2*VY);
end
case 'Jansen'
for k = 1:Np
YABi{k} = sens_monte(Ohandle,ABi{k},1);
SY(k) = 1 - mean((YB - YABi{k}).^2)/(2*VY);
STY(k) = mean( (YA - YABi{k}).^2)/(2*VY);
end
end
%%
%绘图
barh(SY)
set(gca,'ytick',1:Np,'FontSize',10)
title('一阶敏感性指标');
for i=1:Np
if SY(i)>0.3
alignment='right';
color1='white';
else
alignment='left';
color1='blue';
end
text(SY(i),i,[' ' num2str(SY(i),3)],'HorizontalAlignment',alignment,'color',color1,'FontSize',10)
end
