Originally posted by: Mambo12345

Hello everybody,

I have another problem to solve. It is a linear problem but it is pretty big and thats why the function cplexbilp needs more than 30 minutes to solve

options = cplexoptimset('Display', 'iter');
[x,fval] = cplexbilp(M, A, b, Aeq, beq,[],options);

My matrix A is approximately 40.000 x 20.000, so I have about 20.000 variables in my vector x.
I would like to add some constraints that the solver doesnt take to much time. I would like to add for example a constraint where I give him a possible x with a low function value. My idea is that the branch and bound algorithm cuts more branches off because I added a guilty constraint with a very good value. How can I do this ? Do I have to set this as x_0 ??
Or is it possible to put this constraint in A. But I have no idea how to do this.
Thank you very much for your help !

Best regards

 

Did you check the reference documentation and the many examples that ship with CPLEX? Adding a constraint just means to add a row to either A/b or Aeq/beq.

x_0 is the parameter for providing a feasible starting solution. From your description I am not clear what exactly you want to do. If you have a feasible solution then this is best added as starting solution, rather than adding a constraint that the solution must be better than the value provided by that solution.

You may want to take a look at this chapter in the user manual.

Originally posted by: Mambo12345

Thank you so much for your quick answer.
I solved my problem with cplexbilp and afterwards I solved it again and gave him the minimal value as x_0 which I calculated before. This time with the best value as a starting point it took even longer than I did it with no starting point. Thats why I wanted to try to add a row in A.
Lets say for example I found a good solution for example

x1= 0 0 1 0 1 and the function value for that x1 is M*x1=65.
How can I add this constraint to my matrix A and my vector b ?
I hope that you know what I mean.
Does it make sense when I add in my for example 8. row in A my vector M and than set the value b(8) = 65 ?
I am sorry if I ask silly questions, I am just a beginner and for me is this completely new.
Thank you very much !

Originally posted by: EdKlotz

Adding a constraint involving the objective function, which I believe is what you want to do here, is really no different than adding any other constraint.   So yes, as you described in your recent post, if you want to add a constraint of the form M*x1 <= 65, you would append the vector M to your existing A matrix, and you would append the value 65 to your b vector.   I think your previous approach is correct.    And to check your work, make use the Matlab API's ability to export an LP file of your model.

That being said, you can also accomplish essentially the same thing by setting an upper cutoff value of 65 (assuming a minimization problem; otherwise set a lower cutoff value of 65).   This has the advantage of having one fewer constraint in your model, and it also avoids a possible limitation on the information in CPLEX's pseudo cost calculations that we have sometimes encountered with an explicit constraint using the objective function to express the cutoff value.

Originally posted by: Mambo12345

Thank you for your answer. This upper cutoff value is exactly what I want, EdKlotz ! In case of a linear program I just can add the vector M to my Matrix A and and than put the value (for example 65) in the vector b.
But can I give this value also to Matlab/Cplex directly ?
Just for my understanding: Cplex cut every tree off where the upper value is higher then 65 in my case and as lower the value is as faster Cplex calculated the best solution, right ?
I am asking because I also have a quadratic problem. How is it possible to do it in case of cplexmiqp (the ctype is 'B' so binary for every value x).
If I could give an upper cutoff value in case of a quadratic problem this would be perfect for me !
Then I would easily calculate some potential good solutions, take the minimum and give that minimum to a upper cutoff value.
That should reduce the time I guess.
Thank you very much for your help.

Best regards

Originally posted by: EdKlotz

If you have a MIQP, using a cutoff is probably a better idea, because adding an explicit constraint requires adding a quadratic constraint, we potentially turns your MIQP into a MIQCP, which is a harder problem type to solve. 

The upper and lower cutoff values are parameters; since all parameters can be set from the Matlab interface you are using, you can set the cutoff value.  Have a look at the documentation and some of the .m files we provide as examples for the details.

If your model is a binary MIQP that is not performing well, you might want to have a look at the second half of the slide deck I posted at http://orwe-conference.mines.edu/info.html in the tutorials section.   It describes some ways to tighten MIQPs with all binary terms in the quadratic objective.   I don't know if your model has the characteristics needed, but it might be worth a look if good cutoff values are not enough to give you the performance you seek.

Originally posted by: Mambo12345

Wow, this is amazing !!! The slides from the tutorial section are really useful for me !
Absolutely brilliant what you have done there !
But I have one more question to you:
When I have a binary quadratic optimization problem . I can of course solve it with cplexmiqp. But I also can transform it to an linear problem with a transformation y_abcd = x_ac*x_bd. Of course I have to make some more constraints because than I would have constraints to y_abcd and all x_ab. I can add some additional constraints to cut off some branches to make it faster.
When I solve it with cplexmiqp (and set ctype as binary) it is way faster then cplexbilp.
I thought that a corresponding linear problem is faster to solve than the quadratic but not in my case.
Is a QP in general faster than a LP or is it because of my bad implementation ?
I hope that I can learn something more about this and sorry if I ask some silly questions. I am a beginner and do not know that much about optimization.
 

Originally posted by: Mambo12345

I additionally have to say that I have a QAP with around 200 of x_ab and when I do the transformation I have 40.000 y_abcd.
That is probably a reason why cplexbilp take so long to solve it, right ?!
 

Originally posted by: EdKlotz

Yes, you can do the transformation you described, but CPLEX typically will do that for you by default.   So you shouldn't bother spending the time to do it yourself unless you have a really need to do so.   CPLEX has a qtolin parameter that you can use to control this.   You can see whether CPLEX does this linearization for your MIQP by checking the presolved statistics of the run.   If presolve does the linearization, you'll see that the presolved problem is a MILP, not a MIQP.   Here's what I get with a QAP that's only 15x15 by default:

Problem name         : HCP15.lp
Objective sense      : Minimize
Variables            :     225  [Binary: 225,  Qobj: 225]
Objective nonzeros   :       0
Objective Q nonzeros :   50625
Linear constraints   :      30  [Equal: 30]
  Nonzeros           :     450
  RHS nonzeros       :      30

Variables            : Min LB: 0.000000         Max UB: 1.000000       
Objective nonzeros   : Min   : all zero         Max   : all zero       
Objective Q nonzeros : Min   : 767.8828         Max   : 91769.68       
Linear constraints   :
  Nonzeros           : Min   : 1.000000         Max   : 1.000000       
  RHS nonzeros       : Min   : 1.000000         Max   : 1.000000       
CPLEX> mip
Found incumbent of value 2319411.615586 after 0.05 sec. (0.42 ticks)
Tried aggregator 1 time.
MIP Presolve added 50400 rows and 25200 columns.
Reduced MIP has 50430 rows, 25425 columns, and 101250 nonzeros.
Reduced MIP has 25425 binaries, 0 generals, 0 SOSs, and 0 indicators.

So CPLEX did the linearization, and as you pointed out, the transformation adds variables (and also constraints).   You can try turning off the qtolin parameter, in which case CPLEX will instead convexify the objective, and you will see that the presolved model is then a (smaller) MIQP:

CPLEX> s pre qtolin
Present value for indicator to linearize products in the objective involving binary variables (for convex MIQP), or all products of bounded variables (for global QP): -1 (default is -1)
 -1 = automatic
  0 = off
  1 = on
New value for indicator to linearize products in the objective involving binary variables (for convex MIQP), or all products of bounded variables (for global QP): 0

Repairing indefinite Q in the objective.
Reduced MIQP has 30 rows, 225 columns, and 450 nonzeros.
Reduced MIQP has 225 binaries, 0 generals, 0 SOSs, and 0 indicators.
Reduced MIQP objective Q matrix has 50625 nonzeros.

So basically you have a tradeoff here:  solve a larger MILP, or a smaller MIQP.   Since MILPs are typically easier to solve than MIQPs, this can go either way.    CPLEX trieds to make the right decision by default, but sometimes you need to step in and overrule its choice.   A good discussion of this can be found at

http://bob4er.blogspot.com/2015/03/quadratic-optimization-mysteries-part-1.html

Based on my experience (some of which is listed the slides I referenced above), QAPs become challenging just with dimension 15-20,so 200 is challenging, at least regarding proving optimality.  

Originally posted by: Mambo12345

Wow, this is very interesting and good to know. Thank you for all your information !
But your QAP with 15x15 is solved in 0.05 sec ?
I have 12 people and 15 offices and every person has to go in exactly one office. I have a correct matrix A and Aeq but it takes around 300 sec to solve it. Why does it take so long ?
I have to minimize the walking time and I know how often every person has to see another person per day and how long it takes to walk from office to office.
So I have a 12x12 Matrix with the information how often they have to see eachother. A 15x15 matrix with the time from office to office and my solution X is a 12 x 15 matrix with exactly one 1 in every row and one or no 1 in every column. My idea is to formulate the QP as follows

minimize x'Hx .... where x is a vector (180 variables, 180=12*15) and H is a 180x180 matrix.
And if I solve it it takes 5 minutes. I Think it has to be faster or is it usual that it takes so long ?
Is there a faster/better way to formulate the problem ?
Is it faster to have a 15x15 problem and I add 3 people to the office ? (These 3 people would have absolutely no contact to other people, so the walking time is not affected)
It would be very nice if there is a way to solve my problem faster.
Thank you very much for your help.
 

