基準点法によるパレートフロントの計算

clear; close all;

% aspiration levels f1a,f2a
aLevels = [1 3;
           0.5 0.5;
           0.1 0.8
           2 1];

f1 = @(x) x.^2;
f2 = @(x) (x-2).^2;

c = 1;
x = -10:0.1:10;
xopts = [];
for idx=1:length(aLevels)
    aLevel = aLevels(idx,:);
    f1a = aLevel(1);
    f2a = aLevel(2);
    C = (f1(x)-f1a).^2+(f2(x)-f2a).^2 +...
      c*(max(0,f1(x)-f1a).^2 +...
         max(0,f2(x)-f2a).^2);     
    xopt = x(find(C==min(C),1));
    xopts(end+1) = xopt;
end

xopts = unique(sort(xopts));
f1opts = f1(xopts);
f2opts = f2(xopts);
plot(f1opts,f2opts,'o','DisplayName','Pareto Optimal Value');
hold on;

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

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

big;