Detecting problems with numerically challenging models

sstroemer · August 25, 2024, 5:50pm

Consider loading and solving a model using

using JuMP
import Gurobi

mwe = read_from_file("mwe.mps"; format = MOI.FileFormats.FORMAT_MPS)
set_optimizer(mwe, Gurobi.Optimizer)
set_attribute(mwe, "Method", 2)
set_attribute(mwe, "Crossover", 0)

optimize!(mwe)

which shows

primal_status(mwe)       # FEASIBLE_POINT::ResultStatusCode = 1
dual_status(mwe)         # FEASIBLE_POINT::ResultStatusCode = 1
termination_status(mwe)  # OPTIMAL::TerminationStatusCode = 1
has_values(mwe)          # true
has_duals(mwe)           # true

passing all my (current) checks for a “successful” solve.

However …

objective_value(mwe)       # 276039.476076866
dual_objective_value(mwe)  # ERROR: Gurobi Error 10005: Unable to retrieve attribute 'ObjBound'

Further checking the objective function, which is just min θ, shows

objective_function(mwe)         # θ
value(objective_function(mwe))  # 277822.2353388086

As far as I see, this result should not fulfill Gurobi’s default tolerances? What seems to be happening is, that Gurobi returns the primal objective value from iteration 29 (c.f. the log below) - which seems to actually be the value of the dual problem, since it’s smaller than the “dual objective” (indicating it being a max; using PreDual = 0 supports this).

Questions:

Why am I getting this value from objective_value instead of 2.77822235e+05 (c.f. log, iteration 20, dual objective)?
Is there any way to get an “objective value” that I can guarantee to be an upper bound, instead of it possibly (in this case it’s 0.36% lower than the “true” objective obtained from a numerically stable version of the model) being a lower bound? objective_bound does not work here.
Is there any way to catch these kind of instabilities, without parsing the solver log - what’s the most efficient & general approach?

Ad (3): This is a bit complex, because the only visible warning only relates to the large coefficient range - which happens for quite a lot of models a user may construct, without any immediate problems for Gurobi. The Sub-optimal termination print does not happen with default parameters. Running it with

set_attribute(mwe, "NumericFocus", 3)
set_attribute(mwe, "OptimalityTol", 1e-9)
set_attribute(mwe, "FeasibilityTol", 1e-9)
set_attribute(mwe, "BarConvTol", 1e-16)

results in dual_status(mwe) returning UNKNOWN_RESULT_STATUS::ResultStatusCode = 8, but I can’t run every model with 3 as default.

The MWE model file that was used can be found in this gist.

Click to show: Solver Log

Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Ubuntu 24.04 LTS")

CPU model: Intel(R) Core(TM) i9-14900K, instruction set [SSE2|AVX|AVX2]
Thread count: 32 physical cores, 32 logical processors, using up to 32 threads

Optimize a model with 212 rows, 31 columns and 5015 nonzeros
Model fingerprint: 0x40fa2628
Coefficient statistics:
  Matrix range     [2e-10, 1e+07]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e-08, 1e+06]
  RHS range        [2e+05, 2e+12]
Warning: Model contains large matrix coefficient range
Warning: Model contains large rhs
         Consider reformulating model or setting NumericFocus parameter
         to avoid numerical issues.
Presolve removed 6 rows and 11 columns
Presolve time: 0.00s
Presolved: 25 rows, 231 columns, 4998 nonzeros
Ordering time: 0.00s

