How to solve an optimization problem with the objective function is a matrix ?

Hello everyone!

I have no idea how to solve optimization problems.

I have an objective function which is a binary matrix, should this matrix be defined (known) or not?

knowing that I looked for examples but I didn't find the relevant one with my problem.

I would be grateful if you could help me!

[Torsten]

The objective function of a usual optimization problem returns a scalar value that is attempted to be made as small or as large as possible.

So what is this scalar value for your binary matrix ?

[Majid]

sorry! i didn't understand your question. could you please explain more!

Answer 1

[Walter Roberson]

Each element of the matrix is effectively a separate objective, so in order to do this you would need to use a multi-objective optimizer.

You can use Problem Based Optimization https://www.mathworks.com/help/optim/problem-based-approach.html . You would probably create optimization variables that are marked as integers with lower bound 0 and upper bound 1. You would set up whatever constraints are appropriate. You would set up a matrix objective . solve() should then notice that gamultiobj() is needed to minimize the problem.

Note that there is a difference between a problem in which the best location is a binary matrix, compared to a problem in which the objective is a binary matrix. If the inputs are all binary and a matrix, but the output were scalar (for example an economic dispatch cost) then there are a number of different minimizers, but there are only a quite small number that can deal with there being multiple objectives (each element of the output matrix being a different objective.)

[Majid]

@Walter Roberson The problem is that I don't know if the objective function should be already defined or not.

for example i have this objective function :

max M =

in this case u should be known?

[Walter Roberson]

That formula is M = sum(u(1:K, 1:L));

It is not clear exactly what is being varied. Is u(i,j) according to some function of some additional variable(s) and the task is to find the combination of inputs so that the sum is greatest? Or is it the K and L are variable, and the task is to find the sub-rectangle that has the highest sum (imagine for example that u(i,j) contain a mix of positive and negative values and you want to know what the highest score is) ?

M appears to be a scalar.

Is u(i,j) binary? If so then M = nnz(u(1:K, 1:L)) and if K and L span all of u then M = nnz(u)

[Majid]

yes u(i,j) are elements of binary matrix U(K*L).

M is the sum of u(i,j), i will try to find the highest sum but i don't know how to get U which is the optimized result.

it is possible through constraints getting the optimal solution?

[Torsten]

Yes, if there were no constraints, then u(i,j) = 1 for all i and j would be the obvious solution.

[Walter Roberson]

Your objective function would be along the lines of

function M = objective(x)
   calculate u based on x
   M = -sum(M, 'all');
end

You want the negative because you are maximizing, and you can maximize by finding the minimum of the negative of the value.

You do not need any constraints on M or u. You might need constraints on x .

Furthermore, the constraints you put on x might involve computing all or most of u first, and testing whether the result has certain properties (for example, diagonal must be 0, or number of entries in the 4th row must be 2). You can avoid re-calculating u by using a function to calculate it, and memoize and call the function from your objective and also from your nonlinear constraint function.

[Majid]

i have constraint which is related to matrix U, how can i put it and U is not defined.

i'm really wondering how to relate objective function with constraint? should they be linear right?

[Torsten]

Every system of equations has unknowns that are not defined (not known), but follow some equations.

The equations connecting the unknowns lead to values for the unknowns that make the equations true.

This is the same for the elements of your matrix U.

If the constraints (equations or inequalities) are linear in the unknowns, you have a linear optimization problem, otherwise a nonlinear.

[Majid]

@Torsten i didn't get how to know if the constraints linear to the problem or not?

for example i have this constraint :

w-(d(i,j)*u(i,j))>= 0

w is constant

d(i,j) is element of matrix D with same dimention of U.

in this case U is not defined.

[Walter Roberson]

You would use the trial inputs, x, to create a trial u . You would use that trial u in your nonlinear constraint function to signal whether the trial x lead to a u such that the constraints were satisfied. If you use memoize() like I indicated then it will not be necessary to recalculate u (but if you are using genetic algorithm then calculation of the constraints might be postponed relative to the function evaluation, so you might need to have a cache size a bit larger than your population size.)

[Majid]

@Walter Roberson the trial inputs you mean that I choose arbitrarily U?

[Walter Roberson]

In any problem you have a set of model parameters that need to be determined in order to minimize (or maximize) the objective. You calculate various things from those model parameters. You might then have checks to see if the intermediate values are within range -- constraints on the intermediate values. Then if the constraints pass, you continue to calculate a final "cost" for that set of model parameters. The optimization routine then tries other sets of model parameters, getting an associated "cost", and others yet... and at the end returns the set of model parameters that gives the best overall "cost".

