HiGHS: Bug in presolve leads to NaN in solution

This is probably user error, but I’d like to know what I’m doing wrong/how the problem is being set up differently.

In Julia, I get

julia> using JuMP, HiGHS

julia> function run_model(verbose = true)
           model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => verbose))
           @variable(model, v_0 >= -1_000_000_000, Int)
           @variable(model, v_1 >= -1_000_000_000, Int)
           @variable(model, v_2 >= -1_000_000_000, Int)
           @variable(model, v_3 >= -1_000_000_000, Int)
           @variable(model, v_4 >= -1_000_000_000, Int)
           @variable(model, v_5 >= 0, Int)
           @variable(model, v_6 >= 0, Int)
           @variable(model, v_7 >= 0, Int)
           @objective(model, Max, 2v_0 + v_1 - 2v_2 - v_3 - 2v_6 - v_7)
           @constraint(model, -v_3 - v_4 - v_5 - 2v_6 - 3v_7 <= 0)
           @constraint(model, 2v_0 + v_1 - 2v_2 - v_3 - 2v_6 - v_7 <= 1)
           @constraint(model, v_0 - v_2 - v_6 <= 0)
           @constraint(model, v_1 - v_3 - v_7 <= 0)
           optimize!(model)
       end
run_model (generic function with 2 methods)

julia> m = run_model()
Presolving model
3 rows, 6 cols, 12 nonzeros
1 rows, 2 cols, 2 nonzeros
0 rows, 1 cols, 0 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal

Solving report
  Status            Optimal
  Primal bound      -0
  Dual bound        -0
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (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)

It was able to solve the objective, achieving a value of 0. The point of this problem was to check if one of the inequalities in a system was redundant

-v_5 <= 0
-v_6 <= 0
-v_7 <= 0
-v_3 - v_4 - v_5 - 2v_6 - 3v_7 <= 0
v_0 - v_2 - v_6 <= 0
2v_0 + v_1 - 2v_2 - v_3 - 2v_6 - v_7 <= 0
v_1 - v_3 - v_7 <= 0

I.e., is the second to last constraint redundant? To check, I increment the upper bound on the constraint, and then try to maximize it as an objective. If it is able to achieve the incremented upper bound, then the original constraint was necessary. If not, then it is redundant and can be dropped.

A secondary question I have – is there a better way to do this? Particularly from the C++ API?

Note that the redundant constraint here equals 2times the constraint above it plus the constraint below it, but my problems often are not so trivial.

But failing that, this C++ code does not work:

#include <cassert>

#include "Highs.h"

using std::cout;
using std::endl;

