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

ENH: optimize: milp: add function for mixed integer linear programming #31

Closed
wants to merge 5 commits into from

Conversation

mdhaber
Copy link
Collaborator

@mdhaber mdhaber commented Dec 26, 2021

Reference issue

gh-28

What does this implement/fix?

Adds milp function for mixed integer linear programming

Additional information

TODO:

  • Consider adding input validation to Bounds, LinearConstraint and making some arguments optional
  • Consider accepting array of LinearConstraints
  • Is options input validation in _highs_wrapper complete? (I think there is at least some)
  • Decide what fields to return
  • Document options and returns more carefully
  • Debug intermittent segfault in test_milp_6! Haven't been able to reproduce it.
  • Be consistent about l/lb, u/ub
  • Resolve other questions in WIP: HiGHS MIP wrapper #28
  • Open PR on main SciPy issue and solicit feedback on mailing list

Other thoughts:

  • I decided against raising exceptions. If someone solves a problem without try...except and an exception is raised, they can't access the information inside the object. We'll just always return an OptimizeResult but only guarantee the availability of certain attributes.
  • I decided to allow tuples that can be turned into LinearConstraints and Bounds in addition to LinearConstraint/Bound objects. That way the user doesn't need to import anything but milp to use milp.
  • I decided to allow constraints to be a sequence of constraints or a single constraint.

scipy/optimize/_milp.py Outdated Show resolved Hide resolved
@mckib2
Copy link
Owner

mckib2 commented Dec 26, 2021

Debug intermittent segfault in test_milp_6!

I ran this problem 100 times locally in a loop and didn't get a segfault. How often/what OS/platform are you seeing it on?

scipy/optimize/_milp.py Outdated Show resolved Hide resolved
c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options = args_iv

res = _highs_wrapper(c, indptr, indices, data, b_l, b_u,
lb, ub, integrality, options)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When I run locally, it appears that when disp is not explicitly set to False, HiGHS has logging to console set on by default, so we may need to add that option here (or in _highs_wrapper)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had code that set all the options to default values, but then I deleted it. There seemed to be some options input validation in _highs_wrapper, so I figured all that default-setting and such was happening there. If not, LMK where you think all of that should go. Probably should be in one place rather than split between the two.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a default 'log_to_console': False. Does disp map perfectly to log_to_console now? When we were wrapping HiGHS originally, there were a ton of different verbosity options. It looks like that may have changed?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does disp map perfectly to log_to_console now?

For our purposes I think the answer is yes.

It looks like that may have changed?

Correct, message_level went away altogether and there is now a separate dev logging channel in addition to standard information in log_to_console. I have been advised to ignore the dev logging and stick to just log_to_console.

objective function while satisfying the constraints.
fun : float
The optimal value of the objective function ``c @ x``.
simplex_nit : int
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think all simplex_nit/ipm_nit/crossover_nit will always be 0 (or -1 -- I think there might still be upstream thinking to do here) and aren't significant for MIP solutions. We can instead report mip_node_count

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is that from experimentation or did you find some documentation?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought I had asked, but I can't seem to find the thread if I did. In any case, yes -- empirically this what I find happens. The counts are set to -1 at the start of a problem and are never touched by the MIP solver, so stay at -1

@mckib2
Copy link
Owner

mckib2 commented Dec 26, 2021

@mdhaber Added a couple comments, don't need to necessarily address all in this PR. Let me know when you're ready to merge and afterward we can on a PR against the main scipy repo and send a note out to the mailing list

@mdhaber
Copy link
Collaborator Author

mdhaber commented Dec 28, 2021

Let me know when you're ready to merge

I opened this up against your branch but didn't intend to merge it. I think if you open gh-28 on the main repo we can merge that by ourselves because it's pretty minimal (just a new argument). Then we can open this one on the main repo and send a note to mailing list, etc.. Sound good?

@mckib2
Copy link
Owner

mckib2 commented Dec 28, 2021

I opened this up against your branch but didn't intend to merge it. I think if you open gh-28 on the main repo we can merge that by ourselves because it's pretty minimal (just a new argument). Then we can open this one on the main repo and send a note to mailing list, etc.. Sound good?

