Various problems extracting results for some models

[sstroemer]

This is partially related to this previous discourse thread, but now I ran into troubles when using HiGHS.

Expand to show: MWE, code, and full log

MWE file can be found in this gist.

Code:

using JuMP
import HiGHS

mwe = read_from_file("mwe_highs.mps")
set_optimizer(mwe, HiGHS.Optimizer)
set_attribute(mwe, "solver", "ipm")
set_attribute(mwe, "run_crossover", "off")

optimize!(mwe)

solution_summary(mwe)

Full log:

HiGHS log output
Running HiGHS 1.7.2 (git hash: 5ce7a2753): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-09, 1e+07]
  Cost   [8e-02, 2e+00]
  Bound  [1e+03, 1e+06]
  RHS    [1e-13, 4e+11]
WARNING: Problem has excessively large bounds: consider scaling the bounds by 1e+2 or less, or setting option user_bound_scale to -6 or less
Presolving model
4071 rows, 179 cols, 23438 nonzeros  0s
875 rows, 179 cols, 8004 nonzeros  0s
853 rows, 157 cols, 7960 nonzeros  0s
Presolve : Reductions: rows 853(-3350); columns 157(-37); elements 7960(-15816)
Solving the presolved LP
IPX model has 853 rows, 157 columns and 7960 nonzeros
Input
    Number of variables:                                157
    Number of free variables:                           0
    Number of constraints:                              853
    Number of equality constraints:                     0
    Number of matrix entries:                           7960
    Matrix range:                                       [1e-09, 1e+07]
    RHS range:                                          [1e-13, 4e+11]
    Objective range:                                    [4e-02, 2e+00]
    Bounds range:                                       [1e+03, 1e+06]
Preprocessing
    Dualized model:                                     yes
    Number of dense columns:                            0
    Range of scaling factors:                           [4.88e-04, 5.12e+02]
IPX version 1.0
Interior Point Solve
 Iter     P.res    D.res            P.obj           D.obj        mu     Time
   0   1.61e+00 4.06e+08  -3.62747966e+06 -2.73315498e+05  7.27e+08       0s
 Constructing starting basis...
   1   7.45e-01 1.49e+08   2.16503232e+10 -5.82216633e+09  3.36e+08       0s
   2   3.35e-01 1.07e+08   2.63366163e+10 -6.81053553e+09  2.06e+08       0s
   3   1.51e-01 6.50e+07   2.48884472e+10 -7.62696235e+09  1.26e+08       0s
   4   5.94e-02 4.45e+07   1.81062320e+10 -7.24209359e+09  7.43e+07       0s
   5   1.07e-02 6.28e+06   9.16378252e+09 -3.88644637e+09  1.79e+07       0s
   6   9.24e-04 7.96e+05   1.75717528e+09 -7.63840536e+08  2.64e+06       0s
   7   1.46e-05 6.36e+04   1.34956075e+08 -7.48031367e+07  1.92e+05       0s
   8   6.84e-06 9.10e+03   6.44043972e+07 -1.30155029e+07  6.90e+04       0s
   9   9.01e-07 9.10e-03   8.93610184e+06 -2.25615838e+06  9.87e+03       0s
  10   4.00e-07 6.86e-03   6.09280309e+06 -2.24150381e+06  7.32e+03       0s
  11   1.23e-07 3.12e-04   1.90785320e+06 -1.07336378e+06  2.61e+03       0s
  12   3.72e-08 7.05e-05   3.65024020e+05 -6.48372243e+05  8.88e+02       0s
  13   4.57e-09 4.04e-05   1.77437232e+04 -5.72094405e+05  5.16e+02       0s
  14   5.51e-10 2.40e-05  -1.64894005e+05 -4.65095112e+05  2.62e+02       0s
  15   1.78e-10 8.46e-06  -2.26382751e+05 -3.39900358e+05  9.92e+01       0s
  16   1.37e-10 1.67e-06  -2.34368908e+05 -3.07124631e+05  6.36e+01       0s
  17   4.64e-11 4.17e-07  -2.58913594e+05 -2.89300544e+05  2.66e+01       0s
  18   4.96e-12 2.38e-07  -2.68725857e+05 -2.83450094e+05  1.29e+01       0s
  19   9.42e-13 2.38e-07  -2.74756709e+05 -2.79008239e+05  3.72e+00       0s
  20   6.36e-13 1.19e-07  -2.75326767e+05 -2.77990200e+05  2.33e+00       0s
  21   4.23e-13 1.19e-07  -2.75827237e+05 -2.77402153e+05  1.38e+00       0s
  22   8.39e-14 1.79e-07  -2.76758400e+05 -2.77139768e+05  3.33e-01       0s
  23   3.73e-14 1.79e-07  -2.76902232e+05 -2.77081895e+05  1.57e-01       0s
  24   7.77e-15 1.79e-07  -2.77007716e+05 -2.77061901e+05  4.74e-02       0s
  25   1.22e-15 2.98e-07  -2.77037652e+05 -2.77051172e+05  1.18e-02       0s
  26   5.64e-16 2.38e-07  -2.77045541e+05 -2.77048912e+05  2.95e-03       0s
  27   3.89e-16 1.79e-07  -2.77047733e+05 -2.77048357e+05  5.45e-04       0s
  28   3.21e-16 2.38e-07  -2.77048073e+05 -2.77048209e+05  1.19e-04       0s
  29   4.09e-16 1.19e-07  -2.77048176e+05 -2.77048192e+05  1.39e-05       0s
  30*  4.41e-16 1.79e-07  -2.77048188e+05 -2.77048188e+05  3.20e-07       0s
