Pareto set for 2D functions

Hello everyone!

Can someone help me with this problem please?

I represented two functions like these:

fcontour(@(x,y) 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2;
fcontour(@(x,y) 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2;

Then I wrote an optimization problem in order to find their minimum.

Now I would like to find the Pareto set.

I was trying with the following code:

fitness = @simpleMultiObjective;
 nvar = 2;
 [x, fval] = gamultiobj(fitness, nvar);

function z = simpleMultiObjective(xy)
    x = xy(1); y = xy(2);
    z(1) = 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2;    
    z(2) = 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2;
end

But I don't know if I'm on the right way.

Can someone suggest a way to solve this problem, please?

Thank you very much!

Answer 1

[Andreas Apostolatos]

Hi Laura,

You are on the right path, but I would have some recommendations for you accordingly:

1.) When using function 'fcontour' to plot the contours of a function, MATLAB tries to provide the inputs (x,y) of the function as vectors. For this to be possible, the underlying expressions within the function definition should support vector operations. In your definitions above, you used expression 'x.^2' which supports the base to be a vector, thus raising each component of the vector 'x' to 2 while returning a vector containing the results. However, expression 'x*y' cannot be applied on two vectors of the same form (both row or both column vectors). To allow operation 'x*y' to be applied on two vectors elementwise, please consider adding a full stop '.' right before the corresponding operation, namely 'x.*y'. This will render the underlying computations more efficient. In your case, you can change your function definitions as follows,

fcontour(@(x,y) 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^.2 + (-0.02072).*x.*y + 0.5786.*y.^.2);
fcontour(@(x,y) 0.8906 - 9.46e+12.*x +9.46e12.*y + 5.20e+13.*x.^2 + 0.03807.*y.*x + 5.208e+13.*y.^2);

2.) Your objective function is unconstrained and has expressions such as 'x^.2'. The latter expression is in fact equal to 'x^0.2', which leads to imaginary results if 'x' is negative. Therefore, it may be advisable to define 0 as lower bound to your design variables such that the objective function does not attain imaginary values, of course depending on your application and your needs. I have adapted your script to account also for 0 as lower bound, namely,

%% clear output window and workspace
clc, clear;

%% plot the contour of the objective functions
fcontour(@(x, y) 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^.2 + (-0.02072).*x.*y + 0.5786.*y.^.2);
fcontour(@(x, y) 0.8906 - 9.46e+12.*x +9.46e12.*y + 5.20e+13.*x.^2 + 0.03807.*y.*x + 5.208e+13.*y.^2);

%% define optimization options and solve the Pareto front problem
% opts = optimoptions('gamultiobj')';
opts = optimoptions('gamultiobj', 'PopulationSize', 61, 'ParetoFraction', 0.71);
fitness = @simpleMultiObjective;
nvar = 2;
lb = zeros(2, 1);
[x, fval] = gamultiobj(fitness, nvar, [], [], [], [], lb, [], [], opts);

figure
plot(fval(:, 1), fval(:, 2),'o')

%% define the multiobjective function
function z = simpleMultiObjective(xy)
    x = xy(1); y = xy(2);
    z(1) = 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^.2 + (-0.02072).*x.*y + 0.5786.*y.^.2;
    z(2) = 0.8906 - 9.46e+12.*x +9.46e12.*y + 5.20e+13.*x.^2 + 0.03807.*y.*x + 5.208e+13.*y.^2;
end

Please note that I have increased the values for the 'PopulationSize' and 'ParetoFraction' optimization options when using the 'gamultiobj' function in order to obtain a denser Pareto front. For more information, please visit also the following documentation page,

https://www.mathworks.com/help/gads/pareto-front-for-two-objectives.html

In this way, the following Pareto front is produced when running the latter script,

I hope this helps.

Kind regards,

Andreas

[laura bagnale]

Dear Andreas,

Thank you very much for your prompt reply and for your support! This helps me a lot.

I will try it immediately and I will provide a feedback as soon as possible!

Best Regards,

Laura

[laura bagnale]

Dear Andreas,

The program runs well and your explanations were clear, thanks a lot!

I would have another question. In your script you found the Pareto set starting from the two initial functions and setting up the optimization problem. In my previous program it was necessary to find the respective minima of the two functions and then I wanted to find the Pareto set, by considering these minima. I have set up the entire programme in this way

fcontour(@(x,y) 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2, [0 3 0 3]);

p = optimproblem

p.Description = "minimization"

x = optimvar("x", "LowerBound", 0.5, "UpperBound",2);

y = optimvar("y", "LowerBound", 0.5, "UpperBound",2);

p.Objective = 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2

initialPt.x = 0.3;

initialPt.y = 0.5;

hold on

