Introduced LagrangeNewtonSubproblem class that is created in QPMethod…

…/LPMethod and passed on to the QP/LP solver

CMakeLists.txt

@@ -53,6 +53,7 @@ file(GLOB UNO_SOURCE_FILES
    uno/ingredients/inequality_handling_methods/inequality_constrained_methods/*.cpp
    uno/ingredients/inequality_handling_methods/interior_point_methods/*.cpp
    uno/ingredients/subproblem_solvers/*.cpp
+   uno/ingredients/subproblems/*.cpp
    uno/model/*.cpp
    uno/optimization/*.cpp
    uno/options/*.cpp

README.md

@@ -6,7 +6,7 @@
 
 </div>
 
-Uno (Unifying Nonlinear Optimization) is a C++ library that unifies methods for solving nonlinearly constrained optimization problems of the form:
+Uno (Unifying Nonlinear Optimization) is a C++ Lagrange-Newton solver that unifies methods for solving nonlinearly constrained optimization problems of the form:
 
 $$
 \begin{align}
@@ -25,14 +25,15 @@ Oscar Dowson [@odow](https://github.com/odow), David Kiessling [@david0oo](https
 
 ## Unifying nonlinearly constrained nonconvex optimization
 
-We argue that most optimization methods can be broken down into four generic ingredients:
+We argue that most optimization methods can be broken down into the following generic ingredients:
 * a **constraint relaxation strategy**: a systematic way to relax the nonlinear constraints;
-* a **subproblem**: a local model of the (possibly relaxed) problem at the current primal-dual iterate;
+* an **inequality handling method**: a way of handling inequality constraints;
+* a **subproblem**: a local model of the (possibly relaxed) problem at the current primal-dual iterate. Uno currently only supports Lagrange-Newton;
 * a **globalization strategy**: an acceptance test of the trial iterate;
 * a **globalization mechanism**: a recourse action upon rejection of the trial iterate.
 
 <!--
-The following hypergraph illustrates how some of the state-of-the-art solvers can be decomposed in terms of the four ingredients:
+The following hypergraph illustrates how some of the state-of-the-art solvers can be decomposed in terms of the basic ingredients:
 <p align="center">
    <img src="docs/figures/combination_hypergraph.png" alt="Combination hypergraph" width="75%" />
 </p>
@@ -112,7 +113,7 @@ Uno can be installed in Julia via [Uno_jll.jl](https://github.com/JuliaBinaryWra
 ### Combining strategies on the fly
 
 For an overview of the available strategies, type: ```./uno_ampl --strategies```:
+- to pick a constraint relaxation strategy, use the argument: ```constraint_relaxation_strategy=[feasibility_restoration|l1_relaxation]``` 
+- to pick an inequality handling method, use the argument: ```inequality_handling_method=[QP|LP|primal_dual_interior_point]```  
+ - to pick a globalization strategy, use the argument: ```globalization_strategy=[l1_merit|fletcher_filter_method|waechter_filter_method|funnel_method]```  
 - to pick a globalization mechanism, use the argument : ```globalization_mechanism=[LS|TR]```  
-- to pick a constraint relaxation strategy, use the argument: ```constraint_relaxation_strategy=[feasibility_restoration|l1_relaxation]```  
-- to pick a globalization strategy, use the argument: ```globalization_strategy=[l1_merit|fletcher_filter_method|waechter_filter_method|funnel_method]```  
-- to pick a subproblem method, use the argument: ```subproblem=[QP|LP|primal_dual_interior_point]```  

bindings/AMPL/AMPLModel.cpp

@@ -101,8 +101,8 @@ namespace uno {
       fint error_flag = 0;
       // prevent ASL to crash by catching all evaluation errors
       Jmp_buf err_jmp_uno;
-      asl->i.err_jmp_ = &err_jmp_uno;
-      asl->i.err_jmp1_ = &err_jmp_uno;
+      this->asl->i.err_jmp_ = &err_jmp_uno;
+      this->asl->i.err_jmp1_ = &err_jmp_uno;
       if (setjmp(err_jmp_uno.jb)) {
          error_flag = 1;
       }

bindings/AMPL/uno_ampl.cpp

@@ -74,7 +74,7 @@ namespace uno {
       std::cout << "To solve an AMPL model, type ./uno_ampl model.nl -AMPL [option_name=option_value ...]\n";
       std::cout << "To choose a constraint relaxation strategy, use the argument constraint_relaxation_strategy="
                    "[feasibility_restoration|l1_relaxation]\n";
-      std::cout << "To choose a subproblem method, use the argument subproblem=[QP|LP|primal_dual_interior_point]\n";
+      std::cout << "To choose an inequality handling method, use the argument inequality_handling_method=[QP|LP|primal_dual_interior_point]\n";
       std::cout << "To choose a globalization mechanism, use the argument globalization_mechanism=[LS|TR]\n";
       std::cout << "To choose a globalization strategy, use the argument globalization_strategy="
                    "[l1_merit|fletcher_filter_method|waechter_filter_method]\n";

uno/Uno.cpp

@@ -166,7 +166,7 @@ namespace uno {
 
    std::string Uno::get_strategy_combination(const Options& options) {
       return options.get_string("globalization_mechanism") + " " + options.get_string("constraint_relaxation_strategy") + " " +
-                                options.get_string("globalization_strategy") + " " + options.get_string("subproblem");
+                                options.get_string("globalization_strategy") + " " + options.get_string("inequality_handling_method");
 
    }
 

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp

@@ -36,7 +36,7 @@ namespace uno {
                options),
          optimality_problem(std::forward<OptimalityProblem>(optimality_problem)),
          feasibility_problem(std::forward<l1RelaxedProblem>(feasibility_problem)),
-         subproblem_strategy(options.get_string("subproblem")),
+         inequality_handling_method_name(options.get_string("inequality_handling_method")),
          linear_feasibility_tolerance(options.get_double("tolerance")),
          switch_to_optimality_requires_linearized_feasibility(options.get_bool("switch_to_optimality_requires_linearized_feasibility")),
          reference_optimality_primals(optimality_problem.number_variables) {
@@ -71,7 +71,7 @@ namespace uno {
             this->solve_subproblem(statistics, this->optimality_problem, current_iterate, current_iterate.multipliers, direction, warmstart_information);
             if (direction.status == SubproblemStatus::INFEASIBLE) {
                // switch to the feasibility problem, starting from the current direction
-               statistics.set("status", std::string("infeasible " + this->subproblem_strategy));
+               statistics.set("status", std::string("infeasible " + this->inequality_handling_method_name));
                DEBUG << "/!\\ The subproblem is infeasible\n";
                this->switch_to_feasibility_problem(statistics, current_iterate, warmstart_information);
                this->inequality_handling_method->set_initial_point(direction.primals);

uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp

@@ -39,7 +39,7 @@ namespace uno {
       const OptimalityProblem optimality_problem;
       l1RelaxedProblem feasibility_problem;
       Phase current_phase{Phase::OPTIMALITY};
-      const std::string& subproblem_strategy;
+      const std::string& inequality_handling_method_name;
       const double linear_feasibility_tolerance;
       const bool switch_to_optimality_requires_linearized_feasibility;
       ProgressMeasures reference_optimality_progress{};

uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp

@@ -2,35 +2,35 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include <string>
-#include <ingredients/subproblem_solvers/LPSolverFactory.hpp>
 #include "InequalityHandlingMethod.hpp"
 #include "InequalityHandlingMethodFactory.hpp"
-#include "inequality_constrained_methods/QPSubproblem.hpp"
-#include "inequality_constrained_methods/LPSubproblem.hpp"
+#include "inequality_constrained_methods/QPMethod.hpp"
+#include "inequality_constrained_methods/LPMethod.hpp"
 #include "interior_point_methods/PrimalDualInteriorPointMethod.hpp"
+#include "ingredients/subproblem_solvers/LPSolverFactory.hpp"
 #include "ingredients/subproblem_solvers/QPSolverFactory.hpp"
 #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
 #include "options/Options.hpp"
 
 namespace uno {
    std::unique_ptr<InequalityHandlingMethod> InequalityHandlingMethodFactory::create(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros,
          size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options) {
-      const std::string subproblem_strategy = options.get_string("subproblem");
+      const std::string inequality_handling_method = options.get_string("inequality_handling_method");
       // inequality-constrained methods
-      if (subproblem_strategy == "QP") {
-         return std::make_unique<QPSubproblem>(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
+      if (inequality_handling_method == "QP") {
+         return std::make_unique<QPMethod>(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
                number_hessian_nonzeros, options);
       }
-      else if (subproblem_strategy == "LP") {
-         return std::make_unique<LPSubproblem>(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
+      else if (inequality_handling_method == "LP") {
+         return std::make_unique<LPMethod>(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
                options);
       }
       // interior-point method
-      else if (subproblem_strategy == "primal_dual_interior_point") {
+      else if (inequality_handling_method == "primal_dual_interior_point") {
          return std::make_unique<PrimalDualInteriorPointMethod>(number_variables, number_constraints, number_jacobian_nonzeros,
                number_hessian_nonzeros, options);
       }
-      throw std::invalid_argument("Subproblem strategy " + subproblem_strategy + " is not supported");
+      throw std::invalid_argument("Subproblem strategy " + inequality_handling_method + " is not supported");
    }
 
    std::vector<std::string> InequalityHandlingMethodFactory::available_strategies() {

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPSubproblem.cpp → uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.cpp

@@ -1,37 +1,39 @@
 // Copyright (c) 2018-2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
-#include "LPSubproblem.hpp"
-#include "optimization/Direction.hpp"
-#include "optimization/WarmstartInformation.hpp"
+#include "LPMethod.hpp"
 #include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
 #include "ingredients/subproblem_solvers/LPSolver.hpp"
 #include "ingredients/subproblem_solvers/LPSolverFactory.hpp"
+#include "ingredients/subproblems/LagrangeNewtonSubproblem.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/WarmstartInformation.hpp"
 #include "options/Options.hpp"
 
 namespace uno {
-   LPSubproblem::LPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros,
+   LPMethod::LPMethod(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros,
          size_t number_jacobian_nonzeros, const Options& options) :
          InequalityConstrainedMethod("zero", number_variables, number_constraints, 0, false, options),
          solver(LPSolverFactory::create(number_variables, number_constraints,
                number_objective_gradient_nonzeros, number_jacobian_nonzeros, options)) {
    }
 
-   LPSubproblem::~LPSubproblem() { }
+   LPMethod::~LPMethod() { }
 
-   void LPSubproblem::generate_initial_iterate(const OptimizationProblem& /*problem*/, Iterate& /*initial_iterate*/) {
+   void LPMethod::generate_initial_iterate(const OptimizationProblem& /*problem*/, Iterate& /*initial_iterate*/) {
    }
 
-   void LPSubproblem::solve(Statistics& /*statistics*/, const OptimizationProblem& problem, Iterate& current_iterate,
+   void LPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
          const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) {
-      this->solver->solve_LP(problem, current_iterate, this->initial_point, direction, this->trust_region_radius, warmstart_information);
+      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *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++;
       // reset the initial point
       this->initial_point.fill(0.);
    }
 
-   double LPSubproblem::hessian_quadratic_product(const Vector<double>& /*primal_direction*/) const {
+   double LPMethod::hessian_quadratic_product(const Vector<double>& /*primal_direction*/) const {
       return 0.;
    }
 } // namespace

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPSubproblem.hpp → uno/ingredients/inequality_handling_methods/inequality_constrained_methods/LPMethod.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_LPSUBPROBLEM_H
-#define UNO_LPSUBPROBLEM_H
+#ifndef UNO_LPMETHOD_H
+#define UNO_LPMETHOD_H
 
 #include <memory>
 #include "InequalityConstrainedMethod.hpp"
@@ -11,11 +11,11 @@ namespace uno {
    // forward reference
    class LPSolver;
 
-   class LPSubproblem : public InequalityConstrainedMethod {
+   class LPMethod : public InequalityConstrainedMethod {
    public:
-      LPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
+      LPMethod(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
             const Options& options);
-      ~LPSubproblem();
+      ~LPMethod();
 
       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,
@@ -28,4 +28,4 @@ namespace uno {
    };
 } // namespace
 
-#endif // UNO_LPSUBPROBLEM_H
+#endif // UNO_LPMETHOD_H

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPSubproblem.cpp → uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.cpp

@@ -1,20 +1,21 @@
 // Copyright (c) 2018-2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
-#include "QPSubproblem.hpp"
-#include "optimization/Direction.hpp"
-#include "optimization/Iterate.hpp"
-#include "optimization/WarmstartInformation.hpp"
-#include "preprocessing/Preprocessing.hpp"
+#include "QPMethod.hpp"
 #include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
 #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "ingredients/subproblem_solvers/QPSolver.hpp"
 #include "ingredients/subproblem_solvers/QPSolverFactory.hpp"
+#include "ingredients/subproblems/LagrangeNewtonSubproblem.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/Iterate.hpp"
+#include "optimization/WarmstartInformation.hpp"
 #include "options/Options.hpp"
+#include "preprocessing/Preprocessing.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
-   QPSubproblem::QPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros,
+   QPMethod::QPMethod(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros,
          size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options) :
          InequalityConstrainedMethod(options.get_string("hessian_model"), number_variables, number_constraints, number_hessian_nonzeros,
                options.get_string("globalization_mechanism") != "TR" || options.get_bool("convexify_QP"), options),
@@ -26,25 +27,25 @@ namespace uno {
                options)) {
    }
 
-   QPSubproblem::~QPSubproblem() { }
+   QPMethod::~QPMethod() { }
 
-   void QPSubproblem::generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) {
+   void QPMethod::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->solver);
       }
    }
 
