HiGHS reports infeasible for a feasible model

[SupplyChef]

I am using the model below. I get the right answer when setting bigM to 1_000 but I am getting a (wrong) infeasible answer when setting bigM to 1_000_000. I understand a large bigM can create numerical issues but other solvers like Cbc return the right results even with the large bigM. I am wondering if there is something that can be done to better handle this case.

Thanks!

Infeasiblity report:
Presolving model
1 rows, 0 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
  Solution status   infeasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    1 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)

Model:

using HiGHS
function get_orders(arrivals, lead_times, arriving_orders, minimum_order_quantity, maximum_order_quantity)
    m = Model(HiGHS.Optimizer)

    bigM = 1_000_000

    times = 1:length(arrivals)
    @variable(m, orders[times] >= 0)
    @variable(m, has_orders[times], Bin)

    @constraint(m, [t=times], orders[t] <= bigM * has_orders[t])

    if minimum_order_quantity > 0
        @constraint(m, [t=times], orders[t] >= minimum_order_quantity * has_orders[t])
    end
    if !isinf(maximum_order_quantity)
        @constraint(m, [t=times], orders[t] <= maximum_order_quantity)
    end

    @constraint(m, [t=times], arrivals[t] == sum(orders[t1] for t1 in times if t1 + lead_times[t1] == t) + arriving_orders[t])

    optimize!(m)

    return value.(m[:orders])
end

Call: get_orders([1,1,0,3,4], [1,1,1,1,1], [1,0,0,0,0],0,10)

[odow]

Nice find. I'll report this upstream.

As a temporary fix, you should be able to go bigM = min(1_000_000, maximum_order_quantity). That's likely better even if HiGHS fixes the numerical problem.

Here's a simpler example:

julia> using JuMP, HiGHS

julia> function main()
           model = direct_model(HiGHS.Optimizer())
           @variable(model, x, Bin)
           @variable(model, y)
           @constraint(model, y == 1)
           @constraint(model, y <= 1_000_000 * x)
           optimize!(model)
           highs = backend(model)
           Highs_writeModel(highs, "bug.mps")
           return solution_summary(model; verbose = true)
       end
main (generic function with 1 method)

julia> main()
Presolving model
0 rows, 0 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)
                    0 (bound viol.)
                    0 (int. viol.)
                    1 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)
ERROR:   MIP solver claims optimality, but with num/sum/max primal(1/1/1) and dual(-1/inf/inf) infeasibilities
WARNING: There are empty or excessively-long Column names: using constructed names with prefix C
WARNING: There are empty or excessively-long Row names: using constructed names with prefix R
* Solver : HiGHS

* Status
  Termination status : OPTIMAL
  Primal status      : INFEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Result count       : 1
  Has duals          : false
  Message from the solver:
  "kHighsModelStatusOptimal"

* Candidate solution
  Objective value      : 0.00000e+00
  Objective bound      : 0.00000e+00
  Primal solution :
    x : 0.00000e+00
    y : 1.00000e+00

* Work counters
  Solve time (sec)   : 7.32128e-04
  Simplex iterations : -1
  Barrier iterations : -1

julia> using HiGHS_jll

julia> HiGHS_jll.highs() do exe
           run(`$(exe) bug.mps`)
       end
Running HiGHS 1.2.1 [date: , git hash: ]
Copyright (c) 2022 ERGO-Code under MIT licence terms
Number of BV entries in BOUNDS section is 1
MIP  bug has 2 rows; 2 cols; 3 nonzeros; 1 integer variables
Presolving model
0 rows, 0 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)
                    0 (bound viol.)
                    0 (int. viol.)
                    1 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)
ERROR:   MIP solver claims optimality, but with num/sum/max primal(1/1/1) and dual(-1/inf/inf) infeasibilities
Process(`/Users/oscar/.julia/artifacts/5fc32a7aeff18bf5fa92edf45ff0edd6302c4a41/bin/highs bug.mps`, ProcessExited(0))

MPS file:

NAME        
ROWS
 N  COST
 E  R0      
 L  R1      
COLUMNS
    MARK0000  'MARKER'                 'INTORG'
    C0        R1        -1000000
    MARK0001  'MARKER'                 'INTEND'
    C1        R0        1
    C1        R1        1
RHS
    RHS_V     R0        1
BOUNDS
 BV BOUND     C0      
 FR BOUND     C1      
ENDATA

[odow]

Closing in favor of upstream issue https://github.com/ERGO-Code/HiGHS/issues/771.

@SupplyChef you can make your code more robust to these errors sneaking through by doing something like:

optimize!(m)
@assert termination_status(m) == OPTIMAL
@assert primal_status(m) == FEASIBLE_POINT

You should never blindly grab a solution with value without checking the statuses.