diff --git a/uno/ingredients/convexification_strategies/PrimalDualConvexificationStrategy.hpp b/uno/ingredients/convexification_strategies/PrimalDualConvexificationStrategy.hpp index 82fa34d8..96bdd668 100644 --- a/uno/ingredients/convexification_strategies/PrimalDualConvexificationStrategy.hpp +++ b/uno/ingredients/convexification_strategies/PrimalDualConvexificationStrategy.hpp @@ -56,7 +56,6 @@ namespace uno { void PrimalDualConvexificationStrategy::regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, SymmetricMatrix& matrix, size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter) { - DEBUG2 << "Original matrix\n" << matrix << '\n'; this->primal_regularization = ElementType(0.); this->dual_regularization = ElementType(0.); DEBUG << "Testing factorization with regularization factors (0, 0)\n"; diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp index c28803e8..6606888c 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp +++ b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp @@ -1,6 +1,9 @@ // Copyright (c) 2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. +#ifndef UNO_CONVEXIFIEDHESSIAN_H +#define UNO_CONVEXIFIEDHESSIAN_H + #include #include "HessianModel.hpp" @@ -27,4 +30,6 @@ namespace uno { void regularize(Statistics& statistics, SymmetricMatrix& hessian, size_t number_original_variables); }; -} // namespace \ No newline at end of file +} // namespace + +#endif // UNO_CONVEXIFIEDHESSIAN_H \ No newline at end of file diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp index b50c5007..795f1eac 100644 --- a/uno/ingredients/hessian_models/ExactHessian.hpp +++ b/uno/ingredients/hessian_models/ExactHessian.hpp @@ -1,6 +1,9 @@ // Copyright (c) 2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. +#ifndef UNO_EXACTHESSIAN_H +#define UNO_EXACTHESSIAN_H + #include "HessianModel.hpp" namespace uno { @@ -13,4 +16,6 @@ namespace uno { void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector& primal_variables, const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; }; -} // namespace \ No newline at end of file +} // namespace + +#endif // UNO_EXACTHESSIAN_H \ No newline at end of file diff --git a/uno/ingredients/hessian_models/HessianModelFactory.cpp b/uno/ingredients/hessian_models/HessianModelFactory.cpp index 44f97893..360be9db 100644 --- a/uno/ingredients/hessian_models/HessianModelFactory.cpp +++ b/uno/ingredients/hessian_models/HessianModelFactory.cpp @@ -6,6 +6,7 @@ #include "HessianModel.hpp" #include "ConvexifiedHessian.hpp" #include "ExactHessian.hpp" +#include "IdentityHessian.hpp" #include "ZeroHessian.hpp" #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" @@ -23,6 +24,9 @@ namespace uno { else if (hessian_model == "zero") { return std::make_unique(); } + else if (hessian_model == "identity") { + return std::make_unique(); + } throw std::invalid_argument("Hessian model " + hessian_model + " does not exist"); } } // namespace \ No newline at end of file diff --git a/uno/ingredients/hessian_models/IdentityHessian.cpp b/uno/ingredients/hessian_models/IdentityHessian.cpp new file mode 100644 index 00000000..e7eaddb9 --- /dev/null +++ b/uno/ingredients/hessian_models/IdentityHessian.cpp @@ -0,0 +1,24 @@ +// Copyright (c) 2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "IdentityHessian.hpp" +#include "optimization/OptimizationProblem.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "options/Options.hpp" + +namespace uno { + // identity Hessian + IdentityHessian::IdentityHessian(): HessianModel() { + } + + void IdentityHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const { } + + void IdentityHessian::evaluate(Statistics& /*statistics*/, const OptimizationProblem& problem, const Vector& /*primal_variables*/, + const Vector& /*constraint_multipliers*/, SymmetricMatrix& hessian) { + // evaluate Lagrangian Hessian + hessian.set_dimension(problem.number_variables); + for (size_t variable_index: Range(problem.number_variables)) { + hessian.insert(1., variable_index, variable_index); + } + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/hessian_models/IdentityHessian.hpp b/uno/ingredients/hessian_models/IdentityHessian.hpp new file mode 100644 index 00000000..eb2e2bcc --- /dev/null +++ b/uno/ingredients/hessian_models/IdentityHessian.hpp @@ -0,0 +1,21 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_IDENTITYHESSIAN_H +#define UNO_IDENTITYHESSIAN_H + +#include "HessianModel.hpp" + +namespace uno { + // exact Hessian + class IdentityHessian : public HessianModel { + public: + IdentityHessian(); + + void initialize_statistics(Statistics& statistics, const Options& options) const override; + void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector& primal_variables, + const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; + }; +} // namespace + +#endif // UNO_IDENTITYHESSIAN_H \ No newline at end of file diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp index 7032a187..ea02efc0 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.hpp +++ b/uno/ingredients/hessian_models/ZeroHessian.hpp @@ -1,6 +1,9 @@ // Copyright (c) 2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. +#ifndef UNO_ZEROHESSIAN_H +#define UNO_ZEROHESSIAN_H + #include "HessianModel.hpp" namespace uno { @@ -13,4 +16,6 @@ namespace uno { void evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector& primal_variables, const Vector& constraint_multipliers, SymmetricMatrix& hessian) override; }; -} // namespace \ No newline at end of file +} // namespace + +#endif // UNO_ZEROHESSIAN_H \ No newline at end of file diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp index 4743053b..c8c66888 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp @@ -24,7 +24,7 @@ namespace uno { void LPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) { - LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *this->hessian_model, this->trust_region_radius); + LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius); this->solver->solve_LP(statistics, subproblem, this->initial_point, direction, warmstart_information); InequalityConstrainedMethod::compute_dual_displacements(current_multipliers, direction.multipliers); this->number_subproblems_solved++; diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp index bd8d8b15..b4341d4b 100644 --- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp @@ -37,7 +37,7 @@ namespace uno { void QPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) { - LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *this->hessian_model, this->trust_region_radius); + LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius); this->solver->solve_QP(statistics, subproblem, this->initial_point, direction, warmstart_information); InequalityConstrainedMethod::compute_dual_displacements(current_multipliers, direction.multipliers); this->number_subproblems_solved++; diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.cpp index 79e4ea6d..d0410776 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.cpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.cpp @@ -41,7 +41,7 @@ namespace uno { options.get_double("barrier_push_variable_to_interior_k1"), options.get_double("barrier_push_variable_to_interior_k2") }), - // least_square_multiplier_max_norm(options.get_double("least_square_multiplier_max_norm")), + least_square_multiplier_max_norm(options.get_double("least_square_multiplier_max_norm")), damping_factor(options.get_double("barrier_damping_factor")), l1_constraint_violation_coefficient(options.get_double("l1_constraint_violation_coefficient")) { } @@ -100,7 +100,8 @@ namespace uno { // compute least-square multipliers if (problem.is_constrained()) { - this->compute_least_square_multipliers(problem, initial_iterate, initial_iterate.multipliers.constraints); + Preprocessing::compute_least_square_multipliers(problem, *this->linear_solver, initial_iterate, initial_iterate.multipliers, + this->least_square_multiplier_max_norm); } } @@ -140,8 +141,7 @@ namespace uno { // create the barrier reformulation and the subproblem PrimalDualInteriorPointProblem barrier_problem(problem, current_multipliers, this->barrier_parameter()); - LagrangeNewtonSubproblem subproblem(barrier_problem, current_iterate, current_multipliers.constraints, *this->hessian_model, - this->trust_region_radius); + LagrangeNewtonSubproblem subproblem(barrier_problem, current_iterate, current_multipliers, *this->hessian_model, this->trust_region_radius); this->linear_solver->solve_indefinite_system(statistics, subproblem, this->solution, warmstart_information); assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality"); @@ -211,7 +211,8 @@ namespace uno { 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); + Preprocessing::compute_least_square_multipliers(problem, *this->linear_solver, trial_iterate, trial_iterate.multipliers, + this->least_square_multiplier_max_norm); } void PrimalDualInteriorPointMethod::set_auxiliary_measure(const Model& model, Iterate& iterate) { @@ -384,16 +385,6 @@ namespace uno { } } - 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(); - Preprocessing::compute_least_square_multipliers(problem.model, this->augmented_system.matrix, this->augmented_system.rhs, *this->linear_solver, - iterate, constraint_multipliers, this->least_square_multiplier_max_norm); - */ - } - 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()) { diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.hpp index 1512be3e..87687fc9 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/PrimalDualInteriorPointMethod.hpp @@ -58,7 +58,7 @@ namespace uno { double previous_barrier_parameter; const double default_multiplier; const InteriorPointParameters parameters; - // const double least_square_multiplier_max_norm; + const double least_square_multiplier_max_norm; const double damping_factor; // (Section 3.7 in IPOPT paper) const double l1_constraint_violation_coefficient; // (rho in Section 3.3.1 in IPOPT paper) @@ -81,7 +81,6 @@ namespace uno { Vector& direction_primals, Multipliers& direction_multipliers); void compute_bound_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, const Vector& primal_direction, Multipliers& direction_multipliers); - void compute_least_square_multipliers(const OptimizationProblem& problem, Iterate& iterate, Vector& constraint_multipliers); }; } // namespace diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp index ddc0d632..30615d4a 100644 --- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp +++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp @@ -151,16 +151,18 @@ namespace uno { // constraints and Jacobian if (warmstart_information.constraints_changed) { subproblem.evaluate_constraints(this->constraints); + } + if (warmstart_information.constraint_jacobian_changed) { subproblem.evaluate_constraint_jacobian(this->constraint_jacobian); } // Lagrangian Hessian - if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { + if (warmstart_information.objective_changed || warmstart_information.constraints_changed || warmstart_information.constraint_jacobian_changed) { DEBUG << "Evaluating problem Hessian\n"; subproblem.evaluate_hessian(statistics, this->hessian); } - if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { + if (warmstart_information.objective_changed || warmstart_information.constraint_jacobian_changed) { // form the KKT matrix this->augmented_matrix.set_dimension(subproblem.number_variables + subproblem.number_constraints); this->augmented_matrix.reset(); @@ -185,7 +187,7 @@ namespace uno { this->do_symbolic_analysis(this->augmented_matrix); warmstart_information.hessian_sparsity_changed = warmstart_information.jacobian_sparsity_changed = false; } - if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { + if (warmstart_information.objective_changed || warmstart_information.constraint_jacobian_changed) { DEBUG << "Performing numerical factorization of the augmented matrix\n"; this->do_numerical_factorization(this->augmented_matrix); // regularize diff --git a/uno/ingredients/subproblems/LagrangeNewtonSubproblem.cpp b/uno/ingredients/subproblems/LagrangeNewtonSubproblem.cpp index 76d1f32e..46abfcfc 100644 --- a/uno/ingredients/subproblems/LagrangeNewtonSubproblem.cpp +++ b/uno/ingredients/subproblems/LagrangeNewtonSubproblem.cpp @@ -8,7 +8,7 @@ namespace uno { LagrangeNewtonSubproblem::LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate, - const Vector& current_multipliers, HessianModel& hessian_model, double trust_region_radius): + const Multipliers& current_multipliers, HessianModel& hessian_model, double trust_region_radius): number_variables(problem.number_variables), number_constraints(problem.number_constraints), problem(problem), current_iterate(current_iterate), current_multipliers(current_multipliers), hessian_model(hessian_model), trust_region_radius(trust_region_radius) { } @@ -26,7 +26,7 @@ namespace uno { } void LagrangeNewtonSubproblem::evaluate_hessian(Statistics& statistics, SymmetricMatrix& hessian) { - this->hessian_model.evaluate(statistics, this->problem, this->current_iterate.primals, this->current_multipliers, hessian); + this->hessian_model.evaluate(statistics, this->problem, this->current_iterate.primals, this->current_multipliers.constraints, hessian); } void LagrangeNewtonSubproblem::compute_lagrangian_gradient(SparseVector& objective_gradient, RectangularMatrix& jacobian, @@ -40,12 +40,19 @@ namespace uno { // constraints for (size_t constraint_index: Range(this->number_constraints)) { - if (this->current_multipliers[constraint_index] != 0.) { + if (this->current_multipliers.constraints[constraint_index] != 0.) { for (const auto [variable_index, derivative]: jacobian[constraint_index]) { - gradient[variable_index] -= this->current_multipliers[constraint_index] * derivative; + gradient[variable_index] -= this->current_multipliers.constraints[constraint_index] * derivative; } } } + + /* + // bound multipliers + for (size_t variable_index: Range(this->number_variables)) { + gradient[variable_index] -= (this->current_multipliers.lower_bounds[variable_index] + this->current_multipliers.upper_bounds[variable_index]); + } + */ } double LagrangeNewtonSubproblem::dual_regularization_parameter() const { diff --git a/uno/ingredients/subproblems/LagrangeNewtonSubproblem.hpp b/uno/ingredients/subproblems/LagrangeNewtonSubproblem.hpp index 141a3ac2..f6e6f36e 100644 --- a/uno/ingredients/subproblems/LagrangeNewtonSubproblem.hpp +++ b/uno/ingredients/subproblems/LagrangeNewtonSubproblem.hpp @@ -13,6 +13,7 @@ namespace uno { // forward declarations class HessianModel; class Iterate; + class Multipliers; class OptimizationProblem; template class RectangularMatrix; @@ -29,7 +30,7 @@ namespace uno { const size_t number_variables; const size_t number_constraints; - LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate, const Vector& current_multipliers, + LagrangeNewtonSubproblem(const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, HessianModel& hessian_model, double trust_region_radius); void evaluate_objective_gradient(SparseVector& objective_gradient); @@ -48,7 +49,7 @@ namespace uno { protected: const OptimizationProblem& problem; Iterate& current_iterate; - const Vector& current_multipliers; + const Multipliers& current_multipliers; HessianModel& hessian_model; // TODO convexification model double trust_region_radius; diff --git a/uno/optimization/WarmstartInformation.cpp b/uno/optimization/WarmstartInformation.cpp index cdf5dcf3..5d3d266f 100644 --- a/uno/optimization/WarmstartInformation.cpp +++ b/uno/optimization/WarmstartInformation.cpp @@ -17,6 +17,7 @@ namespace uno { void WarmstartInformation::no_changes() { this->objective_changed = false; this->constraints_changed = false; + this->constraint_jacobian_changed = false; this->constraint_bounds_changed = false; this->variable_bounds_changed = false; this->hessian_sparsity_changed = false; @@ -26,6 +27,7 @@ namespace uno { void WarmstartInformation::iterate_changed() { this->objective_changed = true; this->constraints_changed = true; + this->constraint_jacobian_changed = true; this->constraint_bounds_changed = true; this->variable_bounds_changed = true; } @@ -33,6 +35,7 @@ namespace uno { void WarmstartInformation::whole_problem_changed() { this->objective_changed = true; this->constraints_changed = true; + this->constraint_jacobian_changed = true; this->constraint_bounds_changed = true; this->variable_bounds_changed = true; this->hessian_sparsity_changed = true; @@ -42,6 +45,7 @@ namespace uno { void WarmstartInformation::only_objective_changed() { this->objective_changed = true; this->constraints_changed = false; + this->constraint_jacobian_changed = false; this->constraint_bounds_changed = false; this->variable_bounds_changed = false; this->hessian_sparsity_changed = false; diff --git a/uno/optimization/WarmstartInformation.hpp b/uno/optimization/WarmstartInformation.hpp index e429b280..f32dde13 100644 --- a/uno/optimization/WarmstartInformation.hpp +++ b/uno/optimization/WarmstartInformation.hpp @@ -15,6 +15,7 @@ namespace uno { bool objective_changed{true}; bool constraints_changed{true}; + bool constraint_jacobian_changed{true}; bool constraint_bounds_changed{true}; bool variable_bounds_changed{true}; // bool problem_structure_changed{true}; diff --git a/uno/preprocessing/Preprocessing.cpp b/uno/preprocessing/Preprocessing.cpp index 2d548334..eb8123d8 100644 --- a/uno/preprocessing/Preprocessing.cpp +++ b/uno/preprocessing/Preprocessing.cpp @@ -2,80 +2,51 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "Preprocessing.hpp" +#include "ingredients/hessian_models/IdentityHessian.hpp" #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "ingredients/subproblem_solvers/QPSolver.hpp" +#include "ingredients/subproblems/LagrangeNewtonSubproblem.hpp" #include "linear_algebra/SymmetricMatrix.hpp" #include "linear_algebra/RectangularMatrix.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/Iterate.hpp" #include "optimization/WarmstartInformation.hpp" +#include "options/Options.hpp" #include "symbolic/VectorView.hpp" +#include "tools/Statistics.hpp" namespace uno { // compute a least-square approximation of the multipliers by solving a linear system - void Preprocessing::compute_least_square_multipliers(const Model& model, SymmetricMatrix& matrix, Vector& rhs, - DirectSymmetricIndefiniteLinearSolver& linear_solver, Iterate& current_iterate, Vector& multipliers, - double multiplier_max_norm) { - current_iterate.evaluate_objective_gradient(model); - current_iterate.evaluate_constraint_jacobian(model); + void Preprocessing::compute_least_square_multipliers(const OptimizationProblem& problem, DirectSymmetricIndefiniteLinearSolver& linear_solver, + Iterate& current_iterate, Multipliers& multipliers, double multiplier_max_norm) { DEBUG << "Computing least-square multipliers\n"; DEBUG2 << "Current primals: " << current_iterate.primals << '\n'; - /* generate the right-hand side */ - rhs.fill(0.); - // objective gradient - for (const auto [variable_index, derivative]: current_iterate.evaluations.objective_gradient) { - rhs[variable_index] += model.objective_sign * derivative; - } - // variable bound constraints - for (size_t variable_index: Range(model.number_variables)) { - rhs[variable_index] -= current_iterate.multipliers.lower_bounds[variable_index] + current_iterate.multipliers.upper_bounds[variable_index]; - } - DEBUG2 << "RHS for least-square multipliers: "; print_vector(DEBUG2, view(rhs, 0, model.number_variables + model.number_constraints)); - - // if the residuals on the RHS are all 0, the least-square multipliers are all 0 - if (norm_inf(view(rhs, 0, model.number_variables + model.number_constraints)) == 0.) { - multipliers.fill(0.); - DEBUG << "Least-square multipliers are all 0.\n"; - return; - } - - /* build the symmetric matrix */ - matrix.reset(); - // identity block - for (size_t variable_index: Range(model.number_variables)) { - matrix.insert(1., variable_index, variable_index); - matrix.finalize_column(variable_index); - } - // Jacobian of general constraints - for (size_t constraint_index: Range(model.number_constraints)) { - for (const auto [variable_index, derivative]: current_iterate.evaluations.constraint_jacobian[constraint_index]) { - matrix.insert(derivative, variable_index, model.number_variables + constraint_index); - } - matrix.finalize_column(model.number_variables + constraint_index); - } - DEBUG2 << "Matrix for least-square multipliers:\n" << matrix << '\n'; - - /* solve the system */ - Vector solution(matrix.dimension()); - linear_solver.do_symbolic_analysis(matrix); - linear_solver.do_numerical_factorization(matrix); - linear_solver.solve_indefinite_system(matrix, rhs, solution); + // solve the system + IdentityHessian hessian_model{}; + LagrangeNewtonSubproblem subproblem(problem, current_iterate, multipliers, hessian_model, INF); + Vector solution(problem.number_variables + problem.number_constraints); + Options options(false); + options["statistics_print_header_frequency"] = "15"; + Statistics statistics(options); + // do not evaluate the constraints + WarmstartInformation warmstart_information{}; + warmstart_information.constraints_changed = false; + linear_solver.solve_indefinite_system(statistics, subproblem, solution, warmstart_information); // if least-square multipliers too big, discard them. Otherwise, keep them - const auto trial_multipliers = view(solution, model.number_variables, model.number_variables + model.number_constraints); + const auto trial_multipliers = view(-solution, problem.number_variables, problem.number_variables + problem.number_constraints); DEBUG2 << "Trial multipliers: "; print_vector(DEBUG2, trial_multipliers); if (norm_inf(trial_multipliers) <= multiplier_max_norm) { - multipliers = trial_multipliers; + multipliers.constraints = trial_multipliers; } else { DEBUG << "Ignoring the least-square multipliers\n"; } - DEBUG << '\n'; } - size_t count_infeasible_linear_constraints(const Model& model, const std::vector& constraint_values) { + size_t count_infeasible_linear_constraints(const Model& model, const Vector& constraint_values) { size_t infeasible_linear_constraints = 0; for (size_t constraint_index: model.get_linear_constraints()) { if (constraint_values[constraint_index] < model.constraint_lower_bound(constraint_index) || @@ -88,13 +59,13 @@ namespace uno { void Preprocessing::enforce_linear_constraints(const Model& /*model*/, Vector& /*primals*/, Multipliers& /*multipliers*/, QPSolver& /*qp_solver*/) { - WARNING << "Preprocessing::enforce_linear_constraints not implemented yet\n"; /* + WARNING << "Preprocessing::enforce_linear_constraints not implemented yet\n"; const auto& linear_constraints = model.get_linear_constraints(); INFO << "\nPreprocessing phase: the problem has " << linear_constraints.size() << " linear constraints\n"; if (not linear_constraints.empty()) { // evaluate the constraints - std::vector constraints(model.number_constraints); + Vector constraints(model.number_constraints); model.evaluate_constraints(primals, constraints); const size_t infeasible_linear_constraints = count_infeasible_linear_constraints(model, constraints); INFO << "There are " << infeasible_linear_constraints << " infeasible linear constraints at the initial point\n"; diff --git a/uno/preprocessing/Preprocessing.hpp b/uno/preprocessing/Preprocessing.hpp index bd8809f7..b8df59cc 100644 --- a/uno/preprocessing/Preprocessing.hpp +++ b/uno/preprocessing/Preprocessing.hpp @@ -12,6 +12,7 @@ namespace uno { class Iterate; class Model; class Multipliers; + class OptimizationProblem; class QPSolver; template class DirectSymmetricIndefiniteLinearSolver; @@ -22,9 +23,8 @@ namespace uno { class Preprocessing { public: - static void compute_least_square_multipliers(const Model& model, SymmetricMatrix& matrix, Vector& rhs, - DirectSymmetricIndefiniteLinearSolver& linear_solver, Iterate& current_iterate, Vector& multipliers, - double multiplier_max_norm); + static void compute_least_square_multipliers(const OptimizationProblem& problem, DirectSymmetricIndefiniteLinearSolver& linear_solver, + Iterate& current_iterate, Multipliers& multipliers, double multiplier_max_norm); static void enforce_linear_constraints(const Model& model, Vector& primals, Multipliers& multipliers, QPSolver& qp_solver); }; } // namespace