Bug repeatedly resolving an LP with rows deleted in between

[odow]

By the power of git bisect, this commit broke SDDP.jl: b03608c

commit b03608c811f4ea329b60550df86477d2745ca3cf
Author: JAJHall <jajhall@ed.ac.uk>
Date:   Mon Jan 13 17:50:34 2025 +0000

    All unit tests pass

 src/lp_data/HStruct.h         | 14 ++++++++++++--
 src/lp_data/Highs.cpp         |  9 ++-------
 src/lp_data/HighsSolution.cpp | 33 ++++++++++++++++++++-------------
 src/simplex/HApp.h            | 12 ++++++++++++
 4 files changed, 46 insertions(+), 22 deletions(-)

It is not present in v1.9.0, but it is present in v1.10.0.

I'm still chasing a minimal reproducible example. But for now... we have this horrific thing which solves a multistage stochastic program:

using SDDP, HiGHS, Random
function main(; seed)
    Random.seed!(seed)
    S = [0 1 2; 0 0 1; 0 0 0]
    days = [60, 60, 245]
    D = [210, 210, 858]
    q = [4.5, 5.5, 6.5]
    b = [30 75 37.5; 15 37.5 18.25; 7.5 18.75 9.325]
    r = [5, 6, 7]
    graph = SDDP.MarkovianGraph([
        ones(Float64, 1, 1),
        [0.14 0.69 0.17],
        [0.14 0.69 0.17; 0.14 0.69 0.17; 0.14 0.69 0.17],
        [0.14 0.69 0.17; 0.14 0.69 0.17; 0.14 0.69 0.17],
    ])
    model = SDDP.PolicyGraph(
        graph;
        lower_bound = 0.0,
        direct_mode = true,
        optimizer = HiGHS.Optimizer,
    ) do sp, (t, j)
        @variable(sp, 0 <= x_a <= 60, SDDP.State, initial_value = 60)
        @variable(sp, x_b[1:3] >= 0, SDDP.State, initial_value = 0)
        @variable(sp, u[1:5,1:3] >= 0)
        @constraint(sp, x_a.out <= x_a.in)
        if t == 1
            @constraint(sp, [i in 1:3], x_b[i].in == x_b[i].out)
        else
            @expression(sp, e[c in 1:3], x_b[c].in + u[1,c] - u[2,c] - u[3,c] + u[4,c] - u[5,c])
            @constraints(sp, begin
                sum(u[3, :]) >= D[t-1]
                x_b[t-1].out == e[t-1] + x_a.in * b[t-1,j]
                [c = 1:3; c != t - 1], x_b[c].out == e[c]
                sum(x_b[i].out for i in 1:3) <= 3000
                [c in 1:3], u[2, c] <= x_b[c].in
                sum(u[2, :]) <= 3000
            end)
            @stageobjective(
                sp,
                1000 * sum(u[4:5, :]) + 50 * x_a.in +
                sum(0.05 * x_b[i].in * days[t-1] + r[i] * u[1, i] + S[i, t-1] * u[3, i] - q[i] * u[2, i] for i in 1:3)
            )
        end
    end
    SDDP.train(model)
    @assert isapprox(SDDP.calculate_bound(model), 4074.1391; atol = 1e-4)
    return
end

# for seed in 1:100
#     @show seed
#     main(; seed)
# end

main(; seed = 3)

The log file is huge, but here it is:

highs_issue.log

I'll try to find the offending subproblem tomorrow.

[odow]

Here's a better reproducer. I can make this C if needed. Now I know the sequence, I'm trying to write a simpler model from scratch.

julia> highs = Highs_create()
Ptr{Nothing} @0x00007f8e29a3a800

julia> begin
           Highs_readModel(highs, "model2.mps")
           ret = Highs_run(highs)
           Highs_deleteRowsByRange(highs, 10, 10)
           Highs_run(highs)
           num_row = Highs_getNumRow(highs)
           num_col = Highs_getNumCol(highs)
           col_value, col_dual = zeros(num_col), zeros(num_col)
           row_value, row_dual = zeros(num_row), zeros(num_row)
           ret = Highs_getSolution(highs, col_value, col_dual, row_value, row_dual)
       end
Running HiGHS 1.10.0 (git hash: fd8665394e): Copyright (c) 2025 HiGHS under MIT licence terms
LP   model2 has 11 rows; 24 cols; 41 nonzeros
Coefficient ranges:
  Matrix [1e+00, 6e+01]
  Cost   [1e+00, 1e+03]
  Bound  [2e+01, 1e+03]
  RHS    [2e+02, 6e+03]
