From 29685b3191fc93330b0d7a8a82fc87e3c1624dac Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 28 Oct 2024 12:30:44 +0100 Subject: [PATCH] Store the fixed variables in the model. If there are fixed variables, fail in PrimalDualInteriorPointSubproblem for the moment. TO-DO: move the fixed bounds to the set of general constraints --- bindings/AMPL/AMPLModel.cpp | 32 ++++---- bindings/AMPL/AMPLModel.hpp | 5 +- bindings/AMPL/uno_ampl.cpp | 3 +- uno/Uno.cpp | 82 +++++++++---------- uno/Uno.hpp | 2 +- .../PrimalDualInteriorPointSubproblem.cpp | 3 + uno/linear_algebra/Vector.hpp | 4 + uno/model/BoundRelaxedModel.hpp | 1 + .../HomogeneousEqualityConstrainedModel.hpp | 1 + uno/model/Model.hpp | 1 + uno/model/ScaledModel.hpp | 1 + uno/reformulation/OptimizationProblem.hpp | 5 ++ 12 files changed, 80 insertions(+), 60 deletions(-) diff --git a/bindings/AMPL/AMPLModel.cpp b/bindings/AMPL/AMPLModel.cpp index 32fe9b90..b7c4cfa5 100644 --- a/bindings/AMPL/AMPLModel.cpp +++ b/bindings/AMPL/AMPLModel.cpp @@ -13,7 +13,6 @@ #include "tools/Infinity.hpp" #include "options/Options.hpp" #include "symbolic/Concatenation.hpp" -#include "symbolic/ScalarMultiple.hpp" #include "symbolic/UnaryNegation.hpp" #include "Uno.hpp" @@ -27,7 +26,6 @@ namespace uno { int n_discrete = asl->i.nbv_ + asl->i.niv_ + asl->i.nlvbi_ + asl->i.nlvci_ + asl->i.nlvoi_; if (0 < n_discrete) { throw std::runtime_error("Error: " + std::to_string(n_discrete) + " variables are discrete, which Uno cannot handle"); - // asl->i.need_nl_ = 0; } // preallocate initial primal and dual solutions @@ -56,13 +54,14 @@ namespace uno { variable_status(this->number_variables), constraint_type(this->number_constraints), constraint_status(this->number_constraints), - multipliers_with_flipped_sign(this->number_variables, this->number_constraints), + multipliers_with_flipped_sign(this->number_constraints), equality_constraints_collection(this->equality_constraints), inequality_constraints_collection(this->inequality_constraints), lower_bounded_variables_collection(this->lower_bounded_variables), upper_bounded_variables_collection(this->upper_bounded_variables), single_lower_bounded_variables_collection(this->single_lower_bounded_variables), - single_upper_bounded_variables_collection(this->single_upper_bounded_variables) { + single_upper_bounded_variables_collection(this->single_upper_bounded_variables), + fixed_variables_collection(this->fixed_variables) { // evaluate the constraint Jacobian in sparse mode this->asl->i.congrd_mode = 1; @@ -71,6 +70,7 @@ namespace uno { this->upper_bounded_variables.reserve(this->number_variables); this->single_lower_bounded_variables.reserve(this->number_variables); this->single_upper_bounded_variables.reserve(this->number_variables); + this->fixed_variables.reserve(this->number_variables); this->generate_variables(); // constraints @@ -93,6 +93,7 @@ namespace uno { this->variable_upper_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index + 1] : INF; if (this->variable_lower_bounds[variable_index] == this->variable_upper_bounds[variable_index]) { WARNING << "Variable x" << variable_index << " has identical bounds\n"; + this->fixed_variables.emplace_back(variable_index); } } AMPLModel::determine_bounds_types(this->variable_lower_bounds, this->variable_upper_bounds, this->variable_status); @@ -283,16 +284,16 @@ namespace uno { // evaluate the Hessian: store the matrix in a preallocated array this->asl_hessian const int objective_number = -1; // flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g - this->multipliers_with_flipped_sign.constraints = -multipliers; + this->multipliers_with_flipped_sign = -multipliers; if (this->fixed_hessian_sparsity) { (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast(this->asl_hessian.data()), objective_number, &objective_multiplier, - const_cast(this->multipliers_with_flipped_sign.constraints.data())); + const_cast(this->multipliers_with_flipped_sign.data())); } else { double* objective_multiplier_pointer = (objective_multiplier != 0.) ? &objective_multiplier : nullptr; const bool all_zeros_multipliers = are_all_zeros(multipliers); (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast(this->asl_hessian.data()), objective_number, objective_multiplier_pointer, - all_zeros_multipliers ? nullptr : const_cast(this->multipliers_with_flipped_sign.constraints.data())); + all_zeros_multipliers ? nullptr : const_cast(this->multipliers_with_flipped_sign.data())); } // generate the sparsity pattern in the right sparse format @@ -376,24 +377,25 @@ namespace uno { this->asl->p.solve_code_ = 500; } - // flip the signs of the multipliers if we maximize - this->multipliers_with_flipped_sign.constraints = this->objective_sign * iterate.multipliers.constraints; - this->multipliers_with_flipped_sign.lower_bounds = this->objective_sign * iterate.multipliers.lower_bounds; - this->multipliers_with_flipped_sign.upper_bounds = this->objective_sign * iterate.multipliers.upper_bounds; + // flip the signs of the multipliers and the objective if we maximize + iterate.multipliers.constraints *= this->objective_sign; + iterate.multipliers.lower_bounds *= this->objective_sign; + iterate.multipliers.upper_bounds *= this->objective_sign; + iterate.evaluations.objective *= this->objective_sign; // include the bound duals in the .sol file, using suffixes SufDecl lower_bound_suffix{const_cast("lower_bound_duals"), nullptr, ASL_Sufkind_var | ASL_Sufkind_real, 0}; SufDecl upper_bound_suffix{const_cast("upper_bound_duals"), nullptr, ASL_Sufkind_var | ASL_Sufkind_real, 0}; std::array suffixes{lower_bound_suffix, upper_bound_suffix}; - suf_declare_ASL(this->asl, suffixes.data(), 2); - suf_rput_ASL(this->asl, "lower_bound_duals", ASL_Sufkind_var, this->multipliers_with_flipped_sign.lower_bounds.data()); - suf_rput_ASL(this->asl, "upper_bound_duals", ASL_Sufkind_var, this->multipliers_with_flipped_sign.upper_bounds.data()); + suf_declare_ASL(this->asl, suffixes.data(), suffixes.size()); + suf_rput_ASL(this->asl, "lower_bound_duals", ASL_Sufkind_var, iterate.multipliers.lower_bounds.data()); + suf_rput_ASL(this->asl, "upper_bound_duals", ASL_Sufkind_var, iterate.multipliers.upper_bounds.data()); Option_Info option_info{}; option_info.wantsol = 9; // write the solution without printing the message to stdout std::string message = "Uno "; message.append(Uno::current_version()).append(": ").append(status_to_message(termination_status)); - write_sol_ASL(this->asl, message.data(), iterate.primals.data(), this->multipliers_with_flipped_sign.constraints.data(), &option_info); + write_sol_ASL(this->asl, message.data(), iterate.primals.data(), iterate.multipliers.constraints.data(), &option_info); } } diff --git a/bindings/AMPL/AMPLModel.hpp b/bindings/AMPL/AMPLModel.hpp index 37a19222..900d856c 100644 --- a/bindings/AMPL/AMPLModel.hpp +++ b/bindings/AMPL/AMPLModel.hpp @@ -47,6 +47,7 @@ namespace uno { [[nodiscard]] const SparseVector& get_slacks() const override; [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; + [[nodiscard]] const Collection& get_fixed_variables() const override { return this->fixed_variables_collection; } [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override; [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override; @@ -84,7 +85,7 @@ namespace uno { std::vector constraint_status; /*!< Status of the constraints (EQUAL_BOUNDS, BOUNDED_LOWER, BOUNDED_UPPER, BOUNDED_BOTH_SIDES, * UNBOUNDED) */ std::vector linear_constraints; - mutable Multipliers multipliers_with_flipped_sign; + mutable Vector multipliers_with_flipped_sign; // lists of variables and constraints + corresponding collection objects std::vector equality_constraints{}; @@ -100,6 +101,8 @@ namespace uno { CollectionAdapter&> single_lower_bounded_variables_collection; std::vector single_upper_bounded_variables{}; // indices of the single upper-bounded variables CollectionAdapter&> single_upper_bounded_variables_collection; + std::vector fixed_variables; + CollectionAdapter&> fixed_variables_collection; void generate_variables(); void generate_constraints(); diff --git a/bindings/AMPL/uno_ampl.cpp b/bindings/AMPL/uno_ampl.cpp index 6b4d18e2..656705e9 100644 --- a/bindings/AMPL/uno_ampl.cpp +++ b/bindings/AMPL/uno_ampl.cpp @@ -50,8 +50,7 @@ namespace uno { Uno uno = Uno(*globalization_mechanism, options); // solve the instance - Result result = uno.solve(*model, initial_iterate, options); - uno.print_optimization_summary(result); + uno.solve(*model, initial_iterate, options); // std::cout << "memory_allocation_amount = " << memory_allocation_amount << '\n'; } catch (std::exception& exception) { diff --git a/uno/Uno.cpp b/uno/Uno.cpp index ee7fb079..0d502cbb 100644 --- a/uno/Uno.cpp +++ b/uno/Uno.cpp @@ -29,58 +29,58 @@ namespace uno { Level Logger::level = INFO; - Result Uno::solve(const Model& model, Iterate& current_iterate, const Options& options) { + void Uno::solve(const Model& model, Iterate& current_iterate, const Options& options) { Timer timer{}; Statistics statistics = Uno::create_statistics(model, options); - // use the initial primal-dual point to initialize the strategies and generate the initial iterate - this->initialize(statistics, current_iterate, options); - current_iterate.status = TerminationStatus::NOT_OPTIMAL; - // allocate the trial iterate once and for all here - Iterate trial_iterate(current_iterate); - - size_t major_iterations = 0; try { - bool termination = false; - // check for termination - while (not termination) { - major_iterations++; + // use the initial primal-dual point to initialize the strategies and generate the initial iterate + this->initialize(statistics, current_iterate, options); + current_iterate.status = TerminationStatus::NOT_OPTIMAL; + // allocate the trial iterate once and for all here + Iterate trial_iterate(current_iterate); + + size_t major_iterations = 0; + try { + bool termination = false; + // check for termination + while (not termination) { + major_iterations++; + statistics.start_new_line(); + statistics.set("iter", major_iterations); + DEBUG << "### Outer iteration " << major_iterations << '\n'; + + // compute an acceptable iterate by solving a subproblem at the current point + this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate); + termination = this->termination_criteria(trial_iterate.status, major_iterations, timer.get_duration()); + // the trial iterate becomes the current iterate for the next iteration + std::swap(current_iterate, trial_iterate); + } + } + catch (std::exception& exception) { statistics.start_new_line(); - statistics.set("iter", major_iterations); - DEBUG << "### Outer iteration " << major_iterations << '\n'; - - // compute an acceptable iterate by solving a subproblem at the current point - this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate); - termination = this->termination_criteria(trial_iterate.status, major_iterations, timer.get_duration()); - // the trial iterate becomes the current iterate for the next iteration - std::swap(current_iterate, trial_iterate); + statistics.set("status", exception.what()); + if (Logger::level == INFO) statistics.print_current_line(); + DEBUG << exception.what() << '\n'; } + if (Logger::level == INFO) statistics.print_footer(); + + Uno::postprocess_iterate(model, current_iterate, current_iterate.status); + Result result = this->create_result(model, current_iterate, major_iterations, timer); + this->print_optimization_summary(result); } - catch (std::exception& exception) { - statistics.start_new_line(); - statistics.set("status", exception.what()); - if (Logger::level == INFO) statistics.print_current_line(); - DEBUG << exception.what() << '\n'; + catch (const std::exception& e) { + DISCRETE << "An error occurred at the initial iterate: " << e.what() << '\n'; } - if (Logger::level == INFO) statistics.print_footer(); - - Uno::postprocess_iterate(model, current_iterate, current_iterate.status); - return this->create_result(model, current_iterate, major_iterations, timer); } void Uno::initialize(Statistics& statistics, Iterate& current_iterate, const Options& options) { - try { - statistics.start_new_line(); - statistics.set("iter", 0); - statistics.set("status", "initial point"); - this->globalization_mechanism.initialize(statistics, current_iterate, options); - options.print_used(); - if (Logger::level == INFO) statistics.print_current_line(); - } - catch (const std::exception& e) { - DISCRETE << RED << "An error occurred at the initial iterate: " << e.what() << RESET << '\n'; - throw; - } + statistics.start_new_line(); + statistics.set("iter", 0); + statistics.set("status", "initial point"); + this->globalization_mechanism.initialize(statistics, current_iterate, options); + options.print_used(); + if (Logger::level == INFO) statistics.print_current_line(); } Statistics Uno::create_statistics(const Model& model, const Options& options) { diff --git a/uno/Uno.hpp b/uno/Uno.hpp index 03cb2031..577459cf 100644 --- a/uno/Uno.hpp +++ b/uno/Uno.hpp @@ -27,7 +27,7 @@ namespace uno { public: Uno(GlobalizationMechanism& globalization_mechanism, const Options& options); - [[nodiscard]] Result solve(const Model& model, Iterate& initial_iterate, const Options& options); + void solve(const Model& model, Iterate& initial_iterate, const Options& options); static std::string current_version(); static void print_available_strategies(); diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index cf67f14b..9427d78a 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp @@ -58,6 +58,9 @@ namespace uno { if (problem.has_inequality_constraints()) { throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel"); } + if (problem.has_fixed_variables()) { + throw std::runtime_error("The problem has fixed variables. Move them to the set of general constraints."); + } // TODO: enforce linear constraints at initial point //if (options.get_bool("enforce_linear_constraints")) { diff --git a/uno/linear_algebra/Vector.hpp b/uno/linear_algebra/Vector.hpp index 7a5d7f6c..f4008614 100644 --- a/uno/linear_algebra/Vector.hpp +++ b/uno/linear_algebra/Vector.hpp @@ -89,6 +89,10 @@ namespace uno { } } + void operator*=(ElementType factor) { + this->scale(factor); + } + ElementType* data() { return this->vector.data(); } const ElementType* data() const { return this->vector.data(); } diff --git a/uno/model/BoundRelaxedModel.hpp b/uno/model/BoundRelaxedModel.hpp index 14bccea1..105dd849 100644 --- a/uno/model/BoundRelaxedModel.hpp +++ b/uno/model/BoundRelaxedModel.hpp @@ -40,6 +40,7 @@ namespace uno { [[nodiscard]] const SparseVector& get_slacks() const override { return this->model->get_slacks(); } [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); } [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); } + [[nodiscard]] const Collection& get_fixed_variables() const override { return this->model->get_fixed_variables(); } [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override { return this->model->constraint_lower_bound(constraint_index); } [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override { return this->model->constraint_upper_bound(constraint_index); } diff --git a/uno/model/HomogeneousEqualityConstrainedModel.hpp b/uno/model/HomogeneousEqualityConstrainedModel.hpp index d8e122ac..08a29daf 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.hpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.hpp @@ -45,6 +45,7 @@ namespace uno { [[nodiscard]] const Collection& get_equality_constraints() const override { return this->equality_constraints; } [[nodiscard]] const Collection& get_inequality_constraints() const override { return this->inequality_constraints; } [[nodiscard]] const std::vector& get_linear_constraints() const override { return this->model->get_linear_constraints(); } + [[nodiscard]] const Collection& get_fixed_variables() const override { return this->model->get_fixed_variables(); } void initial_primal_point(Vector& x) const override; void initial_dual_point(Vector& multipliers) const override { this->model->initial_dual_point(multipliers); } diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp index ee74a6d9..8e499a4c 100644 --- a/uno/model/Model.hpp +++ b/uno/model/Model.hpp @@ -64,6 +64,7 @@ namespace uno { [[nodiscard]] virtual const SparseVector& get_slacks() const = 0; [[nodiscard]] virtual const Collection& get_single_lower_bounded_variables() const = 0; [[nodiscard]] virtual const Collection& get_single_upper_bounded_variables() const = 0; + [[nodiscard]] virtual const Collection& get_fixed_variables() const = 0; [[nodiscard]] virtual double constraint_lower_bound(size_t constraint_index) const = 0; [[nodiscard]] virtual double constraint_upper_bound(size_t constraint_index) const = 0; diff --git a/uno/model/ScaledModel.hpp b/uno/model/ScaledModel.hpp index c190c551..af9af815 100644 --- a/uno/model/ScaledModel.hpp +++ b/uno/model/ScaledModel.hpp @@ -30,6 +30,7 @@ namespace uno { [[nodiscard]] const SparseVector& get_slacks() const override { return this->model->get_slacks(); } [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); } [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); } + [[nodiscard]] const Collection& get_fixed_variables() const override { return this->model->get_fixed_variables(); } [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override; [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override; diff --git a/uno/reformulation/OptimizationProblem.hpp b/uno/reformulation/OptimizationProblem.hpp index 17454e92..cf20030b 100644 --- a/uno/reformulation/OptimizationProblem.hpp +++ b/uno/reformulation/OptimizationProblem.hpp @@ -34,6 +34,7 @@ namespace uno { [[nodiscard]] bool is_constrained() const; [[nodiscard]] bool has_inequality_constraints() const; + [[nodiscard]] bool has_fixed_variables() const; // function evaluations [[nodiscard]] virtual double get_objective_multiplier() const = 0; @@ -75,6 +76,10 @@ namespace uno { return (not this->model.get_inequality_constraints().is_empty()); } + inline bool OptimizationProblem::has_fixed_variables() const { + return (not this->model.get_fixed_variables().is_empty()); + } + inline size_t OptimizationProblem::get_number_original_variables() const { return this->model.number_variables; }