Skip to content

Commit

Permalink
Interfaced HiGHS as LPSolver + unit test
Browse files Browse the repository at this point in the history
cvanaret committed Nov 5, 2024
1 parent 84ee915 commit 56d2479
Showing 6 changed files with 270 additions and 1 deletion.
11 changes: 11 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -131,6 +131,17 @@ else()
link_to_uno(bqpd ${BQPD})
endif()

# HiGHS
find_package(HIGHS)
if(NOT HIGHS)
message(WARNING "Optional library HiGHS was not found.")
else()
list(APPEND UNO_SOURCE_FILES uno/solvers/HiGHS/HiGHSSolver.cpp)
list(APPEND TESTS_UNO_SOURCE_FILES unotest/HiGHSSolverTests.cpp)
link_to_uno(highs ${HIGHS})
list(APPEND LIBRARIES highs::highs)
endif()

# MUMPS
find_package(MUMPS)
if(NOT MUMPS_LIBRARY)
2 changes: 1 addition & 1 deletion uno/ingredients/subproblems/Direction.cpp
Original file line number Diff line number Diff line change
@@ -28,7 +28,7 @@ namespace uno {
std::string status_to_string(SubproblemStatus status) {
switch (status) {
case SubproblemStatus::OPTIMAL:
return "optimal subproblem";
return "optimal";
case SubproblemStatus::UNBOUNDED_PROBLEM:
return "unbounded subproblem";
case SubproblemStatus::INFEASIBLE:
140 changes: 140 additions & 0 deletions uno/solvers/HiGHS/HiGHSSolver.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
#include <cassert>
#include "HiGHSSolver.hpp"
#include "linear_algebra/RectangularMatrix.hpp"
#include "linear_algebra/SparseVector.hpp"
#include "linear_algebra/Vector.hpp"
#include "ingredients/subproblems/Direction.hpp"
#include "options/Options.hpp"
#include "symbolic/VectorView.hpp"

namespace uno {
HiGHSSolver::HiGHSSolver(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t /*number_hessian_nonzeros*/,
const Options& options): LPSolver(), print_subproblem(options.get_bool("print_subproblem")) {
this->model.lp_.sense_ = ObjSense::kMinimize;
this->model.lp_.offset_ = 0.;
// the linear part of the objective is a dense vector
this->model.lp_.col_cost_.resize(number_variables);
// variable bounds
this->model.lp_.col_lower_.resize(number_variables);
this->model.lp_.col_upper_.resize(number_variables);
// constraint bounds
this->model.lp_.row_lower_.resize(number_constraints);
this->model.lp_.row_upper_.resize(number_constraints);
// constraint matrix in CSR format (each row is a constraint gradient)
this->model.lp_.a_matrix_.format_ = MatrixFormat::kRowwise;
this->model.lp_.a_matrix_.value_.reserve(number_jacobian_nonzeros);
this->model.lp_.a_matrix_.index_.reserve(number_jacobian_nonzeros);
this->model.lp_.a_matrix_.start_.reserve(number_variables + 1);

this->highs_solver.setOptionValue("output_flag", "false");
}

void HiGHSSolver::build_linear_subproblem(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian) {
this->model.lp_.num_col_ = static_cast<HighsInt>(number_variables);
this->model.lp_.num_row_ = static_cast<HighsInt>(number_constraints);

// variable bounds
for (size_t variable_index = 0; variable_index < number_variables; variable_index++) {
this->model.lp_.col_lower_[variable_index] = variables_lower_bounds[variable_index];
this->model.lp_.col_upper_[variable_index] = variables_upper_bounds[variable_index];
// reset the linear part of the objective
this->model.lp_.col_cost_[variable_index] = 0.;
}

// linear part of the objective
for (const auto [variable_index, value]: linear_objective) {
this->model.lp_.col_cost_[variable_index] = value;
}

// constraint bounds
for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) {
this->model.lp_.row_lower_[constraint_index] = constraints_lower_bounds[constraint_index];
this->model.lp_.row_upper_[constraint_index] = constraints_upper_bounds[constraint_index];
}

// constraint matrix
this->model.lp_.a_matrix_.value_.clear();
this->model.lp_.a_matrix_.index_.clear();
this->model.lp_.a_matrix_.start_.clear();

size_t number_nonzeros = 0;
this->model.lp_.a_matrix_.start_.emplace_back(number_nonzeros);
for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) {
for (const auto [variable_index, value]: constraint_jacobian[constraint_index]) {
this->model.lp_.a_matrix_.value_.emplace_back(value);
this->model.lp_.a_matrix_.index_.emplace_back(variable_index);
number_nonzeros++;
}
this->model.lp_.a_matrix_.start_.emplace_back(number_nonzeros);
}

if (this->print_subproblem) {
DEBUG << "Linear objective part: "; print_vector(DEBUG, view(this->model.lp_.col_cost_, 0, number_variables));
DEBUG << "Jacobian:\n";
DEBUG << "J = "; print_vector(DEBUG, this->model.lp_.a_matrix_.value_);
DEBUG << "with column start: "; print_vector(DEBUG, this->model.lp_.a_matrix_.start_);
DEBUG << "and row index: "; print_vector(DEBUG, this->model.lp_.a_matrix_.index_);
for (size_t variable_index = 0; variable_index < number_variables; variable_index++) {
DEBUG << "d" << variable_index << " in [" << this->model.lp_.col_lower_[variable_index] << ", " <<
this->model.lp_.col_upper_[variable_index] << "]\n";
}
for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) {
DEBUG << "linearized c" << constraint_index << " in [" << this->model.lp_.row_lower_[constraint_index] << ", " <<
this->model.lp_.row_upper_[constraint_index]<< "]\n";
}
}
}



void HiGHSSolver::solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints) {
// solve the LP
HighsStatus return_status = this->highs_solver.passModel(this->model);
assert(return_status == HighsStatus::kOk);
return_status = this->highs_solver.run(); // solve
DEBUG << "HiGHS status: " << static_cast<int>(return_status) << '\n';

// if HiGHS could not optimize (e.g. because of indefinite Hessian), return an error
if (return_status == HighsStatus::kError) {
direction.status = SubproblemStatus::ERROR;
return;
}

// TODO check unbounded problems
direction.status = SubproblemStatus::OPTIMAL;
const HighsSolution& solution = this->highs_solver.getSolution();
// read the primal solution and bound dual solution
for (size_t variable_index = 0; variable_index < number_variables; variable_index++) {
direction.primals[variable_index] = solution.col_value[variable_index];
const double bound_multiplier = solution.col_dual[variable_index];
if (0. < bound_multiplier) {
direction.multipliers.lower_bounds[variable_index] = bound_multiplier;
}
else {
direction.multipliers.upper_bounds[variable_index] = bound_multiplier;
}
}
// read the dual solution
for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) {
direction.multipliers.constraints[constraint_index] = solution.row_dual[constraint_index];
}
const HighsInfo& info = this->highs_solver.getInfo();
direction.subproblem_objective = info.objective_function_value;
}

