EDA365电子论坛网

标题: Sobol敏感性分析 [打印本页]

作者: Colbie 时间: 2020-3-18 09:49
标题: Sobol敏感性分析
本帖最后由 Colbie 于 2020-3-18 09:52 编辑

Sobol敏感性分析

function [ Y ] = test1( X)
%UNTITLED3 此处显示有关此函数的摘要
% 此处显示详细说明
x1=X(1) ;x2=X(2) ;x3=X(3) ;
Y=sin(x1)+7*(sin(x2))^2+0.1*x3^4*sin(x1);

end

function [ Y ] = sens_monte(f,M,Nd)
%f 目标函数
% M 参数样本
% Nd 一次采样数目
% 此处显示详细说明
N=size(M,1);
Y = zeros(N,Nd);
for j = 1: ceil(N/100)
parfor i = (j-1)*100+1: (j-1)*100+min(100,N- (j-1)*100)
   simoutarray(i)=f(M(i,: ));
end
end
for i = 1:N
Y(i,: ) = simoutarray(i); % Y  NxNd矩阵
end

% 参数敏感性分析
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

作者: Demyar 时间: 2020-3-18 18:42
Sobol敏感性分析。

作者: NNNei256 时间: 2020-3-19 18:27
Sobol敏感性分析

作者: Demyar 时间: 2020-3-20 18:33
看看楼主的代码。

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)