ε-制約法でパレートフロントを計算する

clear; close all;
clear; close all;
global epsilon;

f1 = @(x) x.^2;
global f2;
f2 = @(x) (x-2).^2;
xopts = [];

A = [];
b = [];
Aeq = [];
beq = [];
x0 = [0];
lb = [-10,-10];
ub = [10,10];

epsilons = linspace(0,4,20);
for idx = 1:length(epsilons);
    epsilon = epsilons(idx);
    xopt = fmincon(f1, x0, A, b, Aeq, beq, lb, ub, @nonlcon);
    xopts(end+1) = xopt;
end

xopts;
f1opts = f1(xopts);
f2opts = f2(xopts);
plot(f1opts,f2opts,'o');
hold on;

x1 = linspace(0,4,20);
y2 = interp1(f2opts, f1opts, x1,'spline');
hold on;
plot(x1,y2);
big;

function [c,ceq] = nonlcon(x)
global f2 epsilon;
c = [f2(x)-epsilon];
ceq = [];
end