void HiGHSSolver::solve_LP(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian, const Vector<double>& /*initial_point*/, Direction& direction,
const WarmstartInformation& /*warmstart_information*/) {
// build the LP in the HiGHS format
this->build_linear_subproblem(number_variables, number_constraints, variables_lower_bounds, variables_upper_bounds, constraints_lower_bounds,
constraints_upper_bounds, linear_objective, constraint_jacobian);

// solve the LP
this->solve_subproblem(direction, number_variables, number_constraints);
}
} // namespace
35 changes: 35 additions & 0 deletions uno/solvers/HiGHS/HiGHSSolver.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
#ifndef UNO_HIGHSSOLVER_H
#define UNO_HIGHSSOLVER_H

#include "solvers/LPSolver.hpp"
#include "Highs.h"

namespace uno {
// forward declaration
class Options;

class HiGHSSolver : public LPSolver {
public:
HiGHSSolver(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros,
const Options& options);

void solve_LP(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian, const Vector<double>& initial_point, Direction& direction,
const WarmstartInformation& warmstart_information) override;

protected:
HighsModel model;
Highs highs_solver;
const bool print_subproblem;

void build_linear_subproblem(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian);
void solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints);
};
} // namespace

#endif // UNO_HIGHSSOLVER_H
11 changes: 11 additions & 0 deletions uno/solvers/LPSolverFactory.hpp
Original file line number Diff line number Diff line change
@@ -13,6 +13,9 @@
#ifdef HAS_BQPD
#include "solvers/BQPD/BQPDSolver.hpp"
#endif
#ifdef HAS_HIGHS
#include "solvers/HiGHS/HiGHSSolver.hpp"
#endif

