From ce7b9571a10d00081afd3523cee8f995aeefc754 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 11 Jun 2025 16:06:24 +0200 Subject: [PATCH 01/11] First draft for LPEQP method based on our 2023 work. TODO: write EQP problem --- CMakeLists.txt | 1 + .../InequalityHandlingMethod.cpp | 12 ++ .../InequalityHandlingMethod.hpp | 3 + .../InequalityHandlingMethodFactory.cpp | 5 + .../EQPProblem.cpp | 10 ++ .../EQPProblem.hpp | 23 +++ .../LPEQPMethod.cpp | 150 ++++++++++++++++++ .../LPEQPMethod.hpp | 64 ++++++++ .../InequalityConstrainedMethod.cpp | 9 -- .../InequalityConstrainedMethod.hpp | 2 - .../subproblem_solvers/QPSolverFactory.hpp | 2 +- ...SymmetricIndefiniteLinearSolverFactory.hpp | 2 +- uno/optimization/Direction.hpp | 22 ++- uno/symbolic/VectorView.hpp | 2 + uno/tools/Infinity.hpp | 4 +- 15 files changed, 295 insertions(+), 16 deletions(-) create mode 100644 uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp create mode 100644 uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp create mode 100644 uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp create mode 100644 uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp create mode 100644 uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index d04d7764a..70be14cfd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -69,6 +69,7 @@ file(GLOB UNO_SOURCE_FILES uno/ingredients/globalization_strategies/switching_methods/funnel_methods/*.cpp uno/ingredients/hessian_models/*.cpp uno/ingredients/inequality_handling_methods/*.cpp + uno/ingredients/inequality_handling_methods/equality_constrained_methods/*.cpp uno/ingredients/inequality_handling_methods/inequality_constrained_methods/*.cpp uno/ingredients/inequality_handling_methods/interior_point_methods/*.cpp uno/ingredients/regularization_strategies/*.cpp diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp new file mode 100644 index 000000000..c742d66e6 --- /dev/null +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp @@ -0,0 +1,12 @@ +#include "InequalityHandlingMethod.hpp" +#include "optimization/Multipliers.hpp" +#include "symbolic/VectorView.hpp" + +namespace uno { + void InequalityHandlingMethod::compute_dual_displacements(const Multipliers& current_multipliers, Multipliers& direction_multipliers) { + // compute dual *displacements* (active-set methods usually compute the new duals, not the displacements) + view(direction_multipliers.constraints, 0, current_multipliers.constraints.size()) -= current_multipliers.constraints; + view(direction_multipliers.lower_bounds, 0, current_multipliers.lower_bounds.size()) -= current_multipliers.lower_bounds; + view(direction_multipliers.upper_bounds, 0, current_multipliers.upper_bounds.size()) -= current_multipliers.upper_bounds; + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp index 8e3675f5d..27493ec39 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp @@ -61,6 +61,9 @@ namespace uno { bool subproblem_definition_changed{false}; [[nodiscard]] virtual std::string get_name() const = 0; + + protected: + static void compute_dual_displacements(const Multipliers& current_multipliers, Multipliers& direction_multipliers); }; } // namespace diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp index 70d5fc340..847e31072 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp @@ -5,6 +5,7 @@ #include "InequalityHandlingMethod.hpp" #include "InequalityHandlingMethodFactory.hpp" #include "inequality_constrained_methods/InequalityConstrainedMethod.hpp" +#include "equality_constrained_methods/LPEQPMethod.hpp" #include "interior_point_methods/PrimalDualInteriorPointMethod.hpp" #include "ingredients/subproblem_solvers/LPSolverFactory.hpp" #include "ingredients/subproblem_solvers/QPSolverFactory.hpp" @@ -19,6 +20,10 @@ namespace uno { if (inequality_handling_method == "inequality_constrained") { return std::make_unique(options); } + // equality-constrained methods + else if (inequality_handling_method == "LPEQP") { + return std::make_unique(options); + } // interior-point method else if (inequality_handling_method == "primal_dual_interior_point") { return std::make_unique(options); diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp new file mode 100644 index 000000000..d4bdfd35f --- /dev/null +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp @@ -0,0 +1,10 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "EQPProblem.hpp" + +namespace uno { + EQPProblem::EQPProblem(const OptimizationProblem& problem, const ActiveSet &active_set): + OptimizationProblem(problem.model), first_reformulation(problem), active_set(active_set) { + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp new file mode 100644 index 000000000..c168baf7d --- /dev/null +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp @@ -0,0 +1,23 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_EQPPROBLEM_H +#define UNO_EQPPROBLEM_H + +#include "optimization/OptimizationProblem.hpp" + +namespace uno { + // forward declaration + class ActiveSet; + + class EQPProblem: public OptimizationProblem { + public: + EQPProblem(const OptimizationProblem& problem, const ActiveSet& active_set); + + protected: + const OptimizationProblem& first_reformulation; + const ActiveSet& active_set; + }; +} // namespace + +#endif // UNO_EQPPROBLEM_H \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp new file mode 100644 index 000000000..18f852fef --- /dev/null +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -0,0 +1,150 @@ +// Copyright (c) 2018-2023 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "LPEQPMethod.hpp" +#include "EQPProblem.hpp" +#include "optimization/Iterate.hpp" +#include "ingredients/constraint_relaxation_strategies/l1RelaxedProblem.hpp" +#include "ingredients/hessian_models/HessianModelFactory.hpp" +#include "ingredients/regularization_strategies/NoRegularization.hpp" +#include "ingredients/regularization_strategies/RegularizationStrategy.hpp" +#include "ingredients/subproblem/Subproblem.hpp" +#include "ingredients/subproblem_solvers/LPSolverFactory.hpp" +#include "ingredients/subproblem_solvers/QPSolverFactory.hpp" +#include "ingredients/subproblem_solvers/SubproblemStatus.hpp" +#include "optimization/Direction.hpp" +#include "options/Options.hpp" +#include "tools/Logger.hpp" + +namespace uno { + LPEQPMethod::LPEQPMethod(const Options& options): + InequalityHandlingMethod(), + enforce_linear_constraints_at_initial_iterate(options.get_bool("enforce_linear_constraints")), + LP_solver(LPSolverFactory::create(options)), + QP_solver(QPSolverFactory::create(options)), + activity_tolerance(options.get_double("TR_activity_tolerance")) { + } + + void LPEQPMethod::initialize(const OptimizationProblem& problem, const HessianModel& hessian_model, + RegularizationStrategy& regularization_strategy) { + this->initial_point.resize(problem.number_variables); + regularization_strategy.initialize_memory(problem, hessian_model); + this->LP_direction = Direction(problem.number_variables, problem.number_constraints); + NoRegularization no_regularization{}; + this->LP_solver->initialize_memory(problem, ZeroHessian(), no_regularization); + this->QP_solver->initialize_memory(problem, hessian_model, regularization_strategy); + } + + void LPEQPMethod::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) { + // do nothing + } + + void LPEQPMethod::generate_initial_iterate(const OptimizationProblem& /*problem*/, Iterate& /*initial_iterate*/) { + /* + if (this->enforce_linear_constraints_at_initial_iterate) { + Preprocessing::enforce_linear_constraints(problem.model, initial_iterate.primals, initial_iterate.multipliers, *this->QP_solver); + } + */ + } + + void LPEQPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, + const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, + RegularizationStrategy& regularization_strategy, double trust_region_radius, + WarmstartInformation& warmstart_information) { + // solve LP subproblem (within trust-region to avoid unboundedness) + if (trust_region_radius == INF) { + trust_region_radius = 100.; // TODO option + } + + Subproblem LP_subproblem{problem, current_iterate, current_multipliers, this->LP_hessian_model, + this->LP_regularization_strategy, trust_region_radius}; + this->solve_LP(statistics, LP_subproblem, current_multipliers, warmstart_information); + DEBUG << "d^*(LP) = " << this->LP_direction << '\n'; + + if (this->LP_direction.status == SubproblemStatus::INFEASIBLE) { + DEBUG << "Infeasible LP, EQP direction will not be computed.\n"; + direction = this->LP_direction; + // reset the initial point + this->initial_point.fill(0.); + return; + } + + // TODO compute Cauchy point + + // reformulate the problem by setting active constraints as equations and inactive bounds as +/-INF + const EQPProblem fixed_active_set_problem(problem, this->LP_direction.active_set); + + // compute EQP direction + Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, + regularization_strategy, trust_region_radius}; + this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, warmstart_information); + DEBUG << "d^*(EQP) = " << direction << '\n'; + // reset the initial point + this->initial_point.fill(0.); + + // TODO compute convex combination of the Cauchy and EQP directions + + } + + void LPEQPMethod::initialize_feasibility_problem(const l1RelaxedProblem& /*problem*/, Iterate& /*current_iterate*/) { + // do nothing + } + + void LPEQPMethod::set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) { + problem.set_elastic_variable_values(current_iterate, [&](Iterate& iterate, size_t /*j*/, size_t elastic_index, double /*jacobian_coefficient*/) { + iterate.primals[elastic_index] = 0.; + iterate.feasibility_multipliers.lower_bounds[elastic_index] = 1.; + iterate.feasibility_multipliers.upper_bounds[elastic_index] = 0.; + }); + } + + double LPEQPMethod::proximal_coefficient() const { + return 0.; + } + + void LPEQPMethod::exit_feasibility_problem(const OptimizationProblem& /*problem*/, Iterate& /*trial_iterate*/) { + // do nothing + } + + // progress measures + double LPEQPMethod::hessian_quadratic_product(const Vector& vector) const { + return this->QP_solver->hessian_quadratic_product(vector); + } + + void LPEQPMethod::set_auxiliary_measure(const OptimizationProblem& /*problem*/, Iterate& iterate) { + iterate.progress.auxiliary = 0.; + } + + double LPEQPMethod::compute_predicted_auxiliary_reduction_model(const OptimizationProblem& /*problem*/, + const Iterate& /*iterate*/, const Vector& /*primal_direction*/, double /*step_length*/) const { + return 0.; + } + + void LPEQPMethod::postprocess_iterate(const OptimizationProblem& /*problem*/, Vector& /*primals*/, Multipliers& /*multipliers*/) { + // do nothing + } + + void LPEQPMethod::set_initial_point(const Vector& point) { + // copy the point into the member + this->initial_point = point; + } + + std::string LPEQPMethod::get_name() const { + return "LP-EQP"; + } + + // protected member functions + void LPEQPMethod::solve_LP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, + const WarmstartInformation& warmstart_information) { + this->LP_solver->solve(statistics, subproblem, this->initial_point, this->LP_direction, warmstart_information); + InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); + this->number_subproblems_solved++; + } + + void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, + Direction& direction, const WarmstartInformation& warmstart_information) { + this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); + InequalityHandlingMethod::compute_dual_displacements(current_multipliers, direction.multipliers); + this->number_subproblems_solved++; + } +} \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp new file mode 100644 index 000000000..bdcb3699f --- /dev/null +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp @@ -0,0 +1,64 @@ +// Copyright (c) 2018-2023 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_LPEQPMETHOD_H +#define UNO_LPEQPMETHOD_H + +#include +#include "ingredients/hessian_models/ZeroHessian.hpp" +#include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp" +#include "ingredients/regularization_strategies/NoRegularization.hpp" +#include "ingredients/subproblem_solvers/QPSolver.hpp" +#include "linear_algebra/Vector.hpp" +#include "optimization/Direction.hpp" + +namespace uno { + class LPEQPMethod : public InequalityHandlingMethod { + public: + explicit LPEQPMethod(const Options& options); + + void initialize(const OptimizationProblem& problem, const HessianModel& hessian_model, + RegularizationStrategy& regularization_strategy) override; + void initialize_statistics(Statistics& statistics, const Options& options) override; + void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; + void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, + const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, + RegularizationStrategy& regularization_strategy, double trust_region_radius, + WarmstartInformation& warmstart_information) override; + + void initialize_feasibility_problem(const l1RelaxedProblem& problem, Iterate& current_iterate) override; + void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) override; + [[nodiscard]] double proximal_coefficient() const override; + void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) override; + + // progress measures + [[nodiscard]] double hessian_quadratic_product(const Vector& vector) const override; + void set_auxiliary_measure(const OptimizationProblem& problem, Iterate& iterate) override; + [[nodiscard]] double compute_predicted_auxiliary_reduction_model(const OptimizationProblem& problem, + const Iterate& iterate, const Vector& primal_direction, double step_length) const override; + + void postprocess_iterate(const OptimizationProblem& problem, Vector& primals, Multipliers& multipliers) override; + + void set_initial_point(const Vector& initial_point) override; + + [[nodiscard]] std::string get_name() const override; + + protected: + Vector initial_point{}; + const bool enforce_linear_constraints_at_initial_iterate; + ZeroHessian LP_hessian_model{}; + NoRegularization LP_regularization_strategy{}; + Direction LP_direction{}; + // pointers to allow polymorphism + const std::unique_ptr LP_solver; + const std::unique_ptr QP_solver; + const double activity_tolerance; + + void solve_LP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, + const WarmstartInformation& warmstart_information); + void solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, + Direction& direction, const WarmstartInformation& warmstart_information); + }; +} // namespace + +#endif // UNO_LPEQPMETHOD_H \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp index ec1509967..4695debed 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp @@ -106,15 +106,6 @@ namespace uno { return evaluation_space.compute_hessian_quadratic_product(vector); } - // compute dual *displacements* - // because of the way we form LPs/QPs, we get the new *multipliers* back from the solver. To get the dual displacements/direction, - // we need to subtract the current multipliers - void InequalityConstrainedMethod::compute_dual_displacements(const Multipliers& current_multipliers, Multipliers& direction_multipliers) { - view(direction_multipliers.constraints, 0, current_multipliers.constraints.size()) -= current_multipliers.constraints; - view(direction_multipliers.lower_bounds, 0, current_multipliers.lower_bounds.size()) -= current_multipliers.lower_bounds; - view(direction_multipliers.upper_bounds, 0, current_multipliers.upper_bounds.size()) -= current_multipliers.upper_bounds; - } - // auxiliary measure is 0 in inequality-constrained methods void InequalityConstrainedMethod::set_auxiliary_measure(const OptimizationProblem& /*problem*/, Iterate& iterate) { iterate.progress.auxiliary = 0.; diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp index af286d314..34ba9f58e 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp @@ -54,8 +54,6 @@ namespace uno { std::unique_ptr solver{}; Vector initial_point{}; const Options& options; // copy of the options for delayed allocation of solver - - static void compute_dual_displacements(const Multipliers& current_multipliers, Multipliers& direction_multipliers); }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/QPSolverFactory.hpp b/uno/ingredients/subproblem_solvers/QPSolverFactory.hpp index df0b2ca27..9d516901d 100644 --- a/uno/ingredients/subproblem_solvers/QPSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/QPSolverFactory.hpp @@ -29,4 +29,4 @@ namespace uno { }; } // namespace -#endif // UNO_QPSOLVERFACTORY_H +#endif // UNO_QPSOLVERFACTORY_H \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp index 1b91bea86..2bbc9d1eb 100644 --- a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp @@ -22,4 +22,4 @@ namespace uno { }; } // namespace -#endif // UNO_LINEARSOLVERFACTORY_H +#endif // UNO_LINEARSOLVERFACTORY_H \ No newline at end of file diff --git a/uno/optimization/Direction.hpp b/uno/optimization/Direction.hpp index 3261e274a..5563d2ae4 100644 --- a/uno/optimization/Direction.hpp +++ b/uno/optimization/Direction.hpp @@ -10,6 +10,25 @@ #include "tools/Infinity.hpp" namespace uno { + struct ActiveConstraints { + std::vector at_lower_bound{}; /*!< List of constraint indices at their lower bound */ + std::vector at_upper_bound{}; /*!< List of constraint indices at their upper bound */ + + explicit ActiveConstraints(size_t capacity) { + this->at_lower_bound.reserve(capacity); + this->at_upper_bound.reserve(capacity); + } + ActiveConstraints() = default; + }; + + struct ActiveSet { + ActiveConstraints constraints{}; /*!< List of general constraints */ + ActiveConstraints bounds{}; /*!< List of bound constraints */ + + ActiveSet(size_t number_variables, size_t number_constraints): constraints(number_constraints), bounds(number_variables) { } + ActiveSet() = default; + }; + class Direction { public: Direction(size_t number_variables, size_t number_constraints); @@ -25,6 +44,7 @@ namespace uno { double norm{INF}; /*!< Norm of \f$x\f$ */ double subproblem_objective{INF}; /*!< Objective value */ + ActiveSet active_set{}; void set_dimensions(size_t new_number_variables, size_t new_number_constraints); void reset(); @@ -33,4 +53,4 @@ namespace uno { }; } // namespace -#endif // UNO_DIRECTION_H +#endif // UNO_DIRECTION_H \ No newline at end of file diff --git a/uno/symbolic/VectorView.hpp b/uno/symbolic/VectorView.hpp index bd5ccc578..ff271b923 100644 --- a/uno/symbolic/VectorView.hpp +++ b/uno/symbolic/VectorView.hpp @@ -4,6 +4,8 @@ #ifndef UNO_VECTORVIEW_H #define UNO_VECTORVIEW_H +#include + namespace uno { // span of an arbitrary container: allocation-free view of a certain length template diff --git a/uno/tools/Infinity.hpp b/uno/tools/Infinity.hpp index 50bce2aa4..db787b586 100644 --- a/uno/tools/Infinity.hpp +++ b/uno/tools/Infinity.hpp @@ -9,7 +9,7 @@ namespace uno { template - const double INF = std::numeric_limits::infinity(); + constexpr double INF = std::numeric_limits::infinity(); template inline bool is_finite(NumericalType value) { @@ -17,4 +17,4 @@ namespace uno { } } // namespace -#endif // UNO_INFINITY_H +#endif // UNO_INFINITY_H \ No newline at end of file From 98278f0f2e94773f8292d5cdb0a6c824a32848c2 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 12 Jun 2025 00:53:07 +0200 Subject: [PATCH 02/11] Minor changes --- .../equality_constrained_methods/LPEQPMethod.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index 18f852fef..f12acf3ac 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -30,8 +30,7 @@ namespace uno { this->initial_point.resize(problem.number_variables); regularization_strategy.initialize_memory(problem, hessian_model); this->LP_direction = Direction(problem.number_variables, problem.number_constraints); - NoRegularization no_regularization{}; - this->LP_solver->initialize_memory(problem, ZeroHessian(), no_regularization); + this->LP_solver->initialize_memory(problem, this->LP_hessian_model, this->LP_regularization_strategy); this->QP_solver->initialize_memory(problem, hessian_model, regularization_strategy); } @@ -78,6 +77,7 @@ namespace uno { Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, regularization_strategy, trust_region_radius}; this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, warmstart_information); + DEBUG << "d^*(EQP) = " << direction << '\n'; // reset the initial point this->initial_point.fill(0.); @@ -140,7 +140,7 @@ namespace uno { InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); this->number_subproblems_solved++; } - + void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); From b30d6f21afb1fbac30061b1ea69bb482d6dd96ea Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 12 Jun 2025 19:09:54 +0200 Subject: [PATCH 03/11] Trying to fix an allocation bug in LPEQP --- .../{EQPProblem.cpp => FixedActiveSetProblem.cpp} | 4 ++-- .../{EQPProblem.hpp => FixedActiveSetProblem.hpp} | 10 +++++----- .../equality_constrained_methods/LPEQPMethod.cpp | 5 ++--- 3 files changed, 9 insertions(+), 10 deletions(-) rename uno/ingredients/inequality_handling_methods/equality_constrained_methods/{EQPProblem.cpp => FixedActiveSetProblem.cpp} (64%) rename uno/ingredients/inequality_handling_methods/equality_constrained_methods/{EQPProblem.hpp => FixedActiveSetProblem.hpp} (59%) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp similarity index 64% rename from uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp rename to uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp index d4bdfd35f..6154a69cd 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp @@ -1,10 +1,10 @@ // Copyright (c) 2018-2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. -#include "EQPProblem.hpp" +#include "FixedActiveSetProblem.hpp" namespace uno { - EQPProblem::EQPProblem(const OptimizationProblem& problem, const ActiveSet &active_set): + FixedActiveSetProblem::FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet &active_set): OptimizationProblem(problem.model), first_reformulation(problem), active_set(active_set) { } } // namespace \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp similarity index 59% rename from uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp rename to uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp index c168baf7d..2399e12a5 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/EQPProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp @@ -1,8 +1,8 @@ // Copyright (c) 2018-2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. -#ifndef UNO_EQPPROBLEM_H -#define UNO_EQPPROBLEM_H +#ifndef UNO_FIXEDACTIVESETPROBLEM_H +#define UNO_FIXEDACTIVESETPROBLEM_H #include "optimization/OptimizationProblem.hpp" @@ -10,9 +10,9 @@ namespace uno { // forward declaration class ActiveSet; - class EQPProblem: public OptimizationProblem { + class FixedActiveSetProblem: public OptimizationProblem { public: - EQPProblem(const OptimizationProblem& problem, const ActiveSet& active_set); + FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet& active_set); protected: const OptimizationProblem& first_reformulation; @@ -20,4 +20,4 @@ namespace uno { }; } // namespace -#endif // UNO_EQPPROBLEM_H \ No newline at end of file +#endif // UNO_FIXEDACTIVESETPROBLEM_H \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index f12acf3ac..e24344f9f 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "LPEQPMethod.hpp" -#include "EQPProblem.hpp" +#include "FixedActiveSetProblem.hpp" #include "optimization/Iterate.hpp" #include "ingredients/constraint_relaxation_strategies/l1RelaxedProblem.hpp" #include "ingredients/hessian_models/HessianModelFactory.hpp" @@ -71,13 +71,12 @@ namespace uno { // TODO compute Cauchy point // reformulate the problem by setting active constraints as equations and inactive bounds as +/-INF - const EQPProblem fixed_active_set_problem(problem, this->LP_direction.active_set); + const FixedActiveSetProblem fixed_active_set_problem(problem, this->LP_direction.active_set); // compute EQP direction Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, regularization_strategy, trust_region_radius}; this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, warmstart_information); - DEBUG << "d^*(EQP) = " << direction << '\n'; // reset the initial point this->initial_point.fill(0.); From 09b46e806bd98a77a6eb80d731ee14f90a7fbda3 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 13 Jun 2025 14:09:12 +0200 Subject: [PATCH 04/11] BQPDSolver: forgot to reset the active set of the direction --- uno/optimization/Direction.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/uno/optimization/Direction.hpp b/uno/optimization/Direction.hpp index 5563d2ae4..9836d2d4d 100644 --- a/uno/optimization/Direction.hpp +++ b/uno/optimization/Direction.hpp @@ -27,6 +27,13 @@ namespace uno { ActiveSet(size_t number_variables, size_t number_constraints): constraints(number_constraints), bounds(number_variables) { } ActiveSet() = default; + + void reset() { + this->constraints.at_lower_bound.clear(); + this->constraints.at_upper_bound.clear(); + this->bounds.at_lower_bound.clear(); + this->bounds.at_upper_bound.clear(); + } }; class Direction { From cf0c70d424345725227899aee4cbc6dc524b0add Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 13 Jun 2025 14:14:06 +0200 Subject: [PATCH 05/11] Fixed variables and constraints in FixedActiveSetProblem + corrected the LP and EQP multipliers --- .../FixedActiveSetProblem.cpp | 61 ++++++++++++++++++- .../FixedActiveSetProblem.hpp | 14 ++++- .../LPEQPMethod.cpp | 26 +++++++- 3 files changed, 97 insertions(+), 4 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp index 6154a69cd..471179f74 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp @@ -2,9 +2,66 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "FixedActiveSetProblem.hpp" +#include "optimization/Direction.hpp" +#include "tools/Infinity.hpp" namespace uno { - FixedActiveSetProblem::FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet &active_set): - OptimizationProblem(problem.model), first_reformulation(problem), active_set(active_set) { + FixedActiveSetProblem::FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet &active_set, + double trust_region_radius): + OptimizationProblem(problem.model), first_reformulation(problem), active_set(active_set), + trust_region_radius(trust_region_radius), + variables_lower_bounds(problem.number_variables), variables_upper_bounds(problem.number_variables), + constraints_lower_bounds(problem.number_constraints), constraints_upper_bounds(problem.number_constraints) { + // initialize all bounds to {-INF,+INF} + for (size_t i: Range(problem.number_variables)) { + this->variables_lower_bounds[i] = -INF; + this->variables_upper_bounds[i] = INF; + } + // set active lower bounds as equality constraints + for (size_t variable_index: this->active_set.bounds.at_lower_bound) { + const double lb = problem.variable_lower_bound(variable_index); + if (-this->trust_region_radius + this->activity_tolerance <= lb) { + this->variables_lower_bounds[variable_index] = this->variables_upper_bounds[variable_index] = lb; + } + } + // set active lower bounds as equality constraints + for (size_t variable_index: this->active_set.bounds.at_upper_bound) { + const double ub = problem.variable_upper_bound(variable_index); + if (ub <= this->trust_region_radius - this->activity_tolerance) { + this->variables_lower_bounds[variable_index] = this->variables_upper_bounds[variable_index] = ub; + } + } + + // initialize all constraint bounds to {-INF,+INF} + for (size_t constraint_index: Range(problem.number_constraints)) { + this->constraints_lower_bounds[constraint_index] = -INF; + this->constraints_upper_bounds[constraint_index] = INF; + } + // set active lower constraint bounds as equality constraints + for (size_t constraint_index: this->active_set.constraints.at_lower_bound) { + this->constraints_lower_bounds[constraint_index] = this->constraints_upper_bounds[constraint_index] = + problem.constraint_lower_bound(constraint_index); + } + // set active upper constraint bounds as equality constraints + for (size_t constraint_index: this->active_set.constraints.at_upper_bound) { + this->constraints_lower_bounds[constraint_index] = this->constraints_upper_bounds[constraint_index] = + problem.constraint_upper_bound(constraint_index); + } + } + + double FixedActiveSetProblem::variable_lower_bound(size_t variable_index) const { + return this->variables_lower_bounds[variable_index]; + } + + double FixedActiveSetProblem::variable_upper_bound(size_t variable_index) const { + return this->variables_upper_bounds[variable_index]; + } + + double FixedActiveSetProblem::constraint_lower_bound(size_t constraint_index) const { + return this->constraints_lower_bounds[constraint_index]; + } + + double FixedActiveSetProblem::constraint_upper_bound(size_t constraint_index) const { + return this->constraints_upper_bounds[constraint_index]; } } // namespace \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp index 2399e12a5..ba26dd6b3 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp @@ -12,11 +12,23 @@ namespace uno { class FixedActiveSetProblem: public OptimizationProblem { public: - FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet& active_set); + FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet& active_set, double trust_region_radius); + + [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const; + [[nodiscard]] virtual double variable_upper_bound(size_t variable_index) const; + [[nodiscard]] virtual double constraint_lower_bound(size_t constraint_index) const; + [[nodiscard]] virtual double constraint_upper_bound(size_t constraint_index) const; protected: const OptimizationProblem& first_reformulation; const ActiveSet& active_set; + double trust_region_radius; + + const double activity_tolerance{1e-6}; // TODO option + Vector variables_lower_bounds; + Vector variables_upper_bounds; + Vector constraints_lower_bounds; + Vector constraints_upper_bounds; }; } // namespace diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index e24344f9f..ad1272cd2 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -50,6 +50,7 @@ namespace uno { const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, double trust_region_radius, WarmstartInformation& warmstart_information) { + //std::cout << "Before solving LP:\n"; warmstart_information.display(); // solve LP subproblem (within trust-region to avoid unboundedness) if (trust_region_radius == INF) { trust_region_radius = 100.; // TODO option @@ -67,11 +68,14 @@ namespace uno { this->initial_point.fill(0.); return; } + else { + this->initial_point = this->LP_direction.primals; + } // TODO compute Cauchy point // reformulate the problem by setting active constraints as equations and inactive bounds as +/-INF - const FixedActiveSetProblem fixed_active_set_problem(problem, this->LP_direction.active_set); + const FixedActiveSetProblem fixed_active_set_problem(problem, this->LP_direction.active_set, trust_region_radius); // compute EQP direction Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, @@ -138,11 +142,31 @@ namespace uno { this->LP_solver->solve(statistics, subproblem, this->initial_point, this->LP_direction, warmstart_information); InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); this->number_subproblems_solved++; + + // reset multipliers for bound constraints active at trust region (except if one of the original bounds is active) + for (size_t variable_index: Range(subproblem.number_variables)) { + if (std::abs(this->LP_direction.primals[variable_index] + subproblem.trust_region_radius) <= this->activity_tolerance) { + this->LP_direction.multipliers.lower_bounds[variable_index] = 0.; + } + if (std::abs(this->LP_direction.primals[variable_index] - subproblem.trust_region_radius) <= this->activity_tolerance) { + this->LP_direction.multipliers.upper_bounds[variable_index] = 0.; + } + } } void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); + + // fix EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) + for (size_t variable_index: this->LP_direction.active_set.bounds.at_lower_bound) { + direction.multipliers.lower_bounds[variable_index] = direction.multipliers.upper_bounds[variable_index]; + direction.multipliers.upper_bounds[variable_index] = 0.; + } + for (size_t variable_index: this->LP_direction.active_set.bounds.at_upper_bound) { + direction.multipliers.upper_bounds[variable_index] = direction.multipliers.lower_bounds[variable_index]; + direction.multipliers.lower_bounds[variable_index] = 0.; + } InequalityHandlingMethod::compute_dual_displacements(current_multipliers, direction.multipliers); this->number_subproblems_solved++; } From a43a23e5d9996246e88d1f35e8b125b2699cfe2b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 13 Jun 2025 14:44:18 +0200 Subject: [PATCH 06/11] Clean up --- .../equality_constrained_methods/LPEQPMethod.cpp | 10 ++++++---- .../equality_constrained_methods/LPEQPMethod.hpp | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index ad1272cd2..21e4d014e 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -93,6 +93,10 @@ namespace uno { // do nothing } + void LPEQPMethod::exit_feasibility_problem(const OptimizationProblem& /*problem*/, Iterate& /*trial_iterate*/) { + // do nothing + } + void LPEQPMethod::set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) { problem.set_elastic_variable_values(current_iterate, [&](Iterate& iterate, size_t /*j*/, size_t elastic_index, double /*jacobian_coefficient*/) { iterate.primals[elastic_index] = 0.; @@ -105,10 +109,6 @@ namespace uno { return 0.; } - void LPEQPMethod::exit_feasibility_problem(const OptimizationProblem& /*problem*/, Iterate& /*trial_iterate*/) { - // do nothing - } - // progress measures double LPEQPMethod::hessian_quadratic_product(const Vector& vector) const { return this->QP_solver->hessian_quadratic_product(vector); @@ -143,6 +143,7 @@ namespace uno { InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); this->number_subproblems_solved++; + // TODO compare radius and original bound // reset multipliers for bound constraints active at trust region (except if one of the original bounds is active) for (size_t variable_index: Range(subproblem.number_variables)) { if (std::abs(this->LP_direction.primals[variable_index] + subproblem.trust_region_radius) <= this->activity_tolerance) { @@ -158,6 +159,7 @@ namespace uno { Direction& direction, const WarmstartInformation& warmstart_information) { this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); + // TODO compare radius and original bound // fix EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) for (size_t variable_index: this->LP_direction.active_set.bounds.at_lower_bound) { direction.multipliers.lower_bounds[variable_index] = direction.multipliers.upper_bounds[variable_index]; diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp index bdcb3699f..465534248 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp @@ -27,9 +27,9 @@ namespace uno { WarmstartInformation& warmstart_information) override; void initialize_feasibility_problem(const l1RelaxedProblem& problem, Iterate& current_iterate) override; + void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) override; void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) override; [[nodiscard]] double proximal_coefficient() const override; - void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) override; // progress measures [[nodiscard]] double hessian_quadratic_product(const Vector& vector) const override; From adb0867cd69fb516459bcbf5cfd0e4ae334ff23c Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 13 Jun 2025 16:53:21 +0200 Subject: [PATCH 07/11] Minor changes --- .../equality_constrained_methods/LPEQPMethod.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index 21e4d014e..eddce576f 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -153,6 +153,7 @@ namespace uno { this->LP_direction.multipliers.upper_bounds[variable_index] = 0.; } } + this->number_subproblems_solved++; } void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, @@ -160,7 +161,7 @@ namespace uno { this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); // TODO compare radius and original bound - // fix EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) + // correct EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) for (size_t variable_index: this->LP_direction.active_set.bounds.at_lower_bound) { direction.multipliers.lower_bounds[variable_index] = direction.multipliers.upper_bounds[variable_index]; direction.multipliers.upper_bounds[variable_index] = 0.; From 2d9afb8c8619f93e564ce80493903e47df8deef3 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 20 Jun 2025 21:53:49 +0200 Subject: [PATCH 08/11] Removed useless comments --- .../equality_constrained_methods/LPEQPMethod.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index eddce576f..85ce2db13 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -143,7 +143,6 @@ namespace uno { InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); this->number_subproblems_solved++; - // TODO compare radius and original bound // reset multipliers for bound constraints active at trust region (except if one of the original bounds is active) for (size_t variable_index: Range(subproblem.number_variables)) { if (std::abs(this->LP_direction.primals[variable_index] + subproblem.trust_region_radius) <= this->activity_tolerance) { @@ -159,8 +158,7 @@ namespace uno { void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); - - // TODO compare radius and original bound + // correct EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) for (size_t variable_index: this->LP_direction.active_set.bounds.at_lower_bound) { direction.multipliers.lower_bounds[variable_index] = direction.multipliers.upper_bounds[variable_index]; From b87420423bb98e6315dbe77936fc60ba4cc974fc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 23 Jun 2025 12:12:45 +0200 Subject: [PATCH 09/11] Added some logger info --- .../FeasibilityRestoration.cpp | 4 ++-- .../equality_constrained_methods/LPEQPMethod.cpp | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index c341270b6..b44a22a43 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -78,7 +78,7 @@ namespace uno { if (this->current_phase == Phase::OPTIMALITY) { statistics.set("phase", "OPT"); try { - DEBUG << "Solving the optimality subproblem\n"; + DEBUG << "Solving the optimality problem\n"; const OptimizationProblem optimality_problem{model}; this->solve_subproblem(statistics, *this->optimality_inequality_handling_method, optimality_problem, current_iterate, direction, *this->optimality_hessian_model, *this->optimality_regularization_strategy, trust_region_radius, @@ -103,7 +103,7 @@ namespace uno { } // solve the feasibility problem (minimize the constraint violation) - DEBUG << "Solving the feasibility subproblem\n"; + DEBUG << "Solving the feasibility problem\n"; statistics.set("phase", "FEAS"); // note: failure of regularization should not happen here, since the feasibility Jacobian has full rank l1RelaxedProblem feasibility_problem{model, 0., this->constraint_violation_coefficient, diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index 85ce2db13..29f5b451d 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -139,6 +139,7 @@ namespace uno { // protected member functions void LPEQPMethod::solve_LP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { + DEBUG << "LP-EQP: solving the LP subproblem\n"; this->LP_solver->solve(statistics, subproblem, this->initial_point, this->LP_direction, warmstart_information); InequalityHandlingMethod::compute_dual_displacements(current_multipliers, this->LP_direction.multipliers); this->number_subproblems_solved++; @@ -157,6 +158,7 @@ namespace uno { void LPEQPMethod::solve_EQP(Statistics& statistics, Subproblem& subproblem, const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { + DEBUG << "LP-EQP: solving the EQP subproblem\n"; this->QP_solver->solve(statistics, subproblem, this->initial_point, direction, warmstart_information); // correct EQP multipliers (the QP solver has no knowledge of the original bounds of fixed variables) From acfa6750a09025b70fbbf92714b687ff8ae77370 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 23 Jun 2025 20:40:27 +0200 Subject: [PATCH 10/11] Pass primal direction to FixedActiveSetProblem to fix active set --- .../TrustRegionStrategy.cpp | 4 ++-- .../FixedActiveSetProblem.cpp | 23 ++++++++++--------- .../FixedActiveSetProblem.hpp | 8 ++++--- .../LPEQPMethod.cpp | 13 +++++++---- .../SymmetricIndefiniteLinearSolver.hpp | 2 +- 5 files changed, 28 insertions(+), 22 deletions(-) diff --git a/uno/ingredients/globalization_mechanisms/TrustRegionStrategy.cpp b/uno/ingredients/globalization_mechanisms/TrustRegionStrategy.cpp index f9c660c0d..8aae72e02 100644 --- a/uno/ingredients/globalization_mechanisms/TrustRegionStrategy.cpp +++ b/uno/ingredients/globalization_mechanisms/TrustRegionStrategy.cpp @@ -116,11 +116,11 @@ namespace uno { assert(0 < this->radius && "The trust-region radius should be positive"); // reset multipliers for bound constraints active at trust region (except if one of the original bounds is active) for (size_t variable_index: Range(model.number_variables)) { - if (std::abs(direction.primals[variable_index] + this->radius) <= this->activity_tolerance && + if (std::abs(direction.primals[variable_index] + this->radius) <= this->activity_tolerance && // d_i = -radius this->activity_tolerance < std::abs(trial_iterate.primals[variable_index] - model.variable_lower_bound(variable_index))) { trial_iterate.multipliers.lower_bounds[variable_index] = 0.; } - if (std::abs(direction.primals[variable_index] - this->radius) <= this->activity_tolerance && + if (std::abs(direction.primals[variable_index] - this->radius) <= this->activity_tolerance && // d_i = radius this->activity_tolerance < std::abs(model.variable_upper_bound(variable_index) - trial_iterate.primals[variable_index])) { trial_iterate.multipliers.upper_bounds[variable_index] = 0.; } diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp index 471179f74..f1179ea98 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.cpp @@ -3,46 +3,47 @@ #include "FixedActiveSetProblem.hpp" #include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" #include "tools/Infinity.hpp" +#include "tools/Logger.hpp" namespace uno { - FixedActiveSetProblem::FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet &active_set, - double trust_region_radius): - OptimizationProblem(problem.model), first_reformulation(problem), active_set(active_set), + FixedActiveSetProblem::FixedActiveSetProblem(const OptimizationProblem& problem, const Vector& direction_primals, + const ActiveSet &active_set, double trust_region_radius): + OptimizationProblem(problem.model), problem(problem), direction_primals(direction_primals), + active_set(active_set), trust_region_radius(trust_region_radius), variables_lower_bounds(problem.number_variables), variables_upper_bounds(problem.number_variables), constraints_lower_bounds(problem.number_constraints), constraints_upper_bounds(problem.number_constraints) { - // initialize all bounds to {-INF,+INF} + // variables: initialize all bounds to {-inf, +inf} for (size_t i: Range(problem.number_variables)) { this->variables_lower_bounds[i] = -INF; this->variables_upper_bounds[i] = INF; } - // set active lower bounds as equality constraints + // variables: set active bounds (except trust region constraint) as fixed variables for (size_t variable_index: this->active_set.bounds.at_lower_bound) { const double lb = problem.variable_lower_bound(variable_index); - if (-this->trust_region_radius + this->activity_tolerance <= lb) { + if (-INF < lb && std::abs(this->direction_primals[variable_index] + this->trust_region_radius) >= this->activity_tolerance) { this->variables_lower_bounds[variable_index] = this->variables_upper_bounds[variable_index] = lb; } } - // set active lower bounds as equality constraints for (size_t variable_index: this->active_set.bounds.at_upper_bound) { const double ub = problem.variable_upper_bound(variable_index); - if (ub <= this->trust_region_radius - this->activity_tolerance) { + if (ub < INF && std::abs(this->direction_primals[variable_index] - this->trust_region_radius) >= this->activity_tolerance) { this->variables_lower_bounds[variable_index] = this->variables_upper_bounds[variable_index] = ub; } } - // initialize all constraint bounds to {-INF,+INF} + // constraints: initialize all constraint bounds to {-inf, +inf} for (size_t constraint_index: Range(problem.number_constraints)) { this->constraints_lower_bounds[constraint_index] = -INF; this->constraints_upper_bounds[constraint_index] = INF; } - // set active lower constraint bounds as equality constraints + // constraints: set active bounds as equality constraints for (size_t constraint_index: this->active_set.constraints.at_lower_bound) { this->constraints_lower_bounds[constraint_index] = this->constraints_upper_bounds[constraint_index] = problem.constraint_lower_bound(constraint_index); } - // set active upper constraint bounds as equality constraints for (size_t constraint_index: this->active_set.constraints.at_upper_bound) { this->constraints_lower_bounds[constraint_index] = this->constraints_upper_bounds[constraint_index] = problem.constraint_upper_bound(constraint_index); diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp index ba26dd6b3..01f92d3d8 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/FixedActiveSetProblem.hpp @@ -8,11 +8,12 @@ namespace uno { // forward declaration - class ActiveSet; + struct ActiveSet; class FixedActiveSetProblem: public OptimizationProblem { public: - FixedActiveSetProblem(const OptimizationProblem& problem, const ActiveSet& active_set, double trust_region_radius); + FixedActiveSetProblem(const OptimizationProblem& problem, const Vector& current_primals, const ActiveSet& active_set, + double trust_region_radius); [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const; [[nodiscard]] virtual double variable_upper_bound(size_t variable_index) const; @@ -20,7 +21,8 @@ namespace uno { [[nodiscard]] virtual double constraint_upper_bound(size_t constraint_index) const; protected: - const OptimizationProblem& first_reformulation; + const OptimizationProblem& problem; + const Vector& direction_primals; const ActiveSet& active_set; double trust_region_radius; diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index 29f5b451d..2417aed51 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -13,6 +13,7 @@ #include "ingredients/subproblem_solvers/QPSolverFactory.hpp" #include "ingredients/subproblem_solvers/SubproblemStatus.hpp" #include "optimization/Direction.hpp" +#include "optimization/WarmstartInformation.hpp" #include "options/Options.hpp" #include "tools/Logger.hpp" @@ -49,7 +50,7 @@ namespace uno { void LPEQPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, double trust_region_radius, - WarmstartInformation& warmstart_information) { + WarmstartInformation& /*warmstart_information*/) { //std::cout << "Before solving LP:\n"; warmstart_information.display(); // solve LP subproblem (within trust-region to avoid unboundedness) if (trust_region_radius == INF) { @@ -58,7 +59,8 @@ namespace uno { Subproblem LP_subproblem{problem, current_iterate, current_multipliers, this->LP_hessian_model, this->LP_regularization_strategy, trust_region_radius}; - this->solve_LP(statistics, LP_subproblem, current_multipliers, warmstart_information); + constexpr WarmstartInformation LP_warmstart_information{}; + this->solve_LP(statistics, LP_subproblem, current_multipliers, LP_warmstart_information); DEBUG << "d^*(LP) = " << this->LP_direction << '\n'; if (this->LP_direction.status == SubproblemStatus::INFEASIBLE) { @@ -75,18 +77,19 @@ namespace uno { // TODO compute Cauchy point // reformulate the problem by setting active constraints as equations and inactive bounds as +/-INF - const FixedActiveSetProblem fixed_active_set_problem(problem, this->LP_direction.active_set, trust_region_radius); + const FixedActiveSetProblem fixed_active_set_problem(problem, this->LP_direction.primals, this->LP_direction.active_set, + trust_region_radius); // compute EQP direction Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, regularization_strategy, trust_region_radius}; - this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, warmstart_information); + constexpr WarmstartInformation EQP_warmstart_information{}; + this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, EQP_warmstart_information); DEBUG << "d^*(EQP) = " << direction << '\n'; // reset the initial point this->initial_point.fill(0.); // TODO compute convex combination of the Cauchy and EQP directions - } void LPEQPMethod::initialize_feasibility_problem(const l1RelaxedProblem& /*problem*/, Iterate& /*current_iterate*/) { diff --git a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp index 20390d524..2966248eb 100644 --- a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp +++ b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp @@ -14,7 +14,7 @@ namespace uno { class Subproblem; template class Vector; - class WarmstartInformation; + struct WarmstartInformation; template class SymmetricIndefiniteLinearSolver { From 7fac25cf542b8f7e876427dc1e94a2a6b1dc3257 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 28 Aug 2025 18:15:20 +0200 Subject: [PATCH 11/11] Fixed rebasing --- .../InequalityHandlingMethod.hpp | 1 + .../LPEQPMethod.cpp | 57 ++++++++++++------- .../LPEQPMethod.hpp | 19 ++++--- .../NoRegularization.hpp | 1 + uno/ingredients/subproblem/Subproblem.cpp | 4 +- uno/ingredients/subproblem/Subproblem.hpp | 2 +- 6 files changed, 55 insertions(+), 29 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp index 27493ec39..454eb98e8 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp @@ -13,6 +13,7 @@ namespace uno { class HessianModel; class Iterate; class l1RelaxedProblem; + class Multipliers; class OptimizationProblem; class Options; template diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp index 2417aed51..ad1436fe7 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.cpp @@ -13,6 +13,7 @@ #include "ingredients/subproblem_solvers/QPSolverFactory.hpp" #include "ingredients/subproblem_solvers/SubproblemStatus.hpp" #include "optimization/Direction.hpp" +#include "optimization/EvaluationSpace.hpp" #include "optimization/WarmstartInformation.hpp" #include "options/Options.hpp" #include "tools/Logger.hpp" @@ -26,13 +27,14 @@ namespace uno { activity_tolerance(options.get_double("TR_activity_tolerance")) { } - void LPEQPMethod::initialize(const OptimizationProblem& problem, const HessianModel& hessian_model, - RegularizationStrategy& regularization_strategy) { + void LPEQPMethod::initialize(const OptimizationProblem& problem, Iterate& current_iterate, + HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, double trust_region_radius) { this->initial_point.resize(problem.number_variables); - regularization_strategy.initialize_memory(problem, hessian_model); this->LP_direction = Direction(problem.number_variables, problem.number_constraints); - this->LP_solver->initialize_memory(problem, this->LP_hessian_model, this->LP_regularization_strategy); - this->QP_solver->initialize_memory(problem, hessian_model, regularization_strategy); + const Subproblem LP_subproblem{problem, current_iterate, this->LP_hessian_model, this->LP_regularization_strategy, trust_region_radius}; + this->LP_solver->initialize_memory(LP_subproblem); + const Subproblem EQP_subproblem{problem, current_iterate, hessian_model, regularization_strategy, trust_region_radius}; + this->QP_solver->initialize_memory(EQP_subproblem); } void LPEQPMethod::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) { @@ -48,19 +50,17 @@ namespace uno { } void LPEQPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, - const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, - RegularizationStrategy& regularization_strategy, double trust_region_radius, - WarmstartInformation& /*warmstart_information*/) { + Direction& direction, HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, + double trust_region_radius, WarmstartInformation& /*warmstart_information*/) { //std::cout << "Before solving LP:\n"; warmstart_information.display(); // solve LP subproblem (within trust-region to avoid unboundedness) if (trust_region_radius == INF) { trust_region_radius = 100.; // TODO option } - Subproblem LP_subproblem{problem, current_iterate, current_multipliers, this->LP_hessian_model, - this->LP_regularization_strategy, trust_region_radius}; + Subproblem LP_subproblem{problem, current_iterate, this->LP_hessian_model, this->LP_regularization_strategy, trust_region_radius}; constexpr WarmstartInformation LP_warmstart_information{}; - this->solve_LP(statistics, LP_subproblem, current_multipliers, LP_warmstart_information); + this->solve_LP(statistics, LP_subproblem, current_iterate.multipliers, LP_warmstart_information); DEBUG << "d^*(LP) = " << this->LP_direction << '\n'; if (this->LP_direction.status == SubproblemStatus::INFEASIBLE) { @@ -81,10 +81,10 @@ namespace uno { trust_region_radius); // compute EQP direction - Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, current_multipliers, hessian_model, + Subproblem EQP_subproblem{fixed_active_set_problem, current_iterate, hessian_model, regularization_strategy, trust_region_radius}; constexpr WarmstartInformation EQP_warmstart_information{}; - this->solve_EQP(statistics, EQP_subproblem, current_multipliers, direction, EQP_warmstart_information); + this->solve_EQP(statistics, EQP_subproblem, current_iterate.multipliers, direction, EQP_warmstart_information); DEBUG << "d^*(EQP) = " << direction << '\n'; // reset the initial point this->initial_point.fill(0.); @@ -103,8 +103,8 @@ namespace uno { void LPEQPMethod::set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) { problem.set_elastic_variable_values(current_iterate, [&](Iterate& iterate, size_t /*j*/, size_t elastic_index, double /*jacobian_coefficient*/) { iterate.primals[elastic_index] = 0.; - iterate.feasibility_multipliers.lower_bounds[elastic_index] = 1.; - iterate.feasibility_multipliers.upper_bounds[elastic_index] = 0.; + iterate.multipliers.lower_bounds[elastic_index] = 1.; + iterate.multipliers.upper_bounds[elastic_index] = 0.; }); } @@ -112,11 +112,30 @@ namespace uno { return 0.; } - // progress measures - double LPEQPMethod::hessian_quadratic_product(const Vector& vector) const { - return this->QP_solver->hessian_quadratic_product(vector); + EvaluationSpace& LPEQPMethod::get_evaluation_space() const { + return this->QP_solver->get_evaluation_space(); + } + + void LPEQPMethod::evaluate_constraint_jacobian(const OptimizationProblem& problem, Iterate& iterate) { + auto& evaluation_space = this->QP_solver->get_evaluation_space(); + evaluation_space.evaluate_constraint_jacobian(problem, iterate); + } + + void LPEQPMethod::compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const { + auto& evaluation_space = this->QP_solver->get_evaluation_space(); + evaluation_space.compute_constraint_jacobian_vector_product(vector, result); } + void LPEQPMethod::compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const { + auto& evaluation_space = this->QP_solver->get_evaluation_space(); + evaluation_space.compute_constraint_jacobian_transposed_vector_product(vector, result); + } + + double LPEQPMethod::compute_hessian_quadratic_product(const Vector& vector) const { + auto& evaluation_space = this->QP_solver->get_evaluation_space(); + return evaluation_space.compute_hessian_quadratic_product(vector); + } + void LPEQPMethod::set_auxiliary_measure(const OptimizationProblem& /*problem*/, Iterate& iterate) { iterate.progress.auxiliary = 0.; } @@ -126,7 +145,7 @@ namespace uno { return 0.; } - void LPEQPMethod::postprocess_iterate(const OptimizationProblem& /*problem*/, Vector& /*primals*/, Multipliers& /*multipliers*/) { + void LPEQPMethod::postprocess_iterate(const OptimizationProblem& /*problem*/, Iterate& /*iterate*/) { // do nothing } diff --git a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp index 465534248..acc7b7fd2 100644 --- a/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/equality_constrained_methods/LPEQPMethod.hpp @@ -17,27 +17,32 @@ namespace uno { public: explicit LPEQPMethod(const Options& options); - void initialize(const OptimizationProblem& problem, const HessianModel& hessian_model, - RegularizationStrategy& regularization_strategy) override; + void initialize(const OptimizationProblem& problem, Iterate& current_iterate, + HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, double trust_region_radius) override; void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, - const Multipliers& current_multipliers, Direction& direction, HessianModel& hessian_model, - RegularizationStrategy& regularization_strategy, double trust_region_radius, - WarmstartInformation& warmstart_information) override; + Direction& direction, HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, + double trust_region_radius, WarmstartInformation& warmstart_information) override; void initialize_feasibility_problem(const l1RelaxedProblem& problem, Iterate& current_iterate) override; void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) override; void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) override; [[nodiscard]] double proximal_coefficient() const override; + // matrix computations + [[nodiscard]] EvaluationSpace& get_evaluation_space() const override; + void evaluate_constraint_jacobian(const OptimizationProblem& problem, Iterate& iterate) override; + void compute_constraint_jacobian_vector_product(const Vector& vector, Vector& result) const override; + void compute_constraint_jacobian_transposed_vector_product(const Vector& vector, Vector& result) const override; + [[nodiscard]] double compute_hessian_quadratic_product(const Vector& vector) const override; + // progress measures - [[nodiscard]] double hessian_quadratic_product(const Vector& vector) const override; void set_auxiliary_measure(const OptimizationProblem& problem, Iterate& iterate) override; [[nodiscard]] double compute_predicted_auxiliary_reduction_model(const OptimizationProblem& problem, const Iterate& iterate, const Vector& primal_direction, double step_length) const override; - void postprocess_iterate(const OptimizationProblem& problem, Vector& primals, Multipliers& multipliers) override; + void postprocess_iterate(const OptimizationProblem& problem, Iterate& iterate) override; void set_initial_point(const Vector& initial_point) override; diff --git a/uno/ingredients/regularization_strategies/NoRegularization.hpp b/uno/ingredients/regularization_strategies/NoRegularization.hpp index 3be3da3bc..03714f0f3 100644 --- a/uno/ingredients/regularization_strategies/NoRegularization.hpp +++ b/uno/ingredients/regularization_strategies/NoRegularization.hpp @@ -5,6 +5,7 @@ #define UNO_NOREGULARIZATION_H #include "RegularizationStrategy.hpp" +#include "ingredients/subproblem/Subproblem.hpp" namespace uno { template diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index 5fa760011..080b159d6 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -13,8 +13,8 @@ namespace uno { Subproblem::Subproblem(const OptimizationProblem& problem, Iterate& current_iterate, HessianModel& hessian_model, RegularizationStrategy& regularization_strategy, double trust_region_radius): number_variables(problem.number_variables), number_constraints(problem.number_constraints), - problem(problem), current_iterate(current_iterate), hessian_model(hessian_model), - regularization_strategy(regularization_strategy), trust_region_radius(trust_region_radius) { + problem(problem), current_iterate(current_iterate), trust_region_radius(trust_region_radius), + hessian_model(hessian_model), regularization_strategy(regularization_strategy) { } void Subproblem::compute_constraint_jacobian_sparsity(int* row_indices, int* column_indices, int solver_indexing, diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index d7e82f0bc..703159812 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -81,11 +81,11 @@ namespace uno { const OptimizationProblem& problem; Iterate& current_iterate; + const double trust_region_radius; protected: HessianModel& hessian_model; RegularizationStrategy& regularization_strategy; - const double trust_region_radius; const ForwardRange empty_set{0}; };