Presolving model
7 rows, 11 cols, 18 nonzeros  0s
Dependent equations search running on 1 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
1 rows, 4 cols, 4 nonzeros  0s
Presolve : Reductions: rows 1(-10); columns 4(-20); elements 4(-37)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
          0     4.1103839839e+03 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : model2
Model status        : Optimal
Objective value     :  4.1103839839e+03
Relative P-D gap    :  2.2126757630e-16
HiGHS run time      :          0.00
LP   model2 has 10 rows; 24 cols; 39 nonzeros
Coefficient ranges:
  Matrix [1e+00, 2e+01]
  Cost   [1e+00, 1e+03]
  Bound  [2e+01, 1e+03]
  RHS    [2e+02, 3e+03]
Solving LP without presolve, or with basis, or unconstrained
Model name          : model2
Model status        : Optimal
Objective value     : -8.2206159000e+02
Relative P-D gap    :  7.7458246212e+00
HiGHS run time      :          0.00
0

julia> begin
           costs, lower, upper = zeros(num_col), zeros(num_col), zeros(num_col)
           Highs_getColsByRange(
               highs,
               0, # from_col,
               num_col - 1, # to_col,
               Ref{Cint}(0), # num_col,
               costs,
               lower,
               upper,
               Ref{Cint}(0), # num_nz,
               C_NULL, # matrix_start,
               C_NULL, # matrix_index,
               C_NULL, # matrix_value,
           )
       end
0

julia> col_value
24-element Vector{Float64}:
    2.27412318
    2.27412318
 1154.473908
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
 1154.473908
    0.0
    0.0
    0.0
    0.0
   -0.0
  272.8947816
    0.0
    0.0
   -0.0
    0.0
    0.0
    0.0
    0.0
    0.0

julia> lower
24-element Vector{Float64}:
   18.19298544
    0.0
 1154.473908
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0

julia> col_value .- lower
24-element Vector{Float64}:
  -15.918862260000001
    2.27412318
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
 1154.473908
    0.0
    0.0
    0.0
    0.0
   -0.0
  272.8947816
    0.0
    0.0
   -0.0
    0.0
    0.0
    0.0
    0.0
    0.0

Note the bug: For x[0], the col_value of 2.2741231800314705 does not satisfy the col_lower of 18.192985440251764

There is also a clue in the log: Relative P-D gap : 7.7458246212e+00

MPS starting point

NAME        model.fail
ROWS
 N  Obj     
 L  r0      
 G  r1      
 E  r2      
 E  r3      
 E  r4      
 L  r5      
 L  r6      
 L  r7      
 L  r8      
 L  r9      
 G  r10     
COLUMNS
    c0        Obj       50
    c0        r0        -1
    c0        r2        -15
    c1        r0        1
    c1        r10       59.00925
    c2        Obj       3
    c2        r3        -1
    c2        r6        -1
    c3        r3        1
    c3        r5        1
    c4        Obj       3
    c4        r2        -1
    c4        r7        -1
    c5        r2        1
    c5        r5        1
    c6        Obj       3
    c6        r4        -1
    c6        r8        -1
    c7        r4        1
    c7        r5        1
    c8        Obj       5
    c8        r3        -1
    c9        Obj       -4.5
    c9        r3        1
    c9        r6        1
    c9        r9        1
    c10       Obj       1
    c10       r1        1
    c10       r3        1
    c11       Obj       1000
    c11       r3        -1
    c12       Obj       1000
    c12       r3        1
    c13       Obj       6
    c13       r2        -1
    c14       Obj       -5.5
    c14       r2        1
    c14       r7        1
    c14       r9        1
    c15       r1        1
    c15       r2        1
    c16       Obj       1000
    c16       r2        -1
    c17       Obj       1000
    c17       r2        1
    c18       Obj       7
    c18       r4        -1
    c19       Obj       -6.5
    c19       r4        1
    c19       r8        1
    c19       r9        1
    c20       r1        1
    c20       r4        1
    c21       Obj       1000
    c21       r4        -1
    c22       Obj       1000
    c22       r4        1
    c23       Obj       1
    c23       r10       1
RHS
    RHS_V     r1        210
    RHS_V     r5        3000
    RHS_V     r9        3000
    RHS_V     r10       6006
BOUNDS
 FX BOUND     c0        18.19298544
 UP BOUND     c1        60
 FX BOUND     c2        1154.473908
 FX BOUND     c4        0
 FX BOUND     c6        0
ENDATA

Debugging process

Writing this down for my future reference:

• Start with Add a macro to trace ccall for debugging jump-dev/HiGHS.jl#258 for both the working and failing builds of HiGHS

