BUG: The status of `scipy.optimize.milp` is not accurate

[nabenabe0928]

Describe your issue.

Why doesn't scipy.optimize.milp give any errors and simply say an optimal solution was found even if one of the bounds was violated?

Reproducing Code Example

import numpy as np
import scipy


integrality = np.array([1] * 5 + [0])
c = np.array([0.0] * 5 + [-1.0])
scipy_bounds = scipy.optimize.Bounds(
    lb=np.zeros(6, dtype=float), ub=np.array([3.0] * 4 + [2.0, 1.0], dtype=float)
)
scipy_constraints = [
    scipy.optimize.LinearConstraint(A=[3.0, 2.0, 2.0, -2.0, 0.0, 1.0], ub=0.0),
    scipy.optimize.LinearConstraint(A=[-3.0, -2.0, -4.0, 2.0, 0.0, 1.0], ub=0.0),
    scipy.optimize.LinearConstraint(A=[1.0, 1.0, 1.0, 1.0, -1.0, 0.0], lb=1.0, ub=1.0),
]
result = scipy.optimize.milp(
    c=c, integrality=integrality, bounds=scipy_bounds, constraints=scipy_constraints
)
print(result.x)

The code above yields [0, 0, 0, 0, -1, 0], which violates the bounds, but it does not raise any errors.
Preferably, the optimization should yield at least one feasible solution, so if you know any good options for it, please let me know 🙏

SciPy/NumPy/Python version and system information

Python 3.11.11
NumPy 2.2.0
Scipy 1.14.1

[nabenabe0928]

The output of disp=True.

Note

Running HiGHS 1.2.0 [date: 2021-07-09, git hash: n/a]
Copyright (c) 2022 ERGO-Code under MIT licence terms
Presolving model
3 rows, 6 cols, 15 nonzeros
2 rows, 5 cols, 10 nonzeros
2 rows, 2 cols, 4 nonzeros
0 rows, 1 cols, 0 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal
WARNING: Untransformed solution with objective 0 is violated by 1 for the original model

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  Solution status   infeasible
                    0 (objective)
                    1 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (total)
                    0.00 (resolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong be.)
                    0 (separation)
                    0 (heuristics)
ERROR:   MIP solver claims optimality, but with num/sum/max primal(1/1/1) infeasibilities

[mdhaber]

What is res.status?

[mdhaber]

res.status is 0 even in SciPy 1.15 (so even with updated HiGHS). This may be a misunderstanding of how we should be reporting the status; it does not seem related to recent wrapper changes.

@jajhall do you happen to know whether we need to look for separate status messages from the MIP solver vs from the underlying LP solver or something?

[jajhall]

Odd. I'll see if I can reproduce it

[mdhaber]

Sorry for the noise @jajhall. I think this was indeed fixed in SciPy 1.15 with the new version of HiGHS.

@nabenabe0928 please upgrade to SciPy 1.15. In SciPy 1.15:

assert np.all(np.asarray(scipy_bounds.residual(result.x)) >= 0)
for constraint in scipy_constraints:
    assert np.all(np.asarray(constraint.residual(result.x)) >= 0)

passes, indicating that the constraint and bounds are satisfied. The result is:

        message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: 0.0
              x: [ 0.000e+00  1.000e+00  0.000e+00  1.000e+00  1.000e+00
                   0.000e+00]
 mip_node_count: 0
 mip_dual_bound: 0.0
        mip_gap: 0.0

[nabenabe0928]

@mdhaber

Thanks, I also confirmed that SciPy 1.15+ resolves the issue:

# box constraints: Yes
# The first constraint: 2 - 2 <= 0 --> 0 <= 0
# The second constraint: -2 + 2 <= 0 --> 0 <= 0
# The third constraint: 1 + 1 - 1 = 1 --> 1 = 1
[0. 1. 0. 1. 1. 0.]

Running HiGHS 1.8.0 (git hash: 222cce7): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 4e+00]
  Cost   [1e+00, 1e+00]
  Bound  [1e+00, 3e+00]
  RHS    [1e+00, 1e+00]
Presolving model
3 rows, 6 cols, 15 nonzeros  0s
2 rows, 5 cols, 10 nonzeros  0s
2 rows, 2 cols, 4 nonzeros  0s
0 rows, 1 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve: Optimal

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               0                  0.00%        0      0      0         0     0.0s

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 0
  Nodes             0
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)