int main() {
  // Create and populate a HighsModel instance for the LP
  //
  // Maximize 2v_0 + v_1 - 2v_2 - v_3 - 2v_6 - v_7
  //
  // -v_5 <= 0
  // -v_6 <= 0
  // -v_7 <= 0
  // -v_3 - v_4 - v_5 - 2v_6 - 3v_7 <= 0
  // v_0 - v_2 - v_6 <= 0
  // 2v_0 + v_1 - 2v_2 - v_3 - 2v_6 - v_7 <= 1
  // v_1 - v_3 - v_7 <= 0
  //
  //
  HighsModel model;

  model.lp_.num_col_ = 8;
  model.lp_.num_row_ = 4;
  model.lp_.sense_ = ObjSense::kMaximize;
  model.lp_.offset_ = 0;
  model.lp_.col_cost_ = {2, 1, -2, -1, 0, 0, -2, -1};
  model.lp_.col_lower_ = {-1e+30, -1e+30, -1e+30, -1e+30, -1e+30, 0, 0, 0};
  model.lp_.col_upper_ = {1e+30, 1e+30, 1e+30, 1e+30,
                          1e+30, 1e+30, 1e+30, 1e+30};
  model.lp_.row_lower_ = {-1e+30, -1e+30, -1e+30, -1e+30};
  model.lp_.row_upper_ = {0, 0, 1, 0};
  //
  // Here the orientation of the matrix is row-wise
  model.lp_.a_matrix_.format_ = MatrixFormat::kRowwise;
  // a_start_ has num_col_1 entries, and the last entry is the number
  // of nonzeros in A, allowing the number of nonzeros in the last
  // column to be defined
  model.lp_.a_matrix_.start_ = {0, 5, 8, 14, 17};
  model.lp_.a_matrix_.index_ = {3, 4, 5, 6, 7, 0, 2, 6, 0,
                                1, 2, 3, 6, 7, 1, 3, 7};
  model.lp_.a_matrix_.value_ = {-1, -1, -1, -2, -3, 1, -1, -1, 2,
                                1,  -2, -1, -2, -1, 1, -1, -1};
  model.lp_.integrality_.resize(model.lp_.num_col_, HighsVarType::kInteger);
  //
  // Create a Highs instance
  Highs highs;
  HighsStatus return_status;
  //
  // Pass the model to HiGHS
  return_status = highs.passModel(model);
  assert(return_status == HighsStatus::kOk);
  //
  // Get a const reference to the LP data in HiGHS
  const HighsLp& lp = highs.getLp();
  //
  // Solve the model
  return_status = highs.run();
  assert(return_status == HighsStatus::kOk);
  //
  // Get the model status
  const HighsModelStatus& model_status = highs.getModelStatus();
  assert(model_status == HighsModelStatus::kOptimal);
  cout << "Model status: " << highs.modelStatusToString(model_status) << endl;
  //
  // Get the solution information
  const HighsInfo& info = highs.getInfo();
  cout << "Simplex iteration count: " << info.simplex_iteration_count << endl;
  cout << "Objective function value: " << info.objective_function_value << endl;
  cout << "Primal  solution status: "
       << highs.solutionStatusToString(info.primal_solution_status) << endl;
  cout << "Dual    solution status: "
       << highs.solutionStatusToString(info.dual_solution_status) << endl;
  cout << "Basis: " << highs.basisValidityToString(info.basis_validity) << endl;
  const bool has_values = info.primal_solution_status;
  const bool has_duals = info.dual_solution_status;
  const bool has_basis = info.basis_validity;
  //
  // Get the solution values and basis
  const HighsSolution& solution = highs.getSolution();
  const HighsBasis& basis = highs.getBasis();
  //
  // Report the primal and solution values and basis
  for (int col = 0; col < lp.num_col_; col++) {
    cout << "Column " << col;
    if (has_values) cout << "; value = " << solution.col_value[col];
    if (has_duals) cout << "; dual = " << solution.col_dual[col];
    if (has_basis)
      cout << "; status: " << highs.basisStatusToString(basis.col_status[col]);
    cout << endl;
  }
  for (int row = 0; row < lp.num_row_; row++) {
    cout << "Row    " << row;
    if (has_values) cout << "; value = " << solution.row_value[row];
    if (has_duals) cout << "; dual = " << solution.row_dual[row];
    if (has_basis)
      cout << "; status: " << highs.basisStatusToString(basis.row_status[row]);
    cout << endl;
  }

  return 0;
}

This is modified from the call_highs_from_cpp example.

When I run this, I get:

Running HiGHS 1.2.2 [date: 2022-05-25, git hash: 6a9d90cf]
Copyright (c) 2022 ERGO-Code under MIT licence terms
Cols:           5 lower bounds exceeding       -1e+20 are treated as -Infinity
Cols:           8 upper bounds exceeding        1e+20 are treated as +Infinity
Rows:           4 lower bounds exceeding       -1e+20 are treated as -Infinity
Presolving model
3 rows, 6 cols, 12 nonzeros
1 rows, 2 cols, 2 nonzeros
0 rows, 1 cols, 0 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal

Solving report
  Status            Optimal
  Primal bound      -0
  Dual bound        -0
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    -nan (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (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 status: Optimal
Simplex iteration count: -1
Objective function value: -nan
Primal  solution status: Feasible
Dual    solution status: None
Basis: Not valid
Column 0; value = -inf
Column 1; value = -inf
Column 2; value = -nan
Column 3; value = -nan
Column 4; value = -inf
Column 5; value = 0
Column 6; value = 0
Column 7; value = 0
Row    0; value = -nan
Row    1; value = -nan
Row    2; value = -nan
Row    3; value = -nan