Asking for Guidance to Weighted-Sum Multiobjective Optimization

Hi,

Frankly speaking I am quite new in multiobjective and also using MATLAB. I have two objective functions which are quadratic.

One I define f1(x) as and the second one f2(x) with the same constraints let say A1x = b1, A2x >= b2 .

Somehow I want to create them as a weighted sum w1 * f1(x) + w2 * f2(x).

But I am not sure whether the code that I have written is correct or not. In the end, I want to find the pareto of the objective function. Could you suggest me what's not correct with the code below?

clc, clear, close all;

% c = [36., 44., 50.50, 64.5, 44.50, 37.];
% c = [51., 42.50, 62., 76.50, 52.50];
% c = [62., 97.50, 95.50, 59.];
% c = [117.5, 87.50, 119.];
c = [187., 165.5];
s = [30, 36, 41, 44, 39, 31];

A1 = kron(eye(length(s)), ones(1,length(c))); % Equality constraints
A2 = kron(ones(1,length(s)), eye(length(c))); % Inequality constraints
b1 = ones(1,length(s))';
b2 = ones(1,length(c))';

%% Set parameter for objective function 1
alpha = 0.7;
q = kron(s, eye(length(c)));
% Positive - negative definite
% Q = 2 * (q' * q);
% Negative Definite
Q = 2 * (q' * q - (max(eig(q' * q)) + alpha) * eye(length(q' * q)));
% Positive Definite
% Q = 2 * (q' * q + (min(eig(q' * q)) + alpha) * eye(length(q' * q)));
% p and r value
p = -2 * c * q;
r = c * c';

%% Objective function 2
m_values = length(c);
% Initialize pairs
pairs = [];

% Generate pairs
for i = 1:(m_values - 1)
    for j = (i + 1):m_values
        pairs = [pairs; [i, j]];
    end
end

% Display pairs
fprintf('For m = %d, pairs are: %s\n', m_values, mat2str(pairs));

% Initialize matrices_B and matrices_P
matrices_B = cell(1, length(pairs));
matrices_P = cell(1, length(pairs));

% Populate matrices_B and matrices_P
for k = 1:pairs
    i = pairs(k, 1);
    j = pairs(k, 2);

    % Initialize matrix_B
    matrix_B = zeros(m_values, m_values);
    matrix_B(i, i) = 1;
    matrix_B(j, j) = 1;
    matrices_B{k} = matrix_B;

    % Initialize matrix_P using Kronecker product
    matrix_P = kron(s' * s, matrix_B);
    matrices_P{k} = matrix_P;
end

% Sum all P_i matrices
sum_matrices_P = sum(cat(3, matrices_P{:}), 3);

% Calculate the final P matrix
% Set P to be negative definite
P = 2 * (sum_matrices_P - (max(eig(sum_matrices_P)) + alpha) * eye(size(sum_matrices_P)));
% P = 2 * (sum_matrices_P );

%% Optimize
options = optimoptions('gamultiobj','PlotFcn','gaplotpareto','UseParallel',true,'UseVectorized',true,'HybridFcn','fgoalattain');
% options = optimoptions('paretosearch','PlotFcn',{'psplotparetof' 'psplotparetox'});
x = optimvar("x",length(c)*length(s),'Type','continuous','LowerBound',0,'UpperBound',1);
w1 = 0.7;
w2 = 1 - w1;
w = [w1, w2]; % Weight vector
prob = optimproblem;
obj1 = 0.5 * x' * Q * x + p * x + r;
obj2 = 0.5 * x' * P * x;

% Combine objectives with the weighted sum
% prob.Objective.obj = w(1) * obj1 + w(2) * obj2;

% Represent the combined objective as two separate objectives
prob.Objective.obj1 =  w(1) * obj1;
prob.Objective.obj2 =  w(2) * obj2;

cons1 = A1 * x == b1;
cons2 = A2 * x >= b2;
% cons3 = obj2 <= 0.3;

prob.Constraints.cons1 = cons1;
prob.Constraints.cons2 = cons2;
% prob.Constraints.cons3 = cons3;

show(prob)
[sol,fval] = solve(prob, Solver="paretosearch");

% Sort the Pareto front solutions
[sorted_obj1, sort_index] = sort(sol.obj1);
sorted_obj2 = sol.obj2(sort_index);

h = paretoplot(sol);
% grid on
% axis auto
h.Marker = "hexagram";
h.MarkerEdgeColor = "b";

% figure;
% scatter(sol.obj1, sol.obj2, Marker="hexagram", MarkerFaceColor="flat", MarkerEdgeColor="flat")
% hold on
% xlabel("obj1");
% ylabel("obj2");
% title("Pareto Front");
% grid on

Thank you!

[Matt J]

What reason is there to think that anything is wrong?

[Muhamad]

Hi,

The reason that I have so far it is because I could not convince with the code that I wrote.

I got the pareto result from paretosearch like a linear line. Also, do I write the objective correctly if I aim to have a weighted sum objective function.

Sorry, if my explanation is hard to understand

Answer 1

[Matt J]

Because your problems have such a simple quadratic form, I might just do a sweep over w1,w2 using quadprog.

w1=linspace(0,1,60);
w2=1-w1;

[Aeq,beq]=deal(A1,b1);
[A,b]=deal(-A2,-b2);

N=width(Aeq);
M=numel(w1);

lb=zeros(N,1);
ub=lb+1;

X=nan(N,M);

for i=1:M
    
   [X(:,i),fval]= quadprog(Q, w1(i)*p, A,b,Aeq,beq,lb,ub);
    
end

[Muhamad]

Hi Matt,

Thank you for your help.

For the code above, it's the way to deal with one Q right?

If I want to include different parameter such P in my code, should I write like this?

[X(:,i),fval]= quadprog(w2(i)*P, [], A,b,Aeq,beq,lb,ub);

What should I modify the code to w1*f1 + w2 * f2 and obtain the pareto of it?

[Matt J]

What I wrote for you is the minimization of w1*f1 + w2 * f2 for a range of w1,w2 between 0 and 1.

Because f1 and f2 share a common Q, the weighting does not affect the second order term of the objective. That is likely why you are seeing a straight-line Pareto front.

[Muhamad]

Although the way I generate Q and P is different, we can say the minimization can be solved using quadprog as in your code, right?

[Matt J]

No, if your f1 and f2 have no linear or Hessian terms in common, the call to quadprog will have the more general form,

   [X(:,i),fval]= quadprog(w1(i)*Q1+w2(i)*Q2, ...
                           w1(i)*p1+w2(i)*p2,...
                           A,b,Aeq,beq,lb,ub);

[Muhamad]

Hi Matt,

Thanks for the help. I think it makes sense now. Maybe this is one last question.

Since the solution that we have from X(:,1) will be 60 based on this line

w1=linspace(0,1,60);

To capture the pareto of it in a plot or using paretosearch, how to do that? do paretoplot is enough or we need to plot like

scatter(X(:,1), X(:,2)) to represent the weighted sum?

[Matt J]

You could do scatter(X(1,:), X(2,:)) or plot(X(:,1), X(:,2)).

[Muhamad]

Thanks a lot!