Summary
    Runtime:                                            0.02s
    Status interior point solve:                        optimal
    Status crossover:                                   not run
    objective value:                                    2.77048188e+05
    interior solution primal residual (abs/rel):        1.83e-04 / 5.17e-16
    interior solution dual residual (abs/rel):          8.86e-10 / 3.53e-10
    interior solution objective gap (abs/rel):          3.76e-04 / 1.36e-09
Ipx: IPM       optimal
Model   status      : Optimal
IPM       iterations: 30
Objective value     :  2.7704818697e+05
HiGHS run time      :          0.02

---

What happens
Running a (yes, unstable) MWE with solver = ipm and run_crossover = off leads to a successful termination of IPX, with the last iteration printing

30*  4.41e-16 1.79e-07  -2.77048188e+05 -2.77048188e+05  3.20e-07       0s

Extracting the dual objective values fails, or goes wrong, as can be seen in the (shortened) solution summary

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver: "kHighsModelStatusOptimal"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 2.77048e+05
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : -2.22082e+05

---

Why does that happen
The following is a list of things that I noticed that might be responsible for this:

1. Dual objective value: This is calculated here relying on MOI.Utilities.get_fallback, which:
   1. seemingly ends up with the wrong sign - but instead it's just
   2. way off, which I assume is due to "unstable" dual results (cleaning up with crossover leads to a proper calculation with the fallback)
2. Objective bound: This is calculated here using an old workaround, that seems to fail here, since Highs_getDoubleInfoValue(model, "mip_dual_bound", p) actually returns 0.0.
3. Relative gap : As with other things, this is probably tuned towards MIPs, and therefore extracts mip_gap, which I just think could be made more consistent by either allowing IPM gaps to be read out with this, or (in the future) providing a separate relative_mip_gap(...) instead of relative_gap(...), etc.

---

Other things

1. JuMP.objective_bound states "Return the best known bound [...]", which (is this correct?) is interpreted as lower bound for minimization and upper bound for maximization problems. Maybe an addition to the docs could help to mention which bound it refers to?
2. As with the discourse thread, I'm kind of unsure how to detect these kind of things properly? The solver reports "optimal", good "residuals", etc.

---

Thoughts