namespace uno {
class LPSolverFactory {
@@ -27,6 +30,11 @@ namespace uno {
return std::make_unique<BQPDSolver>(number_variables, number_constraints, number_objective_gradient_nonzeros,
number_jacobian_nonzeros, 0, BQPDProblemType::LP, options);
}
#endif
#ifdef HAS_HIGHS
if (LP_solver_name == "HiGHS") {
return std::make_unique<HiGHSSolver>(number_variables, number_constraints, number_jacobian_nonzeros, 0, options);
}
#endif
std::string message = "The LP solver ";
message.append(LP_solver_name).append(" is unknown").append("\n").append("The following values are available: ")
@@ -45,6 +53,9 @@ namespace uno {
std::vector<std::string> solvers{};
#ifdef HAS_BQPD
solvers.emplace_back("BQPD");
#endif
#ifdef HAS_HIGHS
solvers.emplace_back("HiGHS");
#endif
return solvers;
}
72 changes: 72 additions & 0 deletions unotest/HiGHSSolverTests.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
// Copyright (c) 2024 Charlie Vanaret
// Licensed under the MIT license. See LICENSE file in the project directory for details.

#include <gtest/gtest.h>
#include "ingredients/subproblems/Direction.hpp"
#include "linear_algebra/RectangularMatrix.hpp"
#include "linear_algebra/SparseVector.hpp"
#include "optimization/WarmstartInformation.hpp"
#include "options/Options.hpp"
#include "solvers/HiGHS/HiGHSSolver.hpp"
#include "tools/Infinity.hpp"

using namespace uno;

TEST(HiGHSSolver, DocumentationLP) {
// https://ergo-code.github.io/HiGHS/stable/interfaces/cpp/library/
// Min f = x_0 + x_1 + 3
// s.t. x_1 <= 7
// 5 <= x_0 + 2x_1 <= 15
// 6 <= 3x_0 + 2x_1
// 0 <= x_0 <= 4; 1 <= x_1

const size_t number_variables = 2;
const size_t number_constraints = 3;
const size_t number_jacobian_nonzeros = 5;
const size_t number_hessian_nonzeros = 0;
Options options(false);
options["print_subproblem"] = "false";
HiGHSSolver highs_solver(number_variables, number_constraints, number_jacobian_nonzeros, number_hessian_nonzeros, options);

// create the LP
SparseVector<double> linear_objective(number_variables);
linear_objective.insert(0, 1.);
linear_objective.insert(1, 1.);

const std::vector<double> variables_lower_bounds{0., 1.};
const std::vector<double> variables_upper_bounds{4., INF<double>};
const std::vector<double> constraints_lower_bounds{-INF<double>, 5., 6.};
const std::vector<double> constraints_upper_bounds{7., 15., INF<double>};

RectangularMatrix<double> constraint_jacobian(number_constraints, number_variables);
constraint_jacobian[0].insert(1, 1.);
constraint_jacobian[1].insert(0, 1.);
constraint_jacobian[1].insert(1, 2.);
constraint_jacobian[2].insert(0, 3.);
constraint_jacobian[2].insert(1, 2.);

Direction direction(number_variables, number_constraints);
WarmstartInformation warmstart_information{};
Vector<double> initial_point{0., 0.};

highs_solver.solve_LP(number_variables, number_constraints, variables_lower_bounds, variables_upper_bounds, constraints_lower_bounds,
constraints_upper_bounds, linear_objective, constraint_jacobian, initial_point, direction, warmstart_information);

const double tolerance = 1e-8;
// check primals
const std::vector<double> primals_reference{0.5, 2.25};
for (size_t index: Range(number_variables)) {
EXPECT_NEAR(direction.primals[index], primals_reference[index], tolerance);
}
// check duals
const std::vector<double> constraint_duals_reference{0., 0.25, 0.25};
for (size_t index: Range(number_constraints)) {
EXPECT_NEAR(direction.multipliers.constraints[index], constraint_duals_reference[index], tolerance);
}
const std::vector<double> lower_bound_duals_reference{0., 0.};
const std::vector<double> upper_bound_duals_reference{0., 0.};
for (size_t index: Range(number_variables)) {
EXPECT_NEAR(direction.multipliers.lower_bounds[index], lower_bound_duals_reference[index], tolerance);
EXPECT_NEAR(direction.multipliers.upper_bounds[index], upper_bound_duals_reference[index], tolerance);
}
}

0 comments on commit 56d2479

Please sign in to comment.