[sol,fval,exitflag,outout] = solve(p,initialPt)

plot(sol.x,sol.y, 'ro','LineWidth',2)

hold on

rnge = [0 3 0 3];

fcontour(@(x,y) 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2, rnge);

p = optimproblem

p.Description = "minimization"

x = optimvar("x", "LowerBound", 0.5, "UpperBound",2);

y = optimvar("y", "LowerBound", 0.5, "UpperBound",2);

p.Objective = 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2

initialPt.x = 0.3;

initialPt.y = 0.5;

hold on

[sol,fval,exitflag,outout] = solve(p,initialPt)

plot(sol.x,sol.y, 'ro','LineWidth',2)

hold off

text(0.6,0.3,'Minimum(f_1(x,y))')

text(2.2,0.4,'Minimum(f_1(x,y))')

fitness = @simpleMultiObjective;

nvar = 2;

[xy, fval] = gamultiobj(fitness, nvar);

x = xy(1); y = xy(2);

x0 = -5:0.5:5;

f1 = 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2;

f2 = 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2;

plot(fval(:, 1), fval(:,2), 'o');

xlabel('f1');

ylabel('f2');

grid;

function z = simpleMultiObjective(xy)

x = xy(1); y = xy(2);

z(1) = 0.9097 + (-0.2332*x) + 0.2645*y + (-1.07)*x^.2 + (-0.02072)*x*y + 0.5786*y^.2;

z(2) = 0.8906 -9.46e+12*x +9.46e12*y + 5.20e+13*x.^2 + 0.03807*y*x + 5.208e+13*y.^2;

end

Is it possible to adapt this situation to your script? Namely, finding the minima of the two initial functions and, after, looking for the Pareto set.

I hope that my question is clear enough.

Sorry but it has been a short time since I started dealing with the Optimization Toolbox.

Best,

Laura

[Andreas Apostolatos]

Hi Laura,

Function 'gamultiobj' utilizes the genetic algorithm for finding the Pareto front which does not employ any starting point. Therefore, the optimization options of function 'gamultiobj' do not provide an API for supplying the algorithm with a starting point.

However, function 'paretosearch' uses the pattern search algorithm for finding the Pareto front, which provides an API for supplying a starting point to the algorithm, see the following documentation pages for more information to the subject,

https://www.mathworks.com/help/gads/paretosearch.html#mw_64c877a5-5afb-48d5-9af5-68ff4dcd986c

https://www.mathworks.com/help/gads/paretosearch-algorithm.html

The corresponding optimization option is called 'InitialPoints' and can accommodate a matrix with 'nvars' columns, where each row represents one initial point while 'nvars' stands for the number of variables. You can then supply the corresponding minima of each objective function that comprise your multiobjective function to the latter optimization option of the 'paretosearch' function, and in this way, you can also skip the specification of lower bounds, since providing the minima of the individual objective functions as starting point suffices for finding the desirable Pareto front.

You can find the adaptation of my previous code snippet that utilizes function 'paretosearch' and optimization option 'InitialPoints' for your convenience below,

%% clear output window and workspace
clc, clear;

%% plot the contour of the objective functions
fcontour(@(x,y) 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^.2 + (-0.02072).*x.*y + 0.5786.*y.^.2);
fcontour(@(x,y) 0.8906 - 9.46e+12.*x +9.46e12.*y + 5.20e+13.*x.^2 + 0.03807.*y.*x + 5.208e+13.*y.^2);

%% find the minimum of objective function 1
p1 = optimproblem;
p1.Description = "minimization";
x1 = optimvar("x", "LowerBound", 0.5, "UpperBound",2);
y1 = optimvar("y", "LowerBound", 0.5, "UpperBound",2);
p1.Objective = 0.9097 + (-0.2332.*x1) + 0.2645.*y1 + (-1.07).*x1.^.2 + (-0.02072).*x1.*y1 + 0.5786.*y1.^.2;
initialPt1.x = 0.3;
initialPt1.y = 0.5;
[sol1, fval1, exitflag1, outout1] = solve(p1, initialPt1);

%% find the minimum of objective function 2
p2 = optimproblem;
p2.Description = "minimization";
x2 = optimvar("x", "LowerBound", 0.5, "UpperBound",2);
y2 = optimvar("y", "LowerBound", 0.5, "UpperBound",2);
p2.Objective = 0.8906 - 9.46e+12.*x2 +9.46e12.*y2 + 5.20e+13.*x2.^2 + 0.03807.*y2.*x2 + 5.208e+13.*y2.^2;
initialPt2.x = 0.3;
initialPt2.y = 0.5;
[sol2, fval2, exitflag2, outout2] = solve(p2, initialPt2);