1. [JuMP.jl] Is there any way to improve "support" for primal-dual solutions with a gap > 0 (e.g., obtained from an ipm)?
2. [HiGHS] It would be really helpful if there was any support to extract more results obtained from IPX.
3. [HiGHS.jl] Depending on (2.), maybe there is another way to implement MOI.DualObjectiveValue?
4. [JuMP.jl] It seems to me that some solvers might report a status (e.g., "optimal"), while knowing that any solution might be "wrong". Indications might be available in the sense of some log messages that "warn" the user, but often an end-user might not see them at all (and tbh they are often issued way too early to be "trusted"). Could an on-demand estimate_solution_quality(...) functionality make sense (not saying I have a definite idea how this could be judged)?

[sstroemer]

[HiGHS]

I've just noticed another weird thing in some logs:

Ipx: IPM       optimal
ERROR:   LP solver claims optimality, but with num/max/sum primal(1/2.38396e-05/2.38396e-05)and dual(1/335143/335143) infeasibilities
Model   status      : Optimal

Note: This is a super small model (IPX model has 2 rows, 7 columns and 8 nonzeros), see MPS at the end of this comment.

1. Maybe - from a solver's point of view: Is it really "good" to terminate with a model status "optimal" for a result that is tagged as being both primal and dual infeasible, and triggers an "ERROR" line in the log?
2. Why does that behaviour change when setting dual_feasibility_tolerance = 1e-10? This triggers running one more iteration, but why does it even terminate with 1e-8 (or 1e-9) early if max_dual_infeasibility clearly shows something larger?
3. What happens "before crossover" internally? Setting run_crossover = on logs that no primal/dual pushes are required, crossover_iteration_count returns 0, which leads me to believe it did not actually perform any pushes? Still, the (really) bad dual results are perfectly cleaned up.

Note ad 3.: I'm interested in understanding if there are some low hanging fruits that can be exploited for a large model, solved with IPX, without running full crossover, but still being able to "clean up a bit".

---

Model MPS file:

Click to expand

NAME
ROWS
 N  OBJ
 L  c1_1
 E  c1
COLUMNS
    x[2112]   c1        2
    x[54531]  c1_1      1
    x[54531]  c1        -1
    x[54531]  OBJ       0.02
    x[61175]  c1_1      -1
    z_pos[2112] c1      2
    z_pos[2112] OBJ     1000000
    z_pos[61175] c1_1   -1
    z_pos[61175] OBJ    1000000
    z_neg[2112] c1      -2
    z_neg[2112] OBJ     1000000
    z_neg[61175] c1_1   1
    z_neg[61175] OBJ    1000000
RHS
    rhs       c1_1      0
    rhs       c1        0
RANGES
BOUNDS
 FX bounds    x[2112]   35932.5417904476
 LO bounds    x[54531]  0
 PL bounds    x[54531]
 FX bounds    x[61175]  0.0011846243591959321
 LO bounds    z_pos[2112] 0
 PL bounds    z_pos[2112]
 LO bounds    z_pos[61175] 0
 PL bounds    z_pos[61175]
 LO bounds    z_neg[2112] 0
 PL bounds    z_neg[2112]
 LO bounds    z_neg[61175] 0
 PL bounds    z_neg[61175]
ENDATA

[sstroemer]

Ad dual_objective_value(...):

This currently does not cover all cases, but maybe it leads somewhere...

1. Calling JuMP.dual_objective_value(model) for a HiGHS-solved model, that falls back to estimating, by manually calculating the value of the dual objective function, is slow. Naively calling it at two locations in my code (that do not know of each other), instead of caching and passing it, means that my iterative solve now spends >90% of the total time in the result extraction.
2. As far as I see, the main culprit is Highs_getRowsByRange, which is used to get the constant of the constraint object. For many special cases, this will be zero anyways (since everything is contained in the ConstraintSet).

So I've tried the following workaround:

function fast_est_dual_obj_val(jump_model::JuMP.Model)
    get_b(set::JuMP.MOI.GreaterThan) = set.lower
    get_b(set::JuMP.MOI.LessThan) = set.upper
    get_b(set::JuMP.MOI.EqualTo) = set.value

    all_con = JuMP.all_constraints(jump_model; include_variable_in_set_constraints = true)

    π = JuMP.dual.(all_con)
    b = get_b.(JuMP.MOI.get.(jump_model, JuMP.MOI.ConstraintSet(), all_con))

    return JuMP.objective_sense(jump_model) == JuMP.MAX_SENSE ? -(π' * b) : (π' * b)
