FIX: optimize.milp: return feasible solution if available on timeout/…

…node limit (#16814)

* FIX: optimize.milp: return feasible solution if available on timeout/node limit 

Co-authored-by: Matt Haberland <mhaberla@calpoly.edu>
Co-authored-by: laenNoCode <36074576+laenNoCode@users.noreply.github.com>

doc/source/tutorial/optimize.rst

@@ -1535,7 +1535,7 @@ Finally, we can solve the transformed problem using :func:`linprog`.
     >>> bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
     >>> result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
     >>> print(result.message)
-    The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is b'At lower/fixed bound')
+    The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is At lower/fixed bound)
 
 The result states that our problem is infeasible, meaning that there is no solution vector that satisfies all the
 constraints. That doesn't necessarily mean we did anything wrong; some problems truly are infeasible.

scipy/optimize/_highs/cython/src/Highs.pxd

@@ -30,7 +30,6 @@ cdef extern from "Highs.h":
         # split up for cython below
         #const HighsModelStatus& getModelStatus(const bool scaled_model = False) const
         const HighsModelStatus & getModelStatus() const
-        const HighsModelStatus & getModelStatus(const bool scaled_model) const
 
         const HighsInfo& getHighsInfo "getInfo" () const
         string modelStatusToString(const HighsModelStatus model_status) const

scipy/optimize/_highs/cython/src/_highs_wrapper.pyx

