Optimization: Capacitated Facility Location Problem

Am tasked to solve the capacitated facility location problem, i.e.,

subject to

I've looked at intlinprog, but I don't understand how to use this function (if it's possible) to solve this type of problem due to the summations.

If it's possible, how would I go about?

Answer 1

[Matt J]

x=optimvar('x',[n,m],'Upper',0);
y=optimvar('y',n,'type','integer','Lower',0,'Upper',1);

objective=sum(f(:).*y(:))+sum(c(:).*x(:));

constraints.colsum=sum(x,1)==d(:).';
constraints.rowsum=sum(x,2)<=s(:).*y(:);

solution=solve( optimproblem('Objective', objective,'Constraints', constraints ) );

[John D'Errico]

Good. I do too often forget about using this formulation for optimizations.

[Sebastian Daneli]

@John D'Errico, @Matt J. Unfortunately, this way of doing it didn’t find a solution. I’m actually using GLPKMEX (https://se.mathworks.com/matlabcentral/fileexchange/75318-glpkmex-gnu-linear-programming-kit-glpk-mex-generator), which uses GLPK, and can be used as a substitute for intlinprog. I am betting on that if I can formulate the problem in such a way that it can be used with intlinprog, GLPKMEX would do the job. To bad that the demo-version of AMPL does not allow me to solve my problems due to their size.

[Matt J]

You can obtain the matrix forms used by intlinprog by doing,

p = prob2struct( optimproblem('Objective', objective,'Constraints', constraints ) )

[Torsten]

According to the OP's formula, the objective

objective=sum(f(:).*y(:))+sum(c(:).*x(:));

is not correct. The inner summation starts at i, not at 1. But I doubt that this is wanted.

[Sebastian Daneli]

@Matt J, @Torsten, I will look into this. Thank you for your help.

[Matt J]

To fix the issue mentioned by Torsten, if indeed we only wish to sum over the upper triangle of x, then this is the required modification.

T=triu(c);

objective=sum(f(:).*y(:))+sum(T(:).*x(:));

[Sebastian Daneli]

My code now look like this:

clear

m=3;                                                % number of locations
n=5;                                                % number of customers

s=[239 225 184];                                    % capacity
d=[92 82 83 69 74];                                 % demand
f=[589 766 886];                                    % fixed cost
c=[14 5 6 24 6; 9 22 26 5 21; 16 11 23 28 24];      % variable cost

x=optimvar('x',[m,n],'Upper',0);
y=optimvar('y',m,'type','integer','Lower',0,'Upper',1);

T=triu(c);
objective=sum(f(:).*y(:))+sum(T(:).*x(:));

constraints.colsum=sum(x,1)==d(:).';
constraints.rowsum=sum(x,2)<=s(:).*y(:);

p=prob2struct(optimproblem('Objective', objective,'Constraints',constraints))

p = struct with fields:
     intcon: [16 17 18]
         lb: [18×1 double]
         ub: [18×1 double]
         x0: []
     solver: 'intlinprog'
      Aineq: [3×18 double]
      bineq: [3×1 double]
        Aeq: [5×18 double]
        beq: [5×1 double]
         f0: 0
          f: [18×1 double]
    options: [1×1 optim.options.Intlinprog]

solution=solve(optimproblem('Objective',objective,'Constraints',constraints));

Solving problem using intlinprog.

No feasible solution found.

Intlinprog stopped because no point satisfies the constraints.

[x,fval,exitflat,output]=intlinprog(p.f,p.intcon,p.Aeq,p.beq,p.Aineq,p.bineq,p.lb,p.ub,p.x0,p.options)

LP:                Optimal objective value is 0.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).

x = 18×1

0 0 0 0 0 0 0 0 0 0

fval = 0

exitflat = 1

output = struct with fields:
        relativegap: 0
        absolutegap: 0
      numfeaspoints: 1
           numnodes: 0
    constrviolation: 0
            message: '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).'

But I'm not getting a reasonable solution.

[Matt J]

The problem is not feasible because the d(j) values that you have supplied are all positive, while xij are constrained to be negative. That makes it impossible to satisfy .

Your second call to intlinprog doesn't show this because you have mixed up the order of some of the inputs:

[x,fval,exitflat,output]=intlinprog(p.f,p.intcon,p.Aineq,p.bineq,p.Aeq,p.beq,p.lb,p.ub)

No feasible solution found.

Intlinprog stopped because no point satisfies the constraints.


x =

     []


fval =

     []

exitflat = -2

output = struct with fields:
        relativegap: []
        absolutegap: []
      numfeaspoints: 0
           numnodes: 0
    constrviolation: []
            message: 'No feasible solution found.↵↵Intlinprog stopped because no point satisfies the constraints.'

Answer 2

[John D'Errico]

First, what is the difference between the summation over i in a sum like this:

sum(f(i)*yi))

and the dot product of two vectors:

f' * y

where f and y are column vectors of the same lengths?

Answer: Nothing, as long as the vectors are real numbers. If they could be complex, then I might need to worry about conjugate tranpose, versus transpose, but that is irrlevant here.

Is that not the first part of your objective? Next, look at the second term in the objective. x is apparently an array, of size n by m. Can you represent that double sum as a linear combination of the elements of x? (Yes.) You will need to use kron to do it.

Next, both x AND y are unknowns in the problem. So you will need to treat x and y as part of the same vector. essentially, you need to pack it all into one long vector.

I won't get more into those questions, since your problem as stated is meaningless. You have a fundamental flaw in the equatinos you wrote as equaity constraints.

You show the sum over i, of the matrix x(i,j). Once you sum over i, i goes away. The result cannot be another matrix d(i,j). Until you resolve that, this problem has no solution.

[Sebastian Daneli]

Thanks for the help. Yes, I made some mistakes when defining the constraints.

[John D'Errico]

Better now. Suppose you have the array x(i,j). Can you convert that set of equality constraints into a problem that a solver like intlinprog would handle? Of course.

INTLINPROG assumes a matrix formulation for the equality constraints, where the objective (x) is a VECTOR. Forget about y for now. Can we write that summation as a set of linear equality constraints?

For example, suppose x was a 2x3 array? So n = 2, m = 3? What would we do?

n = 2;
m = 3;

I'll write x us a symbolic array, so you can visualize what is happening.

x = sym('x',[n,m])

x = 

Now, that sum you want to compute is the sum in MATLAB of

sum(x,1)

ans = 

Do you agree? But can we form that same sum if x were unrolled into a vector? We can see how the elements of x appear, when we treat that matrix as a vector.

xunrolled = x(:)

xunrolled = 

Now, we need to use a variation of what I long ago called the kron trick.

Aeq = kron(eye(m),ones(1,n))

Aeq = 3×6

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

You may not see where this is going immediately, but now try this.

Aeq*xunrolled

ans = 

Do you see we have created the very same sums of elements? That equality constraint now becomes

A*xunrolled == d

where d is a known vector of length m.

Similarly, all of those expressions can be written in terms of matrices and vectors. A matrix times a vector is just a sum of products of elements in the end. So just take each expression in that problem, converting it into what you need in terms of matrices and vectors. Do one piece at a time. Take your time. Get each piece right, then move to the next one.

I'll even give you a hint about how to handle the inequality constraints. That is, what does this do?

Aineq = kron([1 1 1],eye(2))

Aineq = 2×6

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

Again, a variation of the kron trick. kron is terribly useful in these problems. So now, try this:

Aineq*xunrolled

ans = 

Don't forget that you will need to treat the unknowns as a vector that combines x and y, as one long vector. Some of the unknowns will be boolean, so integers from the set {0,1}. But intlinprog can handle a partially integer problem.

[Sebastian Daneli]

Great, I'll look into it. Thanks.

[Sebastian Daneli]

@John D'Errico, I've maneged to get a little bit further, but there are still some things that I don't understand.

For example, how I should express the objective function? Here is an example:

clear 

m = 3;
n = 5;

s = [ 239 225 184 ] ;
d = [ 92 82 83 69 74 ] ;
f = [ 589 766 886 ];
c = [ 14   5   6  24   6  
       9  22  26   5  21  
      16  11  23  28  24 ];

syms x11 x12 x13 x14 x15 x21 x22 x23 x24 x25 x31 x32 x33 x34 x35;
x=[x11 x12 x13 x14 x15; x21 x22 x23 x24 x25; x31 x32 x33 x34 x35];
  
A=[zeros(1,m*n),ones(1,m)]';
Aineq=kron(ones(1,n),eye(m));
Aeq=kron(ones(1,n),eye(m));  

fy=[zeros(1,m*n),f];

The x-matrix is 5 by 3, like the c-matrix. The y-vector is 1 by 3. If i add these into a vector like you say, it becomes a 1 by 18 vector, like

size(fy)

ans = 1×2

1 18

However, the second part of the objective function becomes (if I've done this right):

Aineq*diag(c(:));

which, by using the now symbolic matrix x, becomes:

Aineq*diag(c(:))*x(:)

ans = 

which is a 5 by 3 matrix. Hope you understand my confussion.

[Torsten]

You have c_ij for 1 <= i < = 5 and 1 <= j <= 3

So c in your example should be 5x3, not 3x5.