Introduce the new ingredient LagrangeNewtonSubproblem

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;
       }
@@ -134,7 +134,7 @@ namespace uno {
    }
    */
 
-   void AMPLModel::evaluate_constraints(const Vector<double>& x, std::vector<double>& constraints) const {
+   void AMPLModel::evaluate_constraints(const Vector<double>& x, Vector<double>& constraints) const {
       fint error_flag = 0;
       (*(this->asl)->p.Conval)(this->asl, const_cast<double*>(x.data()), constraints.data(), &error_flag);
       if (0 < error_flag) {

bindings/AMPL/AMPLModel.hpp

@@ -13,8 +13,8 @@
 
 // include AMPL Solver Library (ASL)
 extern "C" {
-#include "asl_pfgh.h"
-#include "getstub.h"
+   #include "asl_pfgh.h"
+   #include "getstub.h"
 }
 
 namespace uno {
@@ -33,7 +33,7 @@ namespace uno {
 
       [[nodiscard]] double evaluate_objective(const Vector<double>& x) const override;
       void evaluate_objective_gradient(const Vector<double>& x, SparseVector<double>& gradient) const override;
-      void evaluate_constraints(const Vector<double>& x, std::vector<double>& constraints) const override;
+      void evaluate_constraints(const Vector<double>& x, Vector<double>& constraints) const override;
       void evaluate_constraint_gradient(const Vector<double>& x, size_t constraint_index, SparseVector<double>& gradient) const override;
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,

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/ConstraintRelaxationStrategy.cpp

@@ -2,16 +2,15 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ConstraintRelaxationStrategy.hpp"
-#include "OptimizationProblem.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp"
 #include "optimization/Direction.hpp"
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.hpp"
-#include "linear_algebra/SymmetricMatrix.hpp"
 #include "model/Model.hpp"
 #include "optimization/Iterate.hpp"
 #include "optimization/Multipliers.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "symbolic/VectorView.hpp"
 #include "symbolic/Expression.hpp"
 #include "options/Options.hpp"

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

@@ -7,8 +7,8 @@
 #include <memory>
 #include "ConstraintRelaxationStrategy.hpp"
 #include "ingredients/globalization_strategies/ProgressMeasures.hpp"
-#include "OptimalityProblem.hpp"
 #include "l1RelaxedProblem.hpp"
+#include "OptimalityProblem.hpp"
 
 namespace uno {
    enum class Phase {FEASIBILITY_RESTORATION = 1, OPTIMALITY = 2};
@@ -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/constraint_relaxation_strategies/OptimalityProblem.cpp

@@ -16,7 +16,7 @@ namespace uno {
       objective_gradient = iterate.evaluations.objective_gradient;
    }
 
-   void OptimalityProblem::evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const {
+   void OptimalityProblem::evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const {
       iterate.evaluate_constraints(this->model);
       constraints = iterate.evaluations.constraints;
    }
@@ -59,7 +59,7 @@ namespace uno {
       }
    }
 
-   double OptimalityProblem::complementarity_error(const Vector<double>& primals, const std::vector<double>& constraints,
+   double OptimalityProblem::complementarity_error(const Vector<double>& primals, const Vector<double>& constraints,
          const Multipliers& multipliers, double shift_value, Norm residual_norm) const {
       // bound constraints
       const Range variables_range = Range(this->model.number_variables);

uno/ingredients/constraint_relaxation_strategies/OptimalityProblem.hpp

@@ -4,7 +4,7 @@
 #ifndef UNO_OPTIMALITYPROBLEM_H
 #define UNO_OPTIMALITYPROBLEM_H
 
-#include "OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 
 namespace uno {
    class OptimalityProblem: public OptimizationProblem {
@@ -13,7 +13,7 @@ namespace uno {
 
       [[nodiscard]] double get_objective_multiplier() const override { return 1.; }
       void evaluate_objective_gradient(Iterate& iterate, SparseVector<double>& objective_gradient) const override;
-      void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
+      void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
 
@@ -32,8 +32,8 @@ namespace uno {
       [[nodiscard]] size_t number_hessian_nonzeros() const override { return this->model.number_hessian_nonzeros(); }
 
       void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate, const Multipliers& multipliers) const override;
-      [[nodiscard]] double complementarity_error(const Vector<double>& primals, const std::vector<double>& constraints,
-            const Multipliers& multipliers, double shift_value, Norm residual_norm) const override;
+      [[nodiscard]] double complementarity_error(const Vector<double>& primals, const Vector<double>& constraints, const Multipliers& multipliers,
+            double shift_value, Norm residual_norm) const override;
    };
 } // namespace
 

uno/ingredients/constraint_relaxation_strategies/l1RelaxedProblem.cpp

@@ -69,7 +69,7 @@ namespace uno {
       }
    }
 
-   void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const {
+   void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const {
       iterate.evaluate_constraints(this->model);
       constraints = iterate.evaluations.constraints;
       // add the contribution of the elastics
@@ -161,8 +161,8 @@ namespace uno {
    }
 
    // complementary slackness error: expression for violated constraints depends on the definition of the relaxed problem
-   double l1RelaxedProblem::complementarity_error(const Vector<double>& primals, const std::vector<double>& constraints,
-         const Multipliers& multipliers, double shift_value, Norm residual_norm) const {
+   double l1RelaxedProblem::complementarity_error(const Vector<double>& primals, const Vector<double>& constraints, const Multipliers& multipliers,
+         double shift_value, Norm residual_norm) const {
       // bound constraints
       const Range variables_range = Range(this->number_variables);
       const VectorExpression bounds_complementarity{variables_range, [&](size_t variable_index) {

uno/ingredients/constraint_relaxation_strategies/l1RelaxedProblem.hpp

@@ -4,7 +4,7 @@
 #ifndef UNO_L1RELAXEDPROBLEM_H
 #define UNO_L1RELAXEDPROBLEM_H
 
-#include "OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "ElasticVariables.hpp"
 #include "symbolic/Concatenation.hpp"
 
@@ -16,7 +16,7 @@ namespace uno {
 
       [[nodiscard]] double get_objective_multiplier() const override;
       void evaluate_objective_gradient(Iterate& iterate, SparseVector<double>& objective_gradient) const override;
-      void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
+      void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
 
@@ -35,8 +35,8 @@ namespace uno {
       [[nodiscard]] size_t number_hessian_nonzeros() const override;
 
       void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate, const Multipliers& multipliers) const override;
-      [[nodiscard]] double complementarity_error(const Vector<double>& primals, const std::vector<double>& constraints,
-            const Multipliers& multipliers, double shift_value, Norm residual_norm) const override;
+      [[nodiscard]] double complementarity_error(const Vector<double>& primals, const Vector<double>& constraints, const Multipliers& multipliers,
+            double shift_value, Norm residual_norm) const override;
 
       // parameterization
       void set_objective_multiplier(double new_objective_multiplier);

uno/ingredients/hessian_models/ConvexifiedHessian.cpp

@@ -2,9 +2,9 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ConvexifiedHessian.hpp"
-#include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "ingredients/hessian_models/UnstableRegularization.hpp"
-#include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
 #include "linear_algebra/SymmetricMatrix.hpp"

uno/ingredients/hessian_models/ExactHessian.cpp

@@ -2,7 +2,7 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ExactHessian.hpp"
-#include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
 

uno/ingredients/hessian_models/ZeroHessian.cpp

@@ -2,7 +2,7 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "ZeroHessian.hpp"
-#include "ingredients/constraint_relaxation_strategies/OptimizationProblem.hpp"
+#include "optimization/OptimizationProblem.hpp"
 #include "linear_algebra/SymmetricMatrix.hpp"
 #include "options/Options.hpp"
 

uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp

@@ -4,32 +4,33 @@
 #include <string>
 #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");
-      // active-set methods
-      if (subproblem_strategy == "QP") {
-         return std::make_unique<QPSubproblem>(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
+      const std::string inequality_handling_method = options.get_string("inequality_handling_method");
+      // inequality-constrained methods
+      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() {
@@ -38,6 +39,9 @@ namespace uno {
          strategies.emplace_back("QP");
          strategies.emplace_back("LP");
       }
+      else if (not LPSolverFactory::available_solvers().empty()) {
+         strategies.emplace_back("LP");
+      }
       if (not SymmetricIndefiniteLinearSolverFactory::available_solvers().empty()) {
          strategies.emplace_back("primal_dual_interior_point");
       }