Barrier statistics:
 Free vars  : 2
 AA' NZ     : 3.000e+02
 Factor NZ  : 3.250e+02
 Factor Ops : 5.525e+03 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0  -7.59084551e+13 -1.00000000e-08  6.16e+04 5.55e+10  3.10e+12     0s
   1  -6.44069412e+13  3.91607307e+12  1.35e+04 7.59e+09  8.76e+11     0s
   2  -5.73215351e+12  3.82831051e+12  8.17e+02 2.40e+08  7.27e+10     0s
   3  -2.08311071e+12  1.28780638e+12  1.36e+02 3.41e+07  1.82e+10     0s
   4  -6.14486535e+11  1.72823434e+09  1.11e-02 4.01e+06  2.51e+09     0s
   5  -6.73161146e+08  1.87472648e+08  1.14e-04 9.42e-05  3.50e+06     0s
   6  -1.07999480e+06  2.21048401e+06  9.50e-06 2.35e-05  1.34e+04     0s
   7  -1.21780429e+05  1.57065824e+06  1.48e-05 2.00e-05  6.90e+03     0s
   8   1.36058977e+04  9.80400032e+05  1.19e-05 5.75e-05  3.95e+03     0s
   9   1.61971020e+05  4.43101449e+05  1.44e-05 1.28e-05  1.16e+03     0s
  10   2.61663385e+05  2.98196976e+05  2.98e-05 4.41e-06  1.64e+02     0s
  11   2.73161379e+05  2.80763734e+05  3.10e-06 1.97e-05  3.08e+01     0s
  12   2.75010006e+05  2.78937021e+05  1.23e-06 1.06e-05  1.22e+01     0s
  13   2.75465037e+05  2.78163248e+05  3.62e-06 3.81e-06  5.80e+00     0s
  14   2.75641249e+05  2.77856133e+05  4.87e-06 5.13e-06  3.22e+00     0s
  15   2.76030806e+05  2.77837142e+05  7.42e-07 5.60e-06  1.59e-01     0s
  16   2.76039271e+05  2.77823529e+05  3.31e-06 3.81e-06  2.23e-02     0s
  17   2.76040424e+05  2.77822719e+05  4.99e-06 1.20e-05  1.05e-02     0s
  18   2.76041510e+05  2.77822250e+05  1.28e-05 7.63e-06  1.28e-03     0s
  19   2.76041593e+05  2.77822266e+05  9.63e-05 6.20e-06  2.67e-04     0s
  20   2.76041610e+05  2.77822236e+05  1.19e-05 5.96e-06  1.61e-05     0s
  21   2.76041620e+05  2.77822236e+05  9.22e-06 6.44e-06  7.39e-07     0s
  22   2.76041609e+05  2.77822234e+05  3.17e-06 5.48e-06  2.00e-07     0s
  23   2.76041619e+05  2.77822236e+05  8.66e-06 6.56e-06  1.64e-09     0s
  24   2.76041621e+05  2.77822237e+05  1.39e-05 6.20e-06  5.37e-10     0s
  25   2.76041610e+05  2.77822234e+05  4.09e-06 1.20e-05  2.79e-10     0s
  26   2.76041614e+05  2.77822236e+05  7.55e-06 2.86e-06  9.26e-11     0s
  27   2.76041617e+05  2.77822237e+05  3.50e-06 6.68e-06  2.11e-13     0s
  28   2.76041613e+05  2.77822238e+05  7.06e-06 5.48e-06  4.26e-14     0s
  29   2.76039476e+05  2.77822235e+05  8.27e-06 6.68e-06  8.85e-19     0s
  30   2.75442817e+05  2.77822237e+05  7.93e-06 8.23e-06  3.69e-19     0s

Barrier solved model in 30 iterations and 0.00 seconds (0.02 work units)
Optimal objective 2.76039476e+05

odow · August 26, 2024, 5:22am

I have no comments. This seems like a good question for Gurobi’s support: cc @torressa

torressa · August 27, 2024, 10:49am

The model has large dual violations (but the dual values are available). You can check this via the attribute DualVio (or even MaxVio):

julia> MOI.get(mwe, Gurobi.ModelAttribute("DualVio"))
0.16951326635776526

julia> MOI.get(mwe, Gurobi.ModelAttribute("MaxVio"))
0.16951326635776526

These are well above the tolerance and are no surprise given the coefficient ranges. Default parameters (not setting Method=2 + Crossover=0) or even just enabling crossover fixes this violation and works as expected.