For example suppose you are need to fit an ellipse, and the model includes the locations of the two foci and the length of the major and minor axes -- four parameters for that part. But you have a constraint that (for example) the two foci have to be within a particular distance of a known curved line. So you would use the trial foci locations passed in in order to calculate distances from the curved line, and your constraint function would test whether those distances are acceptable, and your output function would calculate whatever (such as comparing the area of the convex boundary of the input points to the area of the ellipse -- if you were to mazimize that fraction you would be minimizing the area of the ellipse.)

The constraints do not need to be with respect to anything to do with the objective function; they can be totally arbitrary. Items positioned relative to each other by law rather than by calculated safey index for example.

It happens that sometimes the model parameters are binary inputs. For example the model parameters might be your u matrix, and they might be acting as selectors as to which actions to take or not, with there being tradeoffs that make the choices dependent on each other. "Select exactly 3 of these 7" for example. If you have that kind of situation, then create integer constraints on the inputs, minimum 0, maximum 1, and use the linear inequality and linear equality constraints where feasible. Some systems like this are best handled by intlinprog().

[Majid]

Yes the last model you mentioned is the appropriate to my problem..so i will define intlinprog to find the optimal solution? If it is the case how can i formulate the constraint i mentioned above to linear inequality? Apologize for asking many questions ..i appreciate your patience and time to help me ..

[Majid]

@Walter Roberson @Torsten can i share with you my problem in private , maybe you could more explain it!

[Torsten]

If your matrix u of unknowns is nxm, you would formulate the constraint w-(d(i,j)*u(i,j))>= 0 as

ub = w./D(:)

if all the elements in the matrix D are positive and the objective function as

f = -ones(n*m,1)

[Majid]

@Torsten yes all elements in the matrix D are positive.

so if i use the objective function as @Walter Roberson indicated above , after this i write the line command you mentioned or there is a specific structure to write objective function ? i want to display 'optimal solution found'

[Torsten]

% U is a 12x24 matrix
n = 12;
m = 24;
% Your w-value
w = 12;
% Your distance matrix D, randomly created
D = rand(n,m);
% Objective function is max: sum_i sum_j u_ij
f = -ones(n*m,1);
% u_ij >= 0
lb = zeros(n*m,1);
% u_ij <= w/D
ub = w./D(:);
% u_ij integer
intcon = 1:n*m;
sol = intlinprog(f,intcon,[],[],[],[],lb,ub);

LP:                Optimal objective value is -33234.000000.                                        


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer
within tolerance, options.IntegerTolerance = 1e-05 (the default value).

sol = reshape(sol,n,m)

sol = 12×24

26 105 37 82 367 21 13 9458 16 27 19 61 12 41 22 33 20 20 31 38 113 18 14 30 22 112 15 22 36 270 100 13 26 20 163 13 113 43 33 44 102 52 29 90 12 14 155 264 14 16 13 68 44 135 64 13 60 13 22 346 17 14 29 21 23 14 92 17 17 58 17 34 611 21 13 15 17 36 156 283 77 57 12 25 27 28 15 49 22 27 13 14 35 26 35 16 21 127 40 13 17 17 16 88 12 12 42 17 15 36 25 115 28 67 13 59 408 24 35 18 12 45 12 25 38 21 15 12 17 41 16 20 30 12 13 719 21 12 13 108 16 38 51 34 304 491 256 2767 13 12 60 14 46 116 14 23 101 19 27 15 16 33 16 28 13 28 12 20 17 32 13 14 68 32 68 56 14 20 27 122 24 19 39 14 391 61 900 19 593 36 62 44 36 16 25 29 31 14 91 20 26 82 17 14 33 41 80 14 15 49 27 34 13 12 13 186 1483 18 20 16 23 27 23 42 18 61 138 22 13 38 16 35 23 76 20 128 16 13 107 21

%Check whether the solution satisfies the constraints
any(sol > w./D,'All')

ans = logical
   0

any(sol < 0,'All')

ans = logical
   0

[Majid]

@Torsten what about this error "Intlinprog stopped because the root LP problem is unbounded."

[Torsten]

It means that the value of your objective function can be made as large as you want.

This error usually appears if the constraints don't limit the maximum value of the objective function.

[Majid]

@Torsten so what should i do in this case?knowing that i can't use a large matrix

[Torsten]

What does the error have to do with a large matrix ?

You should try to explain your problem in plain words (not mathematically). Maybe the formulation as an optimization problem becomes clearer this way.

[Majid]

@Torsten okay i will try to explain more :

i have an initial binary matrix F(m*n) , i will choose some of elements of F and generate new matrix U(m*n) respecting the indicated constraint.At the end i will find the sum of matrix U noted (M) and i will try to maximize M.