end

which resulted (for a random model) in the following:

fast_est_dual_obj_val(model)      # 35317.33044169648
JuMP.dual_objective_value(model)  # 35317.33044169651

@b fast_est_dual_obj_val(model)      # 662.679 μs (59947 allocs: 1.512 MiB)
@b JuMP.dual_objective_value(model)  # 10.810 ms (43727 allocs: 1.406 MiB)

So, ~6% of total time are spent in the actual calculation, the rest seems to be repeated calls into HiGHS (maybe it would already be sufficient to just fetch all rows in a single call to Highs_getRowsByRange, instead of doing it row by row?).

[odow]

The Highs_getRowsByRange is #207

I don't know if it is useful for me to respond point-by-point, since it hits many of the same themes as https://discourse.julialang.org/t/detecting-problems-with-numerically-challenging-models/118592.

• You're solving a numerically difficult problem
• You are disabling some features (in the name of performance) which "fix up" the imprecision of the IPM solvers
• The solvers happily return solutions and claim optimality while having primal/dual gaps that are large
• They offer limited ways for novice users to investigate the quality of the solution being returned (I don't think energy system modelers are served by telling them "warning: DualVio is large".

I don't think this is a problem with HiGHS.jl or Gurobi.jl specifically, so much as a general issue that affects all solvers and models in math programming.

You are just noticing it because:

1. You're solving large difficult problems
2. JuMP automatically reports some statistics that other algebraic modeling languages do not, like the dual objective value.

If we didn't report dual objective value in the model summary, you would be none the wiser, and you would happily carry on using your primal solution 😄

Perhaps we should schedule a call to talk through this? It might be easier than posting back-and-forthh.

[odow]

@joaquimg and I have just noticed:

    Matrix range:                                       [1e-09, 1e+07]
    RHS range:                                          [1e-13, 4e+11]

You can't have 16 orders of magnitude in the matrix and 23 orders of magnitude in the RHS. HiGHS really should scream at you that it is going to get the wrong answer. This is a modeling issue, not a solver issue.

[sstroemer]

Hi @odow, sorry late reply because I'm travelling and don't have internet that often. I'll send you an email after my return in October! Please feel free to just close this here, since it may not be something to be resolved here.

Note on your comment about the ranges (and other previous arguments about modeling): These models occur in an automated iterative solve, with silenced models, so with no real control about the models themselves.

My takeaway: Since one cannot guarantee anything about these models, their "quality" should be manually checked before solving (ranges, kappa, etc.) and afterwards (violations, gaps, etc.) since there're "safety net" guarantees from arbitrary solvers.

[odow]

@jajhall this model is worth investigating.

Assuming that I've computed things correctly, the solution violates primal feasibility by a lot.

Also: how do you define the dual objective? What is the convention for the bounds etc?

import HiGHS_jll: libhighs

function libhighs_set_option(highs, option, value)
    typeP = Ref{Cint}()
    ret = @ccall libhighs.Highs_getOptionType(
        highs::Ptr{Cvoid},
        option::Ptr{Cchar},
        typeP::Ptr{Cint},
    )::Cint
    if !iszero(ret)
        # Not an option. Skip.
    elseif typeP[] == Cint(0)  # kHighsOptionTypeBool
        @ccall libhighs.Highs_setBoolOptionValue(
            highs::Ptr{Cvoid},
            option::Ptr{Cchar},
            parse(Cint, value)::Cint,
        )::Cint
    elseif typeP[] == Cint(1)  # kHighsOptionTypeInt
        @ccall libhighs.Highs_setIntOptionValue(
            highs::Ptr{Cvoid},
            option::Ptr{Cchar},
            parse(Cint, value)::Cint,
        )::Cint
    elseif typeP[] == Cint(2)  # kHighsOptionTypeDouble
        @ccall libhighs.Highs_setDoubleOptionValue(
            highs::Ptr{Cvoid},
            option::Ptr{Cchar},
            parse(Float64, value)::Cdouble,
        )::Cint
    else
        @assert typeP[] == Cint(3)  # kHighsOptionTypeString
        @ccall libhighs.Highs_setStringOptionValue(
            highs::Ptr{Cvoid},
            option::Ptr{Cchar},
            value::Ptr{Cchar},
        )::Cint
    end
    return
end

function inner_product(l, u, d)
    if isfinite(l) && !isfinite(u)
        return l * d
    elseif isfinite(u) && !isfinite(l)
        return u * d
    elseif !isfinite(l) && !isfinite(u)
        return 0.0
    elseif d >= 0
        return l * d
    else
        return u * d
    end
end

function libhighs_primal_feasibility(highs)
    m = @ccall libhighs.Highs_getNumRow(highs::Ptr{Cvoid})::Cint
    n = @ccall libhighs.Highs_getNumCol(highs::Ptr{Cvoid})::Cint
    nnz = @ccall libhighs.Highs_getNumNz(highs::Ptr{Cvoid})::Cint
    col_value, col_dual = zeros(n), zeros(n)
    row_value, row_dual = zeros(m), zeros(m)
    @ccall libhighs.Highs_getSolution(
        highs::Ptr{Cvoid},
        col_value::Ptr{Cdouble},
        col_dual::Ptr{Cdouble},
        row_value::Ptr{Cdouble},
        row_dual::Ptr{Cdouble},
    )::Cint
    a_format = Cint(2) # kHighsMatrixFormatRowwise
    q_format = Cint(1) # kHighsHessianFormatTriangular
    num_col, num_row, num_nz = Ref{Cint}(), Ref{Cint}(), Ref{Cint}()
    hessian_num_nz = Ref{Cint}()
    sense = Ref{Cint}()
    offset = Ref{Cdouble}()
    col_cost, col_lower, col_upper = zeros(n), zeros(n), zeros(n)
    row_lower, row_upper = zeros(m), zeros(m)
    a_start, a_index = zeros(Cint, m), zeros(Cint, nnz)
    a_value = zeros(Cdouble, nnz)
    integrality = zeros(Cint, n)
    q_start, q_index, q_value = C_NULL, C_NULL, C_NULL
    @ccall libhighs.Highs_getModel(
        highs::Ptr{Cvoid},
        a_format::Cint,
        q_format::Cint,
        num_col::Ptr{Cint},
        num_row::Ptr{Cint},
        num_nz::Ptr{Cint},
        hessian_num_nz::Ptr{Cint},
        sense::Ptr{Cint},
        offset::Ptr{Cdouble},
        col_cost::Ptr{Cdouble},
        col_lower::Ptr{Cdouble},
        col_upper::Ptr{Cdouble},
        row_lower::Ptr{Cdouble},
        row_upper::Ptr{Cdouble},
        a_start::Ptr{Cint},
        a_index::Ptr{Cint},
        a_value::Ptr{Cdouble},
        q_start::Ptr{Cint},
        q_index::Ptr{Cint},
        q_value::Ptr{Cdouble},
        integrality::Ptr{Cint},
    )::Cint
    primal_infeasibilities = Pair{Int,Float64}[]
    # Column bounds
    for i in 1:n
        li, xi, ui  = col_lower[i], col_value[i], col_upper[i]
        if xi < li
            push!(primal_infeasibilities, i => li - xi)
        end
        if xi > ui
            push!(primal_infeasibilities, i => xi - ui)
        end
    end
    # Row bounds
    for i in 1:m
        li, ui = row_lower[i], row_upper[i]
        lhs = 0.0
        for k in a_start[i]:(get(a_start, i + 1, nnz)-1)
            lhs += col_value[a_index[k+1]+1] * a_value[k+1]
        end
        if lhs < li
            push!(primal_infeasibilities, -i => li - lhs)
        end
        if lhs > ui
            push!(primal_infeasibilities, -i => lhs - ui)
        end
    end
    # Primal objective
    primal_objective = col_cost' * col_value + offset[]
    # Dual objective
    dual_objective = offset[]
    for i in 1:m
        dual_objective += inner_product(row_lower[i], row_upper[i], row_dual[i])
    end
    for i in 1:n
        dual_objective += inner_product(col_lower[i], col_upper[i], col_dual[i])
    end
    println("Manually computed feasibility information")
    println("  primal objective = ", primal_objective)
    println("  dual objective   = ", dual_objective)
    println("  primal infeasibilities")
    println("    count = ", length(primal_infeasibilities))
    println("    max   = ", maximum(last.(primal_infeasibilities)))
    println("    sum   = ", sum(last.(primal_infeasibilities)))
    println("HiGHS claims")
    obj = @ccall(libhighs.Highs_getObjectiveValue(highs::Ptr{Cvoid})::Cdouble)
    println(" primal objective = ", obj)
    pDouble = Ref{Cdouble}()
    @ccall libhighs.Highs_getDoubleInfoValue(
        highs::Ptr{Cvoid},
        "max_primal_infeasibility"::Ptr{Cchar},
        pDouble::Ptr{Cdouble},
    )::Cint
    println("  max_primal_infeasibility   = ", pDouble[])
    @ccall libhighs.Highs_getDoubleInfoValue(
        highs::Ptr{Cvoid},
        "sum_primal_infeasibilities"::Ptr{Cchar},
        pDouble::Ptr{Cdouble},
    )::Cint
    println("  sum_primal_infeasibilities = ", pDouble[])
    return primal_infeasibilities
end

filename = "benchmark/stroemer.mps.gz"
# filename = "test.mps"
highs = @ccall libhighs.Highs_create()::Ptr{Cvoid}
ret = @ccall libhighs.Highs_readModel(
    highs::Ptr{Cvoid},
    filename::Ptr{Cchar},
)::Cint
libhighs_set_option(highs, "solver", "ipm")
libhighs_set_option(highs, "run_crossover", "off")
@ccall libhighs.Highs_run(highs::Ptr{Cvoid})::Cint
primal_infeasibilities = libhighs_primal_feasibility(highs);

julia> include("benchmark/libhighs_utilities.jl");
Running HiGHS 1.7.2 (git hash: 5ce7a2753): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-09, 1e+07]
  Cost   [8e-02, 2e+00]
  Bound  [1e+03, 1e+06]
  RHS    [1e-13, 4e+11]
