From 4505ba5ec3a6b65263a5d4ca4b9d170457cbf413 Mon Sep 17 00:00:00 2001 From: whart222 Date: Tue, 7 May 2024 11:39:22 -0600 Subject: [PATCH] Several changes 1. Adding a few more test cases 2. Adding support for quadratic objectives in Highs 3. Adding Variable::fix() --- lib/coek/coek/api/expression.cpp | 6 + lib/coek/coek/api/expression.hpp | 2 + lib/coek/coek/solvers/gurobi/gurobi.cpp | 308 ++++++++++----------- lib/coek/coek/solvers/highs/coek_highs.hpp | 4 + lib/coek/coek/solvers/highs/highs.cpp | 115 +++++--- lib/coek/test/solver/TestModels.cpp | 109 +++++++- lib/coek/test/solver/test_examples.cpp | 6 +- lib/coek/test/solver/test_gurobi.cpp | 33 ++- lib/coek/test/solver/test_highs.cpp | 35 ++- 9 files changed, 405 insertions(+), 213 deletions(-) diff --git a/lib/coek/coek/api/expression.cpp b/lib/coek/coek/api/expression.cpp index 051a7926..1fcfa77e 100644 --- a/lib/coek/coek/api/expression.cpp +++ b/lib/coek/coek/api/expression.cpp @@ -224,6 +224,12 @@ Variable& Variable::fix(double value) return *this; } +Variable& Variable::fix() +{ + repn->fixed = true; + return *this; +} + Variable& Variable::fixed(bool _flag) { repn->fixed = _flag; diff --git a/lib/coek/coek/api/expression.hpp b/lib/coek/coek/api/expression.hpp index 9a268ed1..84271b82 100644 --- a/lib/coek/coek/api/expression.hpp +++ b/lib/coek/coek/api/expression.hpp @@ -241,6 +241,8 @@ class Variable { Variable& fixed(bool flag); /** Set the variable value and declare the variable fixed. \returns the variable object. */ Variable& fix(double value); + /** Set the variable value and declare the variable fixed. \returns the variable object. */ + Variable& fix(); /** \returns \c true if the variable is fixed */ bool fixed() const; diff --git a/lib/coek/coek/solvers/gurobi/gurobi.cpp b/lib/coek/coek/solvers/gurobi/gurobi.cpp index 414b79cb..6e81339f 100644 --- a/lib/coek/coek/solvers/gurobi/gurobi.cpp +++ b/lib/coek/coek/solvers/gurobi/gurobi.cpp @@ -47,8 +47,8 @@ void add_gurobi_objective(GRBModel* gmodel, Expression& expr, bool sense, std::unordered_map& x, coek::QuadraticExpr& orepn) { orepn.reset(); - // coek::QuadraticExpr orepn; orepn.collect_terms(expr); + if (orepn.linear_coefs.size() + orepn.quadratic_coefs.size() > 0) { GRBLinExpr term1; auto iv = orepn.linear_vars.begin(); @@ -135,128 +135,46 @@ GurobiSolver::~GurobiSolver() delete env; } -void GurobiSolver::collect_results(Model& model, std::shared_ptr& results) +void GurobiSolver::set_gurobi_options() { - try { - int status = gmodel->get(GRB_IntAttr_Status); - if (status == GRB_OPTIMAL) { - results->termination_condition = TerminationCondition::convergence_criteria_satisfied; - results->solution_status = SolutionStatus::optimal; - results->objective_value = gmodel->getObjective().getValue(); - try { - double value = gmodel->get(GRB_DoubleAttr_ObjBound); - results->objective_bound = value; - } - catch (GRBException) { - } - if (not results->objective_bound.has_value()) { - try { - double value = gmodel->get(GRB_DoubleAttr_ObjBoundC); - results->objective_bound = value; - } - catch (GRBException) { - } - } + // 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()) + gmodel->set(it.first, std::to_string(it.second)); +} - // Collect values of Gurobi variables - for (auto& var : model.repn->variables) { - std::shared_ptr& v = var.repn; - if (not v->fixed) { - v->set_value(x[v->index].get(GRB_DoubleAttr_X)); - } - } - } - else if (status == GRB_SUBOPTIMAL) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->solution_status = SolutionStatus::feasible; - results->objective_value = gmodel->getObjective().getValue(); - try { - double value = gmodel->get(GRB_DoubleAttr_ObjBound); - results->objective_bound = value; - } - catch (GRBException) { - } - if (not results->objective_bound.has_value()) { - try { - double value = gmodel->get(GRB_DoubleAttr_ObjBoundC); - results->objective_bound = value; - } - catch (GRBException) { - } - } - results->error_message - = "Unable to satisfy optimality tolerances; a sub-optimal solution is available"; +void GurobiSolver::pre_solve() +{ + results = std::make_shared(); + results->solver_name = "gurobi"; + results->termination_condition = TerminationCondition::error; + results->tic(); - // Collect values of Gurobi variables - for (auto& var : model.repn->variables) { - std::shared_ptr& v = var.repn; - if (not v->fixed) { - v->set_value(x[v->index].get(GRB_DoubleAttr_X)); - } - } - } - else if (status == GRB_INFEASIBLE) { - results->termination_condition = TerminationCondition::proven_infeasible; - results->solution_status = SolutionStatus::infeasible; - } - else if (status == GRB_UNBOUNDED) { - results->termination_condition = TerminationCondition::unbounded; - } - else if (status == GRB_ITERATION_LIMIT) { - results->termination_condition = TerminationCondition::iteration_limit; - } - else if (status == GRB_TIME_LIMIT) { - results->termination_condition = TerminationCondition::time_limit; - } - else if (status == GRB_INTERRUPTED) { - results->termination_condition = TerminationCondition::interrupted; - } - else if (status == GRB_USER_OBJ_LIMIT) { - results->termination_condition = TerminationCondition::objective_limit; - } - else if (status == GRB_WORK_LIMIT) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->error_message - = "Gurobi terminated because the work expended exceeded the value specified in the " - "WorkLimit parameter."; - } - else if (status == GRB_MEM_LIMIT) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->error_message - = "Gurobi terminated because the total amount of allocated memory exceeded the " - "value specified in the SoftMemLimit parameter."; - } - else if (status == GRB_NODE_LIMIT) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->error_message - = "Gurobi terminated because the total number of branch-and-cut nodes explored " - "exceeded the value specified in the NodeLimit parameter."; - } - else if (status == GRB_SOLUTION_LIMIT) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->error_message - = "Gurobi terminated because the number of solutions found reached the value " - "specified in the SolutionLimit parameter."; - } - else if (status == GRB_CUTOFF) { - results->termination_condition = TerminationCondition::other_termination_limit; - results->error_message - = "Gurobi terminated because optimal objective for model was proven to be worse " - "than the value specified in the Cutoff parameter."; - } - else if (status == GRB_NUMERIC) { - results->termination_condition = TerminationCondition::unknown; - results->error_message - = "Gurobi Error: Optimization was terminated due to unrecoverable numerical " - "difficulties."; - } - } - catch (GRBException e) { - results->termination_condition = TerminationCondition::unknown; - results->error_message = "GUROBI Exception: (results) " + e.getMessage(); + 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); } } +void GurobiSolver::post_solve() +{ + delete gmodel; + gmodel = 0; + delete env; + env = 0; + + results->toc(); +} + std::shared_ptr GurobiSolver::solve(Model& model) { pre_solve(); @@ -334,33 +252,6 @@ std::shared_ptr GurobiSolver::solve(Model& model) return results; } -void GurobiSolver::pre_solve() -{ - results = std::make_shared(); - results->solver_name = "gurobi"; - results->termination_condition = TerminationCondition::error; - results->tic(); - - 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); - } -} - -void GurobiSolver::post_solve() -{ - delete gmodel; - gmodel = 0; - delete env; - env = 0; - - results->toc(); -} - #ifdef COEK_WITH_COMPACT_MODEL std::shared_ptr GurobiSolver::solve(CompactModel& compact_model) { @@ -615,17 +506,126 @@ std::shared_ptr GurobiSolver::resolve() return results; } -void GurobiSolver::set_gurobi_options() +void GurobiSolver::collect_results(Model& model, std::shared_ptr& results) { - // 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()) - gmodel->set(it.first, std::to_string(it.second)); + try { + int status = gmodel->get(GRB_IntAttr_Status); + if (status == GRB_OPTIMAL) { + results->termination_condition = TerminationCondition::convergence_criteria_satisfied; + results->solution_status = SolutionStatus::optimal; + results->objective_value = gmodel->getObjective().getValue(); + try { + double value = gmodel->get(GRB_DoubleAttr_ObjBound); + results->objective_bound = value; + } + catch (GRBException) { + } + if (not results->objective_bound.has_value()) { + try { + double value = gmodel->get(GRB_DoubleAttr_ObjBoundC); + results->objective_bound = value; + } + catch (GRBException) { + } + } + + // Collect values of Gurobi variables + for (auto& var : model.repn->variables) { + std::shared_ptr& v = var.repn; + if (not v->fixed) { + v->set_value(x[v->index].get(GRB_DoubleAttr_X)); + } + } + } + else if (status == GRB_SUBOPTIMAL) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->solution_status = SolutionStatus::feasible; + results->objective_value = gmodel->getObjective().getValue(); + try { + double value = gmodel->get(GRB_DoubleAttr_ObjBound); + results->objective_bound = value; + } + catch (GRBException) { + } + if (not results->objective_bound.has_value()) { + try { + double value = gmodel->get(GRB_DoubleAttr_ObjBoundC); + results->objective_bound = value; + } + catch (GRBException) { + } + } + results->error_message + = "Unable to satisfy optimality tolerances; a sub-optimal solution is available"; + + // Collect values of Gurobi variables + for (auto& var : model.repn->variables) { + std::shared_ptr& v = var.repn; + if (not v->fixed) { + v->set_value(x[v->index].get(GRB_DoubleAttr_X)); + } + } + } + else if (status == GRB_INFEASIBLE) { + results->termination_condition = TerminationCondition::proven_infeasible; + results->solution_status = SolutionStatus::infeasible; + } + else if (status == GRB_UNBOUNDED) { + results->termination_condition = TerminationCondition::unbounded; + } + else if (status == GRB_ITERATION_LIMIT) { + results->termination_condition = TerminationCondition::iteration_limit; + } + else if (status == GRB_TIME_LIMIT) { + results->termination_condition = TerminationCondition::time_limit; + } + else if (status == GRB_INTERRUPTED) { + results->termination_condition = TerminationCondition::interrupted; + } + else if (status == GRB_USER_OBJ_LIMIT) { + results->termination_condition = TerminationCondition::objective_limit; + } + else if (status == GRB_WORK_LIMIT) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Gurobi terminated because the work expended exceeded the value specified in the " + "WorkLimit parameter."; + } + else if (status == GRB_MEM_LIMIT) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Gurobi terminated because the total amount of allocated memory exceeded the " + "value specified in the SoftMemLimit parameter."; + } + else if (status == GRB_NODE_LIMIT) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Gurobi terminated because the total number of branch-and-cut nodes explored " + "exceeded the value specified in the NodeLimit parameter."; + } + else if (status == GRB_SOLUTION_LIMIT) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Gurobi terminated because the number of solutions found reached the value " + "specified in the SolutionLimit parameter."; + } + else if (status == GRB_CUTOFF) { + results->termination_condition = TerminationCondition::other_termination_limit; + results->error_message + = "Gurobi terminated because optimal objective for model was proven to be worse " + "than the value specified in the Cutoff parameter."; + } + else if (status == GRB_NUMERIC) { + results->termination_condition = TerminationCondition::unknown; + results->error_message + = "Gurobi Error: Optimization was terminated due to unrecoverable numerical " + "difficulties."; + } + } + catch (GRBException e) { + results->termination_condition = TerminationCondition::unknown; + results->error_message = "GUROBI Exception: (results) " + e.getMessage(); + } } } // namespace coek diff --git a/lib/coek/coek/solvers/highs/coek_highs.hpp b/lib/coek/coek/solvers/highs/coek_highs.hpp index 5e8b0654..07dd93d4 100644 --- a/lib/coek/coek/solvers/highs/coek_highs.hpp +++ b/lib/coek/coek/solvers/highs/coek_highs.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include "Highs.h" #include "coek/solvers/solver_repn.hpp" @@ -10,6 +11,7 @@ class HighsSolver : public SolverRepn { Highs highs; HighsModel hmodel; HighsStatus return_status; + std::shared_ptr results; std::unordered_map x; public: @@ -26,6 +28,8 @@ class HighsSolver : public SolverRepn { protected: void collect_results(Model& model, std::shared_ptr& results); + void pre_solve(); + void post_solve(); }; } // namespace coek diff --git a/lib/coek/coek/solvers/highs/highs.cpp b/lib/coek/coek/solvers/highs/highs.cpp index 14b11ff5..dfbdb530 100644 --- a/lib/coek/coek/solvers/highs/highs.cpp +++ b/lib/coek/coek/solvers/highs/highs.cpp @@ -12,6 +12,7 @@ #include "coek/model/model.hpp" #include "coek/model/model_repn.hpp" #include "coek/util/io_utils.hpp" +#include "coek/util/sequence.hpp" #include "coek_highs.hpp" @@ -52,26 +53,50 @@ void add_objective(HighsModel& model, Expression& expr, bool sense, orepn.collect_terms(expr); model.lp_.col_cost_.resize(x.size()); - if (orepn.linear_coefs.size() + orepn.quadratic_coefs.size() > 0) { + model.lp_.offset_ = orepn.constval; + + if (orepn.linear_coefs.size()) { 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); + if (orepn.quadratic_coefs.size() > 0) { + std::map, double> value; + bool nz = false; + + for (size_t i : indices(orepn.quadratic_coefs)) { + auto& lvar = orepn.quadratic_lvars[i]; + auto& rvar = orepn.quadratic_rvars[i]; + nz = true; + if (x[rvar->index] >= x[lvar->index]) + value[{x[rvar->index], x[lvar->index]}] += orepn.quadratic_coefs[i]; + else + value[{x[lvar->index], x[rvar->index]}] += orepn.quadratic_coefs[i]; + } + + if (nz) { + HighsHessian& hessian = model.hessian_; + size_t prev = 0; + for (auto& it : value) { + auto [i, j] = it.first; + if (j != prev) { + hessian.start_.push_back(hessian.index_.size()); + prev = j; + } + hessian.index_.push_back(i); + if (i == j) + hessian.value_.push_back(2 * it.second); + else + hessian.value_.push_back(it.second); + } + // hessian.dim_ = hessian.index_.size(); + hessian.dim_ = model.lp_.col_lower_.size(); + hessian.start_.push_back(hessian.index_.size()); + hessian.format_ = HessianFormat::kTriangular; } -#endif } if (sense) @@ -131,13 +156,31 @@ void add_constraint(HighsModel& model, Constraint& con, std::unordered_map HighsSolver::solve(Model& coek_model) +void HighsSolver::set_solver_options() { - auto results = std::make_shared(); + 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::pre_solve() +{ + results = std::make_shared(); results->solver_name = "highs"; results->termination_condition = TerminationCondition::error; results->tic(); +} + +void HighsSolver::post_solve() { results->toc(); } +std::shared_ptr HighsSolver::solve(Model& coek_model) +{ + pre_solve(); auto _coek_model = coek_model.repn.get(); hmodel.clear(); @@ -164,12 +207,14 @@ std::shared_ptr HighsSolver::solve(Model& coek_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while creating objectives " + std::string(e.what()); + post_solve(); 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!"; + results->toc(); return results; } @@ -184,6 +229,7 @@ std::shared_ptr HighsSolver::solve(Model& coek_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught exception while creating constraints " + std::string(e.what()); + results->toc(); return results; } @@ -204,12 +250,17 @@ std::shared_ptr HighsSolver::solve(Model& coek_model) 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; + std::cout << "Hessian Dim: " << hmodel.hessian_.dim_ << std::endl; + std::cout << "Hessian Start: " << hmodel.hessian_.start_ << std::endl; + std::cout << "Hessian Row: " << hmodel.hessian_.index_ << std::endl; + std::cout << "Hessian Value: " << hmodel.hessian_.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"; + results->toc(); return results; } @@ -219,12 +270,13 @@ std::shared_ptr HighsSolver::solve(Model& coek_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + results->toc(); return results; } collect_results(coek_model, results); - results->toc(); + post_solve(); return results; } @@ -282,12 +334,14 @@ std::shared_ptr HighsSolver::solve(CompactModel& compact_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while creating objectives " + std::string(e.what()); + results->toc(); 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!"; + results->toc(); return results; } @@ -311,16 +365,18 @@ std::shared_ptr HighsSolver::solve(CompactModel& compact_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught exception while creating constraints " + std::string(e.what()); + results->toc(); return results; } - hmodel.lp_.num_col_ = hmodel.lp_.col_cost_.size(); - hmodel.lp_.num_row_ = hmodel.lp_.row_lower_.size(); + hmodel.lp_.num_col_ = static_cast(hmodel.lp_.col_cost_.size()); + hmodel.lp_.num_row_ = static_cast(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"; + results->toc(); return results; } @@ -330,12 +386,13 @@ std::shared_ptr HighsSolver::solve(CompactModel& compact_model) catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + results->toc(); return results; } collect_results(model, results); - results->toc(); + post_solve(); return results; } #endif @@ -402,12 +459,14 @@ std::shared_ptr HighsSolver::resolve() results->error_message = "Highs Error: Caught highs exception while creating objectives " + std::string(e.what()); + results->toc(); 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!"; + results->toc(); return results; } @@ -422,6 +481,7 @@ std::shared_ptr HighsSolver::resolve() results->error_message = "Highs Error: Caught highs exception while creating constraints " + std::string(e.what()); + results->toc(); return results; } @@ -432,6 +492,7 @@ std::shared_ptr HighsSolver::resolve() catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + results->toc(); return results; } } @@ -466,6 +527,7 @@ std::shared_ptr HighsSolver::resolve() results->termination_condition = TerminationCondition::invalid_model_for_solver; results->error_message = "Error initializing Highs: Cannot optimize models with nonlinear terms."; + results->toc(); return results; break; }; @@ -477,6 +539,7 @@ std::shared_ptr HighsSolver::resolve() catch (const std::exception& e) { results->error_message = "Highs Error: Caught highs exception while optimizing " + std::string(e.what()); + results->toc(); return results; } } @@ -485,22 +548,10 @@ std::shared_ptr HighsSolver::resolve() collect_results(model, results); #endif - results->toc(); + post_solve(); 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 { diff --git a/lib/coek/test/solver/TestModels.cpp b/lib/coek/test/solver/TestModels.cpp index 544a864e..0cfe762b 100644 --- a/lib/coek/test/solver/TestModels.cpp +++ b/lib/coek/test/solver/TestModels.cpp @@ -18,14 +18,6 @@ bool TestModel::check_results(coek::Model& model, std::shared_ptrobjective_value - optimal_objective) > 1e-7) { - std::cout << fmt::format("TEST ERROR: Unexpected objective values ({} != {})", - results->objective_value, optimal_objective) - << std::endl; - std::cout << to_string(*results, 4) << std::endl; - return false; - } - std::vector primal; for (auto& v : model.get_variables()) primal.push_back(v.value()); @@ -49,6 +41,14 @@ bool TestModel::check_results(coek::Model& model, std::shared_ptrobjective_value - optimal_objective) > 1e-7) { + std::cout << fmt::format("TEST ERROR: Unexpected objective values ({} != {})", + results->objective_value, optimal_objective) + << std::endl; + std::cout << to_string(*results, 4) << std::endl; + return false; + } + return true; } @@ -92,15 +92,53 @@ class SimpleLP1 : public TestModel { } }; -class LP_bounds : public TestModel { +class SimpleQP1 : public TestModel { + public: + SimpleQP1() + { + primal_solution = {375, 250}; + optimal_objective = 9531250; + + // Model + model.name("simpleqp1"); + auto& m = model; + auto x = m.add(coek::variable("x").bounds(0, m.inf)); + auto y = m.add(coek::variable("y").bounds(0, m.inf)); + + m.add_objective(50 * x * x + 40 * y * y); + m.add_constraint(2 * x + 3 * y >= 1500); + m.add_constraint(2 * x + y >= 1000); + } +}; + +class SimpleQP2 : public TestModel { + public: + SimpleQP2() + { + primal_solution = {375, 250}; + optimal_objective = 2875000; + + // Model + model.name("simpleqp2"); + auto& m = model; + auto x = m.add(coek::variable("x").bounds(0, m.inf)); + auto y = m.add(coek::variable("y").bounds(0, m.inf)); + + m.add_objective(10 * x * x + 5 * x * y + 4 * y * x + 10 * y * y); + m.add_constraint(2 * x + 3 * y >= 1500); + m.add_constraint(2 * x + y >= 1000); + } +}; + +class LP_bounds1 : public TestModel { public: - LP_bounds() + LP_bounds1() { primal_solution = {1, -1, 1, -1, 1, 1, -1, 1, -1, 1}; optimal_objective = 10.0; // Model - model.name("lp_bounds"); + model.name("lp_bounds1"); auto a = model.add_variable().value(2.0); auto b = model.add_variable().value(-2.0); auto c = model.add_variable().value(2.0); @@ -128,6 +166,45 @@ class LP_bounds : public TestModel { } }; +class LP_bounds2 : public TestModel { + public: + LP_bounds2() + { + primal_solution = {0, -2, 0, -2, 0, 0, -2, 0, -2, 0, 1}; + optimal_objective = 8.0; + + // Model + model.name("lp_bounds2"); + auto a = model.add_variable().value(2.0); + auto b = model.add_variable().value(-2.0); + auto c = model.add_variable().value(2.0); + auto d = model.add_variable().value(-2.0); + auto e = model.add_variable().value(2.0); + + auto aa = model.add_variable().value(2.0); + auto bb = model.add_variable().value(-2.0); + auto cc = model.add_variable().value(2.0); + auto dd = model.add_variable().value(-2.0); + auto ee = model.add_variable().value(2.0); + + auto x = model.add_variable().value(1.0); + x.fix(); + + model.add_objective(a - b + c - d + e + aa - bb + cc - dd + ee); + + model.add(a + x >= 1); + model.add(b + x <= -1); + model.add(aa + x - 1 >= 0); + model.add(bb + x + 1 <= 0); + model.add(inequality(1, c + x, 3)); + model.add(inequality(-3, d + x, -1)); + model.add(inequality(0, cc + x - 1, 2)); + model.add(inequality(-2, dd + x + 1, 0)); + model.add(e + x == 1); + model.add(ee + x - 1 == 0); + } +}; + class QP_bounds : public TestModel { public: QP_bounds() @@ -171,10 +248,16 @@ std::shared_ptr model(const std::string& name) return std::make_shared(); if (name == "simplelp1") return std::make_shared(); + if (name == "simpleqp1") + return std::make_shared(); + if (name == "simpleqp2") + return std::make_shared(); if (name == "qp_bounds") return std::make_shared(); - if (name == "lp_bounds") - return std::make_shared(); + if (name == "lp_bounds1") + return std::make_shared(); + if (name == "lp_bounds2") + return std::make_shared(); return std::make_shared(); } diff --git a/lib/coek/test/solver/test_examples.cpp b/lib/coek/test/solver/test_examples.cpp index 624a80b3..44c1b549 100644 --- a/lib/coek/test/solver/test_examples.cpp +++ b/lib/coek/test/solver/test_examples.cpp @@ -89,7 +89,7 @@ std::initializer_list adnames{ }; } // namespace -TEST_CASE("ipopt_examples", "[solver][ipopt]") +TEST_CASE("ipopt_examples", "[solvers][ipopt]") { INFO("TEST_CASE ipopt_examples"); @@ -302,7 +302,7 @@ TEST_CASE("ipopt_examples", "[solver][ipopt]") #endif } -TEST_CASE("gurobi_examples", "[solver][gurobi]") +TEST_CASE("gurobi_examples", "[solvers][gurobi]") { coek::Solver solver("gurobi"); solver.set_option("OutputFlag", 0); @@ -329,7 +329,7 @@ TEST_CASE("gurobi_examples", "[solver][gurobi]") } } -TEST_CASE("highs_examples", "[solver][highs]") +TEST_CASE("highs_examples", "[solvers][highs]") { coek::Solver solver("highs"); solver.set_option("output_flag", false); diff --git a/lib/coek/test/solver/test_gurobi.cpp b/lib/coek/test/solver/test_gurobi.cpp index e57e5ad1..490f3544 100644 --- a/lib/coek/test/solver/test_gurobi.cpp +++ b/lib/coek/test/solver/test_gurobi.cpp @@ -11,7 +11,7 @@ // Gurobi Checks // -TEST_CASE("gurobi_checks", "[solver][gurobi]") +TEST_CASE("gurobi_checks", "[solvers][gurobi]") { coek::Solver solver("gurobi"); solver.set_option("OutputFlag", 0); @@ -25,7 +25,6 @@ TEST_CASE("gurobi_checks", "[solver][gurobi]") auto res = solver.solve(m); REQUIRE(res->termination_condition == coek::TerminationCondition::empty_model); } - SECTION("simplelp1") { auto test = test::model("simplelp1"); @@ -34,11 +33,35 @@ TEST_CASE("gurobi_checks", "[solver][gurobi]") auto res = solver.solve(m); REQUIRE(test->check_results(m, res) == true); } - SECTION("lp_bounds") + SECTION("simpleqp1") + { + auto test = test::model("simpleqp1"); + auto m = test->model; + REQUIRE(m.name() == "simpleqp1"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("simpleqp2") + { + auto test = test::model("simpleqp2"); + auto m = test->model; + REQUIRE(m.name() == "simpleqp2"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("lp_bounds1") + { + auto test = test::model("lp_bounds1"); + auto m = test->model; + REQUIRE(m.name() == "lp_bounds1"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("lp_bounds2") { - auto test = test::model("lp_bounds"); + auto test = test::model("lp_bounds2"); auto m = test->model; - REQUIRE(m.name() == "lp_bounds"); + REQUIRE(m.name() == "lp_bounds2"); auto res = solver.solve(m); REQUIRE(test->check_results(m, res) == true); } diff --git a/lib/coek/test/solver/test_highs.cpp b/lib/coek/test/solver/test_highs.cpp index 8fce4d7f..d2541f2c 100644 --- a/lib/coek/test/solver/test_highs.cpp +++ b/lib/coek/test/solver/test_highs.cpp @@ -9,7 +9,7 @@ // Highs Checks // -TEST_CASE("highs_checks", "[solver][highs]") +TEST_CASE("highs_checks", "[solvers][highs]") { coek::Solver solver("highs"); solver.set_option("output_flag", false); @@ -32,11 +32,35 @@ TEST_CASE("highs_checks", "[solver][highs]") auto res = solver.solve(m); REQUIRE(test->check_results(m, res) == true); } - SECTION("lp_bounds") + SECTION("simpleqp1") { - auto test = test::model("lp_bounds"); + auto test = test::model("simpleqp1"); auto m = test->model; - REQUIRE(m.name() == "lp_bounds"); + REQUIRE(m.name() == "simpleqp1"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("simpleqp2") + { + auto test = test::model("simpleqp2"); + auto m = test->model; + REQUIRE(m.name() == "simpleqp2"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("lp_bounds1") + { + auto test = test::model("lp_bounds1"); + auto m = test->model; + REQUIRE(m.name() == "lp_bounds1"); + auto res = solver.solve(m); + REQUIRE(test->check_results(m, res) == true); + } + SECTION("lp_bounds2") + { + auto test = test::model("lp_bounds2"); + auto m = test->model; + REQUIRE(m.name() == "lp_bounds2"); auto res = solver.solve(m); REQUIRE(test->check_results(m, res) == true); } @@ -46,8 +70,7 @@ TEST_CASE("highs_checks", "[solver][highs]") auto m = test->model; REQUIRE(m.name() == "qp_bounds"); auto res = solver.solve(m); - REQUIRE(test->check_results(m, res) - == false); // Expected failure (Our HIGHS interface doesn't yet support QIPs) + REQUIRE(test->check_results(m, res) == true); } } else {