満足化トレードオフ法（STOM）によるパレートフロントの計算

その１

clear; close all;

aLevels = [2 1;
           1 3;
           0.4 0.5];

f1 = @(x) x.^2;
f2 = @(x) (x-2).^2;
x=-10:0.01:10;
f1i = min(f1(x));
f2i = min(f2(x));
c1 = 2;
c2 = 2;
f1a = f1i + c1;
f2a = f2i + c2;

xopts = [];
w1s = [];
w2s = [];

for idx=1:length(aLevels)
aLevel = aLevels(idx,:);
f1a = aLevel(1);
f2a = aLevel(2);

w1 = 1/(f1a-f1i);
w2 = 1/(f2a-f2i);
w1s(end+1) = w1;
w2s(end+1) = w2;

C = max(w1*(f1(x)-f1a),w2*(f2(x)-f2a));
xopt = x(find(C==min(C)));
xopts(end+1) = xopt;
f1opt = f1(xopt);
f2opt = f2(xopt);

end

xopts = unique(sort(xopts));
f1opts = f1(xopts);
f2opts = f2(xopts);
ws = [w1s;w2s]';

plot(f1opts,f2opts,'o','DisplayName','Pareto Optimal Value');
hold on;

x1=0:0.1:2;
y1=interp1(f2opts,f1opts,x1,'spline');
plot(x1,y1,'-','DisplayName','Pareto front');

plot(aLevels(:,1),aLevels(:,2),'x','DisplayName','Aspiration Level');
legend();

plot(f1i,f2i,'xk','DisplayName','Ideal point');

big;