Multi-objective optimization

Driver:

ActiveLearning

Download script: multi_objective_optimization.m

Many optimization task do not consist of only a single objective \(f(\mathbf{x})\) to be minimized in some domain \(\mathbf{x}\in \mathcal{X}\subset \mathbb{R}^d\), but rather in a set of objectives \(\mathbf{f}(\mathbf{x}) = [f_1(\mathbf{x}), f_2(\mathbf{x}),\dots,f_k(\mathbf{x})]^{\rm T}\) to be minimized together.

A solution \(\mathbf{a} = \mathbf{f}(\mathbf{x}_1)\) is said to dominate another solution \(\mathbf{b} = \mathbf{f}(\mathbf{x}_2)\) if \(a_i \le b_i\) for all \(i \in [1, \dots, k]\) and there exists at least on entry \(j \in [1, \dots, k]\) such that \(a_j < b_j\).

A solution is called Pareto optimal if it is not dominated by any other solution. The set of Pareto optimal solutions is called the Pareto front. The aim of a multi-objective optimization is to find the Pareto front for the set of objectives \(\mathbf{f}(\mathbf{x})\).

We consider an analytic multi-objective optimization problem with \(k=2\) objectives:

\[ \begin{align}\begin{aligned}\text{Minimize } f_1(x_1,x_2) &=& x_1\\\text{Minimize } f_2(x_1,x_2) &=& \frac{g(x_2)}{x_1}\end{aligned}\end{align} \]

It can be shown, that the problem has Pareto-optimal solutions \((x_1, x_2^*)\), where \(x_2^*\) is the globally minimum solution of \(g(x_2)\) and \(x_1\) can take any value [DEB1999].

As an example, we consider the multi-modal function

\[g(x_2) = 20 + (x_2^2 - 10\cos(2\pi x_2))\]

with local minima approximately located at integer values of \(x_2\) and the global minimum at \(g(x_2^* = 0)=10\). That is, the Pareto front is defined as \(f_2 = \frac{10}{f_1}\).

[DEB1999]

Deb, K., “Multi-objective genetic algorithms: Problem difficulties and construction of test problems”. Evolutionary computation, 7(3), pp.205-230 (1999)

jcm_optimizer_path = '<JCM_OPTIMIZER_PATH>';
addpath(fullfile(jcm_optimizer_path, 'interface', 'matlab'));
 
server = jcmoptimizer.Server();
client = jcmoptimizer.Client(server.port); 

% Definition of the search domain
design_space = { ...
    struct('name', 'x1', 'type', 'continuous', 'domain', [1, 4]), ... 
    struct('name', 'x2', 'type', 'continuous', 'domain', [1, 4]), ... 
};

 % Creation of the study object with study_id 'multi_objective_optimization'
study = client.create_study( ...
    'design_space', design_space, ...
    'driver','ActiveLearning',...
    'name','Multi-objective optimization',...
    'study_id', 'multi_objective_optimization');
%Lower and upper reference point for hypervolume definition (see figure)
lower_ref=[0, 0];
upper_ref=[5, 50];
 
study.configure( ...
    'max_iter', 50, ...
    'surrogates', { ...
        ...A multi-output Gaussian process that learns the dependence of
        ...the model on the design parameters.
        struct('type', "GP", 'name', "model_vector", 'output_dim', 2, ...
             'correlate_outputs', false)
    }, ...
    'variables', { ...
        ...Selectors for the values of f1 and f2
        struct('type', "SingleSelector", 'name', "f1", ...
             'input', "model_vector", 'select_by_name', "model_vector0"), ...
        struct('type', "SingleSelector", 'name', "f2", ...
              'input', "model_vector", 'select_by_name', "model_vector1") ...
    }, ...
    'objectives', { ...
        ... Multi-minimization objective for f1, f2
        struct('type', "MultiMinimizer", 'variables', ["f1", "f2"], 'name', "objective", ...
             'lower_reference', lower_ref, 'upper_reference', upper_ref) ...
    } ...
);

% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    pause(2); % make objective expensive
    observation = study.new_observation();
    %determine list of the two objective values
    g = 20 + (sample.x2^2 - 10*cos(2*pi*sample.x2));
    observation.add([sample.x1, g/sample.x1]);
    
end  

% Run the minimization
study.set_evaluator(@evaluate);
study.run(); 


fig = figure('Position', [0, 0, 1000, 500]);

%Plot all observations of (f1, f2)
f1 = study.driver.get_observed_values('object_type', "variable", 'name', "f1");
f2 = study.driver.get_observed_values('object_type', "variable", 'name', "f2");
plot(f1.means, f2.means, "o", 'DisplayName', "Donimated observations");
hold on;
 
%Plot the estimate of the Paretro front, i.e. nondominated observations
pf = study.driver.get_state("pareto_front");
nondom_f1 = pf.f1;
nondom_f2 = pf.f2;
plot(nondom_f1, nondom_f2, "o", 'DisplayName', "Pareto front estimate");
 
%Plot analytic Pareto front
X = linspace(1, 4, 100);
plot(X, 10./X, "k--", 'DisplayName', "Analytic Pareto front");
 
%Plot nondominated yypervolume
plot(lower_ref(1), lower_ref(2), "rx", 'DisplayName', "Lower reference");
plot(upper_ref(1), upper_ref(2), "bx", 'DisplayName', "Upper reference");
X = [lower_ref(1), repelem(nondom_f1, 2), upper_ref(1)];
nondom_lower = repelem(lower_ref(2), length(X));
nondom_f2_aug = [upper_ref(2), nondom_f2];
nondom_upper = repelem(nondom_f2_aug, 2);
patch([X, fliplr(X)], [nondom_lower, fliplr(nondom_upper)], 'b', ...
             'FaceAlpha', 0.2, 'EdgeAlpha', 0.2, ...
             'DisplayName', "Nondominated hypervolume");

xlabel("f1");
ylabel("f2");
legend('Location', "north");
grid();
saveas(fig, "multi_objective_minimization.png");

Result of the multi-objective minimization in the \((f_1, f_2)\)-plane.

The blue dots show all observations, whereas the orange dots show nondominated solution, i.e. an estimate of the Pareto front. The dashed line shows the true analytic Pareto front within the domain of \((x_1, x_2)\).

The progress of the study is shown in terms of the nondominated hypervolume (blue shaded area) that is defined by the nondominating observations as well as a lower and upper reference point (red and blue cross).