-   void QPSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
+   void QPMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,
          Direction& direction, WarmstartInformation& warmstart_information) {
-      this->solver->solve_QP(statistics, problem, current_iterate, current_multipliers.constraints, this->initial_point, direction,
-            *this->hessian_model, this->trust_region_radius, warmstart_information);
+      LagrangeNewtonSubproblem subproblem(problem, current_iterate, current_multipliers.constraints, *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++;
       // reset the initial point
       this->initial_point.fill(0.);
    }
 
-   double QPSubproblem::hessian_quadratic_product(const Vector<double>& primal_direction) const {
+   double QPMethod::hessian_quadratic_product(const Vector<double>& primal_direction) const {
       return this->solver->hessian_quadratic_product(primal_direction);
    }
 } // namespace

uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPSubproblem.hpp → uno/ingredients/inequality_handling_methods/inequality_constrained_methods/QPMethod.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_QPSUBPROBLEM_H
-#define UNO_QPSUBPROBLEM_H
+#ifndef UNO_QPMETHOD_H
+#define UNO_QPMETHOD_H
 
 #include <memory>
 #include "InequalityConstrainedMethod.hpp"
@@ -11,11 +11,11 @@ namespace uno {
    // forward reference
    class QPSolver;
 
-   class QPSubproblem : public InequalityConstrainedMethod {
+   class QPMethod : public InequalityConstrainedMethod {
    public:
-      QPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
+      QPMethod(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
             size_t number_hessian_nonzeros, const Options& options);
-      ~QPSubproblem();
+      ~QPMethod();
 
       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,
@@ -29,4 +29,4 @@ namespace uno {
    };
 } // namespace
 
-#endif // UNO_QPSUBPROBLEM_H
+#endif // UNO_QPMETHOD_H