Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incremental Solving not supported? #7

Closed
drlagos opened this issue Aug 21, 2019 · 16 comments

Comments

@drlagos
Copy link

commented Aug 21, 2019

Hello, I've been trying to learn how to do boolean modelling in Mixed Integer Programming, and I would like to get all solutions to the problem. Here is a simple example model where blocks with ports need to be connected, with some symmetry breaking on the blocks and ports (should have 3 solutions):

import mip

mm = mip.Model(solver_name="cbc")
mm.verbose = False

# Blocks
bus = mm.add_var("bus", 0, 1, var_type=mip.BINARY)
battery_1 = mm.add_var("battery_1", 0, 1, var_type=mip.BINARY)
battery_2 = mm.add_var("battery_2", 0, 1, var_type=mip.BINARY)

# Ports
bus_p1 = mm.add_var("bus_p1", 0, 1, var_type=mip.BINARY)
bus_p2 = mm.add_var("bus_p2", 0, 1, var_type=mip.BINARY)
battery_1_p = mm.add_var("battery_1_p", 0, 1, var_type=mip.BINARY)
battery_2_p = mm.add_var("battery_2_p", 0, 1, var_type=mip.BINARY)

# Constraints
mm.add_constr(battery_1 + battery_2 >= 0)
mm.add_constr(battery_1 + battery_2 <= 2)
mm.add_constr(bus == 1)
mm.add_constr(battery_1 >= battery_2)
mm.add_constr(battery_1 == battery_1_p)
mm.add_constr(battery_2 == battery_2_p)
mm.add_constr(bus >= bus_p1)
mm.add_constr(bus >= bus_p2)
mm.add_constr(bus_p1 <= bus_p2)
mm.add_constr(bus_p1 + bus_p2 == battery_1_p + battery_2_p)

count = 0
for i in range(0, 10):
    status = mm.optimize()
    if status == status.FEASIBLE or status == status.OPTIMAL:
        print(f"{bus}: {bus.x}")
        print(f"{battery_1}: {battery_1.x}")
        print(f"{battery_2}: {battery_2.x}")
    else:
        print("INFEASIBLE")
        break
    count += 1
    # or(a != a.x, b != b.x, c != by)
    flip = [1 - v if v.x else v for v in (bus, battery_1, battery_2)]
    mm.add_constr(mip.xsum(flip) >= 1)
print(f"CBC (Python-MIP) found {count} solutions")

However this code only gives me 1 solution, 10 times.

If I do the same with google Or-Tools, it works as expected:

from ortools.linear_solver import pywraplp as or_tools