For the matrix F is not mentioned in my problem because it already existed in the constraint. Elements of D related to each "1" existed in matrix F. otherwise 0.

[Torsten]

Sorry, but I can't make sense of your description.

[Torsten]

I don't see X appearing anywhere in your problem formulation (your objective function does not depend on X explicitly).

And any special properties of D and E^pi ?

[Majid]

D and E_pi are matrix with positive values.

I want to optimize U if X is known. the two constraints depends to X. it means i used X to get them.

[Torsten]

But X (and thus the constraints) are fixed and don't change with U ?

Then your problem doesn't need an optimizer:

U = floor(min(w./D(:),(E_total-E_zi)./E_pi(:),2))

[Majid]

yes, but U is undefined, U should be a binary matrix.

[Torsten]

ok, then

W = min(w./D(:),(E_total-E_zi)./E_pi(:),2);
U = zeros(size(W));
U(W>=1) = 1.0;

[Majid]

@Torsten why the problem doesn't need an optimizer ?

from the code you mentioned, it seems like i select the elements of matrix U, but in real i want to find it through constraints. i don't know if you understand me or not?

[Torsten]

Let U be unknown.

The constraints say

U <= w./D(:)

and

U <= (E_total-E_zi)./E_pi(:)

Thus

U <= min(w./D(:),(E_total-E_zi)./E_pi(:),2)

Now we come to solve for U.

You want to maximize sum U where U is binary, thus you must take

U = 0 at the positions where min(w./D(:),(E_total-E_zi)./E_pi(:),2) < 1 

and you can take

U = 1 at the positions where min(w./D(:),(E_total-E_zi)./E_pi(:),2) >=1.

And that's the solution for U that takes into account your constraints.

[Majid]

@Torsten i'm sorry, it seems like i failed to explain what i need.

the matrix U should be generated first, i already tell you that i have intial matrix F (binary matrix also) from it i will selected elements satisfied the constraints.

[Torsten]

If U is generated at the beginning, what do you want to optimize ?

What is the relationship between U and F ?

I don't understand

i have an initial binary matrix F(m*n) , i will choose some of elements of F and generate new matrix U(m*n)

If F and U are both m*n, then F and U must be identical if the elements of U are chosen from F, don't they ?

[Majid]

F is binary matrix , for each "1" i calculated D et E_pi.

Now i want to generate matrix U with elements "1" and "0" , each "1" refers to constraint satisfied.

all this with optimization.

If F and U are both m*n, then F and U must be identical if the elements of U are chosen from F, don't they ? yes in the ideal case

[Torsten]

But that's exactly what the code does.

U is a copy of the matrix F where the elements that are 1 in F and for which

min(w./D(:),(E_total-E_zi)./E_pi(:),2) < 1

are replaced by 0.

[Majid]

@Torsten yes i know this, but how can i do this with optimization? for example can i use CPLEX optimizer in matlab?

[Torsten]

First state the optimization problem properly. Then in the next step we can look for a suitable optimizer.

It might be the case that I don't understand the problem correctly. But up to now, I don't see the necessity to use an optimizer.

[Majid]

@Torsten i really appreciate your time and patience.

the optimization is needed because i can't select elements through equations.

I have to linearize the problem and the constraints so that I can solve the problem.

[Torsten]

Then please write down the complete optimization problem and I will give you the code.

Why did you delete the problem formulation you already included ? Was it wrong ? Was it incomplete ?

[Torsten]

As said, for this problem you don't need optimization.

But if you insist, here is an implementation:

% U is a 12x24 matrix
n = 12;
m = 24;
% Your constant values
w = 12;
E_zi = 4;
E_total = 4.1;
% Your distance matrix D, randomly created
D = rand(n,m);
% Your matrix E_pi, randomly created
E_pi = rand(n,m);
% Objective function is max: sum_i sum_j u_ij
f = -ones(n*m,1);
% u_ij <= min(w/D,E_total-E_zi)/E_pi)
Aineq = eye(n*m);
bineq = min(w./D(:),(E_total-E_zi)./E_pi(:));
% u_ij >= 0
lb = zeros(n*m,1);
% u_ij <= 1
ub = ones(n*m,1);
% u_ij integer
intcon = 1:n*m;
sol = intlinprog(f,intcon,Aineq,bineq,[],[],lb,ub);

LP:                Optimal objective value is -24.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer
within tolerance, options.IntegerTolerance = 1e-05 (the default value).

sol = reshape(sol,n,m)

sol = 12×24

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0

any(sol < 0,'All')

ans = logical
   0

any(sol > 1,'All')

ans = logical
   0

any(sol > w./D,'All')

ans = logical
   0

any(sol > (E_total-E_zi)./E_pi,'All')

