diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b1cf83b..7704ea44 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/uno/ingredients/subproblems/Direction.cpp b/uno/ingredients/subproblems/Direction.cpp index 3ab6b453..c58a74a6 100644 --- a/uno/ingredients/subproblems/Direction.cpp +++ b/uno/ingredients/subproblems/Direction.cpp @@ -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: diff --git a/uno/solvers/HiGHS/HiGHSSolver.cpp b/uno/solvers/HiGHS/HiGHSSolver.cpp new file mode 100644 index 00000000..a6c769c0 --- /dev/null +++ b/uno/solvers/HiGHS/HiGHSSolver.cpp @@ -0,0 +1,140 @@ +#include +#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& variables_lower_bounds, + const std::vector& variables_upper_bounds, const std::vector& constraints_lower_bounds, + const std::vector& constraints_upper_bounds, const SparseVector& linear_objective, + const RectangularMatrix& constraint_jacobian) { + this->model.lp_.num_col_ = static_cast(number_variables); + this->model.lp_.num_row_ = static_cast(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(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& variables_lower_bounds, + const std::vector& variables_upper_bounds, const std::vector& constraints_lower_bounds, + const std::vector& constraints_upper_bounds, const SparseVector& linear_objective, + const RectangularMatrix& constraint_jacobian, const Vector& /*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 \ No newline at end of file diff --git a/uno/solvers/HiGHS/HiGHSSolver.hpp b/uno/solvers/HiGHS/HiGHSSolver.hpp new file mode 100644 index 00000000..5f78285b --- /dev/null +++ b/uno/solvers/HiGHS/HiGHSSolver.hpp @@ -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& variables_lower_bounds, + const std::vector& variables_upper_bounds, const std::vector& constraints_lower_bounds, + const std::vector& constraints_upper_bounds, const SparseVector& linear_objective, + const RectangularMatrix& constraint_jacobian, const Vector& 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& variables_lower_bounds, + const std::vector& variables_upper_bounds, const std::vector& constraints_lower_bounds, + const std::vector& constraints_upper_bounds, const SparseVector& linear_objective, + const RectangularMatrix& constraint_jacobian); + void solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints); + }; +} // namespace + +#endif // UNO_HIGHSSOLVER_H \ No newline at end of file diff --git a/uno/solvers/LPSolverFactory.hpp b/uno/solvers/LPSolverFactory.hpp index 548a6347..8613ec27 100644 --- a/uno/solvers/LPSolverFactory.hpp +++ b/uno/solvers/LPSolverFactory.hpp @@ -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(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(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 solvers{}; #ifdef HAS_BQPD solvers.emplace_back("BQPD"); +#endif +#ifdef HAS_HIGHS + solvers.emplace_back("HiGHS"); #endif return solvers; } diff --git a/unotest/HiGHSSolverTests.cpp b/unotest/HiGHSSolverTests.cpp new file mode 100644 index 00000000..d21304d5 --- /dev/null +++ b/unotest/HiGHSSolverTests.cpp @@ -0,0 +1,72 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#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 linear_objective(number_variables); + linear_objective.insert(0, 1.); + linear_objective.insert(1, 1.); + + const std::vector variables_lower_bounds{0., 1.}; + const std::vector variables_upper_bounds{4., INF}; + const std::vector constraints_lower_bounds{-INF, 5., 6.}; + const std::vector constraints_upper_bounds{7., 15., INF}; + + RectangularMatrix 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 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 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 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 lower_bound_duals_reference{0., 0.}; + const std::vector 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); + } +} \ No newline at end of file