@@ -12,10 +12,14 @@ from libcpp.map cimport map as cppmap
 from libcpp.cast cimport reinterpret_cast
 
 from .HConst cimport (
+    HIGHS_CONST_INF,
+
     HighsModelStatus,
     HighsModelStatusNOTSET,
     HighsModelStatusMODEL_ERROR,
     HighsModelStatusOPTIMAL,
+    HighsModelStatusREACHED_TIME_LIMIT,
+    HighsModelStatusREACHED_ITERATION_LIMIT,
 
     HighsOptionTypeBOOL,
     HighsOptionTypeINT,
@@ -123,6 +127,7 @@ cdef apply_options(dict options, Highs & highs):
             'ipm_iteration_limit',
             'keep_n_rows',
             'max_threads',
+            'mip_max_nodes',
             'highs_debug_level',
             'min_threads',
             'simplex_crash_strategy',
@@ -653,21 +658,36 @@ def _highs_wrapper(
 
     # Extract what we need from the solution
     cdef HighsModelStatus model_status = highs.getModelStatus()
-    cdef HighsModelStatus scaled_model_status = highs.getModelStatus(True)
-    cdef HighsModelStatus unscaled_model_status = model_status
 
     # We might need an info object if we can look up the solution and a place to put solution
     cdef HighsInfo info = highs.getHighsInfo() # it should always be safe to get the info object
     cdef HighsSolution solution
     cdef HighsBasis basis
     cdef double[:, ::1] marg_bnds = np.zeros((2, numcol))  # marg_bnds[0, :]: lower
 
-    # If the status is bad, don't look up the solution
-    if model_status != HighsModelStatusOPTIMAL:
+    # Failure modes:
+    #     LP: if we have anything other than an Optimal status, it
+    #         is unsafe (and unhelpful) to read any results
+    #    MIP: has a non-Optimal status or has timed out/reached max iterations
+    #             1) If not Optimal/TimedOut/MaxIter status, there is no solution
+    #             2) If TimedOut/MaxIter status, there may be a feasible solution.
+    #                if the objective function value is not Infinity, then the
+    #                current solution is feasible and can be returned.  Else, there
+    #                is no solution.
+    mipFailCondition = model_status not in {
+        HighsModelStatusOPTIMAL,
+        HighsModelStatusREACHED_TIME_LIMIT,
+        HighsModelStatusREACHED_ITERATION_LIMIT,
+    } or (model_status in {
+        HighsModelStatusREACHED_TIME_LIMIT,
+        HighsModelStatusREACHED_ITERATION_LIMIT,
+    } and (info.objective_function_value == HIGHS_CONST_INF))
+    lpFailCondition = model_status != HighsModelStatusOPTIMAL
+    if (highs.getLp().isMip() and mipFailCondition) or (not highs.getLp().isMip() and lpFailCondition):
         return {
             'status': <int> model_status,
             'message': f'model_status is {highs.modelStatusToString(model_status).decode()}; '
-                       f'primal_status is {utilBasisStatusToString(<HighsBasisStatus> info.primal_solution_status)}',
+                       f'primal_status is {utilBasisStatusToString(<HighsBasisStatus> info.primal_solution_status).decode()}',
             'simplex_nit': info.simplex_iteration_count,
             'ipm_nit': info.ipm_iteration_count,
             'fun': None,
@@ -689,7 +709,6 @@ def _highs_wrapper(
         res = {
             'status': <int> model_status,
             'message': highs.modelStatusToString(model_status).decode(),
-            'unscaled_status': <int> unscaled_model_status,
 
             # Primal solution
             'x': [solution.col_value[ii] for ii in range(numcol)],

scipy/optimize/_linprog_highs.py

@@ -115,6 +115,7 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
                    primal_feasibility_tolerance=None,
                    ipm_optimality_tolerance=None,
                    simplex_dual_edge_weight_strategy=None,
+                   mip_max_nodes=None,
                    **unknown_options):
     r"""
     Solve the following linear programming problem using one of the HiGHS
@@ -179,6 +180,10 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
         Curently, using ``None`` always selects ``'steepest-devex'``, but this
         may change as new options become available.
 
+    mip_max_nodes : int
+        The maximum number of nodes allotted to solve the problem; default is
+        the largest possible value for a ``HighsInt`` on the platform.
+        Ignored if not using the MIP solver.
     unknown_options : dict
         Optional arguments not used by this particular solver. If
         ``unknown_options`` is non-empty, a warning is issued listing all
@@ -339,6 +344,7 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
         'dual_feasibility_tolerance': dual_feasibility_tolerance,
         'ipm_optimality_tolerance': ipm_optimality_tolerance,
         'log_to_console': disp,
+        'mip_max_nodes': mip_max_nodes,
         'output_flag': disp,
         'primal_feasibility_tolerance': primal_feasibility_tolerance,
         'simplex_dual_edge_weight_strategy':

scipy/optimize/_milp.py

@@ -130,13 +130,14 @@ def _milp_iv(c, integrality, bounds, constraints, options):
 
     # options IV
     options = options or {}
-    supported_options = {'disp', 'presolve', 'time_limit'}
+    supported_options = {'disp', 'presolve', 'time_limit', 'node_limit'}
     unsupported_options = set(options).difference(supported_options)
     if unsupported_options:
         message = (f"Unrecognized options detected: {unsupported_options}. "
                    "These will be passed to HiGHS verbatim.")
         warnings.warn(message, RuntimeWarning, stacklevel=3)
-    options_iv = {'log_to_console': options.get("disp", False)}
+    options_iv = {'log_to_console': options.pop("disp", False),
+                  'mip_max_nodes': options.pop("node_limit", None)}
     options_iv.update(options)
 
     return c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options_iv
@@ -225,6 +226,9 @@ def milp(c, *, integrality=None, bounds=None, constraints=None, options=None):
         disp : bool (default: ``False``)
             Set to ``True`` if indicators of optimization status are to be
             printed to the console during optimization.
+        node_limit : int, optional
+            The maximum number of nodes (linear program relaxations) to solve
+            before stopping. Default is no maximum number of nodes.
         presolve : bool (default: ``True``)
             Presolve attempts to identify trivial infeasibilities,
             identify trivial unboundedness, and simplify the problem before

scipy/optimize/tests/test_milp.py

@@ -273,6 +273,48 @@ def test_infeasible_prob_16609():
     np.testing.assert_equal(res.status, 2)
 
 
+_msg_time = "Time limit reached. (HiGHS Status 13:"
+_msg_iter = "Iteration limit reached. (HiGHS Status 14:"
+
+
+@pytest.mark.skipif(np.intp(0).itemsize < 8,
+                    reason="Unhandled 32-bit GCC FP bug")
+@pytest.mark.slow
+@pytest.mark.timeout(360)
+@pytest.mark.parametrize(["options", "msg"], [({"time_limit": 10}, _msg_time),
+                                              ({"node_limit": 1}, _msg_iter)])
+def test_milp_timeout_16545(options, msg):
+    # Ensure solution is not thrown away if MILP solver times out
+    # -- see gh-16545
+    rng = np.random.default_rng(5123833489170494244)
+    A = rng.integers(0, 5, size=(100, 100))
+    b_lb = np.full(100, fill_value=-np.inf)
+    b_ub = np.full(100, fill_value=25)
+    constraints = LinearConstraint(A, b_lb, b_ub)
+    variable_lb = np.zeros(100)
+    variable_ub = np.ones(100)
+    variable_bounds = Bounds(variable_lb, variable_ub)
+    integrality = np.ones(100)
+    c_vector = -np.ones(100)
+    res = milp(
+        c_vector,
+        integrality=integrality,
+        bounds=variable_bounds,
+        constraints=constraints,
+        options=options,
+    )
+
+    assert res.message.startswith(msg)
+    assert res["x"] is not None
+
+    # ensure solution is feasible
+    x = res["x"]
+    tol = 1e-8  # sometimes needed due to finite numerical precision
+    assert np.all(b_lb - tol <= A @ x) and np.all(A @ x <= b_ub + tol)
+    assert np.all(variable_lb - tol <= x) and np.all(x <= variable_ub + tol)
+    assert np.allclose(x, np.round(x))
+
+
 def test_three_constraints_16878():
     # `milp` failed when exactly three constraints were passed
     # Ensure that this is no longer the case.