WARNING: Problem has excessively large bounds: consider scaling the bounds by 1e+2 or less, or setting option user_bound_scale to -6 or less
Presolving model
4071 rows, 179 cols, 23438 nonzeros  0s
868 rows, 179 cols, 7917 nonzeros  0s
846 rows, 157 cols, 7873 nonzeros  0s
Presolve : Reductions: rows 846(-3357); columns 157(-37); elements 7873(-15903)
Solving the presolved LP
IPX model has 846 rows, 157 columns and 7873 nonzeros
Input
    Number of variables:                                157
    Number of free variables:                           0
    Number of constraints:                              846
    Number of equality constraints:                     0
    Number of matrix entries:                           7873
    Matrix range:                                       [1e-09, 1e+07]
    RHS range:                                          [1e-13, 4e+11]
    Objective range:                                    [4e-02, 2e+00]
    Bounds range:                                       [1e+03, 1e+06]
Preprocessing
    Dualized model:                                     yes
    Number of dense columns:                            0
    Range of scaling factors:                           [4.88e-04, 5.12e+02]
IPX version 1.0
Interior Point Solve
 Iter     P.res    D.res            P.obj           D.obj        mu     Time
   0   1.61e+00 4.03e+08  -3.66929576e+06 -2.80164908e+05  7.23e+08       0s
 Constructing starting basis...
   1   7.55e-01 1.42e+08   2.11231148e+10 -5.88836059e+09  3.31e+08       0s
   2   3.84e-01 1.03e+08   2.52449514e+10 -6.77601358e+09  2.15e+08       0s
   3   1.62e-01 6.26e+07   2.45424246e+10 -7.58869432e+09  1.27e+08       0s
   4   6.54e-02 4.27e+07   1.78431446e+10 -7.17364101e+09  7.46e+07       0s
   5   1.24e-02 6.96e+06   9.76090911e+09 -4.31990229e+09  2.00e+07       0s
   6   6.89e-04 6.63e+05   1.58326057e+09 -7.13738837e+08  2.39e+06       0s
   7   7.76e-05 4.13e+04   2.46404688e+08 -5.04335915e+07  2.76e+05       0s
   8   9.91e-06 5.17e+03   3.06155218e+07 -8.00019998e+06  3.52e+04       0s
   9   1.59e-06 3.00e+02   6.51271708e+06 -3.11647157e+06  8.62e+03       0s
  10   6.94e-07 2.03e+02   4.42252902e+06 -2.61475272e+06  6.25e+03       0s
  11   1.95e-07 1.15e+01   1.16247901e+06 -9.33907161e+05  1.86e+03       0s
  12   7.09e-08 3.29e+00   2.33052140e+05 -6.30830702e+05  7.64e+02       0s
  13   1.10e-08 1.49e+00  -6.37331850e+04 -5.16162486e+05  3.98e+02       0s
  14   1.89e-09 6.45e-01  -1.87443027e+05 -3.91529128e+05  1.80e+02       0s
  15   5.11e-10 2.73e-01  -2.43048280e+05 -3.25901641e+05  7.29e+01       0s
  16   4.22e-10 6.53e-02  -2.45993831e+05 -3.10931335e+05  5.71e+01       0s
  17   2.36e-10 9.89e-03  -2.56473399e+05 -2.87552996e+05  2.73e+01       0s
  18   5.85e-11 2.38e-03  -2.70031405e+05 -2.81458762e+05  1.01e+01       0s
  19   2.62e-11 9.68e-04  -2.72995455e+05 -2.79622587e+05  5.83e+00       0s
  20   7.84e-12 3.31e-04  -2.75624088e+05 -2.77953994e+05  2.05e+00       0s
  21   3.97e-12 1.07e-04  -2.76191853e+05 -2.77430268e+05  1.09e+00       0s
  22   1.25e-12 4.19e-05  -2.76727983e+05 -2.77201765e+05  4.17e-01       0s
  23   3.88e-13 9.30e-06  -2.76925394e+05 -2.77086691e+05  1.42e-01       0s
  24   2.21e-13 3.34e-06  -2.76971166e+05 -2.77066709e+05  8.40e-02       0s
  25   3.74e-14 6.56e-07  -2.77032496e+05 -2.77051948e+05  1.71e-02       0s
  26   1.47e-14 3.58e-07  -2.77041263e+05 -2.77049236e+05  7.01e-03       0s
  27   6.79e-15 2.38e-07  -2.77044889e+05 -2.77048574e+05  3.24e-03       0s
  28   8.88e-16 1.79e-07  -2.77047774e+05 -2.77048240e+05  4.09e-04       0s
  29   5.45e-16 2.38e-07  -2.77048099e+05 -2.77048191e+05  8.16e-05       0s
  30   5.26e-16 1.79e-07  -2.77048185e+05 -2.77048188e+05  2.61e-06       0s
  31*  5.18e-16 2.38e-07  -2.77048188e+05 -2.77048188e+05  7.97e-08       0s