sstroemer · August 30, 2024, 12:29pm

Hi @torressa, thanks for checking this out and taking the time. Preface: I know, this is a MWE for a numerically challenging model. Solving such models is sometimes necessary. Other parameters are not feasible for a “real” model: Crossover=0 actually worsens any solution, since I’m explicitly not interested in “extreme” basis solutions. Even for those where it would be okay to activate crossover, Method=2 is the only way to actually get through the solve.

Therefore, let me rephrase my questions / points that I cannot wrap my head around:

Can you confirm that my guess is correct, that with default settings (PreDual = -1) it chooses to solve the dual in this MWE?

pre 2.)
BarConvTol states “The barrier solver terminates when the relative difference between the primal and dual objective values is less than the specified tolerance (with a GRB_OPTIMAL status).” Similar conditions also exist for other parameters, e.g., OptimalityTol.

Why does the MWE result in status code 2 (OPTIMAL), instead of 13 (SUBOPTIMAL), which states “Unable to satisfy optimality tolerances; a sub-optimal solution is available.” The rel. gap I see from the log is greater than 6.4e-3, while BarConvTol defaults to 1e-8.

pre 3.) When using the barrier, I can only expect the true objective to be between the lower and upper bound. Therefore, depending on whether I’m asking “How much is my objective at most/least?” I need to rely on the upper/lower bound.

How can I know which objective bound Gurobi returns as its “objective” (ObjVal)? For a parameter like PreDual I actually want to trust the heuristic a lot of times. If that however prevents me from extracting proper objective values …
Is there any way to query the second bound, that is being logged^[1]?

Test Cases

From the outside it seems like that what happens is:

ObjVal accesses the “primal” objective value, as the log also shows it.
For PreDual=0 this seems to return the incumbent primal objective value.
For PreDual=1, the primal being logged is actually the “primal objective value of the dual, that is getting solved”. ObjVal again accesses this, which then means - from the point of view of the initial problem to be solved - this is actually the incumbent dual objective value.

That means:

I have no idea whether the returned value is a lower or an upper bound on the true objective.
Something like objective_value(mwe) ≈ value(objective_function(mwe)), which a user might expect to hold (with some tol.), does not.

`PreDual=0`

set_attribute(mwe, "PreDual", 0)
set_attribute(mwe, "BarConvTol", 1e-2)

# Last iteration in log:
#   13   2.76998764e+05  2.76994994e+05  0.00e+00 2.97e-08  2.19e-02     0s

MOI.get(mwe, Gurobi.ModelAttribute("ObjVal"))  # 276998.7644756779
objective_value(mwe)                           # 276998.7644756779
value(objective_function(mwe))                 # 276998.7644756779

`PreDual=1`

set_attribute(mwe, "PreDual", 1)
set_attribute(mwe, "BarConvTol", 1e-2)

# Last iteration in log:
#   14   2.75641249e+05  2.77856133e+05  4.87e-06 5.13e-06  3.22e+00     0s

MOI.get(mwe, Gurobi.ModelAttribute("ObjVal"))  # 275641.248734008
objective_value(mwe)                           # 275641.248734008
value(objective_function(mwe))                 # 277856.1332178363

I was tempted to query ObjBound, since that seems to be what JuMP.dual_objective_value(...) seems to use for Gurobi.jl. Looking at the documentation makes me suspect that this may only be guaranteed to “work” for MILPs? @odow Might this be worth changing, or throwing a “better” error than letting it run into Gurobi Error 10005? Currently both dual_objective_value and objective_bound may “unexpectedly” fail (even for simpler models). ↩︎

torressa · September 2, 2024, 4:28pm

Can you confirm that my guess is correct, that with default settings (PreDual = -1) it chooses to solve the dual in this MWE?

I looks like it, yes.

In general, though, I would not trust the solution values of any variables with such high violations (check MaxVio) and low barconvtol.