# Detecting problems with numerically challenging models

**URL:** <https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592>\
**Category:** Optimization (Mathematical)\
**Tags:** gurobi\
**Created:** [August 25, 2024, 5:50pm UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592 "2024-08-25T17:50:51Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![sstroemer](https://avatars.discourse-cdn.com/v4/letter/s/a88e4f/32.png) [@sstroemer](https://discourse.julialang.org/u/sstroemer)\
**Post date:** [August 25, 2024, 5:50pm UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592/1 "2024-08-25T17:50:51Z")

</div>

Consider loading and solving a model using

```julia
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

```julia
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 …

```julia
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

```julia
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:

1. Why am I getting this value from `objective_value` instead of `2.77822235e+05` (c.f. log, iteration 20, dual objective)?
2. 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.
3. 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

```julia
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](https://gist.github.com/sstroemer/9bb90661228d17be1f2a1992d2ffb481).

* * *

> **Click to show: Solver Log**
>
> ```console
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 26, 2024, 5:22am UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592/2 "2024-08-26T05:22:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![torressa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torressa/32/202204_2.png) [@torressa](https://discourse.julialang.org/u/torressa)\
**Post date:** [August 27, 2024, 10:49am UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592/3 "2024-08-27T10:49:41Z")

</div>

The model has large dual violations (but the dual values are available). You can check this via the attribute [`DualVio`](https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/quality.html#dualvio) (or even [`MaxVio`](https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/quality.html#maxvio)):

```julia
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.

---

<div class="post-metadata">

**Author:** ![sstroemer](https://avatars.discourse-cdn.com/v4/letter/s/a88e4f/32.png) [@sstroemer](https://discourse.julialang.org/u/sstroemer)\
**Post date:** [August 30, 2024, 12:29pm UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592/4 "2024-08-30T12:29:27Z")

</div>

> [@torressa](#):
>
> The model has large dual violations (but the dual values are available). You can check this via the attribute [`DualVio`](https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/quality.html#dualvio) (or even [`MaxVio`](https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/quality.html#maxvio)):
> 
> […]
> 
> 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.

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:

1. 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](https://www.gurobi.com/documentation/current/refman/barconvtol.html) 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`.

1. Why does the MWE result in status code `2` (`OPTIMAL`), instead of `13` (`SUBOPTIMAL`), [which states](https://docs.gurobi.com/projects/optimizer/en/current/reference/numericcodes/statuscodes.html#secstatuscodes) _“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.

1. 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 …
2. 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:

1. `ObjVal` accesses the “primal” objective value, as the log also shows it.
2. For `PreDual=0` this seems to return the incumbent primal objective value.
3. 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:

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

#### `PreDual=0`

```julia
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`

```julia
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

```

* * *

* * *

1. 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](https://www.gurobi.com/documentation/11.0/refman/objbound.html) 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).

---

<div class="post-metadata">

**Author:** ![torressa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torressa/32/202204_2.png) [@torressa](https://discourse.julialang.org/u/torressa)\
**Post date:** [September 2, 2024, 4:28pm UTC](https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592/5 "2024-09-02T16:28:47Z")

</div>

> 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.