Summary
    Runtime:                                            0.04s
    Status interior point solve:                        optimal
    Status crossover:                                   not run
    objective value:                                    2.77048188e+05
    interior solution primal residual (abs/rel):        2.44e-04 / 6.89e-16
    interior solution dual residual (abs/rel):          2.48e-09 / 9.88e-10
    interior solution objective gap (abs/rel):          1.62e-04 / 5.84e-10
Ipx: IPM       optimal
Model   status      : Optimal
IPM       iterations: 31
Objective value     :  2.7704818743e+05
HiGHS run time      :          0.06
Manually computed feasibility information
  primal objective = 277048.1874263953
  dual objective   = -231222.881764544
  primal infeasibilities
    count = 63
    max   = 68.27635192871094
    sum   = 418.93546406958944
HiGHS claims
 primal objective = 277048.1874263954
  max_primal_infeasibility   = 0.0
  sum_primal_infeasibilities = 0.0

[jajhall]

I need the model to investigate this

[odow]

Sorry it was linked at the top but not super clearly https://gist.github.com/sstroemer/07406d9422993aab2ee4d2fbd3ee8a8c

[jajhall]

I note that

• Gurobi gives an optimal objective value of 2.770481876e+05, warning: of unscaled primal violation = 2.04475e-05 and residual = 0.00390625
• HiGHS gives an optimal objective value of 2.7704818743e+05