om = or_tools.Solver("cbc", problem_type=or_tools.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# Instances
bus = om.BoolVar("bus")
battery_1 = om.BoolVar("battery_1")
battery_2 = om.BoolVar("battery_2")

# Ports
bus_p1 = om.BoolVar("bus_p1")
bus_p2 = om.BoolVar("bus_p2")
battery_1_p = om.BoolVar("battery_1_p")
battery_2_p = om.BoolVar("battery_2_p")

# Constraints
om.Add(battery_1 + battery_2 >= 0)
om.Add(battery_1 + battery_2 <= 2)
om.Add(bus == 1)
om.Add(battery_1 >= battery_2)
om.Add(battery_1 == battery_1_p)
om.Add(battery_2 == battery_2_p)
om.Add(bus >= bus_p1)
om.Add(bus >= bus_p2)
om.Add(bus_p1 <= bus_p2)
om.Add(bus_p1 + bus_p2 == battery_1_p + battery_2_p)

count = 0
# Limit number of solutions
for i in range(0, 10):
    status = om.Solve()
    if status == om.FEASIBLE or status == om.OPTIMAL:
        print(f"{bus}: {bus.solution_value()}")
        print(f"{battery_1}: {battery_1.solution_value()}")
        print(f"{battery_2}: {battery_2.solution_value()}")
    else:
        print("INFEASIBLE")
        break
    count += 1
    # or(a != a.x, b != b.x, c != by)
    flip = [1 - v if v.solution_value() else v for v in (bus, battery_1, battery_2)]
    om.Add(sum(flip) >= 1)
print(f"CBC (Or-Tools) found {count} solutions")

This gives the expected 3 solutions.

Is there anything I'm missing, like a different way to restrict solutions, or is python-mip just not intended to be used like this?

@drlagos

This comment has been minimized.

Copy link
Author

commented Aug 21, 2019

Just tested with PuLP, it also works correctly with 3 solutions (albeit very, very slow).

@tuliotoffolo

This comment has been minimized.

Copy link
Contributor

commented Aug 21, 2019

@drlagos

This comment has been minimized.

Copy link
Author

commented Aug 21, 2019

Great, thanks! I'm quite interested in python-mip since it lets me tell CBC to aim for feasible solutions or aim for optimal solutions (the ability to do both is what I'm interested in, hence investigating MIP solvers, OptSMT is still quite new and there's only 1 solver that is both good at it and has a decent license).

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Aug 21, 2019

Hi,

There are two bugs here:

  1. CBC storing updated bounds on variables after optimizing (we'll fix it)
  2. your code has also a dangerous operation:

if v.x tests for zero. as computation in linear programming solvers use limited accuracy it is often the case that instead zero you got 0.000000000001 (which is also zero for practical purposes)
it would be better to use here:

abs(v.x) <= 1e-7

@drlagos

This comment has been minimized.

Copy link
Author

commented Aug 21, 2019

I did not realize integer variables were actually handled as reals rather than integers.
Thank you for the warning! You probably just saved me some frustrating nights debugging that issue.

@tkralphs

This comment has been minimized.

Copy link
Member

commented Aug 21, 2019

@h-g-s Although integer variables are handled as reals inside the solver, the value of any variable that is intended to be an integer should be handed to the user as an integer. I believe Cbc itself should do this rounding before handing the user the solution, but if not Cbc, it seems your interface should do it. I don't think users should need to be aware of these details.

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Aug 21, 2019

@tkralphs true, I'll check it in CBC

@spoorendonk

This comment has been minimized.

Copy link
Contributor

commented Aug 21, 2019

@h-g-s I think that 1. (the storing of updated bounds) was also the problem I had in #3 and the reason why the pricing problem could not be re-optimized. Perhaps the cutting stock example should be reverted (except for the model bug) when this issue is fixed.

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Aug 22, 2019

@spoorendonk I think that the CBC fix will not take too long...

@aopt

This comment has been minimized.

Copy link

commented Aug 26, 2019

@h-g-s Although integer variables are handled as reals inside the solver, the value of any variable that is intended to be an integer should be handed to the user as an integer. I believe Cbc itself should do this rounding before handing the user the solution, but if not Cbc, it seems your interface should do it. I don't think users should need to be aware of these details.

Rounding can lead to larger infeasibilities. I think this is why the solver would not do this automatically. Safer is to use a tighter integer feasibility tolerance (Cplex allows epint=0). Of course that can lead to longer solution times.

@tkralphs

This comment has been minimized.

Copy link
Member

commented Aug 26, 2019

The infeasibility that the user probably cares most about is the one computed with the the rounded values for the integer variables. From a practical standpoint, it doesn't matter (to the user) that the infeasibility of the given continuous solution is within a tolerance (for the linear constraints) if rounding the values of the integer variables causes it to exceed that tolerance. For practical reasons, it's difficult to ensure that the rounded solution satisfy a given feasibility tolerance, but that shouldn't prevent the interface from handing the user a solution that is fully integer feasible. Whether you do the rounding on behalf of the user or not, there is some potential for confusion, but rounding seems to eliminate the most likely source.

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Sep 6, 2019

1.4.0 released with (hopefully) a fix for the MIP reoptimization issue in CBC. Could you try it ???

@drlagos

This comment has been minimized.

Copy link
Author

commented Sep 7, 2019

I have good news and bad news.

The good news is that the bug is fixed!

The bad news is that there is a new one, setting verbosity to 0 doesn't do anything and I get a spam of messages from CBC. Though maybe i'm just doing something wrong.

my_model.verbosity = 0

Should be this right?

max_solutions also appears to be ignored, I get this from CBC:

No match for max_solutions - ? for list of commands

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Sep 7, 2019

@h-g-s

This comment has been minimized.

Copy link
Contributor

commented Sep 7, 2019

@drlagos

This comment has been minimized.

Copy link
Author

commented Sep 7, 2019

Verbosity works now :) just max_solutions still missing.
I will close the issue since not having max_solutions does not really affect me.

Excellent work!

@drlagos drlagos closed this Sep 7, 2019

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
6 participants
You can’t perform that action at this time.