HiGHS
HiGHS copied to clipboard
ERGO-Code
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
Your models are not equivalent because the lower bounds are different. In JuMP, if you want a large bound, just omit it. In C++, use highs.getInfinity() for the value that HiGHS treats as infinity.
I also don't know if manually poking the fields in .lp is the way to go (even if it's documented as the example). The NaN might come from the fact that something wasn't set right.
You can see how the C API calls the C++ API here: https://github.com/ERGO-Code/HiGHS/blob/abaf8907a25ed40fcd2b33eb26d22b3c79ef5dbc/src/interfaces/highs_c_api.cpp#L82-L113 I think that's probably what you need.
p.s. what is this for???
Okay, here is an updated script:
#include <cassert>
#include <vector>
#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
//
//
// Create a Highs instance
Highs highs;
const HighsInt num_col = 8;
const HighsInt num_row = 4;
const HighsInt sense = HighsInt(ObjSense::kMaximize);
const HighsInt offset = 0;
std::vector<double> col_cost = {2, 1, -2, -1, 0, 0, -2, -1};
const auto Inf = highs.getInfinity();
std::vector<double> col_lower = {-Inf, -Inf, -Inf, -Inf, -Inf, 0, 0, 0};
std::vector<double> col_upper = {Inf, Inf, Inf, Inf,
Inf, Inf, Inf, Inf};
std::vector<double> row_lower = {-Inf, -Inf, -Inf, -Inf};
std::vector<double> row_upper = {0, 0, 1, 0};
const HighsInt a_format = HighsInt(MatrixFormat::kRowwise);
std::vector<HighsInt> a_start = {0, 5, 8, 14, 17};
std::vector<HighsInt> a_index = {3, 4, 5, 6, 7, 0, 2, 6, 0,
1, 2, 3, 6, 7, 1, 3, 7};
std::vector<double> a_value = {-1, -1, -1, -2, -3, 1, -1, -1, 2,
1, -2, -1, -2, -1, 1, -1, -1};
std::vector<HighsInt> integrality(num_col, HighsInt(HighsVarType::kInteger));
//
// Pass the model to HiGHS
HighsStatus return_status = highs.passModel(
num_col, num_row, a_value.size(), a_format, sense, offset,
col_cost.data(), col_lower.data(), col_upper.data(), row_lower.data(),
row_upper.data(), a_start.data(), a_index.data(), a_value.data(),
integrality.data());;
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;
}
Output still NaN:
Running HiGHS 1.2.2 [date: 2022-05-25, git hash: 6a9d90cf]
Copyright (c) 2022 ERGO-Code under MIT licence terms
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
p.s. what is this for???
I need/want to remove redundant constraints while performing Fourier-Motzkin elimination. This can prevent the exponential explosion in number of constraints as we reduce the dimension.
I wrote an algorithm that does this and can handle symbolic bounds, but I also need to run it with integer bounds. I'm not 100% confident in the algorithm's results -- I've fixed a ton of bugs and the answers are starting to look correct -- and I wanted something to validate it against. Plus, because I need to eliminate a lot of integer constraints in the actual program, it'd be worth benchmarking Highs on those problems so that I can use it instead when it's faster (and stick with my algorithm only when I need symbolic bound support).
My program will eventually need an ILP solver anyway, so using it there wouldn't be an extra dependency.
I don't have any background in integer optimization, so odds are a lot of things could be done better / I am open to suggestions. So in reference to the XY_problem/wanting to avoid doing things sub-optimally, a more high level overview:
I am modeling loop nests (i.e. for (size_t i = ...)).
Loop bounds themselves are constraints; when reordering affine loops (e.g. swapping the inner most and outer most loops in a nest), I use Fourier-Motzkin to eliminate all loops interior to a given loop to get its bounds.
Aside from avoiding the exponential explosion in the number of constraints, we also want to remove any redundant constraints to make the loop as simple as possible for the sake of generating better code.
Additionally, dependencies between accesses to memory can be represented as systems of constraints ("polyehdra"). If one of these dependence polyhedra is empty, then there is no dependency. I.e., there are no combinations of loop iterations between the two memory accesses where they interact with the same address.
I use Fourier-Motzkin to keep eliminating variables until there are none left (non-empty), or a constraint proving it is empty (e.g. 1 <= 0).
Assuming non-empty dependencies, we are not allowed to reorder loop iterations that violate these dependencies.
That is, if there is a dependence between memory access x and memory access y, and x happened before y in the original program, we cannot reorder loop iterations in a manner that lets y happen before x.
So we need to place constraints on legal schedules.
Farkas' lemma allows us to turn dependency polyehdra into constraints on affine schedules, but adds a ton of farkas multipliers that we can eliminate with Fourier-Motzkin.
I'm planning on first trying to code ad-hoc adjustments to legalize an optimal schedule, but if this fails for a particular code being optimized, I'll need to run an ILP optimizer to find a set of legal transformations. This requires another set of constraints that must also be produced via Farkas + Fourier-Motzkin.
The ILP solver must be run many times:
- for each loop level in the nest...
- we add a constraint to the existing constraints to guarantee the new solution is not a linear combination of older solutions
- we minimize multiple objectives to minimize their lexicographical order. This means minimizing the first, adding constraints to pin this objective's value, solve the next...
Anyway, to reduce the complexity of the future steps, I prune redundant constraints aggressively.
Okay, so I ran this with Gurobi as comparison:
julia> using JuMP, HiGHS, Gurobi
julia> function run_model(optimizer)
model = Model(optimizer)
@variable(model, v_0, Int)
@variable(model, v_1, Int)
@variable(model, v_2, Int)
@variable(model, v_3, Int)
@variable(model, v_4, 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)
return solution_summary(model; verbose = true)
end
run_model (generic function with 2 methods)
julia> run_model(Gurobi.Optimizer)
Set parameter Username
Academic license - for non-commercial use only - expires 2022-06-19
Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (mac64[x86])
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 4 rows, 8 columns and 17 nonzeros
Model fingerprint: 0xe6b080f0
Variable types: 0 continuous, 8 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 3e+00]
Objective range [1e+00, 2e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 4 rows and 8 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)
Solution count 1: -0
No other solutions better than -0
Optimal solution found (tolerance 1.00e-04)
Best objective -0.000000000000e+00, best bound -0.000000000000e+00, gap 0.0000%
User-callback calls 182, time in user-callback 0.00 sec
* Solver : Gurobi
* Status
Termination status : OPTIMAL
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Result count : 1
Has duals : false
Message from the solver:
"Model was solved to optimality (subject to tolerances), and an optimal solution is available."
* Candidate solution
Objective value : -0.00000e+00
Objective bound : -0.00000e+00
Relative gap : 0.00000e+00
Dual objective value : -0.00000e+00
Primal solution :
v_0 : 0.00000e+00
v_1 : 0.00000e+00
v_2 : 0.00000e+00
v_3 : 0.00000e+00
v_4 : 0.00000e+00
v_5 : 0.00000e+00
v_6 : -0.00000e+00
v_7 : -0.00000e+00
* Work counters
Solve time (sec) : 6.02961e-04
Barrier iterations : 0
Node count : 0
julia> run_model(HiGHS.Optimizer)
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)
* Solver : HiGHS
* Status
Termination status : OPTIMAL
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Result count : 1
Has duals : false
Message from the solver:
"kHighsModelStatusOptimal"
* Candidate solution
Objective value : NaN
Objective bound : NaN
Relative gap : 0.00000e+00
Primal solution :
v_0 : NaN
v_1 : NaN
v_2 : NaN
v_3 : NaN
v_4 : 0.00000e+00
v_5 : 0.00000e+00
v_6 : 0.00000e+00
v_7 : 0.00000e+00
* Work counters
Solve time (sec) : 5.34238e-04
Simplex iterations : -1
Barrier iterations : -1
Node count : 0
So even your original case was only working because of the large bounds and thus was returning an incorrect solution (one reason not to use fake bounds).
If we turn off presolve, we get:
julia> run_model(optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "off"))
Presolve is switched off
Objective function is integral with scale 1
Solving MIP model with:
4 rows
8 cols (0 binary, 8 integer, 0 implied int., 0 continuous)
17 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% inf -inf inf 0 0 0 0 0.0s
T 0 0 0 0.00% inf -0 Large 0 0 0 0 0.0s
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 1
LP iterations 0 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
* Solver : HiGHS
* Status
Termination status : OPTIMAL
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Result count : 1
Has duals : false
Message from the solver:
"kHighsModelStatusOptimal"
* Candidate solution
Objective value : 0.00000e+00
Objective bound : 0.00000e+00
Relative gap : 0.00000e+00
Primal solution :
v_0 : 0.00000e+00
v_1 : 0.00000e+00
v_2 : -0.00000e+00
v_3 : -0.00000e+00
v_4 : 0.00000e+00
v_5 : 0.00000e+00
v_6 : 0.00000e+00
v_7 : 0.00000e+00
* Work counters
Solve time (sec) : 1.15691e-03
Simplex iterations : -1
Barrier iterations : -1
Node count : 1
So this looks like a bug in presolve, and probably does some / 0 that introduces a NaN into the constraint matrix.
Thanks for your contributions to the discussion @odow Let's see what @lgottwald makes of it
Any comments on the approach toward passModel?
If the C-api of passModel fills out internal data structures, it'd be faster to just fill them directly, e.g. like in my first C++ version, instead of filling temporary objects just to copy them to a destination.
Also, auto complete makes it look like there's an API for updating existing models? Documentation seems rather sparse. This would be useful, as I need to do many solves from modifications of earlier models. E.g., for minimizing a series of lexicographically ordered objectives. Each lexicographically ordered objective is itself a series of solves, adding each found solution as new equality constraints (i.e. that the cost function equals the optimal objective). Here, it would also probably be useful to be able to use the previous solution as an initial condition, as that'd be a feasible point on the boundary, which I imagine the simplex algorithm would like as a starting point.
The C++ method Highs::passModel uses std::move to avoid copies, but I guessed this wasn't possible in C when passing pointers.
As for the hot starting, yes, this comes automatically with HiGHS when you add rows or columns after solving an LP with simplex
There is a bug in postsolve. I did not handle the rare case of duplicate/parallel columns where the columns are free and integer. In the postsolve case I just used the infinite bounds without checking. For LP this usually is no issue as duplicate columns with infinite bounds usually imply a dominated column which is handled differently. Only when there is integrality the two columns are replaced by one merged column and postsolve needs to split the assigned value into the original columns.
And for setting up the model in C++ using the version that fills the HighsLp can avoid copying, but the above code does not. It still needs to call
highs.passModel(std::move(model));
instead of
highs.passModel(model);
A few comments:
I don't have any background in integer optimization, so odds are a lot of things could be done better / I am open to suggestions.
I came up with what is probably a far better way for removing extra variables: adding augment variables to turn linear inequalities into equalities, similar to the simplex algorithm. It's easy to eliminate variables from equalities, and inequalities can be recovered in the end, needing only a little bit of cleaning to remove redundancies (which I may use HiGHS for, on top of needing it for the bigger ILPs). The initial implementation seems to work well so far.
For LP this usually is no issue as duplicate columns with infinite bounds usually imply a dominated column which is handled differently.
Okay. As a somewhat related problem where one column is a linear combination of two others, if we have the system
x + y + z <= 10
2x + y - z <= 0
-x + y + 5z <= 30
here, the third constraint equals 3 times the first minus 2 times the second. However, the third constraint is not redundant, as the solution [-9, 17, 1] satisfies the first two constraints while violating the third:
julia> [ 1 1 1;
2 1 -1
-1 1 5] * [-9, 17, 1]
3-element Vector{Int64}:
9
-2
31
Oddly, HiGHS fails to find this solution while claiming is result is optimal:
julia> function run_model(verbose = true)
model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => verbose))
@variable(model, x, Int)
@variable(model, y, Int)
@variable(model, z, Int)
@objective(model, Max, -x + y + 5z)
@constraint(model, x + y + z <= 10)
@constraint(model, 2x + y - z <= 0)
@constraint(model, -x + y + 5z <= 31)
# print(model)
optimize!(model)
return model
end
run_model (generic function with 2 methods)
julia> m = run_model()
Presolving model
3 rows, 3 cols, 9 nonzeros
3 rows, 3 cols, 9 nonzeros
Objective function is integral with scale 1
Solving MIP model with:
3 rows
3 cols (0 binary, 3 integer, 0 implied int., 0 continuous)
9 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% inf -inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 31 30 3.33% 0 0 0 1 0.0s
Solving report
Status Optimal
Primal bound 30
Dual bound 30
Gap 0% (tolerance: 0.01%)
Solution status feasible
30 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 1 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
A JuMP Model
Maximization problem with:
Variables: 3
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 3 constraints
`VariableRef`-in-`MathOptInterface.Integer`: 3 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: HiGHS
Names registered in the model: x, y, z
julia> solution_summary(m, verbose=true)
* Solver : HiGHS
* Status
Termination status : OPTIMAL
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Result count : 1
Has duals : false
Message from the solver:
"kHighsModelStatusOptimal"
* Candidate solution
Objective value : 3.00000e+01
Objective bound : 3.00000e+01
Relative gap : 0.00000e+00
Primal solution :
x : 0.00000e+00
y : 0.00000e+00
z : 6.00000e+00
* Work counters
Solve time (sec) : 1.65844e-03
Simplex iterations : -1
Barrier iterations : -1
Node count : 1
It finds [0, 0, 6] as the optimal solution, even though [-9, 17, 1] is larger and also satisfies all the constraints.
Perhaps to @odow's point about adding random constraints:
julia> function run_model_lb(verbose = true)
model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => verbose))
@variable(model, x >= -10, Int)
@variable(model, y >= -10, Int)
@variable(model, z >= -10, Int)
@objective(model, Max, -x + y + 5z)
@constraint(model, x + y + z <= 10)
@constraint(model, 2x + y - z <= 0)
@constraint(model, -x + y + 5z <= 31)
# print(model)
optimize!(model)
return model
end
run_model_lb (generic function with 2 methods)
julia> m_lb = run_model_lb()
Presolving model
3 rows, 3 cols, 9 nonzeros
3 rows, 3 cols, 9 nonzeros
Objective function is integral with scale 1
Solving MIP model with:
3 rows
3 cols (0 binary, 3 integer, 0 implied int., 0 continuous)
9 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 90 -inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 31 29 6.90% 0 0 0 3 0.0s
C 0 0 0 0.00% 31 30 3.33% 6 2 0 5 0.0s
T 0 0 0 0.00% 31 31 0.00% 6 2 0 5 0.0s
Solving report
Status Optimal
Primal bound 31
Dual bound 31
Gap 0% (tolerance: 0.01%)
Solution status feasible
31 (objective)
0 (bound viol.)
1.06581410364e-14 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 5 (total)
0 (strong br.)
2 (separation)
0 (heuristics)
A JuMP Model
Maximization problem with:
Variables: 3
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 3 constraints
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 3 constraints
`VariableRef`-in-`MathOptInterface.Integer`: 3 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: HiGHS
Names registered in the model: x, y, z
julia> solution_summary(m_lb, verbose=true)
* Solver : HiGHS
* Status
Termination status : OPTIMAL
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Result count : 1
Has duals : false
Message from the solver:
"kHighsModelStatusOptimal"
* Candidate solution
Objective value : 3.10000e+01
Objective bound : 3.10000e+01
Relative gap : 0.00000e+00
Primal solution :
x : -9.00000e+00
y : 1.70000e+01
z : 1.00000e+00
* Work counters
Solve time (sec) : 2.33388e-03
Simplex iterations : -1
Barrier iterations : -1
Node count : 1
adding bounds on the variables, but leaving the objective and constraints the same, we now find a better objective (the solution [-9, 17, 1])!
Should I file a separate issue for this bug?
I also could not reproduce this. The implied lower bounds explain this.
Another note: You should set the relative gap tolerance, mip_rel_gap to 0.0 if you use HiGHS to determine redundancy. Otherwise it may happen that you declare a constraint redundant when it is not redundant in absolute tolerances but only in relative tolerance of 1e-4. E.g. a constraint x + y + z <= 10000 where x + y + z can attain 10001 might still be declared redundant since the solution with value 10000 has a maximum gap of 1e-4 which means optimal intolerances. This is not guaranteed of course, since it might happen that HiGHS finds the solution with value 10001 before proving a sufficient dual bound to terminate.
If your objective function stays in a reasonable scale this will usually be ok though and the default tolerances match those of gurobi's default tolerances.
The implied lower bounds explain this.
JuMP doesn't have implied lower bounds. The variables are free with -Inf, Inf bounds.
Could you write the model as MPS please then, so I can compare with the same model (.lp can change order columns/rows)
Here's what I get with v1.2.1:
shell> cat problem.mps
NAME
ROWS
N COST
L R0
L R1
L R2
COLUMNS
MARK0000 'MARKER' 'INTORG'
C0 COST 1
C0 R0 1
C0 R1 2
C0 R2 -1
C1 COST -1
C1 R0 1
C1 R1 1
C1 R2 1
C2 COST -5
C2 R0 1
C2 R1 -1
C2 R2 5
RHS
RHS_V R0 10
RHS_V R2 31
BOUNDS
FR BOUND C0
FR BOUND C1
FR BOUND C2
ENDATA
julia> HiGHS_jll.highs() do exe
run(`$exe --solution_file highs.out problem.mps`)
end
Running HiGHS 1.2.1 [date: , git hash: ]
Copyright (c) 2022 ERGO-Code under MIT licence terms
MIP problem has 3 rows; 3 cols; 9 nonzeros; 3 integer variables
Presolving model
3 rows, 3 cols, 9 nonzeros
3 rows, 3 cols, 9 nonzeros
Objective function is integral with scale 1
Solving MIP model with:
3 rows
3 cols (0 binary, 3 integer, 0 implied int., 0 continuous)
9 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% -inf inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% -31 -30 3.33% 0 0 0 1 0.0s
Solving report
Status Optimal
Primal bound -30
Dual bound -30
Gap 0% (tolerance: 0.01%)
Solution status feasible
-30 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 1 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
Process(`/Users/oscar/.julia/artifacts/b22cc32ff7b7fd1398043f1f6276b38d00b63864/bin/highs --solution_file highs.out problem.mps`, ProcessExited(0))
shell> cat highs.out
Model status
Optimal
# Primal solution values
Feasible
Objective -30
# Columns 3
C0 0
C1 0
C2 6
# Rows 3
R0 6
R1 -6
R2 30
# Dual solution values
None
# Basis
HiGHS v1
None
JuMP doesn't assume a lower bound of 0. What happens if you turn presolve off?
On Fri, 27 May 2022, 6:36 PM Julian Hall, @.***> wrote:
Oddly, HiGHS fails to find this solution while claiming is result is optimal:
It finds [0, 0, 6] as the optimal solution, even though [-9, 17, 1] is larger and also satisfies all the constraints.
It looks very much as if your variable definition implies (as us usually the case in modelling systems) that variables are non-negative unless a lower bound is specified.
@variable(model, x >= -10, Int) @variable(model, y >= -10, Int) @variable(model, z >= -10, Int)Now you've allowed variables to be negative, HiGHS (unsurprisingly) gets the optimal solution you expected earlier
Solving report Status Optimal Primal bound 31 Dual bound 31 Gap 0% (tolerance: 0.01%) Solution status feasible 31 (objective) 0 (bound viol.) 1.06581410364e-14 (int. viol.) 0 (row viol.) Timing 0.00 (total) 0.00 (presolve) 0.00 (postsolve) Nodes 1 LP iterations 5 (total)
adding bounds on the variables, but leaving the objective and constraints the same, we now find a better objective (the solution [-9, 17, 1])!
Should I file a separate issue for this bug?
No bug that I can see
— Reply to this email directly, view it on GitHub https://github.com/ERGO-Code/HiGHS/issues/888#issuecomment-1139328225, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6MQJMS77XRMWC2CFTQQNLVMBUPLANCNFSM5W7L6PNA . You are receiving this because you were mentioned.Message ID: @.***>