【MATLAB第67期】# 源码分享 | 基于MATLAB的morris全局敏感性分析

【MATLAB第67期】# 源码分享 | 基于MATLAB的morris全局敏感性分析

一、代码展示

clear all
npoint=100;%在分位数超空间中要采样的点数(计算次数iter=npoint*(nfac+1)
nfac=20;%研究函数的不确定因素数量
[mu, order] = morris_sa1(@(x)test_function(x), nfac, npoint)for t=1:size(mu,2)
W(1,t)=mu(1,t)/sum(mu);
end
figure()
bar(W)
grid on
xlabel('Feature index')
ylabel('Feature weight')%输出
%1）μ：
%每个因素，按降序排列。
%2）order：各因素的指标。考虑修正最后的因素，
%那些“mu”非常接近于零的数。function [ mu, order ] = morris_sa1( studied_function, nfac, npoint )%
%输出
%1）μ：
%每个因素，按降序排列。
%2）顺序：各因素的指标。考虑修正最后的因素，
%那些“mu”非常接近于零的人。
delta=1/npoint;
mini=delta/2;
maxi=mini+delta*(npoint-1);coord = 0:npoint-1;
for i=1:nfacpoints(1:npoint,i) = coord(randperm(length(coord)));
end
points = points/(npoint-1)*(maxi-mini)+mini;for i=1:npointtable_outputs(i,1) = studied_function(points(i,:)); % Output at the sampled point.for j=1:nfacif points(i,j) < 0.5 % If the coordinate is smaller than 0.5, a positive variation is appliedtable_outputs(i,1+j) = studied_function([points(i,1:j-1) points(i,j)+0.5 points(i,j+1:nfac)]); % Output after the variation of the j-th factor.table_ee(i,j) = (table_outputs(i,1+j)-table_outputs(i,1))/0.5; % Elementary effect of the j-th factor.else % If the coordinate if larger than 0.5, a negative variation is appliedtable_outputs(i,1+j) = studied_function([points(i,1:j-1) points(i,j)-0.5 points(i,j+1:nfac)]);table_ee(i,j) = (table_outputs(i,1+j)-table_outputs(i,1))/(-0.5);endend
end% Estimation of the factors influence with the average of the absolute
% values of the elementary effects
for j=1:nfacmu_temp(j) = mean(abs(table_ee(:,j)));
end[mu, order] = sort(mu_temp,'descend'); % Ordering.end

目标函数

a = [100 0 100 100 100 100 1 10 0 0 9 0 100 100 4 100 100 7 100 2];
alpha = [1 4 1 1 1 1 0.4 3 0.8 0.7 2 1.3 1 1 0.3 1 1 1.5 1 0.6];
delt = [0.2942 0.2560 0.3004 0.5150 0.7723 0.4567 0.8390 0.1369 0.1558 0.4356 0.0257 0.3248 0.0718 0.9155 0.6877 0.5548 0.5835 0.8083 0.6309 0.8071];for i=1:20y(i) = ((1+alpha(i))*abs(2*(X(i)+delt(i)-fix(X(i)+delt(i)))-1)^alpha(i)+a(i))/(1+a(i));
end

二、效果展示

三、代码获取

后台私信回复“67期”即可获取下载链接。