diff --git a/build_dev.sh b/build_dev.sh index 7bb61915..0f5d130d 100755 --- a/build_dev.sh +++ b/build_dev.sh @@ -57,7 +57,7 @@ else . ${SPACK_HOME}/share/spack/setup-env.sh spack env create coekenv spack env activate coekenv - spack add asl cppad fmt rapidjson catch2 + spack add asl cppad fmt rapidjson catch2 highs spack install spack env deactivate fi @@ -67,5 +67,5 @@ echo "" \rm -Rf build mkdir build cd build -cmake -DCMAKE_PREFIX_PATH=${SPACK_HOME}/var/spack/environments/coekenv/.spack-env/view -Dwith_python=${with_python} -Dwith_gurobi=$with_gurobi -Dwith_cppad=ON -Dwith_fmtlib=ON -Dwith_rapidjson=ON -Dwith_catch2=ON -Dwith_tests=ON -Dwith_asl=ON -Dwith_openmp=OFF .. +cmake -DCMAKE_PREFIX_PATH=${SPACK_HOME}/var/spack/environments/coekenv/.spack-env/view -Dwith_python=${with_python} -Dwith_gurobi=$with_gurobi -Dwith_highs=ON -Dwith_cppad=ON -Dwith_fmtlib=ON -Dwith_rapidjson=ON -Dwith_catch2=ON -Dwith_tests=ON -Dwith_asl=ON -Dwith_openmp=OFF .. make -j20 diff --git a/lib/coek/CMakeLists.txt b/lib/coek/CMakeLists.txt index 32296856..0d141c56 100644 --- a/lib/coek/CMakeLists.txt +++ b/lib/coek/CMakeLists.txt @@ -28,6 +28,7 @@ option(with_asl "Use the ASL autograd library" OFF) # Solver Dependencies option(with_gurobi "Use the Gurobi solver" OFF) set(GUROBI_HOME "" CACHE FILEPATH "Set the path to gurobi") +option(with_highs "Use the Highs solver" OFF) ##################### Other Config ##################### @@ -42,8 +43,13 @@ endif() MESSAGE("-- With Gurobi Solver: ${with_gurobi}") if(with_gurobi) - find_package(GUROBI) - list(APPEND solver_flags "-DWITH_GUROBI") + find_package(GUROBI REQUIRED) +endif() + +MESSAGE("-- With Highs Solver: ${with_highs}") +if(with_highs) + find_package(HIGHS REQUIRED) + find_package(Threads REQUIRED) endif() diff --git a/lib/coek/coek/CMakeLists.txt b/lib/coek/coek/CMakeLists.txt index 3d70bcff..c098994e 100644 --- a/lib/coek/coek/CMakeLists.txt +++ b/lib/coek/coek/CMakeLists.txt @@ -34,6 +34,7 @@ SET(sources solvers/solver_results.cpp solvers/solver.cpp solvers/solver_repn.cpp + solvers/create_solver.cpp solvers/testsolver.cpp autograd/autograd.cpp abstract/expr_rule.cpp @@ -154,16 +155,19 @@ endif() if(with_gurobi) list(APPEND sources solvers/gurobi/gurobi.cpp) - #if(NOT "$ENV{GUROBI_HOME}" STREQUAL "") - # MESSAGE("-- GUROBI_HOME $ENV{GUROBI_HOME}") - # list(APPEND coek_link_directories $ENV{GUROBI_HOME}/lib) - #else() - # MESSAGE("-- GUROBI_HOME Not provided") - #endif() list(APPEND coek_compile_options -DWITH_GUROBI) list(APPEND coek_link_libraries ${GUROBI_CXX_LIBRARY} ${GUROBI_LIBRARY}) list(APPEND coek_include_directories ${GUROBI_INCLUDE_DIRS}) -#gurobi_g++5.2 gurobi90) +endif() + +# +# HIGHS LIBRARY +# +if(with_highs) + list(APPEND sources + solvers/highs/highs.cpp) + list(APPEND coek_compile_options -DWITH_HIGHS) + list(APPEND coek_link_libraries highs::highs) endif() # diff --git a/lib/coek/coek/api/expression_visitor.hpp b/lib/coek/coek/api/expression_visitor.hpp index a0e5583a..a6d26d6f 100644 --- a/lib/coek/coek/api/expression_visitor.hpp +++ b/lib/coek/coek/api/expression_visitor.hpp @@ -5,6 +5,18 @@ namespace coek { +class ParameterTerm; +class IndexParameterTerm; +class VariableTerm; +class BaseExpressionTerm; +class SubExpressionTerm; + +typedef std::shared_ptr ParameterRepn; +typedef std::shared_ptr IndexParameterRepn; +typedef std::shared_ptr VariableRepn; +typedef std::shared_ptr ExpressionRepn; +typedef std::shared_ptr SubExpressionRepn; + class VariableTerm; class Expression; class Objective; diff --git a/lib/coek/coek/autograd/asl_repn.cpp b/lib/coek/coek/autograd/asl_repn.cpp index a5070e81..bc423519 100644 --- a/lib/coek/coek/autograd/asl_repn.cpp +++ b/lib/coek/coek/autograd/asl_repn.cpp @@ -467,10 +467,10 @@ void ASL_Repn::call_hesset() nnz_lag_h = static_cast(sphsetup(-1, coeff_obj, mult_supplied, uptri)); } -bool ASL_Repn::get_option(const std::string& option, std::string& value) const +bool ASL_Repn::get_option(const std::string& option, bool& value) const { - if (option == "temp_directory") { - value = temp_directory; + if (option == "remove_nl_file") { + value = remove_nl_file; return true; } return false; @@ -485,16 +485,31 @@ bool ASL_Repn::get_option(const std::string& option, int& value) const return false; } -void ASL_Repn::set_option(const std::string& option, const std::string value) +bool ASL_Repn::get_option(const std::string& option, std::string& value) const { - if (option == "temp_directory") - temp_directory = value; + if (option == "temp_directory") { + value = temp_directory; + return true; + } + return false; } -void ASL_Repn::set_option(const std::string& option, int value) +void ASL_Repn::set_option(const std::string& option, bool value) { if (option == "remove_nl_file") remove_nl_file = value; } +void ASL_Repn::set_option(const std::string& option, int value) +{ + if (option == "remove_nl_file") + remove_nl_file = (value == 1); +} + +void ASL_Repn::set_option(const std::string& option, const std::string value) +{ + if (option == "temp_directory") + temp_directory = value; +} + } // namespace coek diff --git a/lib/coek/coek/autograd/asl_repn.hpp b/lib/coek/coek/autograd/asl_repn.hpp index 8e78f9d7..e8f2cff6 100644 --- a/lib/coek/coek/autograd/asl_repn.hpp +++ b/lib/coek/coek/autograd/asl_repn.hpp @@ -22,8 +22,8 @@ class ASL_Repn : public NLPModelRepn { ASL_pfgh* asl_; // The temporary file directory used to write the NL file std::string temp_directory; - // If set to 1, then the NL file is removed - int remove_nl_file = 1; + // If true, then the NL file is removed + bool remove_nl_file = true; // ---------------------------------------------------------------------- // Problem information @@ -103,10 +103,12 @@ class ASL_Repn : public NLPModelRepn { void compute_J(std::vector& J); public: - bool get_option(const std::string& option, std::string& value) const; + bool get_option(const std::string& option, bool& value) const; bool get_option(const std::string& option, int& value) const; - void set_option(const std::string& option, const std::string value); + bool get_option(const std::string& option, std::string& value) const; + void set_option(const std::string& option, bool value); void set_option(const std::string& option, int value); + void set_option(const std::string& option, const std::string value); protected: void* nerror_; diff --git a/lib/coek/coek/autograd/autograd.hpp b/lib/coek/coek/autograd/autograd.hpp index 1605bbde..b637ac47 100644 --- a/lib/coek/coek/autograd/autograd.hpp +++ b/lib/coek/coek/autograd/autograd.hpp @@ -66,6 +66,17 @@ class NLPModelRepn { void find_used_variables(); public: + /** Get the value of an boolean option + * + * The option value is returned by reference if it has + * a value. + * + * \param option the option name + * \param value a boolean value that is passed by reference + * + * \returns \c true if the option is found + */ + virtual bool get_option(const std::string& /*option*/, bool& /*value*/) const { return false; } /** Get the value of an integer option * * The option value is returned by reference if it has @@ -106,6 +117,12 @@ class NLPModelRepn { return false; } + /** Set a boolean option + * + * \param option the option name + * \param value the boolean value + */ + virtual void set_option(const std::string& /*option*/, bool /*value*/) {} /** Set an integer option * * \param option the option name @@ -123,7 +140,17 @@ class NLPModelRepn { * \param option the option name * \param value the string value */ - virtual void set_option(const std::string& /*option*/, const std::string /*value*/) {} + virtual void set_option(const std::string& /*option*/, const std::string& /*value*/) {} + + /** Set a string option + * + * \param option the option name + * \param value the string value + */ + void set_option(const std::string& option, const char* value) + { + this->set_option(option, std::string(value)); + } }; NLPModelRepn* create_NLPModelRepn(Model& model, const std::string& name); diff --git a/lib/coek/coek/autograd/cppad_repn.cpp b/lib/coek/coek/autograd/cppad_repn.cpp index 8a7d9b61..d7c6a9e4 100644 --- a/lib/coek/coek/autograd/cppad_repn.cpp +++ b/lib/coek/coek/autograd/cppad_repn.cpp @@ -610,6 +610,19 @@ void CppAD_Repn::reset(void) set_variables(currx); } +bool CppAD_Repn::get_option(const std::string& option, bool& value) const +{ + if (option == "sparse_JH") { + value = sparse_JH; + return true; + } + else if (option == "simplify_expressions") { + value = simplify_expressions; + return true; + } + return false; +} + bool CppAD_Repn::get_option(const std::string& option, int& value) const { if (option == "sparse_JH") { @@ -623,6 +636,14 @@ bool CppAD_Repn::get_option(const std::string& option, int& value) const return false; } +void CppAD_Repn::set_option(const std::string& option, bool value) +{ + if (option == "sparse_JH") + sparse_JH = value; + else if (option == "simplify_expressions") + simplify_expressions = value; +} + void CppAD_Repn::set_option(const std::string& option, int value) { if (option == "sparse_JH") diff --git a/lib/coek/coek/autograd/cppad_repn.hpp b/lib/coek/coek/autograd/cppad_repn.hpp index 58a0d093..27526af9 100644 --- a/lib/coek/coek/autograd/cppad_repn.hpp +++ b/lib/coek/coek/autograd/cppad_repn.hpp @@ -127,8 +127,10 @@ class CppAD_Repn : public NLPModelRepn { void compute_J(std::vector& J); public: + bool get_option(const std::string& option, bool& value) const; bool get_option(const std::string& option, int& value) const; + void set_option(const std::string& option, bool value); void set_option(const std::string& option, int value); public: diff --git a/lib/coek/coek/model/nlp_model.cpp b/lib/coek/coek/model/nlp_model.cpp index aff9219c..1da4c32e 100644 --- a/lib/coek/coek/model/nlp_model.cpp +++ b/lib/coek/coek/model/nlp_model.cpp @@ -23,6 +23,8 @@ void NLPModel::initialize(Model& model, std::string type) std::shared_ptr tmp(create_NLPModelRepn(model, type)); repn = tmp; + for (const auto& ioption : boolean_options()) + repn->set_option(ioption.first, ioption.second); for (const auto& ioption : integer_options()) repn->set_option(ioption.first, ioption.second); for (const auto& doption : double_options()) diff --git a/lib/coek/coek/solvers/create_solver.cpp b/lib/coek/coek/solvers/create_solver.cpp new file mode 100644 index 00000000..6e0aa9e9 --- /dev/null +++ b/lib/coek/coek/solvers/create_solver.cpp @@ -0,0 +1,49 @@ +#include "coek/solvers/solver_repn.hpp" + +namespace coek { + +SolverRepn* create_coektest_solver(); +NLPSolverRepn* create_ipopt_solver(); +#ifdef WITH_GUROBI +SolverRepn* create_gurobi_solver(); +#endif +#ifdef WITH_HIGHS +SolverRepn* create_highs_solver(); +#endif + +SolverRepn* create_solver(std::string& name, OptionCache& options) +{ + if (name == "test") + return create_coektest_solver(); + +#ifdef WITH_GUROBI + if (name == "gurobi") { + auto tmp = create_gurobi_solver(); + tmp->set_options(options); + return tmp; + } +#endif + +#ifdef WITH_HIGHS + if (name == "highs") { + auto tmp = create_highs_solver(); + tmp->set_options(options); + return tmp; + } +#endif + + return 0; +} + +NLPSolverRepn* create_nlpsolver(std::string& name, OptionCache& options) +{ + if (name == "ipopt") { + auto tmp = create_ipopt_solver(); + tmp->set_options(options); + return tmp; + } + + return 0; +} + +} // namespace coek diff --git a/lib/coek/coek/solvers/gurobi/gurobi.cpp b/lib/coek/coek/solvers/gurobi/gurobi.cpp index f0e855f8..d4eb0455 100644 --- a/lib/coek/coek/solvers/gurobi/gurobi.cpp +++ b/lib/coek/coek/solvers/gurobi/gurobi.cpp @@ -18,6 +18,8 @@ namespace coek { +SolverRepn* create_gurobi_solver() { return new GurobiSolver(); } + namespace { auto add_gurobi_variable(GRBModel* gmodel, double lb, double ub, const Variable& eval) @@ -340,9 +342,12 @@ std::shared_ptr GurobiSolver::solve(Model& model) } #ifdef COEK_WITH_COMPACT_MODEL -int GurobiSolver::solve(CompactModel& compact_model) +std::shared_ptr GurobiSolver::solve(CompactModel& compact_model) { - std::cout << "STARTING GUROBI" << std::endl << std::flush; + auto results = std::make_shared(); + results->solver_name = "gurobi"; + results->termination_condition = TerminationCondition::error; + results->tic(); env = new GRBEnv(true); auto it = integer_options().find("OutputFlag"); @@ -351,8 +356,6 @@ int GurobiSolver::solve(CompactModel& compact_model) env->start(); gmodel = new GRBModel(*env); - std::cout << "BUILDING GUROBI MODEL" << std::endl << std::flush; - // Add Gurobi variables for (auto& val : compact_model.repn->variables) { if (auto eval = std::get_if(&val)) { @@ -385,7 +388,7 @@ int GurobiSolver::solve(CompactModel& compact_model) else { auto& seq = std::get(val); for (auto jt = seq.begin(); jt != seq.end(); ++jt) { - model.repn->objectives.push_back(*jt); + // model.repn->objectives.push_back(*jt); Expression tmp = jt->expr(); add_gurobi_objective(gmodel, tmp, jt->sense(), x, orepn); nobj++; @@ -394,15 +397,15 @@ int GurobiSolver::solve(CompactModel& compact_model) } } catch (GRBException e) { - std::cerr << "GUROBI Exception: (objective) " << e.getMessage() << std::endl; - throw; + results->error_message + = "Gurobi Error: Caught gurobi exception while creating objectives " + e.getMessage(); + return results; } if (nobj > 1) { - // // TODO - is this an error? - // - std::cerr << "Error initializing Gurobi: More than one objective defined!" << std::endl; - return -1; + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message = "Error initializing Gurobi: More than one objective defined!"; + return results; } // Add Gurobi constraints @@ -422,47 +425,30 @@ int GurobiSolver::solve(CompactModel& compact_model) } } catch (GRBException e) { - std::cerr << "GUROBI Exception: (constraint) " << e.getMessage() << std::endl; - throw; + results->error_message + = "Gurobi Error: Caught gurobi exception while creating constraints " + e.getMessage(); + return results; } - std::cout << "OPTIMIZING GUROBI MODEL" << std::endl << std::flush; - - set_gurobi_options(); try { + set_gurobi_options(); gmodel->optimize(); - - int status = gmodel->get(GRB_IntAttr_Status); - if (status == GRB_OPTIMAL) { - // TODO: Are there other conditions where the variables have valid values? - // TODO: If we do not update the COEK variable values, should we set them to NAN? - // TODO: We need to cache the optimization status in COEK somewhere - // TODO: Is there a string description of the solver status? - -# if 0 - -WEH - What is the 'results object' for compact models? This is not defined yet. - - // Collect values of Gurobi variables - for (auto& var : model.variables) { - var.set_value(x[var.index].get(GRB_DoubleAttr_X)); - } -# endif - } } catch (GRBException e) { - std::cerr << "GUROBI Exception: (solver) " << e.getMessage() << std::endl; - // TODO: We should raise a CoekException object, to ensure that COEK can manage exceptions - // in a uniform manner. - throw; + results->error_message + = "Gurobi Error: Caught gurobi exception while optimizing " + e.getMessage(); + return results; } + collect_results(model, results); + delete gmodel; gmodel = 0; delete env; env = 0; - return 0; + results->toc(); + return results; } #endif @@ -620,6 +606,8 @@ void GurobiSolver::set_gurobi_options() // All options are converted to strings for Gurobi for (auto& it : string_options()) gmodel->set(it.first, it.second); + for (auto& it : boolean_options()) + gmodel->set(it.first, std::to_string(it.second)); for (auto& it : integer_options()) gmodel->set(it.first, std::to_string(it.second)); for (auto& it : double_options()) diff --git a/lib/coek/coek/solvers/highs/coek_highs.hpp b/lib/coek/coek/solvers/highs/coek_highs.hpp new file mode 100644 index 00000000..5e8b0654 --- /dev/null +++ b/lib/coek/coek/solvers/highs/coek_highs.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include "Highs.h" +#include "coek/solvers/solver_repn.hpp" + +namespace coek { + +class HighsSolver : public SolverRepn { + public: + Highs highs; + HighsModel hmodel; + HighsStatus return_status; + std::unordered_map x; + + public: + HighsSolver() : SolverRepn() {} + + std::shared_ptr resolve(); + std::shared_ptr solve(Model&); +#ifdef COEK_WITH_COMPACT_MODEL + std::shared_ptr solve(CompactModel&); +#endif + + public: + void set_solver_options(); + + protected: + void collect_results(Model& model, std::shared_ptr& results); +}; + +} // namespace coek diff --git a/lib/coek/coek/solvers/highs/highs.cpp b/lib/coek/coek/solvers/highs/highs.cpp new file mode 100644 index 00000000..b3cde6bb --- /dev/null +++ b/lib/coek/coek/solvers/highs/highs.cpp @@ -0,0 +1,614 @@ +#include +#include +#include +#include + +#include "../../ast/value_terms.hpp" +#include "coek/api/constraint.hpp" +#include "coek/api/expression.hpp" +#include "coek/api/objective.hpp" +#include "coek/compact/constraint_sequence.hpp" +#include "coek/compact/expression_sequence.hpp" +#include "coek/compact/objective_sequence.hpp" +#include "coek/model/model.hpp" +#include "coek/model/model_repn.hpp" +#include "coek/util/io_utils.hpp" + +#include "coek_highs.hpp" + +namespace coek { + +SolverRepn* create_highs_solver() { return new HighsSolver(); } + +namespace { + +bool add_variable(HighsModel& model, std::unordered_map& x, double lower, + double upper, const Variable& var) +{ + x[var.id()] = model.lp_.col_lower_.size(); + + if (var.is_binary() or var.is_integer()) { + model.lp_.col_lower_.push_back(lower); + model.lp_.col_upper_.push_back(upper); + model.lp_.integrality_.push_back(HighsVarType::kInteger); + return false; + } + + else { + if (upper >= 1e19) + upper = 1.0e30; + if (lower >= 1e19) + lower = -1.0e30; + model.lp_.col_lower_.push_back(lower); + model.lp_.col_upper_.push_back(upper); + model.lp_.integrality_.push_back(HighsVarType::kContinuous); + return true; + } +} + +void add_objective(HighsModel& model, Expression& expr, bool sense, + std::unordered_map& x, coek::QuadraticExpr& orepn) +{ + orepn.reset(); + orepn.collect_terms(expr); + model.lp_.col_cost_.resize(x.size()); + + if (orepn.linear_coefs.size() + orepn.quadratic_coefs.size() > 0) { + auto it = orepn.linear_coefs.begin(); + for (auto& var : orepn.linear_vars) { + model.lp_.col_cost_[x[var->index]] = *it; + ++it; + } + model.lp_.offset_ = orepn.constval; + +#if 0 + if (orepn.quadratic_coefs.size() == 0) + gmodel->setObjective(term1); + else { + GRBQuadExpr quadexpr; + for (size_t i = 0; i < orepn.quadratic_coefs.size(); i++) + quadexpr.addTerm(orepn.quadratic_coefs[i], x[orepn.quadratic_lvars[i]->index], + x[orepn.quadratic_rvars[i]->index]); + quadexpr.add(term1); + gmodel->setObjective(quadexpr); + } +#endif + } + + if (sense) + model.lp_.sense_ = ObjSense::kMinimize; + else + model.lp_.sense_ = ObjSense::kMaximize; +} + +void add_constraint(HighsModel& model, Constraint& con, std::unordered_map& x, + coek::QuadraticExpr& repn) +{ + repn.reset(); + repn.collect_terms(con); + + if (repn.quadratic_coefs.size() > 0) + throw std::runtime_error("Highs cannot handle nonlinear expressions in constraints"); + + if (repn.linear_coefs.size() == 0) + return; + + if (model.lp_.a_matrix_.start_.size() == 1) { + model.lp_.a_matrix_.start_.push_back(repn.linear_coefs.size()); + } + else { + HighsInt tmp = model.lp_.a_matrix_.start_.back(); + tmp += static_cast(repn.linear_coefs.size()); + model.lp_.a_matrix_.start_.push_back(tmp); + } + // Row coefs + auto it = repn.linear_coefs.begin(); + for (auto& var : repn.linear_vars) { + model.lp_.a_matrix_.index_.push_back(static_cast(x[var->index])); + model.lp_.a_matrix_.value_.push_back(*it); + ++it; + } + + // Row lower/upper + double lower = -1.0e30; + double upper = 1.0e30; + if (con.is_inequality()) { + if (con.lower().repn) { + double tmp = con.lower().value(); + if (tmp > -COEK_INFINITY) + lower = -repn.constval + tmp; + } + if (con.upper().repn) { + double tmp = con.upper().value(); + if (tmp < COEK_INFINITY) + upper = -repn.constval + tmp; + } + } + else + lower = upper = -repn.constval + con.lower().value(); + model.lp_.row_lower_.push_back(lower); + model.lp_.row_upper_.push_back(upper); +} + +} // namespace + +std::shared_ptr HighsSolver::solve(Model& coek_model) +{ + auto results = std::make_shared(); + results->solver_name = "highs"; + results->termination_condition = TerminationCondition::error; + results->tic(); + + auto _coek_model = coek_model.repn.get(); + assert(_coek_model->objectives.size() == 1); + + hmodel.clear(); + + // Add variables + bool continuous = true; + for (auto& var : _coek_model->variables) { + if (not var.fixed()) + continuous = continuous and add_variable(hmodel, x, var.lower(), var.upper(), var); + } + if (continuous) + hmodel.lp_.integrality_.resize(0); + + // Add objective + unsigned int nobj = 0; + try { + coek::QuadraticExpr orepn; + for (auto& obj : _coek_model->objectives) { + Expression tmp = obj.expr(); + add_objective(hmodel, tmp, obj.sense(), x, orepn); + nobj++; + } + } + catch (const std::exception& e) { + results->error_message = "Highs Error: Caught highs exception while creating objectives " + + std::string(e.what()); + return results; + } + if (nobj > 1) { + // TODO - is this an error? + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message = "Error initializing Highs: More than one objective defined!"; + return results; + } + + // Add constraints + hmodel.lp_.a_matrix_.format_ = MatrixFormat::kRowwise; + try { + coek::QuadraticExpr repn; + for (auto& con : _coek_model->constraints) { + add_constraint(hmodel, con, x, repn); + } + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught exception while creating constraints " + std::string(e.what()); + return results; + } + + hmodel.lp_.num_col_ = hmodel.lp_.col_cost_.size(); + hmodel.lp_.num_row_ = hmodel.lp_.row_lower_.size(); + +#if 0 + std::cout << "Ncol: " << hmodel.lp_.num_col_ << std::endl; + std::cout << "Nrow: " << hmodel.lp_.num_row_ << std::endl; + std::cout << "ColCost: " << hmodel.lp_.col_cost_ << std::endl; + std::cout << "Offset: " << hmodel.lp_.offset_ << std::endl; + std::cout << "Integrality: "; + for (auto& val : hmodel.lp_.integrality_) + std::cout << (int)val << " "; + std::cout << std::endl; + std::cout << "Lower: " << hmodel.lp_.row_lower_ << std::endl; + std::cout << "Upper: " << hmodel.lp_.row_upper_ << std::endl; + std::cout << "Start: " << hmodel.lp_.a_matrix_.start_ << std::endl; + std::cout << "Index: " << hmodel.lp_.a_matrix_.index_ << std::endl; + std::cout << "Value: " << hmodel.lp_.a_matrix_.value_ << std::endl; +#endif + + set_solver_options(); + return_status = highs.passModel(hmodel); + if (return_status != HighsStatus::kOk) { + results->error_message = "Highs Error: Error initializing model"; + return results; + } + + try { + return_status = highs.run(); + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + return results; + } + + collect_results(coek_model, results); + + results->toc(); + return results; +} + +#ifdef COEK_WITH_COMPACT_MODEL +std::shared_ptr HighsSolver::solve(CompactModel& compact_model) +{ + auto results = std::make_shared(); + results->solver_name = "highs"; + results->termination_condition = TerminationCondition::error; + results->tic(); + + // Add variables + bool continuous = true; + for (auto& val : compact_model.repn->variables) { + if (auto eval = std::get_if(&val)) { + Expression lb = eval->lower_expression().expand(); + auto lb_ = lb.value(); + Expression ub = eval->upper_expression().expand(); + auto ub_ = ub.value(); + continuous = continuous and add_variable(hmodel, x, lb_, ub_, *eval); + model.add(*eval); + } + else { + auto& seq = std::get(val); + for (auto jt = seq.begin(); jt != seq.end(); ++jt) { + continuous = continuous and add_variable(hmodel, x, jt->lower(), jt->upper(), *jt); + model.add(*jt); + } + } + } + if (continuous) + hmodel.lp_.integrality_.resize(0); + + // Add objective + unsigned int nobj = 0; + try { + coek::QuadraticExpr orepn; + for (auto& val : compact_model.repn->objectives) { + if (auto eval = std::get_if(&val)) { + Expression tmp = eval->expr().expand(); + add_objective(hmodel, tmp, eval->sense(), x, orepn); + nobj++; + } + else { + auto& seq = std::get(val); + for (auto jt = seq.begin(); jt != seq.end(); ++jt) { + // model.repn->objectives.push_back(*jt); + Expression tmp = jt->expr(); + add_objective(hmodel, tmp, jt->sense(), x, orepn); + nobj++; + } + } + } + } + catch (const std::exception& e) { + results->error_message = "Highs Error: Caught highs exception while creating objectives " + + std::string(e.what()); + return results; + } + if (nobj > 1) { + // TODO - is this an error? + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message = "Error initializing Highs: More than one objective defined!"; + return results; + } + + // Add constraints + hmodel.lp_.a_matrix_.format_ = MatrixFormat::kRowwise; + try { + coek::QuadraticExpr repn; + for (auto& val : compact_model.repn->constraints) { + if (auto cval = std::get_if(&val)) { + Constraint c = cval->expand(); + add_constraint(hmodel, c, x, repn); + } + else { + auto& seq = std::get(val); + for (auto jt = seq.begin(); jt != seq.end(); ++jt) { + add_constraint(hmodel, *jt, x, repn); + } + } + } + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught exception while creating constraints " + std::string(e.what()); + return results; + } + + hmodel.lp_.num_col_ = hmodel.lp_.col_cost_.size(); + hmodel.lp_.num_row_ = hmodel.lp_.row_lower_.size(); + + set_solver_options(); + return_status = highs.passModel(hmodel); + if (return_status != HighsStatus::kOk) { + results->error_message = "Highs Error: Error initializing model"; + return results; + } + + try { + return_status = highs.run(); + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + return results; + } + + collect_results(model, results); + + results->toc(); + return results; +} +#endif + +std::shared_ptr HighsSolver::resolve() +{ + auto results = std::make_shared(); + results->solver_name = "highs"; + results->termination_condition = TerminationCondition::error; + results->tic(); + +#if 0 + auto _model = model.repn.get(); + + if (initial_solve()) { + env = new GRBEnv(true); + auto it = integer_options().find("OutputFlag"); + if (it != integer_options().end()) + env->set(GRB_IntParam_OutputFlag, it->second); + env->start(); + gmodel = new GRBModel(*env); + + assert(_model->objectives.size() == 1); + + // Add variables + for (auto& var : _model->variables) { + std::shared_ptr v = var.repn; + if (not v->fixed) { + double lb = v->lb->eval(); + double ub = v->ub->eval(); + if (v->binary) + x[v->index] = gmodel->addVar(lb, ub, 0, GRB_BINARY); + else if (v->integer) + x[v->index] = gmodel->addVar(lb, ub, 0, GRB_INTEGER); + else { + if (ub >= 1e19) { + if (lb <= -1e19) + x[v->index] + = gmodel->addVar(-GRB_INFINITY, GRB_INFINITY, 0, GRB_CONTINUOUS); + else + x[v->index] = gmodel->addVar(lb, GRB_INFINITY, 0, GRB_CONTINUOUS); + } + else { + if (lb <= -1e19) + x[v->index] = gmodel->addVar(-GRB_INFINITY, ub, 0, GRB_CONTINUOUS); + else + x[v->index] = gmodel->addVar(lb, ub, 0, GRB_CONTINUOUS); + } + } + } + } + + gmodel->update(); + + // Add objective + int nobj = 0; + try { + coek::QuadraticExpr orepn; + for (auto& obj : model.repn->objectives) { + Expression tmp = obj.expr(); + add_objective(gmodel, tmp, obj.sense(), x, orepn); + nobj++; + } + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while creating objectives " + + std::string(e.what()); + return results; + } + if (nobj > 1) { + // TODO - is this an error? + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message = "Error initializing Highs: More than one objective defined!"; + return results; + } + + // Add constraints + try { + coek::QuadraticExpr repn; + for (auto& con : model.repn->constraints) { + add_constraint(gmodel, con, x, repn); + } + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while creating constraints " + + std::string(e.what()); + return results; + } + + set_solver_options(); + try { + gmodel->optimize(); + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + return results; + } + } + + else { + for (auto it = updated_coefs.begin(); it != updated_coefs.end(); ++it) { + size_t i = std::get<0>(*it); + size_t where = std::get<1>(*it); + size_t j = std::get<2>(*it); + + switch (where) { + case 0: // Constant Value + if (i > 0) + gmodel->getConstr(i - 1).set(GRB_DoubleAttr_RHS, -repn[i].constval->eval()); + else + gmodel->set(GRB_DoubleAttr_ObjCon, repn[0].constval->eval()); + break; + + case 1: // Linear Coef + if (i > 0) + gmodel->chgCoeff(gmodel->getConstr(i - 1), x[repn[i].linear_vars[j]->index], + repn[i].linear_coefs[j]->eval()); + else + x[repn[0].linear_vars[j]->index].set(GRB_DoubleAttr_Obj, + repn[0].linear_coefs[j]->eval()); + break; + + case 2: // Quadratic Coef + break; + + case 3: // Nonlinear terms + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message + = "Error initializing Highs: Cannot optimize models with nonlinear terms."; + return results; + break; + }; + } + + try { + gmodel->optimize(); + } + catch (const std::exception& e) { + results->error_message + = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + return results; + } + } + + // Collect values of solver variables + collect_results(model, results); + +#endif + results->toc(); + return results; +} + +void HighsSolver::set_solver_options() +{ + for (auto& it : string_options()) + highs.setOptionValue(it.first, it.second); + for (auto& it : boolean_options()) + highs.setOptionValue(it.first, it.second); + for (auto& it : integer_options()) + highs.setOptionValue(it.first, it.second); + for (auto& it : double_options()) + highs.setOptionValue(it.first, it.second); +} + +void HighsSolver::collect_results(Model& model, std::shared_ptr& results) +{ + try { + const auto& model_status = highs.getModelStatus(); + + if (model_status == HighsModelStatus::kOptimal) { + results->termination_condition = TerminationCondition::convergence_criteria_satisfied; + results->solution_status = SolutionStatus::optimal; + const auto& solution = highs.getSolution(); + const auto& info = highs.getInfo(); + + // Objective + results->objective_value = info.objective_function_value; + if (highs.getModel().lp_.isMip()) + results->objective_value = info.mip_dual_bound; + else + results->objective_bound = results->objective_value; + // highs.getOptionValue("objective_bound", results->objective_value); + + // Collect values of variables + const bool has_primal = info.primal_solution_status; + if (has_primal) { + for (auto& var : model.repn->variables) { + if (not var.fixed()) { + var.value(solution.col_value[x[var.id()]]); + } + } + + // TODO: collect duals and basis if the user requests this information + } + } + else if (model_status == HighsModelStatus::kNotset) { + results->termination_condition = TerminationCondition::empty_model; + } + else if (model_status == HighsModelStatus::kLoadError) { + results->termination_condition = TerminationCondition::invalid_model_for_solver; + results->error_message = "Highs terminated with an error loading the model."; + } + else if (model_status == HighsModelStatus::kModelError) { + results->termination_condition = TerminationCondition::error; + results->error_message = "Highs terminated with a model error."; + } + else if (model_status == HighsModelStatus::kPresolveError) { + results->termination_condition = TerminationCondition::error; + results->error_message = "Highs terminated with a pre-solve error."; + } + else if (model_status == HighsModelStatus::kSolveError) { + results->termination_condition = TerminationCondition::error; + results->error_message = "Highs terminated with a solve error."; + } + else if (model_status == HighsModelStatus::kPostsolveError) { + results->termination_condition = TerminationCondition::error; + results->error_message = "Highs terminated with a post-solve error."; + } + else if (model_status == HighsModelStatus::kModelEmpty) { + results->termination_condition = TerminationCondition::empty_model; + } + else if (model_status == HighsModelStatus::kInfeasible) { + results->termination_condition = TerminationCondition::proven_infeasible; + results->solution_status = SolutionStatus::infeasible; + } + else if (model_status == HighsModelStatus::kUnboundedOrInfeasible) { + results->termination_condition = TerminationCondition::infeasible_or_unbounded; + } + else if (model_status == HighsModelStatus::kUnbounded) { + results->termination_condition = TerminationCondition::unbounded; + } + else if (model_status == HighsModelStatus::kObjectiveBound) { + results->termination_condition = TerminationCondition::objective_limit; + results->error_message + = "Highs terminated because optimal objective for model was proven to be worse " + "than the value specified in the objective_bound option."; + } + else if (model_status == HighsModelStatus::kObjectiveTarget) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Highs terminated because optimal objective for model was proven to be worse " + "than the value specified in the objective_target option."; + } + else if (model_status == HighsModelStatus::kTimeLimit) { + results->termination_condition = TerminationCondition::time_limit; + } + else if (model_status == HighsModelStatus::kIterationLimit) { + results->termination_condition = TerminationCondition::iteration_limit; + } + else if (model_status == HighsModelStatus::kSolutionLimit) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Highs terminated because the number of solutions found reached the value " + "specified in the SolutionLimit parameter."; + } + else if (model_status == HighsModelStatus::kInterrupt) { + results->termination_condition = TerminationCondition::interrupted; + } + else if (model_status == HighsModelStatus::kUnknown) { + results->termination_condition = TerminationCondition::unknown; + results->error_message + = "Highs Error: Optimization was terminated due to unknown difficulties."; + } + } + catch (const std::exception& e) { + results->termination_condition = TerminationCondition::unknown; + results->error_message = "HIGHS Exception: (results) " + std::string(e.what()); + } +} + +} // namespace coek diff --git a/lib/coek/coek/solvers/ipopt/ipopt_capi.cpp b/lib/coek/coek/solvers/ipopt/ipopt_capi.cpp index b58e8ada..daff3349 100644 --- a/lib/coek/coek/solvers/ipopt/ipopt_capi.cpp +++ b/lib/coek/coek/solvers/ipopt/ipopt_capi.cpp @@ -576,11 +576,13 @@ class IpoptSolverRepn_CAPI : public IpoptSolverRepn { void set_start_from_last_x(bool flag) { nlp->start_from_last_x = flag; } void set_options(const std::map& string_options, + const std::map& boolean_options, const std::map& integer_options, const std::map& double_options); }; void IpoptSolverRepn_CAPI::set_options(const std::map& string_options, + const std::map& boolean_options, const std::map& integer_options, const std::map& double_options) { @@ -589,6 +591,10 @@ void IpoptSolverRepn_CAPI::set_options(const std::map& char* tmp2 = const_cast(it->second.c_str()); (*AddIpoptStrOption_func_ptr)(nlp->app, tmp1, tmp2); } + for (auto it = boolean_options.begin(); it != boolean_options.end(); ++it) { + char* tmp1 = const_cast(it->first.c_str()); + (*AddIpoptIntOption_func_ptr)(nlp->app, tmp1, it->second); + } for (auto it = integer_options.begin(); it != integer_options.end(); ++it) { char* tmp1 = const_cast(it->first.c_str()); (*AddIpoptIntOption_func_ptr)(nlp->app, tmp1, it->second); diff --git a/lib/coek/coek/solvers/ipopt/ipopt_solver.cpp b/lib/coek/coek/solvers/ipopt/ipopt_solver.cpp index 1158bb86..2b922a1b 100644 --- a/lib/coek/coek/solvers/ipopt/ipopt_solver.cpp +++ b/lib/coek/coek/solvers/ipopt/ipopt_solver.cpp @@ -9,6 +9,8 @@ namespace coek { +NLPSolverRepn* create_ipopt_solver() { return new IpoptSolver(); } + std::shared_ptr IpoptSolver::solve(NLPModel& _model) { if (not available_) { @@ -29,7 +31,7 @@ std::shared_ptr IpoptSolver::solve(NLPModel& _model) "time."; return res; } - repn->set_options(string_options(), integer_options(), double_options()); + repn->set_options(string_options(), boolean_options(), integer_options(), double_options()); return repn->perform_solve(); } @@ -46,7 +48,7 @@ std::shared_ptr IpoptSolver::resolve_exec() if (not initial_solve()) model->reset(); - repn->set_options(string_options(), integer_options(), double_options()); + repn->set_options(string_options(), boolean_options(), integer_options(), double_options()); auto it = string_options().find("warm_start_init_point"); if ((it != string_options().end()) and (it->second == "yes")) repn->set_start_from_last_x(true); diff --git a/lib/coek/coek/solvers/ipopt/ipopt_solver.hpp b/lib/coek/coek/solvers/ipopt/ipopt_solver.hpp index 61444af7..bbc6d83f 100644 --- a/lib/coek/coek/solvers/ipopt/ipopt_solver.hpp +++ b/lib/coek/coek/solvers/ipopt/ipopt_solver.hpp @@ -14,6 +14,7 @@ class IpoptSolverRepn { virtual std::shared_ptr perform_solve() = 0; virtual void set_start_from_last_x(bool flag) = 0; virtual void set_options(const std::map& string_options, + const std::map& boolean_options, const std::map& integer_options, const std::map& double_options) = 0; diff --git a/lib/coek/coek/solvers/solver_repn.cpp b/lib/coek/coek/solvers/solver_repn.cpp index 8f14c314..41d5e37f 100644 --- a/lib/coek/coek/solvers/solver_repn.cpp +++ b/lib/coek/coek/solvers/solver_repn.cpp @@ -9,11 +9,16 @@ #include "coek/model/model_repn.hpp" #include "coek/model/nlp_model.hpp" #include "coek/solvers/solver_repn.hpp" +/* #include "coek/solvers/ipopt/ipopt_solver.hpp" #include "testsolver.hpp" #ifdef WITH_GUROBI # include "coek/solvers/gurobi/coek_gurobi.hpp" #endif +#ifdef WITH_GUROBI +# include "coek/solvers/gurobi/coek_highs.hpp" +#endif +*/ namespace coek { @@ -52,6 +57,7 @@ void SolverCache::reset_cache() pcache.clear(); } +/* SolverRepn* create_solver(std::string& name, OptionCache& options) { if (name == "test") @@ -78,6 +84,7 @@ NLPSolverRepn* create_nlpsolver(std::string& name, OptionCache& options) return 0; } +*/ std::shared_ptr NLPSolverRepn::resolve(bool reset_nlpmodel) { @@ -187,7 +194,7 @@ void SolverRepn::load(CompactModel& _model) load(model); } -int SolverRepn::solve(CompactModel& _model) +std::shared_ptr SolverRepn::solve(CompactModel& _model) { model = _model.expand(); return solve(model); diff --git a/lib/coek/coek/solvers/solver_results.hpp b/lib/coek/coek/solvers/solver_results.hpp index 40402377..2d3a0a59 100644 --- a/lib/coek/coek/solvers/solver_results.hpp +++ b/lib/coek/coek/solvers/solver_results.hpp @@ -194,10 +194,16 @@ std::string to_string(const SolverResults& res); namespace std { -inline ostream& operator<<(const coek::SolverResults& res, ostream& os) +inline std::ostream& operator<<(std::ostream& os, const coek::SolverResults& res) { os << coek::to_string(res); return os; } +inline ostream& operator<<(ostream& os, const std::shared_ptr& res) +{ + os << coek::to_string(*res); + return os; +} + } // namespace std diff --git a/lib/coek/coek/solvers/testsolver.cpp b/lib/coek/coek/solvers/testsolver.cpp index 3a5912b8..415dd97a 100644 --- a/lib/coek/coek/solvers/testsolver.cpp +++ b/lib/coek/coek/solvers/testsolver.cpp @@ -11,6 +11,8 @@ namespace coek { +SolverRepn* create_coektest_solver() { return new TestSolver(); } + std::shared_ptr TestSolver::solve(Model& model) { assert(initial_solve()); diff --git a/lib/coek/coek/util/option_cache.cpp b/lib/coek/coek/util/option_cache.cpp index a96db3c7..abb7878f 100644 --- a/lib/coek/coek/util/option_cache.cpp +++ b/lib/coek/coek/util/option_cache.cpp @@ -9,6 +9,7 @@ namespace coek { class OptionCacheRepn { public: std::map string_options; + std::map boolean_options; std::map integer_options; std::map double_options; }; @@ -29,6 +30,13 @@ const std::map& OptionCache::string_options() const return options->string_options; } +std::map& OptionCache::boolean_options() { return options->boolean_options; } + +const std::map& OptionCache::boolean_options() const +{ + return options->boolean_options; +} + std::map& OptionCache::integer_options() { return options->integer_options; } const std::map& OptionCache::integer_options() const @@ -43,6 +51,15 @@ const std::map& OptionCache::double_options() const return options->double_options; } +bool OptionCache::get_option(const std::string& option, bool& value) const +{ + auto it = options->boolean_options.find(option); + if (it == options->boolean_options.end()) + return false; + value = it->second; + return true; +} + bool OptionCache::get_option(const std::string& option, int& value) const { auto it = options->integer_options.find(option); @@ -70,6 +87,11 @@ bool OptionCache::get_option(const std::string& option, std::string& value) cons return true; } +void OptionCache::set_option(const std::string& option, bool value) +{ + options->boolean_options[option] = value; +} + void OptionCache::set_option(const std::string& option, int value) { options->integer_options[option] = value; @@ -80,7 +102,12 @@ void OptionCache::set_option(const std::string& option, double value) options->double_options[option] = value; } -void OptionCache::set_option(const std::string& option, const std::string value) +void OptionCache::set_option(const std::string& option, const std::string& value) +{ + options->string_options[option] = value; +} + +void OptionCache::set_option(const std::string& option, const char* value) { options->string_options[option] = value; } diff --git a/lib/coek/coek/util/option_cache.hpp b/lib/coek/coek/util/option_cache.hpp index df7ee1f0..5ebd21cb 100644 --- a/lib/coek/coek/util/option_cache.hpp +++ b/lib/coek/coek/util/option_cache.hpp @@ -17,11 +17,28 @@ class OptionCache { std::map& string_options(); const std::map& string_options() const; + + std::map& boolean_options(); + const std::map& boolean_options() const; + std::map& integer_options(); const std::map& integer_options() const; + std::map& double_options(); const std::map& double_options() const; + /** Get the value of a boolean option + * + * The option value is returned by reference if it has + * a value. + * + * \param option the option name + * \param value a boolean value that is passed by reference + * + * \returns \c true if the option is found + */ + virtual bool get_option(const std::string& option, bool& value) const; + /** Get the value of an integer option * * The option value is returned by reference if it has @@ -58,6 +75,13 @@ class OptionCache { */ virtual bool get_option(const std::string& option, std::string& value) const; + /** Set a boolean option + * + * \param option the option name + * \param value the boolean value + */ + virtual void set_option(const std::string& option, bool value); + /** Set an integer option * * \param option the option name @@ -77,7 +101,14 @@ class OptionCache { * \param option the option name * \param value the string value */ - virtual void set_option(const std::string& option, const std::string value); + virtual void set_option(const std::string& option, const std::string& value); + + /** Set a string option + * + * \param option the option name + * \param value the string value + */ + virtual void set_option(const std::string& option, const char* value); }; } // namespace coek diff --git a/lib/coek/test/test_examples.cpp b/lib/coek/test/test_examples.cpp index 32676d97..30c4f758 100644 --- a/lib/coek/test/test_examples.cpp +++ b/lib/coek/test/test_examples.cpp @@ -6,6 +6,7 @@ #include "catch2/generators/catch_generators.hpp" #include "coek/ast/base_terms.hpp" #include "coek/coek.hpp" +#include "coek/solvers/solver_results.hpp" // // EXAMPLES @@ -327,3 +328,30 @@ TEST_CASE("gurobi_examples", "[smoke]") } } } + +TEST_CASE("highs_examples", "[smoke]") +{ + coek::Solver solver("highs"); + if (solver.available()) { + SECTION("simplelp1") + { + auto m = simplelp1(); + solver.set_option("output_flag", false); + auto res = solver.solve(m); + // std::cout << res << std::endl; + REQUIRE(coek::check_optimal_termination(res)); + + REQUIRE(res->objective_value == Catch::Approx(28750.0)); + REQUIRE(res->objective_bound == Catch::Approx(28750.0)); + check(m.get_variables(), simplelp1_soln); + } + } + else { + SECTION("simplelp1") + { + auto m = simplelp1(); + auto res = solver.solve(m); + REQUIRE(res->termination_condition == coek::TerminationCondition::solver_not_available); + } + } +} diff --git a/lib/coek/test/test_indexed.cpp b/lib/coek/test/test_indexed.cpp index aac8d3e7..b81cd7c1 100644 --- a/lib/coek/test/test_indexed.cpp +++ b/lib/coek/test/test_indexed.cpp @@ -1,6 +1,9 @@ #include #include "catch2/catch_test_macros.hpp" +#include "catch2/catch_approx.hpp" +#include "catch2/matchers/catch_matchers.hpp" + #include "coek/ast/base_terms.hpp" #include "coek/ast/value_terms.hpp" #include "coek/coek.hpp"