ans = logical
   0

[Majid]

@Torsten after doing sol = reshape(sol,n,m) , i don't get the needed matrix, i want to get "1" only when the constraints are satisfied.

LP:                Optimal objective value is -9894.000000.                                         


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).


sol =

  Columns 1 through 20

     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1

[Torsten]

You get U(i,j) = 1 if w./D(i,j) is >=1 and (E_total-E_zi)./E_pi(i,j) >=1

You get U(i,j) = 0 if w./D(i,j) is >=0 and (E_total-E_zi)./E_pi(i,j) >=0, but one of them is < 1.

You get an infeasible problem if one of the constraints is < 0.

[Majid]

@Torsten but in the matrix D and E_pi not all elements are values, there are zeros also, i won't consider them.

[Torsten]

Now this is exactly the optimization problem you posted.

If something has to be changed, you must change it in the formulation of your problem.

rng('Default')
% U is a 12x24 matrix
n = 12;
m = 24;
% Your constant values
w = 12;
E_zi = 4;
E_total = 4.1;
% Your distance matrix D, randomly created
D = rand(n,m);
% Your matrix E_pi, randomly created
E_pi = rand(n,m);
% Objective function is max: sum_i sum_j u_ij
f = -ones(n*m,1);
% D_ij * u_ij <= w
% E_pi_ij * u_ij <= E_total - E_zi
Aineq = zeros(2*n*m,n*m);
bineq = zeros(2*n*m,1);
icount = 0;
for i = 1:n
  for j = 1:m
    icount = icount + 1;
    Aineq(icount,(i-1)*m+j) = D(i,j);
    bineq(icount) = w;
    icount = icount + 1;
    Aineq(icount,(i-1)*m+j) = E_pi(i,j);
    bineq(icount) = E_total - E_zi;
  end
end
% u_ij >= 0
lb = zeros(n*m,1);
% u_ij <= 1
ub = ones(n*m,1);
% u_ij integer
intcon = 1:n*m;
sol = intlinprog(f,intcon,Aineq,bineq,[],[],lb,ub);

LP:                Optimal objective value is -24.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer
within tolerance, options.IntegerTolerance = 1e-05 (the default value).

sol = reshape(sol,m,n).'

sol = 12×24

0 1 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

any(sol < 0,'All')

ans = logical
   0

any(sol > 1,'All')

ans = logical
   0

any(sol.*D > w,'All')

ans = logical
   0

any(sol.*E_pi > E_total-E_zi,'All')

ans = logical
   0

[Majid]

@Torsten it was the same result, it seems that zeros are considered too.

Optimal objective value is -9887.000000. and i have 9960 elements in my matrix.

[Torsten]

ok, then U(i,j) = 1 for 9887 positions and U(i,j) = 0 for 73 positions.

[Majid]

@Torsten yes but i need a condition to don't consider zeros in both matrix.

[Torsten]

You can fix U(i,j) as 0 if D(i,j) = 0 and/or E_pi(i,j) = 0. I don't know what makes sense for your task.

[Majid]

@Torsten it means after getting optimal solution , i will add this condition?

[Torsten]

U(i,j) = 0 is also a constraint on the U-matrix that can be set for the optimizer.

So you can set it before calling "intlinprog" within the problem formulation.

But all these optimizations are not necessary - I don't understand why you prefer such a difficult solution for such a simple problem.

[Majid]

@Torsten it is not simple but i can't explain it clearly. when i add condition before calling 'intlinprog" nothing was change. i should add constraints for both matrix D and E_pi to consider only values > 0.

how can i do this ?

[Torsten]

You mean U(i,j) should only be constrained by D resp. E_pi if D(i,j) resp.E_pi(i,j) > 0 ?

[Majid]

@Torsten exactly

[Torsten]

Then replace

for i = 1:n
  for j = 1:m
    icount = icount + 1;
    Aineq(icount,(i-1)*m+j) = D(i,j);
    bineq(icount) = w;
    icount = icount + 1;
    Aineq(icount,(i-1)*m+j) = E_pi(i,j);
    bineq(icount) = E_total - E_zi;
  end
end

by

for i = 1:n
  for j = 1:m
    if D(i,j) > 0
      icount = icount + 1;
      Aineq(icount,(i-1)*m+j) = D(i,j);
      bineq(icount) = w;
    end
    if E_pi(i,j) > 0
      icount = icount + 1;
      Aineq(icount,(i-1)*m+j) = E_pi(i,j);
      bineq(icount) = E_total - E_zi;
    end
  end
end
Aineq = Aineq(1:icount,:);
bineq = bineq(1:icount);

[Majid]

@Torsten nothing was change, i didn't understand where is the problem?