BUG: 1.9.0rc1: `OptimizeResult` not populated when `optimize.milp` runs into `time_limit` even though a feasible solution was found.

[woutervbk]

Describe your issue.

When optimize.milp terminates due to the time limit being reached, the OptimizeResult object is not populated, even though one or more feasible solutions have been found.

Example HiGHS log below shows that a feasible solution has been found (see solution status) when the 1 second time limit is reached. HiGHS log.txt

but the corresponding OptimizeResult object (see "Error message" section) is more or less empty and the only thing that is propagated correctly is the message, saying that the time limit was reached.

Reproducing Code Example

import numpy as np
from random import randint
from scipy.optimize import milp, LinearConstraint, Bounds

A = np.zeros((100, 100))
for c in range(100):
    for v in range(100):
        A[c, v] = randint(0, 5)
b_lb = [-1*np.inf for c in range(100)]
b_ub = [25 for c in range(100)]
constraints = LinearConstraint(A, b_lb, b_ub)
variable_lb = [0 for v in range(100)]
variable_ub = [1 for v in range(100)]
variable_bounds = Bounds(variable_lb, variable_ub)
integrality = [1 for v in range(100)]
c_vector = 100 * [-1]
res = milp(
    c_vector,
    integrality=integrality,
    bounds=variable_bounds,
    constraints=constraints,
    options={"time_limit": 1, "disp": True}
)

Error message

{'fun': None,
 'message': 'Time limit reached. (HiGHS Status 13: model_status is Time limit '
            "reached; primal_status is b'At upper bound')",
 'mip_dual_bound': None,
 'mip_gap': None,
 'mip_node_count': None,
 'status': 1,
 'success': False,
 'x': None}

SciPy/NumPy/Python version information

1.9.0rc1 1.23.0 sys.version_info(major=3, minor=9, micro=13, releaselevel='final', serial=0)

[tupui]

Thank you for reporting @woutervbk. It should indeed return something in this case.

@mckib2 (I think it's you 😅), I can see that a time limit would be considered as a HighsStatusError (scipy/optimize/_highs/cython/src/_highs_wrapper.pyx:L664). ~Commenting this section with OP's code it does work just fine to return a result when there is a time limit.~ EDIT: the object is still empty. Is there a reason this was not made possible?

Presolving model
100 rows, 100 cols, 8375 nonzeros
100 rows, 100 cols, 8375 nonzeros
Objective function is integral with scale 1
Solving MIP model with:
   100 rows
   100 cols (100 binary, 0 integer, 0 implied int., 0 continuous)
   8375 nonzeros
( 0.0s) Starting symmetry detection
( 0.0s) No symmetry present
Solving root node LP relaxation
        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time
 S       0       0         0   0.00%   -inf            0               9999.00%        0      0      0         0     0.0s
 L       0       0         0   0.00%   -9.925402033    -8                24.07%     7237     28      0       817     1.0s
Solving report
  Status            Time limit reached
  Primal bound      -8
  Dual bound        -9.92540203302
  Solution status   feasible
                    -8 (objective)
                    0 (bound viol.)
                    3.88578058619e-15 (int. viol.)
                    0 (row viol.)
  Timing            1.00 (total)
                    0.02 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     2769 (total)
                    0 (strong br.)
                    669 (separation)
                    1952 (heuristics)

[tupui]

Actually it's more L682 with HighsModelStatusOPTIMAL. Forcing this to be optimal in this case, the returned object has the solution.

            fun: -7.999999999999588
        message: 'Time limit reached. (HiGHS Status 13: Time limit reached)'
 mip_dual_bound: -9.921725263529199
        mip_gap: 24.021565794121372
 mip_node_count: 2
         status: 1
        success: False
              x: array([0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 1., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 1., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
       0., 1., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.])

[woutervbk]

That looks better! But should success not also be True in this case?

[tupui]

Maybe, I don't know this optimizer. Here I just did some quick tweaks.

[woutervbk]

I would say that the industry standard is to consider a MIP optimisation to be successful if

1. a feasible solution has been found
2. and the solver terminated due to:
   • the optimal solution was found (0% mip gap)
   • the desired MIP gap being achieved (> 0%)
   • the time limit being reached

[mckib2]

This behavior is probably a result of this upstream conversation about LP solvers:

I've given some thought to whether to return a dual feasible solution if the simplex solver in HiGHS is interrupted by time or iterations (or otherwise), or when the LP is found to be infeasible, but felt that it wasn't worth it.

This was before HiGHS had a mature MIP solver available, so some tweaking of _highs_wrapper may be warranted to make sure a feasible MIP solution is available if it exists after timeout

[woutervbk]

I agree with the above statement. Due to the nature of MIP problems, many of these problems are not solved to 100% proven optimality, but rather terminated if the mip gap is deemed to be below a certain acceptable threshold.

Btw, there currently does not seem to be any way to set such an absolute or relative mip gap termination criterium in the options that are passed to HiGHS. This would be a very useful addition (much more useful than f.e. termination after max number of iterations).

[woutervbk]

@tupui, will this be fixed in the upcoming 1.9.0 release? Could we add it to the milestone?

[tupui]

@woutervbk I am afraid not. We are already close to the release with RC3 and I did not see any work started there.

[woutervbk]

Okay gotcha. Can we then add it to the 1.9.1 milestone, alt. 1.10.0? Solving this issue would make milp() a lot more useable in practice and would make its behaviour closer to commercial MILP solvers (e.g. gurobi, CPLEX, ..).

[tupui]

@woutervbk as I said above, we need someone to do the actual work. Would you like to try make a PR?

[mdhaber]

@mckib2 and I plan to do this shortly.