diff --git a/uno/ingredients/subproblems/SubproblemFactory.cpp b/uno/ingredients/subproblems/SubproblemFactory.cpp index 9820c522..61ea10cd 100644 --- a/uno/ingredients/subproblems/SubproblemFactory.cpp +++ b/uno/ingredients/subproblems/SubproblemFactory.cpp @@ -6,7 +6,7 @@ #include "SubproblemFactory.hpp" #include "ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp" #include "ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp" -#include "ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp" +#include "ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp" #include "solvers/QPSolverFactory.hpp" #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp" #include "options/Options.hpp" @@ -26,7 +26,7 @@ namespace uno { } // interior-point method else if (subproblem_strategy == "primal_dual_interior_point") { - return std::make_unique(number_variables, number_constraints, number_jacobian_nonzeros, + return std::make_unique(number_variables, number_constraints, number_jacobian_nonzeros, number_hessian_nonzeros, options); } throw std::invalid_argument("Subproblem strategy " + subproblem_strategy + " is not supported"); diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp similarity index 75% rename from uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp rename to uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp index 06a07436..eb5148d3 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp @@ -2,10 +2,10 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include -#include "PrimalDualInteriorPointSubproblem.hpp" +#include "PrimalDualInteriorPointMethod.hpp" +#include "PrimalDualInteriorPointProblem.hpp" #include "optimization/Direction.hpp" #include "optimization/Iterate.hpp" -#include "ingredients/hessian_models/HessianModelFactory.hpp" #include "linear_algebra/SparseStorageFactory.hpp" #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp" @@ -16,7 +16,7 @@ #include "tools/Infinity.hpp" namespace uno { - PrimalDualInteriorPointSubproblem::PrimalDualInteriorPointSubproblem(size_t number_variables, size_t number_constraints, + PrimalDualInteriorPointMethod::PrimalDualInteriorPointMethod(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options): Subproblem("exact", number_variables, number_hessian_nonzeros, false, options), objective_gradient(2 * number_variables), // original variables + barrier terms @@ -50,12 +50,12 @@ namespace uno { l1_constraint_violation_coefficient(options.get_double("l1_constraint_violation_coefficient")) { } - void PrimalDualInteriorPointSubproblem::initialize_statistics(Statistics& statistics, const Options& options) { + void PrimalDualInteriorPointMethod::initialize_statistics(Statistics& statistics, const Options& options) { statistics.add_column("regulariz", Statistics::double_width - 4, options.get_int("statistics_regularization_column_order")); statistics.add_column("barrier", Statistics::double_width - 5, options.get_int("statistics_barrier_parameter_column_order")); } - void PrimalDualInteriorPointSubproblem::generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) { + void PrimalDualInteriorPointMethod::generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) { if (problem.has_inequality_constraints()) { throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel"); } @@ -72,7 +72,7 @@ namespace uno { initial_iterate.set_number_variables(problem.number_variables); // make the initial point strictly feasible wrt the bounds for (size_t variable_index: Range(problem.number_variables)) { - initial_iterate.primals[variable_index] = PrimalDualInteriorPointSubproblem::push_variable_to_interior(initial_iterate.primals[variable_index], + initial_iterate.primals[variable_index] = PrimalDualInteriorPointMethod::push_variable_to_interior(initial_iterate.primals[variable_index], problem.variable_lower_bound(variable_index), problem.variable_upper_bound(variable_index)); } @@ -83,7 +83,7 @@ namespace uno { // set the slacks to the constraint values for (const auto [constraint_index, slack_index]: problem.model.get_slacks()) { - initial_iterate.primals[slack_index] = PrimalDualInteriorPointSubproblem::push_variable_to_interior(initial_iterate.evaluations.constraints[constraint_index], + initial_iterate.primals[slack_index] = PrimalDualInteriorPointMethod::push_variable_to_interior(initial_iterate.evaluations.constraints[constraint_index], problem.variable_lower_bound(slack_index), problem.variable_upper_bound(slack_index)); } // since the slacks have been set, the function evaluations should also be updated @@ -106,11 +106,11 @@ namespace uno { } } - double PrimalDualInteriorPointSubproblem::barrier_parameter() const { + double PrimalDualInteriorPointMethod::barrier_parameter() const { return this->barrier_parameter_update_strategy.get_barrier_parameter(); } - double PrimalDualInteriorPointSubproblem::push_variable_to_interior(double variable_value, double lower_bound, double upper_bound) const { + double PrimalDualInteriorPointMethod::push_variable_to_interior(double variable_value, double lower_bound, double upper_bound) const { const double range = upper_bound - lower_bound; const double perturbation_lb = std::min(this->parameters.push_variable_to_interior_k1 * std::max(1., std::abs(lower_bound)), this->parameters.push_variable_to_interior_k2 * range); @@ -121,62 +121,24 @@ namespace uno { return variable_value; } - void PrimalDualInteriorPointSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, - const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { - // barrier Lagrangian Hessian - if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { - // original Lagrangian Hessian - this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints); - - // diagonal barrier terms (grouped by variable) - for (size_t variable_index: Range(problem.number_variables)) { - double diagonal_barrier_term = 0.; - if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded - const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index); - diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound; - } - if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded - const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index); - diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound; - } - this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index); - } - } - + void PrimalDualInteriorPointMethod::evaluate_functions(Statistics& statistics, const PrimalDualInteriorPointProblem& barrier_problem, + Iterate& current_iterate, const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { // barrier objective gradient if (warmstart_information.objective_changed) { - // original objective gradient - problem.evaluate_objective_gradient(current_iterate, this->objective_gradient); - - // barrier terms - for (size_t variable_index: Range(problem.number_variables)) { - double barrier_term = 0.; - if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded - barrier_term += -this->barrier_parameter()/(current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index)); - // damping - if (not is_finite(problem.variable_upper_bound(variable_index))) { - barrier_term += this->damping_factor * this->barrier_parameter(); - } - } - if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded - barrier_term += -this->barrier_parameter()/(current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index)); - // damping - if (not is_finite(problem.variable_lower_bound(variable_index))) { - barrier_term -= this->damping_factor * this->barrier_parameter(); - } - } - this->objective_gradient.insert(variable_index, barrier_term); - } + barrier_problem.evaluate_objective_gradient(current_iterate, this->objective_gradient); } - // constraints and Jacobian if (warmstart_information.constraints_changed) { - problem.evaluate_constraints(current_iterate, this->constraints); - problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian); + barrier_problem.evaluate_constraints(current_iterate, this->constraints); + barrier_problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian); + } + // barrier Lagrangian Hessian + if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { + this->hessian_model->evaluate(statistics, barrier_problem, current_iterate.primals, current_multipliers.constraints); } } - void PrimalDualInteriorPointSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, + void PrimalDualInteriorPointMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { if (problem.has_inequality_constraints()) { throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel"); @@ -195,11 +157,14 @@ namespace uno { } statistics.set("barrier", this->barrier_parameter()); + // create a primal-dual interior-point reformulation of the problem + PrimalDualInteriorPointProblem barrier_problem(problem, current_multipliers, this->barrier_parameter(), this->damping_factor); + // evaluate the functions at the current iterate - this->evaluate_functions(statistics, problem, current_iterate, current_multipliers, warmstart_information); + this->evaluate_functions(statistics, barrier_problem, current_iterate, current_multipliers, warmstart_information); // compute the primal-dual solution - this->assemble_augmented_system(statistics, problem, current_multipliers); + this->assemble_augmented_system(statistics, current_multipliers, barrier_problem.number_variables, barrier_problem.number_constraints); this->augmented_system.solve(*this->linear_solver); assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality"); this->number_subproblems_solved++; @@ -208,24 +173,24 @@ namespace uno { direction.subproblem_objective = this->evaluate_subproblem_objective(direction); } - void PrimalDualInteriorPointSubproblem::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, - const Multipliers& current_multipliers) { + void PrimalDualInteriorPointMethod::assemble_augmented_system(Statistics& statistics, const Multipliers& current_multipliers, + size_t number_variables, size_t number_constraints) { // assemble, factorize and regularize the augmented matrix - this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints); - this->augmented_system.factorize_matrix(problem.model, *this->linear_solver); + this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, number_variables, number_constraints); + this->augmented_system.factorize_matrix(*this->linear_solver); const double dual_regularization_parameter = std::pow(this->barrier_parameter(), this->parameters.regularization_exponent); - this->augmented_system.regularize_matrix(statistics, problem.model, *this->linear_solver, problem.number_variables, problem.number_constraints, + this->augmented_system.regularize_matrix(statistics, *this->linear_solver, number_variables, number_constraints, dual_regularization_parameter); // check the inertia [[maybe_unused]] auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = this->linear_solver->get_inertia(); - assert(number_pos_eigenvalues == problem.number_variables && number_neg_eigenvalues == problem.number_constraints && number_zero_eigenvalues == 0); + assert(number_pos_eigenvalues == number_variables && number_neg_eigenvalues == number_constraints && number_zero_eigenvalues == 0); // rhs - this->assemble_augmented_rhs(problem, current_multipliers); + this->assemble_augmented_rhs(current_multipliers, number_variables, number_constraints); } - void PrimalDualInteriorPointSubproblem::initialize_feasibility_problem(const l1RelaxedProblem& /*problem*/, Iterate& current_iterate) { + void PrimalDualInteriorPointMethod::initialize_feasibility_problem(const l1RelaxedProblem& /*problem*/, Iterate& current_iterate) { this->solving_feasibility_problem = true; this->first_feasibility_iteration = true; this->subproblem_definition_changed = true; @@ -248,7 +213,7 @@ namespace uno { } // set the elastic variables of the current iterate - void PrimalDualInteriorPointSubproblem::set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) { + void PrimalDualInteriorPointMethod::set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) { DEBUG << "IPM: setting the elastic variables and their duals\n"; // c(x) - p + n = 0 // analytical expression for p and n: @@ -273,18 +238,18 @@ namespace uno { problem.set_elastic_variable_values(current_iterate, elastic_setting_function); } - double PrimalDualInteriorPointSubproblem::proximal_coefficient(const Iterate& /*current_iterate*/) const { + double PrimalDualInteriorPointMethod::proximal_coefficient(const Iterate& /*current_iterate*/) const { return std::sqrt(this->barrier_parameter()); } - void PrimalDualInteriorPointSubproblem::exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) { + void PrimalDualInteriorPointMethod::exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) { assert(this->solving_feasibility_problem && "The barrier subproblem did not know it was solving the feasibility problem."); this->barrier_parameter_update_strategy.set_barrier_parameter(this->previous_barrier_parameter); this->solving_feasibility_problem = false; this->compute_least_square_multipliers(problem, trial_iterate, trial_iterate.multipliers.constraints); } - void PrimalDualInteriorPointSubproblem::set_auxiliary_measure(const Model& model, Iterate& iterate) { + void PrimalDualInteriorPointMethod::set_auxiliary_measure(const Model& model, Iterate& iterate) { // auxiliary measure: barrier terms double barrier_terms = 0.; for (const size_t variable_index: model.get_lower_bounded_variables()) { @@ -305,14 +270,14 @@ namespace uno { iterate.progress.auxiliary = barrier_terms; } - double PrimalDualInteriorPointSubproblem::compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate, + double PrimalDualInteriorPointMethod::compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate, const Vector& primal_direction, double step_length) const { const double directional_derivative = this->compute_barrier_term_directional_derivative(model, current_iterate, primal_direction); return step_length * (-directional_derivative); // }, "α*(μ*X^{-1} e^T d)"}; } - double PrimalDualInteriorPointSubproblem::compute_barrier_term_directional_derivative(const Model& model, const Iterate& current_iterate, + double PrimalDualInteriorPointMethod::compute_barrier_term_directional_derivative(const Model& model, const Iterate& current_iterate, const Vector& primal_direction) const { double directional_derivative = 0.; for (const size_t variable_index: model.get_lower_bounded_variables()) { @@ -333,7 +298,7 @@ namespace uno { return directional_derivative; } - void PrimalDualInteriorPointSubproblem::update_barrier_parameter(const OptimizationProblem& problem, const Iterate& current_iterate, + void PrimalDualInteriorPointMethod::update_barrier_parameter(const OptimizationProblem& problem, const Iterate& current_iterate, const Multipliers& current_multipliers, const DualResiduals& residuals) { const bool barrier_parameter_updated = this->barrier_parameter_update_strategy.update_barrier_parameter(problem, current_iterate, current_multipliers, residuals); @@ -342,7 +307,7 @@ namespace uno { } // Section 3.9 in IPOPT paper - bool PrimalDualInteriorPointSubproblem::is_small_step(const OptimizationProblem& problem, const Vector& current_primals, + bool PrimalDualInteriorPointMethod::is_small_step(const OptimizationProblem& problem, const Vector& current_primals, const Vector& direction_primals) const { const Range variables_range = Range(problem.number_variables); const VectorExpression relative_direction_size{variables_range, [&](size_t variable_index) { @@ -352,14 +317,14 @@ namespace uno { return (norm_inf(relative_direction_size) <= this->parameters.small_direction_factor * machine_epsilon); } - double PrimalDualInteriorPointSubproblem::evaluate_subproblem_objective(const Direction& direction) const { + double PrimalDualInteriorPointMethod::evaluate_subproblem_objective(const Direction& direction) const { const double linear_term = dot(direction.primals, this->objective_gradient); const double quadratic_term = this->hessian_model->hessian.quadratic_product(direction.primals, direction.primals) / 2.; return linear_term + quadratic_term; } // TODO use a single function for primal and dual fraction-to-boundary rules - double PrimalDualInteriorPointSubproblem::primal_fraction_to_boundary(const OptimizationProblem& problem, const Vector& current_primals, + double PrimalDualInteriorPointMethod::primal_fraction_to_boundary(const OptimizationProblem& problem, const Vector& current_primals, const Vector& primal_direction, double tau) { double step_length = 1.; for (const size_t variable_index: problem.get_lower_bounded_variables()) { @@ -382,7 +347,7 @@ namespace uno { return step_length; } - double PrimalDualInteriorPointSubproblem::dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers, + double PrimalDualInteriorPointMethod::dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers, Multipliers& direction_multipliers, double tau) { double step_length = 1.; for (const size_t variable_index: problem.get_lower_bounded_variables()) { @@ -406,7 +371,8 @@ namespace uno { } // generate the right-hand side - void PrimalDualInteriorPointSubproblem::assemble_augmented_rhs(const OptimizationProblem& problem, const Multipliers& current_multipliers) { + void PrimalDualInteriorPointMethod::assemble_augmented_rhs(const Multipliers& current_multipliers, size_t number_variables, + size_t number_constraints) { this->augmented_system.rhs.fill(0.); // objective gradient @@ -415,7 +381,7 @@ namespace uno { } // constraint: evaluations and gradients - for (size_t constraint_index: Range(problem.number_constraints)) { + for (size_t constraint_index: Range(number_constraints)) { // Lagrangian if (current_multipliers.constraints[constraint_index] != 0.) { for (const auto [variable_index, derivative]: this->constraint_jacobian[constraint_index]) { @@ -423,12 +389,12 @@ namespace uno { } } // constraints - this->augmented_system.rhs[problem.number_variables + constraint_index] = -this->constraints[constraint_index]; + this->augmented_system.rhs[number_variables + constraint_index] = -this->constraints[constraint_index]; } - DEBUG2 << "RHS: "; print_vector(DEBUG2, view(this->augmented_system.rhs, 0, problem.number_variables + problem.number_constraints)); DEBUG << '\n'; + DEBUG2 << "RHS: "; print_vector(DEBUG2, view(this->augmented_system.rhs, 0, number_variables + number_constraints)); DEBUG << '\n'; } - void PrimalDualInteriorPointSubproblem::assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, + void PrimalDualInteriorPointMethod::assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, Vector& direction_primals, Multipliers& direction_multipliers) { // form the primal-dual direction direction_primals = view(this->augmented_system.solution, 0, problem.number_variables); @@ -438,15 +404,15 @@ namespace uno { this->compute_bound_dual_direction(problem, current_primals, current_multipliers, direction_primals, direction_multipliers); // determine if the direction is a "small direction" (Section 3.9 of the Ipopt paper) TODO - const bool is_small_step = PrimalDualInteriorPointSubproblem::is_small_step(problem, current_primals, direction_primals); + const bool is_small_step = PrimalDualInteriorPointMethod::is_small_step(problem, current_primals, direction_primals); if (is_small_step) { DEBUG << "This is a small step\n"; } // "fraction-to-boundary" rule for primal variables and constraints multipliers const double tau = std::max(this->parameters.tau_min, 1. - this->barrier_parameter()); - const double primal_step_length = PrimalDualInteriorPointSubproblem::primal_fraction_to_boundary(problem, current_primals, direction_primals, tau); - const double bound_dual_step_length = PrimalDualInteriorPointSubproblem::dual_fraction_to_boundary(problem, current_multipliers, direction_multipliers, tau); + const double primal_step_length = PrimalDualInteriorPointMethod::primal_fraction_to_boundary(problem, current_primals, direction_primals, tau); + const double bound_dual_step_length = PrimalDualInteriorPointMethod::dual_fraction_to_boundary(problem, current_multipliers, direction_multipliers, tau); DEBUG << "Fraction-to-boundary rules:\n"; DEBUG << "primal step length = " << primal_step_length << '\n'; DEBUG << "bound dual step length = " << bound_dual_step_length << "\n\n"; @@ -457,7 +423,7 @@ namespace uno { direction_multipliers.upper_bounds.scale(bound_dual_step_length); } - void PrimalDualInteriorPointSubproblem::compute_bound_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, + void PrimalDualInteriorPointMethod::compute_bound_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, const Vector& primal_direction, Multipliers& direction_multipliers) { direction_multipliers.lower_bounds.fill(0.); direction_multipliers.upper_bounds.fill(0.); @@ -475,7 +441,7 @@ namespace uno { } } - void PrimalDualInteriorPointSubproblem::compute_least_square_multipliers(const OptimizationProblem& problem, Iterate& iterate, + void PrimalDualInteriorPointMethod::compute_least_square_multipliers(const OptimizationProblem& problem, Iterate& iterate, Vector& constraint_multipliers) { this->augmented_system.matrix.set_dimension(problem.number_variables + problem.number_constraints); this->augmented_system.matrix.reset(); @@ -483,7 +449,7 @@ namespace uno { iterate, constraint_multipliers, this->least_square_multiplier_max_norm); } - void PrimalDualInteriorPointSubproblem::postprocess_iterate(const OptimizationProblem& problem, Iterate& iterate) { + void PrimalDualInteriorPointMethod::postprocess_iterate(const OptimizationProblem& problem, Iterate& iterate) { // rescale the bound multipliers (Eq. 16 in Ipopt paper) for (const size_t variable_index: problem.get_lower_bounded_variables()) { const double coefficient = this->barrier_parameter() / (iterate.primals[variable_index] - problem.variable_lower_bound(variable_index)); @@ -522,7 +488,7 @@ namespace uno { } } - void PrimalDualInteriorPointSubproblem::set_initial_point(const Vector& /*point*/) { + void PrimalDualInteriorPointMethod::set_initial_point(const Vector& /*point*/) { // do nothing } } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp similarity index 85% rename from uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp rename to uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp index bef8aefc..292bad34 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.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_INFEASIBLEINTERIORPOINTSUBPROBLEM_H -#define UNO_INFEASIBLEINTERIORPOINTSUBPROBLEM_H +#ifndef UNO_INFEASIBLEINTERIORPOINTMETHOD_H +#define UNO_INFEASIBLEINTERIORPOINTMETHOD_H #include "ingredients/subproblems/Subproblem.hpp" #include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp" @@ -13,6 +13,7 @@ namespace uno { template class DirectSymmetricIndefiniteLinearSolver; class DualResiduals; + class PrimalDualInteriorPointProblem; struct InteriorPointParameters { double tau_min; @@ -23,9 +24,9 @@ namespace uno { double push_variable_to_interior_k2; }; - class PrimalDualInteriorPointSubproblem : public Subproblem { + class PrimalDualInteriorPointMethod : public Subproblem { public: - PrimalDualInteriorPointSubproblem(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, + PrimalDualInteriorPointMethod(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options); void initialize_statistics(Statistics& statistics, const Options& options) override; @@ -67,8 +68,8 @@ namespace uno { [[nodiscard]] double barrier_parameter() const; [[nodiscard]] double push_variable_to_interior(double variable_value, double lower_bound, double upper_bound) const; - void evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - const WarmstartInformation& warmstart_information); + void evaluate_functions(Statistics& statistics, const PrimalDualInteriorPointProblem& barrier_problem, Iterate& current_iterate, + const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information); void update_barrier_parameter(const OptimizationProblem& problem, const Iterate& current_iterate, const Multipliers& current_multipliers, const DualResiduals& residuals); [[nodiscard]] bool is_small_step(const OptimizationProblem& problem, const Vector& current_primals, const Vector& direction_primals) const; @@ -79,8 +80,8 @@ namespace uno { const Vector& primal_direction, double tau); [[nodiscard]] static double dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers, Multipliers& direction_multipliers, double tau); - void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers); - void assemble_augmented_rhs(const OptimizationProblem& problem, const Multipliers& current_multipliers); + void assemble_augmented_system(Statistics& statistics, const Multipliers& current_multipliers, size_t number_variables, size_t number_constraints); + void assemble_augmented_rhs(const Multipliers& current_multipliers, size_t number_variables, size_t number_constraints); void assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, Vector& direction_primals, Multipliers& direction_multipliers); void compute_bound_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, @@ -89,4 +90,4 @@ namespace uno { }; } // namespace -#endif // UNO_INFEASIBLEINTERIORPOINTSUBPROBLEM_H \ No newline at end of file +#endif // UNO_INFEASIBLEINTERIORPOINTMETHOD_H \ No newline at end of file diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.cpp new file mode 100644 index 00000000..50d12fad --- /dev/null +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.cpp @@ -0,0 +1,159 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "PrimalDualInteriorPointProblem.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "optimization/Iterate.hpp" +#include "tools/Infinity.hpp" + +namespace uno { + PrimalDualInteriorPointProblem::PrimalDualInteriorPointProblem(const OptimizationProblem& problem, const Multipliers& current_multipliers, + double barrier_parameter, double damping_factor) : + OptimizationProblem(problem.model, problem.number_variables, problem.number_constraints), + problem(problem), current_multipliers(current_multipliers), barrier_parameter(barrier_parameter), damping_factor(damping_factor) { } + + double PrimalDualInteriorPointProblem::get_objective_multiplier() const { + return this->problem.get_objective_multiplier(); + } + + void PrimalDualInteriorPointProblem::evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const { + this->problem.evaluate_objective_gradient(iterate, objective_gradient); + + // barrier terms + for (size_t variable_index: Range(this->problem.number_variables)) { + double barrier_term = 0.; + if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded + barrier_term += -this->barrier_parameter / (iterate.primals[variable_index] - problem.variable_lower_bound(variable_index)); + // damping + if (not is_finite(problem.variable_upper_bound(variable_index))) { + barrier_term += this->damping_factor * this->barrier_parameter; + } + } + if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded + barrier_term += -this->barrier_parameter / (iterate.primals[variable_index] - problem.variable_upper_bound(variable_index)); + // damping + if (not is_finite(problem.variable_lower_bound(variable_index))) { + barrier_term -= this->damping_factor * this->barrier_parameter; + } + } + objective_gradient.insert(variable_index, barrier_term); + } + } + + void PrimalDualInteriorPointProblem::evaluate_constraints(Iterate& iterate, std::vector& constraints) const { + problem.evaluate_constraints(iterate, constraints); + } + + void PrimalDualInteriorPointProblem::evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const { + problem.evaluate_constraint_jacobian(iterate, constraint_jacobian); + } + + void PrimalDualInteriorPointProblem::evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, + SymmetricMatrix& hessian) const { + // original Lagrangian Hessian + this->problem.evaluate_lagrangian_hessian(x, multipliers, hessian); + + // barrier terms + for (size_t variable_index: Range(this->problem.number_variables)) { + double diagonal_barrier_term = 0.; + if (is_finite(this->problem.variable_lower_bound(variable_index))) { // lower bounded + const double distance_to_bound = x[variable_index] - problem.variable_lower_bound(variable_index); + diagonal_barrier_term += this->current_multipliers.lower_bounds[variable_index] / distance_to_bound; + } + if (is_finite(this->problem.variable_upper_bound(variable_index))) { // upper bounded + const double distance_to_bound = x[variable_index] - problem.variable_upper_bound(variable_index); + diagonal_barrier_term += this->current_multipliers.upper_bounds[variable_index] / distance_to_bound; + } + hessian.insert(diagonal_barrier_term, variable_index, variable_index); + } + } + + double PrimalDualInteriorPointProblem::variable_lower_bound(size_t /*variable_index*/) const { + return -INF; + } + + double PrimalDualInteriorPointProblem::variable_upper_bound(size_t /*variable_index*/) const { + return INF; + } + + double PrimalDualInteriorPointProblem::constraint_lower_bound(size_t /*constraint_index*/) const { + return 0.; + } + + double PrimalDualInteriorPointProblem::constraint_upper_bound(size_t /*constraint_index*/) const { + return 0.; + } + + const Collection& PrimalDualInteriorPointProblem::get_lower_bounded_variables() const { + return this->problem.get_lower_bounded_variables(); + } + + const Collection& PrimalDualInteriorPointProblem::get_upper_bounded_variables() const { + return this->problem.get_upper_bounded_variables(); + } + + const Collection& PrimalDualInteriorPointProblem::get_single_lower_bounded_variables() const { + return this->problem.get_single_lower_bounded_variables(); + } + + const Collection& PrimalDualInteriorPointProblem::get_single_upper_bounded_variables() const { + return this->problem.get_single_upper_bounded_variables(); + } + + size_t PrimalDualInteriorPointProblem::number_objective_gradient_nonzeros() const { + size_t number_nonzeros = this->problem.number_objective_gradient_nonzeros(); + // barrier contribution + for (size_t variable_index: Range(this->problem.number_variables)) { + if (is_finite(this->problem.variable_lower_bound(variable_index)) || is_finite(this->problem.variable_upper_bound(variable_index))) { + number_nonzeros++; + } + } + return number_nonzeros; + } + + size_t PrimalDualInteriorPointProblem::number_jacobian_nonzeros() const { + return this->problem.number_jacobian_nonzeros(); + } + + size_t PrimalDualInteriorPointProblem::number_hessian_nonzeros() const { + size_t number_zeros = this->problem.number_hessian_nonzeros(); + // barrier contribution + for (size_t variable_index: Range(this->problem.number_variables)) { + if (is_finite(this->problem.variable_lower_bound(variable_index)) || is_finite(this->problem.variable_upper_bound(variable_index))) { + number_zeros++; + } + } + return number_zeros; + } + + void PrimalDualInteriorPointProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, + const Multipliers& multipliers) const { + this->problem.evaluate_lagrangian_gradient(lagrangian_gradient, iterate, multipliers); + + // barrier terms + for (size_t variable_index: Range(this->problem.number_variables)) { + double barrier_term = 0.; + if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded + barrier_term += -this->barrier_parameter / (iterate.primals[variable_index] - problem.variable_lower_bound(variable_index)); + // damping + if (not is_finite(problem.variable_upper_bound(variable_index))) { + barrier_term += this->damping_factor * this->barrier_parameter; + } + } + if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded + barrier_term += -this->barrier_parameter / (iterate.primals[variable_index] - problem.variable_upper_bound(variable_index)); + // damping + if (not is_finite(problem.variable_lower_bound(variable_index))) { + barrier_term -= this->damping_factor * this->barrier_parameter; + } + } + // the objective contribution of the Lagrangian gradient may be scaled. Barrier terms go into the constraint contribution + lagrangian_gradient.constraints_contribution[variable_index] += barrier_term; + } + } + + double PrimalDualInteriorPointProblem::complementarity_error(const Vector& primals, const std::vector& constraints, + const Multipliers& multipliers, double shift_value, Norm residual_norm) const { + return this->problem.complementarity_error(primals, constraints, multipliers, shift_value, residual_norm); + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.hpp new file mode 100644 index 00000000..4bc26895 --- /dev/null +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointProblem.hpp @@ -0,0 +1,49 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_PRIMALDUALINTERIORPOINTPROBLEM_H +#define UNO_PRIMALDUALINTERIORPOINTPROBLEM_H + +#include "reformulation/OptimizationProblem.hpp" + +namespace uno { + class PrimalDualInteriorPointProblem : public OptimizationProblem { + public: + PrimalDualInteriorPointProblem(const OptimizationProblem& problem, const Multipliers& current_multipliers, double barrier_parameter, + double damping_factor); + + // function evaluations + [[nodiscard]] double get_objective_multiplier() const override; + void evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const override; + void evaluate_constraints(Iterate& iterate, std::vector& constraints) const override; + void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const override; + void evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, + SymmetricMatrix& hessian) const override; + + [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; + [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; + [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override; + [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override; + [[nodiscard]] const Collection& get_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_upper_bounded_variables() const override; + [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; + + [[nodiscard]] size_t number_objective_gradient_nonzeros() const override; + [[nodiscard]] size_t number_jacobian_nonzeros() const override; + [[nodiscard]] size_t number_hessian_nonzeros() const override; + + void evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, + const Multipliers& multipliers) const override; + [[nodiscard]] double complementarity_error(const Vector& primals, const std::vector& constraints, + const Multipliers& multipliers, double shift_value, Norm residual_norm) const override; + + protected: + const OptimizationProblem& problem; + const Multipliers& current_multipliers; + const double barrier_parameter; + const double damping_factor; + }; +} // namespace + +#endif // UNO_PRIMALDUALINTERIORPOINTPROBLEM_H diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 69176875..d1b33a52 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -9,7 +9,6 @@ #include "SparseStorageFactory.hpp" #include "RectangularMatrix.hpp" #include "ingredients/hessian_models/UnstableRegularization.hpp" -#include "model/Model.hpp" #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "options/Options.hpp" #include "tools/Statistics.hpp" @@ -26,8 +25,8 @@ namespace uno { const Options& options); void assemble_matrix(const SymmetricMatrix& hessian, const RectangularMatrix& constraint_jacobian, size_t number_variables, size_t number_constraints); - void factorize_matrix(const Model& model, DirectSymmetricIndefiniteLinearSolver& linear_solver); - void regularize_matrix(Statistics& statistics, const Model& model, DirectSymmetricIndefiniteLinearSolver& linear_solver, + void factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver); + void regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter); void solve(DirectSymmetricIndefiniteLinearSolver& linear_solver); // [[nodiscard]] T get_primal_regularization() const; @@ -89,8 +88,7 @@ namespace uno { } template - void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& /*model*/, - DirectSymmetricIndefiniteLinearSolver& linear_solver) { + void SymmetricIndefiniteLinearSystem::factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver) { if (true || this->number_factorizations == 0) { linear_solver.do_symbolic_analysis(this->matrix); } @@ -99,7 +97,7 @@ namespace uno { } template - void SymmetricIndefiniteLinearSystem::regularize_matrix(Statistics& statistics, const Model& model, + void SymmetricIndefiniteLinearSystem::regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter) { DEBUG2 << "Original matrix\n" << this->matrix << '\n'; @@ -141,7 +139,7 @@ namespace uno { while (not good_inertia) { DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n"; DEBUG2 << this->matrix << '\n'; - this->factorize_matrix(model, linear_solver); + this->factorize_matrix(linear_solver); number_attempts++; if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) {