diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index e9b6104e..6192f184 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/uno/solvers/BQPD/BQPDSolver.cpp @@ -80,7 +80,7 @@ namespace uno { if (warmstart_information.objective_changed || warmstart_information.constraints_changed || warmstart_information.problem_changed) { this->save_hessian(subproblem); } - this->solve_subproblem(initial_point, direction, warmstart_information, subproblem); + this->solve_subproblem(subproblem, initial_point, direction, warmstart_information); } void BQPDSolver::solve_LP(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, @@ -89,7 +89,7 @@ namespace uno { DEBUG << "LP:\n"; } this->set_up_subproblem(warmstart_information, subproblem); - this->solve_subproblem(initial_point, direction, warmstart_information, subproblem); + this->solve_subproblem(subproblem, initial_point, direction, warmstart_information); } void BQPDSolver::set_up_subproblem(const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem) { @@ -151,8 +151,8 @@ namespace uno { WSC.ll += sizeof(intptr_t); } - void BQPDSolver::solve_subproblem(const Vector& initial_point, Direction& direction, - const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem) { + void BQPDSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, + const WarmstartInformation& warmstart_information) { direction.primals = initial_point; const int n = static_cast(subproblem.number_variables); const int m = static_cast(subproblem.number_constraints); @@ -293,11 +293,7 @@ namespace uno { } } // namespace -void hessian_vector_product(int *n, const double x[], const double /*ws*/[], const int lws[], double v[]) { - for (int i = 0; i < *n; i++) { - v[i] = 0.; - } - +void hessian_vector_product([[maybe_unused]] int *n, const double x[], const double /*ws*/[], const int lws[], double v[]) { /* int footer_start = lws[0]; for (int i = 0; i < *n; i++) { @@ -317,13 +313,19 @@ void hessian_vector_product(int *n, const double x[], const double /*ws*/[], con std::copy(reinterpret_cast(lws), reinterpret_cast(lws) + sizeof(intptr_t), reinterpret_cast(&pointer_to_subproblem)); const uno::LagrangeNewtonSubproblem* subproblem = reinterpret_cast(pointer_to_subproblem); - uno::Vector my_x(*n); - for (size_t variable_index: uno::Range(*n)) { - my_x[variable_index] = x[variable_index]; + assert(static_cast(*n) == subproblem->number_variables && "BQPD and the subproblem do not have the same number of variables"); + for (size_t i = 0; i < subproblem->number_variables; i++) { + v[i] = 0.; + } + // convert x[] and v[] into Vector + // TODO improve that + uno::Vector primals(subproblem->number_variables); + for (size_t variable_index: uno::Range(subproblem->number_variables)) { + primals[variable_index] = x[variable_index]; } - uno::Vector result(*n); - subproblem->compute_hessian_vector_product(my_x, result); - for (size_t variable_index: uno::Range(*n)) { + uno::Vector result(subproblem->number_variables); + subproblem->compute_hessian_vector_product(primals, result); + for (size_t variable_index: uno::Range(subproblem->number_variables)) { v[variable_index] = result[variable_index]; } } \ No newline at end of file diff --git a/uno/solvers/BQPD/BQPDSolver.hpp b/uno/solvers/BQPD/BQPDSolver.hpp index 2d4d8644..a4db6336 100644 --- a/uno/solvers/BQPD/BQPDSolver.hpp +++ b/uno/solvers/BQPD/BQPDSolver.hpp @@ -84,8 +84,8 @@ namespace uno { RectangularMatrix constraint_jacobian; /*!< Sparse Jacobian of the constraints */ void set_up_subproblem(const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem); - void solve_subproblem(const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information, - const LagrangeNewtonSubproblem& subproblem); + void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, + const WarmstartInformation& warmstart_information); void categorize_constraints(size_t number_variables, Direction& direction); void save_hessian(const LagrangeNewtonSubproblem& subproblem); void save_gradients_to_local_format(size_t number_constraints, const SparseVector& linear_objective, diff --git a/uno/solvers/HiGHS/HiGHSSolver.cpp b/uno/solvers/HiGHS/HiGHSSolver.cpp index fe92ffe7..dd20ba41 100644 --- a/uno/solvers/HiGHS/HiGHSSolver.cpp +++ b/uno/solvers/HiGHS/HiGHSSolver.cpp @@ -1,15 +1,18 @@ #include #include "HiGHSSolver.hpp" +#include "ingredients/local_models/LagrangeNewtonSubproblem.hpp" #include "linear_algebra/RectangularMatrix.hpp" #include "linear_algebra/SparseVector.hpp" #include "linear_algebra/Vector.hpp" #include "optimization/Direction.hpp" +#include "optimization/WarmstartInformation.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")) { + const Options& options): LPSolver(), print_subproblem(options.get_bool("print_subproblem")), + linear_objective(number_variables), constraints(number_constraints), constraint_jacobian(number_constraints, number_variables) { this->model.lp_.sense_ = ObjSense::kMinimize; this->model.lp_.offset_ = 0.; // the linear part of the objective is a dense vector @@ -29,66 +32,78 @@ namespace uno { 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); + void HiGHSSolver::solve_LP(const LagrangeNewtonSubproblem& subproblem, const Vector& /*initial_point*/, Direction& direction, + const WarmstartInformation& warmstart_information) { + this->set_up_subproblem(subproblem, warmstart_information); + this->solve_subproblem(subproblem, direction); + } + + // build the LP in the HiGHS format + void HiGHSSolver::set_up_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information) { + this->model.lp_.num_col_ = static_cast(subproblem.number_variables); + this->model.lp_.num_row_ = static_cast(subproblem.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.; + if (warmstart_information.variable_bounds_changed) { + subproblem.set_direction_bounds(this->model.lp_.col_lower_, this->model.lp_.col_upper_); } // linear part of the objective - for (const auto [variable_index, value]: linear_objective) { - this->model.lp_.col_cost_[variable_index] = value; + if (warmstart_information.objective_changed) { + for (size_t variable_index: Range(subproblem.number_variables)) { + this->model.lp_.col_cost_[variable_index] = 0.; + } + subproblem.evaluate_objective_gradient(this->linear_objective); + for (const auto [variable_index, value]: this->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]; + if (warmstart_information.constraints_changed) { + subproblem.evaluate_constraints(this->constraints); + } + if (warmstart_information.constraint_bounds_changed) { + subproblem.set_constraint_bounds(this->constraints, this->model.lp_.row_lower_, this->model.lp_.row_upper_); } // 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++; - } + if (warmstart_information.constraints_changed) { + // TODO evaluate directly into this->model.lp_.a_matrix_ + subproblem.evaluate_constraint_jacobian(this->constraint_jacobian); + 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 < subproblem.number_constraints; constraint_index++) { + for (const auto [variable_index, value]: this->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 << "Linear objective part: "; print_vector(DEBUG, view(this->model.lp_.col_cost_, 0, subproblem.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++) { + for (size_t variable_index = 0; variable_index < subproblem.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++) { + for (size_t constraint_index = 0; constraint_index < subproblem.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) { + void HiGHSSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction) { // solve the LP HighsStatus return_status = this->highs_solver.passModel(this->model); assert(return_status == HighsStatus::kOk); @@ -116,7 +131,7 @@ namespace uno { 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++) { + for (size_t variable_index: Range(subproblem.number_variables)) { direction.primals[variable_index] = solution.col_value[variable_index]; const double bound_multiplier = solution.col_dual[variable_index]; if (0. < bound_multiplier) { @@ -127,23 +142,10 @@ namespace uno { } } // read the dual solution - for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) { + for (size_t constraint_index: Range(subproblem.number_constraints)) { 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*/, const OptimizationProblem& problem) { - // 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 index 4986528a..68f4a31a 100644 --- a/uno/solvers/HiGHS/HiGHSSolver.hpp +++ b/uno/solvers/HiGHS/HiGHSSolver.hpp @@ -3,6 +3,9 @@ #include "solvers/LPSolver.hpp" #include "Highs.h" +#include "linear_algebra/RectangularMatrix.hpp" +#include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/Vector.hpp" namespace uno { // forward declaration @@ -13,22 +16,20 @@ namespace uno { 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, const OptimizationProblem& problem) override; + void solve_LP(const LagrangeNewtonSubproblem& subproblem, 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, const OptimizationProblem& problem); - void solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints); + SparseVector linear_objective; + Vector constraints; + RectangularMatrix constraint_jacobian; + + void set_up_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information); + void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction); }; } // namespace