diff --git a/src/gen/libhighs.jl b/src/gen/libhighs.jl
index 384d5ba..075a42e 100644
--- a/src/gen/libhighs.jl
+++ b/src/gen/libhighs.jl
@@ -43,12 +43,16 @@ macro trace(def)
             $fn_body
         else
             open(log_file, "a") do $io
+                println($io, "# === INPUTS === ")
                 $print_args
                 println($io, "ret = ", $call)
             end
             result = (() -> $fn_body)()
             open(log_file, "a") do $io
-                println($io, "# ret = ", $_trace_format_value(result), "\n")
+                println($io, "# === OUTPUTS === ")
+                println($io, "# ret = ", $_trace_format_value(result))
+                $print_args
+                println($io)
             end
             result
         end

• Chunk the log into separate function calls
• Split the logs into separate files for each subproblem that was solved
• Compare the outputs to find a model and function call that differ
• Drop all unnecessary function calls
• Concatenate into a runnable script
• Manually simplify

function get_model(fn_call)
    regex = if occursin("Highs_create", fn_call)
        r"ret = Ptr\{Nothing\} (\@0x[0-9a-f]+)"
    else
        r"highs = Ptr\{Nothing\} (\@0x[0-9a-f]+)"
    end
    m = match(regex, fn_call)
    if m === nothing
        @show fn_call
        error()
    end
    return String(m[1])
end

function split_into_models(filename)
    ordered_models = String[]
    models = Dict{String,Vector{String}}()
    for fn_call in function_calls(filename)
        if occursin("Highs_versionMajor", fn_call)
            continue
        elseif occursin("Highs_versionMinor", fn_call)
            continue
        end
        model = get_model(fn_call)
        if !haskey(models, model)
            models[model] = String[]
            push!(ordered_models, model)
        end
        push!(models[model], fn_call)
    end
    return ordered_models, models
end
ordered_models_pass, models_pass = split_into_models("/tmp/highs_pass.log");
ordered_models_fail, models_fail = split_into_models("/tmp/highs_fail.log");

for (m_pass, m_fail) in zip(ordered_models_pass, ordered_models_fail)
    @warn "m_pass, m_fail = $m_pass, $m_fail"
    fn_pass = models_pass[m_pass]
    fn_fail = models_fail[m_fail]
    for (i, (pass, fail)) in enumerate(zip(fn_pass, fn_fail))
        if preprocess(pass) != preprocess(fail)
            @warn "i = $i"
            println()
            @info "PASS"
            println(pass)
            println()
            @info "FAIL"
            println(fail)
            error()
        end
    end
end

function filter_rows(fn_calls)
    filter!(fn_calls) do fn_call
        return !any([
            "Highs_getNumCols",
            "Highs_getNumRows",
            "Highs_getOptionType",
            "Highs_zeroAllClocks",
            "Highs_getObjectiveSense",
            "Highs_getObjectiveOffset",
            "Highs_getRowsByRange",
            "Highs_setCallback",
            "Highs_startCallback",
            "Highs_stopCallback",
            "Highs_getModelStatus",
            "Highs_setBoolOptionValue",
            "Highs_getIntInfoValue",
            "Highs_getBasis",
            "Highs_getObjectiveValue",
        ]) do needle
            return occursin(needle, fn_call)
        end
    end
    map(fn_calls) do fn_call
        return split(fn_call, "# === OUTPUTS ==="; keepempty = false)[1]
    end
end


m_pass, m_fail = "@0x00007fc6dff77e00", "@0x00007fdf8c9a1e00"
fn_pass = filter_rows(models_pass[m_pass][1:352])
fn_fail = filter_rows(models_fail[m_fail][1:352])
write("/tmp/model.pass.log", join(fn_pass, "\n\n"))
write("/tmp/model.fail.log", join(fn_fail, "\n\n"))

[odow]

Okaaay! I have a 3 variable, 2 constraint model that gets the wrong answer:

HiGHS v1.10

julia> highs = Highs_create()
Ptr{Nothing} @0x00007f8e29b19600

julia> Highs_addCol(highs, 0.0, 2.0, 2.0, 0, C_NULL, C_NULL)
Running HiGHS 1.10.0 (git hash: fd8665394e): Copyright (c) 2025 HiGHS under MIT licence terms
0

julia> Highs_addCol(highs, 0.0, -Inf, Inf, 0, C_NULL, C_NULL)
0

julia> Highs_addCol(highs, 0.0, -Inf, Inf, 0, C_NULL, C_NULL)
0

julia> Highs_addRow(highs, 0.0, 0.0, 2, Cint[0, 2], [2.0, -1.0])
0

julia> Highs_addRow(highs, 10.0, Inf, 1, Cint[1], [6.0])
0