%% find the Pareto front using the paretosearch algorithm
x0(1, :) = [sol1.x sol1.y];
x0(2, :) = [sol2.x sol2.y];
opts = optimoptions('paretosearch', 'InitialPoints', x0);
fitness = @simpleMultiObjective;
nvar = 2;
[x, fval] = paretosearch(fitness, nvar,[],[],[],[],[],[],[],opts);

figure
plot(fval(:,1), fval(:,2),'o')
xlabel('Objective 1')
ylabel('Objective 2')
title('Pareto front')

%% define the multiobjective function
function z = simpleMultiObjective(xy)
    x = xy(1); y = xy(2);
    z(1) = 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^.2 + (-0.02072).*x.*y + 0.5786.*y.^.2;
    z(2) = 0.8906 - 9.46e+12.*x +9.46e12.*y + 5.20e+13.*x.^2 + 0.03807.*y.*x + 5.208e+13.*y.^2;
end

There is also a short video that introduces the latter topic, see the link below,

https://www.mathworks.com/videos/pareto-sets-for-multiobjective-optimization-1535999427864.html

I hope that this information helps you further.

Kind regards,

Andreas

Answer 2

[laura bagnale]

Hi Andreas,

it is perfect for me, thank you very much for your support!

Your answers were clear and very helpful!

I really want to thank you a lot!

Kind Regards,

Laura!

[Andreas Apostolatos]

Hi Laura,

I am pleased that I could help!

I wish you a lot of success ahead.

Kind regards,

Andreas

[laura bagnale]

Hi Andreas,

I modified the last code a bit with one different equation (the second one). But an error occurs:

fcontour(@(x,y) 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^2 + (-0.02072).*x.*y + 0.5786.*y.^2);

fcontour(@(x,y) 0.9076 + 0.007993.*x + 0.02571.*y + 0.0382.*x.^2 + 0.04141.*y.*x -0.6286.*y.^2);

%% find the minimum of objective function 1

p1 = optimproblem;

p1.Description = "minimization";

x1 = optimvar("x", "LowerBound", 0.5, "UpperBound",2);

y1 = optimvar("y", "LowerBound", 0.5, "UpperBound",2);

p1.Objective = 0.9097 + (-0.2332.*x1) + 0.2645.*y1 + (-1.07).*x1.^.2 + (-0.02072).*x1.*y1 + 0.5786.*y1.^.2;

initialPt1.x = 0.3;

initialPt1.y = 0.3;

[sol1, fval1, exitflag1, outout1] = solve(p1, initialPt1);

%% find the minimum of objective function 2

p2 = optimproblem;

p2.Description = "minimization";

x2 = optimvar("x", "LowerBound", 0.5, "UpperBound",2);

y2 = optimvar("y", "LowerBound", 0.5, "UpperBound",2);

p2.Objective = 0.9076 + 0.007993.*x2 + 0.02571.*y2 + 0.0382.*x2.^.2 + 0.04141.*y2.*x2 -0.6286.*y2.^2.;

initialPt2.x = 0.3;

initialPt2.y = 0.3;

[sol2, fval2, exitflag2, outout2] = solve(p2, initialPt2);

%% find the Pareto front using the paretosearch algorithm

x0(1, :) = [sol1.x sol1.y];

x0(2, :) = [sol2.x sol2.y];

opts = optimoptions('paretosearch', 'InitialPoints', x0);

fitness = @simpleMultiObjective;

nvar = 2;

[x, fval] = paretosearch(fitness, nvar,[],[],[],[],[],[],[],opts);

figure

plot(fval(:,2), fval(:,1),'o')

xlabel('Objective 1')

ylabel('Objective 2')

title('Pareto front')

%% define the multiobjective function

function z = simpleMultiObjective(xy)

x = xy(1); y = xy(2);

z(1) = 0.9097 + (-0.2332.*x) + 0.2645.*y + (-1.07).*x.^2 + (-0.02072).*x.*y + 0.5786.*y.^2;

z(2) = 0.9076 + 0.007993.*x + 0.02571.*y + 0.0382.*x.^2 + 0.04141.*y.*x -0.6286.*y.^2;

end

In plot(fval(:,2), fval(:,1),'o') it apperas: Index in position 2 exceeds array bounds.

Do you thin that you could help me once again, Please?

Thanks a lot!

Laura

[laura bagnale]

Dear Andreas,

sorry if I am disturbing you. It seems that I solved the problem by replacing z(1) and z(2) with corrected form in the syntax. In particular I replaced x.^2 with x.^.2 and for y.^2 too. But I don't understand what does it mean exactly.

I try to looking for it in the matlab documentation but if you could clarify this, you would be doing me a great favour.

If you can't, no problem I will find it! :)

Thank you very much!

Laura