I'll keep looking at why you're getting the anomalies that you observe

[odow]

The problem isn't the objective value. It's the feasibility of the actual primal (and dual) solution.

[jajhall]

The problem isn't the objective value. It's the feasibility of the actual primal (and dual) solution.

Sure, but to get the objective function so close to that of Gurobi and get such large infeasibilities looks odd.

[odow]

If I run with crossover, I get

Manually computed feasibility information
  primal objective = 277048.1875806275
  dual objective   = 277048.1875806275
  primal infeasibilities
    count = 200
    max   = 6.103515625e-5
    sum   = 0.00020435247858031306
HiGHS claims
 primal objective = 277048.18758062756
  max_primal_infeasibility   = 0.0
  sum_primal_infeasibilities = 0.0

That looks much better.

[odow]

I think the primal feasibilities are just because of the large relative magnitude of the LHS:
abs(lhs) = 1.4074251393645157e11. So we get an absolute primal infeasibility of 1e2, but that's small relatively.

[jajhall]

Thanks, it looks as if crossover is moving on the facet of points with objective value 277048.1875806275

I think the primal feasibilities are just because of the large relative magnitude of the LHS:
abs(lhs) = 1.4074251393645157e11. So we get an absolute primal infeasibility of 1e2, but that's small relatively.