Yep, sounds good. I'll close gh-28 with latest doc changes (I also removed rounding as discussed with Julian) and I'll open a PR against the main repo shortly

scipy/optimize/_milp.py Outdated Show resolved Hide resolved
scipy/optimize/_milp.py Outdated Show resolved Hide resolved
scipy/optimize/_milp.py Outdated Show resolved Hide resolved
@mdhaber
Copy link
Collaborator Author

mdhaber commented Jan 9, 2022

@mckib2 could you let me know about the things above? Once I get the documentation and your thoughts on which of the options we should include in the documentation, I'll submit the PR to SciPy and send a post to the mailing list.

@mckib2
Copy link
Owner

mckib2 commented Jan 9, 2022

@mckib2 could you let me know about the things above? Once I get the documentation and your thoughts on which of the options we should include in the documentation, I'll submit the PR to SciPy and send a post to the mailing list.

Did I respond to all the things you needed info about?

@mdhaber
Copy link
Collaborator Author

mdhaber commented Jan 9, 2022

Is this resolved?

Could you also consider merging scipygh-15297 if you agree that Andy can add what he wants later? That way I can open mdhaber#67 on SciPy and hopefully we can get that merged, because milp relies for some input standardization.

I would also appreciate your thoughts on what options and results we should document:

Options
disp (as discussed above), time_limit, presolve are definitely good
Reviewing gh-28, I still don't think we've heard back about which other options definitely have an effect, and I don't want to include options unless we're sure they do something. OK if I only document the three options above for now? It looks like I do need to do input validation so that we return an OptimizeWaning for unknown options. I will allow the other highs options to be passed in without warning, but I won't document them.

Results
We'll guarantee that there will be success, status, and message. Should we convert status to the corresponding number used in linprog, or leave it as HiGHS outputs it? Which status numbers are considered success?
We'll mention x, fun, mip_dual_bound, mip_gap, and max_integrality_violation may be available, depending on the status.
We won't return all the other things that _highs_wrapper currently returns, e.g. slack, unscaled_status, lambda, marg_bnds.
Once scipygh-15297 and mdhaber#67 merge, we can mention how those can be used to check constraint/bound residuals.

@mckib2
Copy link
Owner

mckib2 commented Jan 10, 2022

Is this resolved?

yes

Options

If we pull the latest HiGHS in I think we'll get a dual gap tolerance we can add in addition to those you listed. Starting out with a small set of options is not necessarily a bad thing -- we'll get to add in over time only the ones that users need instead of throwing irrelevant ones at everyone

Should we convert status to the corresponding number used in linprog

Good question -- I think it might be better to convert them to linprog-style so they are comparable. If there's a different status code depending on whether you're running linprog(integrality is not Noe) or milp that might cause some confusion.

Which status numbers are considered success?

Might be a good question for Julian -- my inclination is to say if there's a feasible solution available in the result object, then it's at least a partial success.

We won't return all the other things that _highs_wrapper currently returns, e.g. slack, unscaled_status, lambda, marg_bnds.

Sounds good

@mckib2
Copy link
Owner

mckib2 commented Jan 10, 2022

FYI, this seems related to status reporting: ERGO-Code/HiGHS#665

@mdhaber
Copy link
Collaborator Author

mdhaber commented Jan 23, 2022

  • I added input validation for options. Unrecognized options will be passed to HiGHS verbatim.
  • I documented the three supported options.
  • I assembled the desired OptimizeResult
  • I documented the attributes of an OptimizeResult.
  • I hooked up 'disp' to 'log_to_console'
  • I confirmed that 'presolve' works by looking at the 'log_to_console' output.
  • I confirmed that 'time_limit' works by running a problem that takes a long time to solve fully with a time limit of 1. I get the following output:
 crossover_nit: -1
           fun: None
       ipm_nit: -1
       message: "model_status is Time limit reached; primal_status is b'At lower/fixed bound'"
   simplex_nit: -1
        status: 13

Seems to me we should at least get the nits mip_node_count, no?

_highs_wrapper is returning x as a list instead of, e.g., an array. Is that expected?

@mdhaber mdhaber closed this Jan 23, 2022
@mdhaber
Copy link
Collaborator Author

mdhaber commented Jan 23, 2022

Opened on main repo.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants