Find all feasible solutions for a series of MIP starts

[Morgan Rodgers]

I am using Gurobi to find all feasible solutions to a (binary) MIP problem (I have a junk F(x) == 0 objective function I don't care about). 

At the moment I am using 

model.setParam('MIPFocus', 1)

model.setParam('PoolSolutions', 2000000000)

model.setParam('PoolSearchMode', 2)

There are a couple of pieces that I want to vary:

sum of all variables == size         # for a list of 5 or 6 possible integer values of size

(fix set S_0 of variables all == 0)(fix set S_1 of variables all == 1)            # for a big list (hundreds? thousands?) of pairs of sets <S_0, S_1>

Currently I am using constraints to enforce size, S_0, and S_1 conditions; then I solve the model; then I modify those constraints and solve again.

It seems there has to be a better way using MIP starts; I'm hoping I can make size and integer variable, 

make an MIP start for each triple <size, S_0,S_1>, and then run the model.

BUT I need to make sure I find ALL feasible solutions using that partial MIP start, and I certainly don't want Gurobi to start finding solutions that don't use the start.

What settings do I need to have to make sure this happens?

[Daniel Espinoza]

Hi Morgan,

I don't think that using mip starts will help in your case, as you are truly enumerating all feasible solutions satisfying a pattern.

The only thing to be careful about is that the solution pool could (in principle) overflow (but I guess you would run out of memory before). Also, you should disable presolve, as it only guarantees to preserve `an` optimal solution.

Hope it helps.

Best,

Daniel

[Morgan Rodgers]

Hi Daniel,

I'm not sure if I'm understanding about the mip starts:  Some geometric considerations of the problem let me say that if there is a feasible solution, it is "equivalent" to one that satisfies one of the mip starts in my list. Since these are very difficult models, they are more likely to finish if I am able to exploit this and run for each of the particular patterns in question.  Are you saying that mip starts are not the best way to do this? I am certainly having more success when I fix each variable to the required values using constraints. But it feels tricky to automate my python script to set these constraints, run the model, then modify all the constraints and run again. I was hoping it was easier programmatically to just replace a partial mip start with another.

(Though I suppose that aside from the single "size" variable that is very easy to modify, I am just manually implementing "orbital fixing" which I can set Gurobi to do for me?)

There is not much danger of the solution pool overflowing, or at least I should be able to tell if that happens. In most cases I will be lucky if there is a single feasible solution.

I did not know this about the presolve, thank you for pointing it out.  I assumed the presolve just removed redundant variables and constraints. 

In a related question: I noticed that when Gurobi is listing all of the nodes it is processing (under "Current Node/Obj"), it sometimes says a number, sometimes says infeasible, and occasionally says "cutoff". What does this mean exactly? Is it abandoning this node without finding all of the feasible solutions?

[Daniel Espinoza]

Hi Morgan,

In Gurobi, mipstart only help the solver a first feasible solution, this in turn is used to use heuristics (as RINS) that require a feasible solution to further improve the primal bound, which in turn lets the solver close the gap faster, and find a solution within the user-required optimality tolerances. In your case, you said your objective function is just zero, and furthermore, you are asking for all solutions (through poolsearchmode=2), so all this will not be relevant.

When you set a bunch of variables to zero or one, and add some constraints, what you are doing is slicing the feasible region, and (maybe) allowing the solver to further reduce the fractional region through cuts. That may help in getting the solutions satisfying that set of constraints quicker, and enumerate overall all solutions, but it has nothing to do with mipstarts.

Now, Orbital branching goes exactly in the oposite direction of finding all solutions, it tries to use only one solution representative for a symmetry class, so I think there is a misunderstanding there. The Orbital Branching/fixing paper is very clear on what it does, I would recommend you to take a look; but as I said, it goes in a very different direction from what you say you want to accomplish.

Regarding what the log-file means, I would suggest you to reed this:

http://www.gurobi.com/documentation/8.0/refman/mip_logging.html

Best,

Daniel

[Morgan Rodgers]

Hi Daniel,

Thanks again for your reply.

I have downloaded the Orbital branching paper due to Ostrowski et al and will look over the details.

When I say I am looking for "all" solutions, up to equivalence is the important part. As long as I find at least one representative for each equivalence class, I can use the symmetry group in my CAS to construct all of them if necessary (though I am more interested in reduction to exactly one from each class). For example, when I am fixing variables, I am doing this based on group orbits; this slices the feasible region and so I am no longer finding all solutions in the space, but by fixing variables in a way that corresponds to a representative of each orbit (on eg triples of variables), I guarantee that I find all solutions up to equivalence.

But I do need to make absolutely certain that Gurobi finds at least one example from each class; if I am going to publish a paper saying that I have used Gurobi to classify all examples of a particular type of rare object up to isomorphism, it needs to be clear that I ensured that Gurobi did not miss some examples because it already had an optimal solution.

I have read that page on the log several times, I am not sure it contains the information I am looking for. I understand that that column is telling me about the objective of the LP relaxation at that node. In my model, I see either "-", "0.00000", "infeasible", or "cutoff".  It makes sense that "0.00000" means that the LP relaxation has feasible solutions that of course optimize the objective function, and that "infeasible" means that there are not even feasible solutions at that node with continuous variables. I don't know how exactly to interpret the "-" or the "cutoff", the latter in particular since the documentation seems to say this means the objective value at that node is much worse than the optimal value; but given that the objective function is constant, this can only mean the node is infeasible?

Best,

Morgan

[Daniel Espinoza]

Hi Morgan,

I'll try to clarify inline.

I have downloaded the Orbital branching paper due to Ostrowski et al and will look over the details.

It's a nice paper! 

 

When I say I am looking for "all" solutions, up to equivalence is the important part. As long as I find at least one representative for each equivalence class, I can use the symmetry group in my CAS to construct all of them if necessary (though I am more interested in reduction to exactly one from each class). For example, when I am fixing variables, I am doing this based on group orbits; this slices the feasible region and so I am no longer finding all solutions in the space, but by fixing variables in a way that corresponds to a representative of each orbit (on eg triples of variables), I guarantee that I find all solutions up to equivalence.

If that is the case, you only need to find ONE feasible solution for each assignment of variables, and then poolsearchmode=2 is not needed, and also a partialmipstart might help

 

But I do need to make absolutely certain that Gurobi finds at least one example from each class; if I am going to publish a paper saying that I have used Gurobi to classify all examples of a particular type of rare object up to isomorphism, it needs to be clear that I ensured that Gurobi did not miss some examples because it already had an optimal solution.

Well.... that IS a tricky thing to say.... Gurobi makes statements of correctness `within tolerances`, this is due to the fact that almost all computations are carried out in floating point arithmetic.

Now, of you don't have weird numbers in the formulation, everything should be fine....

You can always claim that the feasible clases are OK (after all, you do have a solution). Now, for the infeasible cases, what you could do is to compute an IIS, and show that indeed is an IIS. That would be an actual proof that the process is correct.

 

I have read that page on the log several times, I am not sure it contains the information I am looking for. I understand that that column is telling me about the objective of the LP relaxation at that node. In my model, I see either "-", "0.00000", "infeasible", or "cutoff".  It makes sense that "0.00000" means that the LP relaxation has feasible solutions that of course optimize the objective function, and that "infeasible" means that there are not even feasible solutions at that node with continuous variables. I don't know how exactly to interpret the "-" or the "cutoff", the latter in particular since the documentation seems to say this means the objective value at that node is much worse than the optimal value; but given that the objective function is constant, this can only mean the node is infeasible?

It could be that the node (during phase 1) is deemed infeasible (and depending on where is coming) show in the log as cutoff, but yes, it means that the node is infeasible. Now, I am not sure what you mean by `-`, can you post some lines of the log?

Hope this helps.

Best,

Daniel

[Morgan Rodgers]

If that is the case, you only need to find ONE feasible solution for each assignment of variables, and then poolsearchmode=2 is not needed, and also a partialmipstart might help

Unfortunately it is not quite that simple. While I need just one feasible solution for each equivalence class of solutions, each set of start values I use could potentially lead to several nonequivalent solutions.  I am convinced that the partialmipstart is not a good tool for what I am trying, I had misunderstood how they were used. Fortunately I think I can either make a copy of my model and fix my variables in the copy, or else name these constraints in a way that they are easy to programmatically remove and replace.

(Basically I am working with a large graph, looking for special vertex sets; if I know my special set will contain a k-clique, I can look at all of the k-cliques of the graph and take one representative from each orbit under the automorphism group to use as a start. But I do not know that a k-clique will uniquely determine the set I am looking for.)

I'm confident I am okay working within the tolerances (definitely no weird values in my model), I just wanted to make sure I had the settings correct so that I am not ignoring part of the search tree.

I have attached some of my output; in this example I had not imposed a size restriction so the zero vector was a trivial optimal solution.

Presolve time: 0.01s

Presolved: 41 rows, 81 columns, 1026 nonzeros

Variable types: 0 continuous, 81 integer (81 binary)

Root relaxation: objective 0.000000e+00, 74 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       0.0000000    0.00000  0.00%     -    0s

Optimal solution found at node 0 - now completing solution pool...

     0     0          -    0         0.00000    0.00000  0.00%     -    0s

     0     0          -    0         0.00000    0.00000  0.00%     -    0s

     0     2          -    0         0.00000    0.00000  0.00%     -    0s

 86174  3444          -   33         0.00000    0.00000  0.00%   1.9    5s

 143417  5281          -   41         0.00000    0.00000  0.00%   1.8   10s

 186948  5974          -   31         0.00000    0.00000  0.00%   1.8   15s

 224275  6404          -   38         0.00000    0.00000  0.00%   1.8   20s

 256517  6258     cutoff   30         0.00000    0.00000  0.00%   1.8   25s

 284810  6892     cutoff   46         0.00000    0.00000  0.00%   1.8   30s

 310984  7442 infeasible   35         0.00000    0.00000  0.00%   1.8   35s

 335080  7527 infeasible   44         0.00000    0.00000  0.00%   1.7   40s

 358430  7791          -   35         0.00000    0.00000  0.00%   1.7   45s

 379985  7925          -   35         0.00000    0.00000  0.00%   1.7   50s

 401791  8604          -   44         0.00000    0.00000  0.00%   1.7   55s

 422596  8650 infeasible   40         0.00000    0.00000  0.00%   1.7   60s