I was looking at the terms in the row activities that you're computing in cases where you find a large infeasibility, and they don't get bigger than 1e12, so there shouldn't be problems adding them up.

I have to go to ZIB now, but will carry on looking at this later with VScode

[odow]

Perhaps the dual objective is just an issue because of how we map the duals to the bound.

If I have a double-sided constraint, then we either:

• look at the basis to see which side was binding
• use a heuristic of "if >=0, the lower side is binding, otherwise the top side is"

The heuristic can fail if the dual is slightly dual infeasible at say, `-1e-9, but we should still multiply by the lower bound.

Perhaps we also need to instead look at the row value.

[jajhall]

If IPM is run without crossover then there's no concept of "binding". The same is true for PDLP, which is why HighsInfo now has the following (with values for this model)

Max complementarity violation [type: double]
max_complementarity_violation = 1.340212873e-05

Sum of complementarity violations [type: double]
sum_complementarity_violations = 0.0003667661438

[odow]

So how should I compute the dual objective if I have two finite bounds and a single dual for the row?

[jajhall]

When I'm checking for KKT errors without a basis, I look to see whether the primal value is above or below the midpoint of [l, u], and choose that bound, so I guess you should do that when you compute the dual objective.

[odow]

@sstroemer I've fixed the dual objective value issue:

[odow]

I've also fixed the ObjectiveBound and RelativeGap to report the primal-dual gap if the model is an LP in #229.

Julian is also working on fixing the underlying imprecision in the IPM solver, which relates to how HiGHS refactors Ax >= b into Ax - s = 0; s >= b. Because of the numerics of the problem, errors in Ax - s = 0 propagated into the feasibility tolerance that weren't being checked for.

I'm closing this issue because I don't think there is anything left to do in HiGHS.jl. I'll keep the model around in the open-energy-modeling-benchmarks repo.

Edit: oh, and I've made a note about JuMP-level solution quality checker in jump-dev/JuMP.jl#3664

[jajhall]

Primal and dual residual errors are now checked and corrected when IPM or PDLP claim to solve LPs to optimality.

Running sstroemer/mwe.mps with IPM and without presolve or crossover, optimality is no longer claimed.

Fix is in ERGO-Code/HiGHS#1978