If you think the solve is too slow then please at least share the full log output for the solve, otherwise nobody will be able to tell what is up.

Even better, share your full model (preferably exported to SAV format) then we can tell whether the issue may be in the model.

Originally posted by: Mambo12345

Elapsed time = 6410.86 sec. (807832.56 ticks, tree = 2.75 MB, solutions = 15)
 277891    16    infeasible           6777.0000     6773.0000  9241049    0.06%
 277892    19        cutoff           6777.0000     6773.0000  9241044    0.06
Cover cuts applied:  2
Zero-half cuts applied:  76
Gomory fractional cuts applied:  4

Root node processing (before b&c):
  Real time             =    3.97 sec. (897.48 ticks)
Parallel b&c, 4 threads:
  Real time             = 6407.02 sec. (806940.21 ticks)
  Sync time (average)   =  440.90 sec.
  Wait time (average)   =    1.74 sec.
                          ------------
Total (root+branch&cut) = 6410.98 sec. (807837.69 ticks)

This is what I got. I hope that you mean that with the log output.
I dont know how to export it into a SAV model thats why I attached my model.

I will also post my code for that problem:

[Aeq, A, beq, b] = Constraints(T,Y,c'); %this generates my constraints
H=quadraticmatrix(T,Y); %this generates my matrix
opt = cplexoptimset('cplex');
options.Display = 'iter';
Cplex.Param.mip.cuts.bqp = 3;
Cplex.Param.emphasis.numerical = 1;
Cplex.Param.mip.cuts.gomory =2
 

f = zeros(1,size(Y,1)*size(T,1))
ctype = [];
      for i=1:length(A),
     ctype = [ctype;'B'];
      end
ctype = ctype';

[x,fval,exitflag,output] = cplexmiqp(H, f, A, b', Aeq, beq', [], [], [], [], [], ctype ,[],options); %this solves my problem
lsg2= reshape(x,[size(T,1),size(Y,1)]); %this is to see my solution as a matrix
lsg2=lsg2'
min2=fval

__________________________________

Now my functions:

function H =quadraticmatrix(T,Y)

 merge = cell(size(Y));
for k1 = 1:size(Y,1)
for k2 = 1:size(Y,1)
merge{k1,k2} = Y(k1,k2) * T;
end
end
merge = cell2mat(merge);
H=merge+merge' % this is to get a symmetric matrix

end

 

_________________________
This is my other function to generate the constraints:

function [Aeq, A, beq, b] = Constraints(T,Y,c)
 

b = ones(1,size(T,1));
beq = ones(1,size(Y,1));

Aeq = zeros(size(Y,1),size(Y,1)*size(T,1));
A = zeros(size(T,1),size(Y,1)*size(T,1));

counter=0;
counter2=0;
counter3=0;

grenze1 = 8;
grenze2 = 14; 
grenze3 = 20; 
grenze4 = 26; 

%this is to generate Aeq

%the people are just allowed to sit in some areas, that depends on the vector c

%if Y is a 10x10 matrix than c hast do be a vector with 10 elements
for j = 0:size(Y,1) - 1
    
    if c(j+1) == 0
   
        Aeq(j+1,size(T,1)*j+1:size(T,1)*j+grenze1)=1;
        Aeq(j+1,size(T,1)*j+grenze2+1:size(T,1)*j+grenze3)=1;
    else if c(j+1) == 1
     
        Aeq(j+1,size(T,1)*j+1:size(T,1)*j+grenze1)=1;
        Aeq(j+1,size(T,1)*j+grenze2+1:size(T,1)*j+grenze3)=1;
        Aeq(j+1,size(T,1)*j+grenze3+1:size(T,1)*j+grenze4)=1;
        else
         
        Aeq(j+1,size(T,1)*j+grenze1+1:size(T,1)*j+grenze2)=1;
        Aeq(j+1,size(T,1)*j+grenze3+1:size(T,1)*j+grenze4)=1;  
            
        end
    end
end

%this is to generate my matrix A:

for j = 0:size(T,1)-1
for i = 1:size(T,1):size(Y,1)*size(T,1)  
  A(j+1,i+j) = 1;
end
end
end

In my example now my matrix T is of size 26x26 (I have 26 offices) and I have 11 people, so Y is 11x11and c is 11x1 with the information in which area the people are allowed to sit.

I hope that helps to understand my model.
Thank you very much for your help.