julia> Highs_run(highs)
LP   has 2 rows; 3 cols; 3 nonzeros
Coefficient ranges:
  Matrix [1e+00, 6e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [1e+01, 1e+01]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-2); columns 0(-3); elements 0(-3) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
0

julia> Highs_deleteRowsByRange(highs, 1, 1)
0

julia> Highs_run(highs)
LP   has 1 rows; 3 cols; 2 nonzeros
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [0e+00, 0e+00]
Solving LP without presolve, or with basis, or unconstrained
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
0

julia> col_value = zeros(3);

julia> Highs_getSolution(highs, col_value, C_NULL, C_NULL, C_NULL)
0

julia> col_value[1]  # Should be 2.0
1.0

HiGHS v1.9

julia> highs = Highs_create()
Ptr{Nothing} @0x00007fc5311aca00

julia> Highs_addCol(highs, 0.0, 2.0, 2.0, 0, C_NULL, C_NULL)
Running HiGHS 1.9.0 (git hash: 66f735e60): Copyright (c) 2024 HiGHS under MIT licence terms
0

julia> Highs_addCol(highs, 0.0, -Inf, Inf, 0, C_NULL, C_NULL)
0

julia> Highs_addCol(highs, 0.0, -Inf, Inf, 0, C_NULL, C_NULL)
0

julia> Highs_addRow(highs, 0.0, 0.0, 2, Cint[0, 2], [2.0, -1.0])
0

julia> Highs_addRow(highs, 10.0, Inf, 1, Cint[1], [6.0])
0

julia> Highs_run(highs)
Coefficient ranges:
  Matrix [1e+00, 6e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [1e+01, 1e+01]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-2); columns 0(-3); elements 0(-3) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
0

julia> Highs_deleteRowsByRange(highs, 1, 1)
0

julia> Highs_run(highs)
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [0e+00, 0e+00]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-1); columns 0(-3); elements 0(-2) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
0

julia> col_value = zeros(3);

julia> Highs_getSolution(highs, col_value, C_NULL, C_NULL, C_NULL)
0

julia> col_value[1]  # Should be 2.0
2.0

[odow]

Here's a C test

void test_delete_row_resolve_with_basis() {
    void* highs = Highs_create();
    int ret;
    double INF = Highs_getInfinity(highs);
    ret = Highs_addCol(highs, 0.0, 2.0, 2.0, 0, NULL, NULL);
    ret = Highs_addCol(highs, 0.0, -INF, INF, 0, NULL, NULL);
    ret = Highs_addCol(highs, 0.0, -INF, INF, 0, NULL, NULL);
    int index_1[2] = {0, 2};
    double value_1[2] = {2.0, -1.0};
    ret = Highs_addRow(highs, 0.0, 0.0, 2, index_1, value_1);
    int index_2[1] = {1};
    double value_2[1] = {6.0};
    ret = Highs_addRow(highs, 10.0, INF, 1, index_2, value_2);
    Highs_run(highs);
    double col_value[3] = {0.0, 0.0, 0.0};
    Highs_getSolution(highs, col_value, NULL, NULL, NULL);
    assertDoubleValuesEqual("col_value[0]", col_value[0], 2.0);
    ret = Highs_deleteRowsByRange(highs, 1, 1);
    assert(ret == 0);
    ret = Highs_run(highs);
    assert(ret == 0);
    ret = Highs_getSolution(highs, col_value, NULL, NULL, NULL);
    assert(ret == 0);
    assertDoubleValuesEqual("col_value[0]", col_value[0], 2.0);
    Highs_destroy(highs);
}

v1.10.0

Running HiGHS 1.10.0 (git hash: fa215aaa7): Copyright (c) 2025 HiGHS under MIT licence terms
LP   has 2 rows; 3 cols; 3 nonzeros
Coefficient ranges:
  Matrix [1e+00, 6e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [1e+01, 1e+01]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-2); columns 0(-3); elements 0(-3) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
LP   has 1 rows; 3 cols; 2 nonzeros
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [0e+00, 0e+00]
Solving LP without presolve, or with basis, or unconstrained
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Value col_value[0] = 1 differs from 2 by 1 but should be equal
Assertion failed: (1 == 0), function assertDoubleValuesEqual, file main.c, line 16.
./main.sh: line 15: 35239 Abort trap: 6           ./main

v1.9.0

Running HiGHS 1.9.0 (git hash: 66f735e60): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 6e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [1e+01, 1e+01]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-2); columns 0(-3); elements 0(-3) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [0e+00, 0e+00]
  Bound  [2e+00, 2e+00]
  RHS    [0e+00, 0e+00]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-1); columns 0(-3); elements 0(-2) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00

[jajhall]

Thanks

[jajhall]

Closed by #2268

[odow]

Thanks!