Pareto set for 2D functions

7 views (last 30 days)
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!

Accepted Answer

Andreas Apostolatos
Andreas Apostolatos on 20 May 2021
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,
In this way, the following Pareto front is produced when running the latter script,
I hope this helps.
Kind regards,
Andreas
  3 Comments
laura bagnale
laura bagnale on 21 May 2021
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
Andreas Apostolatos on 21 May 2021
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,
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,
I hope that this information helps you further.
Kind regards,
Andreas

Sign in to comment.

More Answers (1)

laura bagnale
laura bagnale on 24 May 2021
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!
  3 Comments
laura bagnale
laura bagnale on 26 May 2021
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
laura bagnale on 26 May 2021
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!