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

news2024/11/24 11:59:19

【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:nfac
    points(1:npoint,i) = coord(randperm(length(coord)));
end
points = points/(npoint-1)*(maxi-mini)+mini;


for i=1:npoint
    table_outputs(i,1) = studied_function(points(i,:)); % Output at the sampled point.
    for j=1:nfac
        if points(i,j) < 0.5 % If the coordinate is smaller than 0.5, a positive variation is applied
            table_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 applied
            table_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);
        end
    end
end

% Estimation of the factors influence with the average of the absolute
% values of the elementary effects
for j=1:nfac
    mu_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:20
    y(i) = ((1+alpha(i))*abs(2*(X(i)+delt(i)-fix(X(i)+delt(i)))-1)^alpha(i)+a(i))/(1+a(i));
end