From f118d9203a10b668ab25ab4156281a71eed507ef Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Mon, 3 Nov 2025 19:56:07 +0100
Subject: [PATCH 01/10] Added new class l1Relaxation (for now, based on
 UnconstrainedStrategy)

---
 .../ConstraintRelaxationStrategyFactory.cpp   |   4 +
 .../l1Relaxation.cpp                          | 107 ++++++++++++++++++
 .../l1Relaxation.hpp                          |  49 ++++++++
 3 files changed, 160 insertions(+)
 create mode 100644 uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
 create mode 100644 uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp

diff --git a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.cpp b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.cpp
index f15ba9721..1318e6a59 100644
--- a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.cpp
@@ -4,6 +4,7 @@
 #include <string>
 #include "ConstraintRelaxationStrategyFactory.hpp"
 #include "FeasibilityRestoration.hpp"
+#include "l1Relaxation.hpp"
 #include "UnconstrainedStrategy.hpp"
 #include "options/Options.hpp"
 #include "tools/Logger.hpp"
@@ -20,6 +21,9 @@ namespace uno {
       if (constraint_relaxation_type == "feasibility_restoration") {
          return std::make_unique<FeasibilityRestoration>(options);
       }
+      else if (constraint_relaxation_type == "l1_relaxation") {
+         return std::make_unique<l1Relaxation>(options);
+      }
       throw std::invalid_argument("ConstraintRelaxationStrategy " + constraint_relaxation_type + " is not supported");
    }
 } // namespace
\ No newline at end of file
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
new file mode 100644
index 000000000..9eb522f70
--- /dev/null
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -0,0 +1,107 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include <cassert>
+#include "l1Relaxation.hpp"
+#include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
+#include "ingredients/hessian_models/HessianModelFactory.hpp"
+#include "ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.hpp"
+#include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategyFactory.hpp"
+#include "optimization/Direction.hpp"
+#include "optimization/Iterate.hpp"
+#include "optimization/OptimizationProblem.hpp"
+#include "optimization/WarmstartInformation.hpp"
+#include "options/Options.hpp"
+#include "symbolic/VectorView.hpp"
+#include "tools/Logger.hpp"
+
+namespace uno {
+   l1Relaxation::l1Relaxation(const Options& options) :
+         ConstraintRelaxationStrategy(options),
+         inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
+         inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
+         constraint_violation_coefficient(options.get_double(("l1_constraint_violation_coefficient"))) {
+   }
+
+   void l1Relaxation::initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate,
+         Direction& direction, double trust_region_radius, const Options& options) {
+      this->relaxed_problem = std::make_unique<const l1RelaxedProblem>(model, this->penalty_parameter,
+         this->constraint_violation_coefficient);
+      assert(this->relaxed_problem != nullptr);
+      this->feasibility_problem = std::make_unique<const l1RelaxedProblem>(model, 0., 1.);
+      assert(this->feasibility_problem != nullptr);
+
+      // Hessian model
+      this->hessian_model = HessianModelFactory::create(model, options);
+
+      // memory allocation
+      this->inequality_handling_method->initialize(*this->relaxed_problem, initial_iterate, *this->hessian_model,
+         *this->inertia_correction_strategy, trust_region_radius);
+      direction = Direction(this->relaxed_problem->number_variables, this->relaxed_problem->number_constraints);
+
+      // statistics
+      this->inertia_correction_strategy->initialize_statistics(statistics, options);
+      this->inequality_handling_method->initialize_statistics(statistics, options);
+
+      // initial iterate
+      this->inequality_handling_method->generate_initial_iterate(initial_iterate);
+      initial_iterate.evaluate_objective_gradient(model);
+      initial_iterate.evaluate_constraints(model);
+      this->inequality_handling_method->evaluate_constraint_jacobian(initial_iterate);
+      this->relaxed_problem->evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, *this->inequality_handling_method,
+         initial_iterate);
+      this->compute_primal_dual_residuals(*this->relaxed_problem, initial_iterate);
+   }
+
+   void l1Relaxation::compute_feasible_direction(Statistics& statistics, GlobalizationStrategy& /*globalization_strategy*/,
+         Iterate& current_iterate, Direction& direction, double trust_region_radius, WarmstartInformation& warmstart_information) {
+      direction.reset();
+      DEBUG << "Solving the subproblem\n";
+      direction.set_dimensions(this->relaxed_problem->number_variables, this->relaxed_problem->number_constraints);
+      this->inequality_handling_method->solve(statistics, current_iterate, direction, *this->hessian_model,
+         *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
+      direction.norm = norm_inf(view(direction.primals, 0, this->relaxed_problem->get_number_original_variables()));
+      DEBUG3 << direction << '\n';
+      warmstart_information.no_changes();
+   }
+
+   bool l1Relaxation::solving_feasibility_problem() const {
+      return false;
+   }
+
+   void l1Relaxation::switch_to_feasibility_problem(Statistics& /*statistics*/, GlobalizationStrategy& /*globalization_strategy*/,
+         Iterate& /*current_iterate*/, double /*trust_region_radius*/,
+         WarmstartInformation& /*warmstart_information*/) {
+      throw std::runtime_error("The problem is unconstrained, switching to the feasibility problem should not happen");
+   }
+
+   bool l1Relaxation::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         double trust_region_radius, const Model& model, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, WarmstartInformation& warmstart_information, UserCallbacks& user_callbacks) {
+      const bool accept_iterate = this->inequality_handling_method->is_iterate_acceptable(statistics, globalization_strategy,
+         *this->hessian_model, *this->inertia_correction_strategy, trust_region_radius, current_iterate, trial_iterate,
+         direction, step_length, user_callbacks);
+      trial_iterate.status = this->check_termination(model, trial_iterate);
+      warmstart_information.no_changes();
+      return accept_iterate;
+   }
+
+   SolutionStatus l1Relaxation::check_termination(const Model& model, Iterate& iterate) {
+      iterate.evaluate_objective_gradient(model);
+      iterate.evaluate_constraints(model);
+
+      this->relaxed_problem->evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, *this->inequality_handling_method,
+         iterate);
+      ConstraintRelaxationStrategy::compute_primal_dual_residuals(*this->relaxed_problem, iterate);
+      return ConstraintRelaxationStrategy::check_termination(*this->relaxed_problem, iterate);
+   }
+
+   std::string l1Relaxation::get_name() const {
+      return this->inequality_handling_method->get_name() + " with " + this->hessian_model->name + " Hessian and " +
+         this->inertia_correction_strategy->get_name() + " regularization";
+   }
+
+   size_t l1Relaxation::get_number_subproblems_solved() const {
+      return this->inequality_handling_method->number_subproblems_solved;
+   }
+} // namespace
\ No newline at end of file
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
new file mode 100644
index 000000000..31dba79f1
--- /dev/null
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
@@ -0,0 +1,49 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_L1RELAXATION_H
+#define UNO_L1RELAXATION_H
+
+#include <memory>
+#include "ConstraintRelaxationStrategy.hpp"
+#include "l1RelaxedProblem.hpp"
+#include "ingredients/hessian_models/HessianModel.hpp"
+#include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
+
+namespace uno {
+   class l1Relaxation : public ConstraintRelaxationStrategy {
+   public:
+      explicit l1Relaxation(const Options& options);
+      ~l1Relaxation() override = default;
+
+      void initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate, Direction& direction,
+         double trust_region_radius, const Options& options) override;
+
+      // direction computation
+      void compute_feasible_direction(Statistics& statistics, GlobalizationStrategy& globalization_strategy, Iterate& current_iterate,
+         Direction& direction, double trust_region_radius, WarmstartInformation& warmstart_information) override;
+      [[nodiscard]] bool solving_feasibility_problem() const override;
+      void switch_to_feasibility_problem(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         Iterate& current_iterate, double trust_region_radius, WarmstartInformation& warmstart_information) override;
+
+      // trial iterate acceptance
+      [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         double trust_region_radius, const Model& model, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, WarmstartInformation& warmstart_information, UserCallbacks& user_callbacks) override;
+      [[nodiscard]] SolutionStatus check_termination(const Model& model, Iterate& iterate) override;
+
+      [[nodiscard]] std::string get_name() const override;
+      [[nodiscard]] size_t get_number_subproblems_solved() const override;
+
+   private:
+      std::unique_ptr<const l1RelaxedProblem> relaxed_problem{};
+      std::unique_ptr<const l1RelaxedProblem> feasibility_problem{};
+      std::unique_ptr<InequalityHandlingMethod> inequality_handling_method;
+      std::unique_ptr<HessianModel> hessian_model{};
+      std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
+      const double constraint_violation_coefficient;
+      double penalty_parameter{1.}; // TODO
+   };
+} // namespace
+
+#endif //UNO_L1RELAXATION_H
\ No newline at end of file

From 7db3677edc93bd63ba9eb491ecb97a1ee26d0f00 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Wed, 5 Nov 2025 12:25:33 +0100
Subject: [PATCH 02/10] Copied most of the 2024 SQuID implementation

---
 .../l1Relaxation.cpp                          | 263 +++++++++++++++++-
 .../l1Relaxation.hpp                          |  32 ++-
 .../relaxed_problems/l1RelaxedProblem.cpp     |  12 +-
 .../relaxed_problems/l1RelaxedProblem.hpp     |   9 +-
 .../InequalityHandlingMethod.hpp              |   1 -
 .../InequalityConstrainedMethod.cpp           |   5 -
 .../InequalityConstrainedMethod.hpp           |   1 -
 .../InteriorPointMethod.hpp                   |   7 -
 uno/optimization/OptimizationProblem.hpp      |   2 +-
 9 files changed, 296 insertions(+), 36 deletions(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index 9eb522f70..289c0ad6d 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -8,25 +8,44 @@
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.hpp"
 #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategyFactory.hpp"
 #include "optimization/Direction.hpp"
+#include "optimization/EvaluationSpace.hpp"
 #include "optimization/Iterate.hpp"
 #include "optimization/OptimizationProblem.hpp"
 #include "optimization/WarmstartInformation.hpp"
 #include "options/Options.hpp"
+#include "symbolic/Sum.hpp"
+#include "symbolic/ScalarMultiple.hpp"
 #include "symbolic/VectorView.hpp"
 #include "tools/Logger.hpp"
+#include "tools/Statistics.hpp"
+
+/*
+ * A Sequential Quadratic Optimization Algorithm with Rapid Infeasibility Detection
+ * James V. Burke, Frank E. Curtis, and Hao Wang
+ * https://epubs.siam.org/doi/abs/10.1137/120880045
+ */
 
 namespace uno {
    l1Relaxation::l1Relaxation(const Options& options) :
          ConstraintRelaxationStrategy(options),
          inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
          inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
-         constraint_violation_coefficient(options.get_double(("l1_constraint_violation_coefficient"))) {
+         penalty_parameter(0.1), // TODO add option
+         tolerance(options.get_double("primal_tolerance")),
+         parameters({
+            1e-2, /* beta */
+            0.1, /* theta */
+            10., /* kappa_rho */
+            10., /* kappa_lambda */
+            1e-2, /* epsilon */
+            1. - 1e-18, /* omega */
+            0.5 /* delta */
+         }) {
    }
 
    void l1Relaxation::initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate,
          Direction& direction, double trust_region_radius, const Options& options) {
-      this->relaxed_problem = std::make_unique<const l1RelaxedProblem>(model, this->penalty_parameter,
-         this->constraint_violation_coefficient);
+      this->relaxed_problem = std::make_unique<l1RelaxedProblem>(model, this->penalty_parameter, 1.);
       assert(this->relaxed_problem != nullptr);
       this->feasibility_problem = std::make_unique<const l1RelaxedProblem>(model, 0., 1.);
       assert(this->feasibility_problem != nullptr);
@@ -42,9 +61,22 @@ namespace uno {
       // statistics
       this->inertia_correction_strategy->initialize_statistics(statistics, options);
       this->inequality_handling_method->initialize_statistics(statistics, options);
+      statistics.add_column("penalty param.", Statistics::double_width, options.get_int("statistics_penalty_parameter_column_order"));
+      statistics.set("penalty param.", this->penalty_parameter);
 
       // initial iterate
+      initial_iterate.set_number_variables(this->feasibility_problem->number_variables);
+      initial_iterate.residuals.lagrangian_gradient.resize(this->feasibility_problem->number_variables);
+      initial_iterate.multipliers.lower_bounds.resize(this->feasibility_problem->number_variables);
+      initial_iterate.multipliers.upper_bounds.resize(this->feasibility_problem->number_variables);
       this->inequality_handling_method->generate_initial_iterate(initial_iterate);
+      this->inequality_handling_method->initialize_feasibility_problem(initial_iterate);
+      this->inequality_handling_method->set_elastic_variable_values(*this->feasibility_problem, initial_iterate);
+      // initialize the feasibility multipliers
+      this->feasibility_multipliers.constraints.resize(this->feasibility_problem->number_constraints);
+      this->feasibility_multipliers.lower_bounds.resize(this->feasibility_problem->number_variables);
+      this->feasibility_multipliers.upper_bounds.resize(this->feasibility_problem->number_variables);
+
       initial_iterate.evaluate_objective_gradient(model);
       initial_iterate.evaluate_constraints(model);
       this->inequality_handling_method->evaluate_constraint_jacobian(initial_iterate);
@@ -55,24 +87,104 @@ namespace uno {
 
    void l1Relaxation::compute_feasible_direction(Statistics& statistics, GlobalizationStrategy& /*globalization_strategy*/,
          Iterate& current_iterate, Direction& direction, double trust_region_radius, WarmstartInformation& warmstart_information) {
+      statistics.set("penalty param.", this->penalty_parameter);
       direction.reset();
-      DEBUG << "Solving the subproblem\n";
-      direction.set_dimensions(this->relaxed_problem->number_variables, this->relaxed_problem->number_constraints);
-      this->inequality_handling_method->solve(statistics, current_iterate, direction, *this->hessian_model,
+
+      // solve feasibility problem
+      DEBUG << "Solving the l1 feasibility problem\n";
+      Direction feasibility_direction(this->feasibility_problem->number_variables, this->feasibility_problem->number_constraints);
+      std::swap(current_iterate.multipliers, this->feasibility_multipliers);
+      this->feasibility_inequality_handling_method->solve(statistics, current_iterate, feasibility_direction, *this->hessian_model,
+         *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
+      assert(direction.status == SubproblemStatus::OPTIMAL && "The feasibility subproblem was not solved to optimality");
+      std::swap(current_iterate.multipliers, this->feasibility_multipliers);
+      // assemble multipliers for feasibility problem
+      Multipliers trial_feasibility_multipliers(current_iterate.number_variables, current_iterate.number_constraints);
+      trial_feasibility_multipliers.constraints = this->feasibility_multipliers.constraints + feasibility_direction.multipliers.constraints;
+      trial_feasibility_multipliers.lower_bounds = this->feasibility_multipliers.lower_bounds + feasibility_direction.multipliers.lower_bounds;
+      trial_feasibility_multipliers.upper_bounds = this->feasibility_multipliers.upper_bounds + feasibility_direction.multipliers.upper_bounds;
+      const double feasibility_dual_error = this->Rinf(current_iterate, trial_feasibility_multipliers);
+      DEBUG << "Infeasibility dual error = " << feasibility_dual_error << '\n';
+
+      // update penalty parameter and duals
+      // test equation (3.12)
+      if (this->tolerance < this->v(current_iterate) &&
+            this->delta_l(feasibility_direction, current_iterate) <= this->parameters.theta * this->v(current_iterate)) {
+         // update penalty parameter according to (3.13)
+         this->penalty_parameter = std::min(this->penalty_parameter, this->parameters.kappa_rho * feasibility_dual_error * feasibility_dual_error);
+         DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+
+         // update current multipliers according to (3.14)
+         const double multipliers_distance = norm_2(
+               current_iterate.multipliers.constraints + (-1) * this->feasibility_multipliers.constraints);
+         const double alpha = std::min(1., this->parameters.kappa_lambda * feasibility_dual_error * feasibility_dual_error / multipliers_distance);
+         current_iterate.multipliers.constraints = alpha * current_iterate.multipliers.constraints +
+            (1. - alpha) * this->feasibility_multipliers.constraints;
+         DEBUG2 << "Updated multipliers: " << current_iterate.multipliers.constraints << '\n';
+      }
+
+      // solve l1 relaxed problem
+      DEBUG << "\nSolving the regular l1 relaxed problem\n";
+      this->relaxed_problem->set_objective_multiplier(this->penalty_parameter);
+      //this->subproblem->set_initial_point(feasibility_direction.primals);
+      warmstart_information.only_objective_changed();
+      Direction optimality_direction(direction.number_variables, direction.number_constraints);
+      this->inequality_handling_method->solve(statistics, current_iterate, optimality_direction, *this->hessian_model,
          *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
-      direction.norm = norm_inf(view(direction.primals, 0, this->relaxed_problem->get_number_original_variables()));
+      assert(direction.status == SubproblemStatus::OPTIMAL && "The l1 relaxed subproblem was not solved to optimality");
+      // assemble multipliers for l1 relaxed problem
+      Multipliers trial_multipliers(current_iterate.number_variables, current_iterate.number_constraints);
+      trial_multipliers.constraints = current_iterate.multipliers.constraints + optimality_direction.multipliers.constraints;
+      trial_multipliers.lower_bounds = current_iterate.multipliers.lower_bounds + optimality_direction.multipliers.lower_bounds;
+      trial_multipliers.upper_bounds = current_iterate.multipliers.upper_bounds + optimality_direction.multipliers.upper_bounds;
+      const double optimality_dual_error = this->Ropt(current_iterate, this->penalty_parameter, trial_multipliers);
+      DEBUG << "Optimality dual error = " << optimality_dual_error << '\n';
+
+      // interpolate between two directions
+      DEBUG << "\nInterpolating between the two directions:\n";
+      DEBUG2 << "Feasibility direction: ";
+      print_vector(DEBUG2, view(feasibility_direction.primals, 0, this->relaxed_problem->model.number_variables));
+      DEBUG2 << "Optimality direction:  ";
+      print_vector(DEBUG2, view(optimality_direction.primals, 0, this->relaxed_problem->model.number_variables));
+      const double w = this->compute_w(feasibility_direction, optimality_direction, current_iterate);
+      const auto interpolated_direction = w * feasibility_direction.primals + (1 - w) * optimality_direction.primals;
+      for (size_t variable_index: Range(direction.number_variables)) {
+         direction.primals[variable_index] = interpolated_direction[variable_index];
+      }
+
+      // update the penalty parameter (3.18)
+      const double multipliers_inf_norm = norm_inf(trial_multipliers.constraints, trial_multipliers.lower_bounds, trial_multipliers.upper_bounds);
+      if (this->penalty_parameter * multipliers_inf_norm > 1.) {
+         this->penalty_parameter = std::min(this->parameters.delta * this->penalty_parameter, (1. - this->parameters.epsilon) / multipliers_inf_norm);
+         DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+      }
+      // update the penalty parameter (3.19)
+      const double delta_m = this->delta_m(direction, current_iterate, this->penalty_parameter);
+      const double delta_l = this->delta_l(direction, current_iterate);
+      if (delta_m >= this->parameters.epsilon * delta_l && w >= this->parameters.omega) {
+         this->penalty_parameter *= this->parameters.delta;
+      }
+      else if (delta_m < this->parameters.epsilon * delta_l) {
+         const double zeta = this->compute_zeta(direction, current_iterate);
+         this->penalty_parameter = std::min(this->parameters.delta * this->penalty_parameter, zeta);
+      }
+      DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+
+      // construct the final direction
+      direction.multipliers = optimality_direction.multipliers;
+      direction.norm = norm_inf(view(direction.primals, 0, this->relaxed_problem->model.number_variables));
       DEBUG3 << direction << '\n';
       warmstart_information.no_changes();
    }
 
    bool l1Relaxation::solving_feasibility_problem() const {
-      return false;
+      return (this->penalty_parameter == 0.);
    }
 
    void l1Relaxation::switch_to_feasibility_problem(Statistics& /*statistics*/, GlobalizationStrategy& /*globalization_strategy*/,
          Iterate& /*current_iterate*/, double /*trust_region_radius*/,
          WarmstartInformation& /*warmstart_information*/) {
-      throw std::runtime_error("The problem is unconstrained, switching to the feasibility problem should not happen");
+      throw std::runtime_error("l1Relaxation::switch_to_feasibility_problem is not implemented\n");
    }
 
    bool l1Relaxation::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
@@ -81,6 +193,9 @@ namespace uno {
       const bool accept_iterate = this->inequality_handling_method->is_iterate_acceptable(statistics, globalization_strategy,
          *this->hessian_model, *this->inertia_correction_strategy, trust_region_radius, current_iterate, trial_iterate,
          direction, step_length, user_callbacks);
+      if (accept_iterate) {
+         this->check_exact_relaxation(trial_iterate);
+      }
       trial_iterate.status = this->check_termination(model, trial_iterate);
       warmstart_information.no_changes();
       return accept_iterate;
@@ -104,4 +219,134 @@ namespace uno {
    size_t l1Relaxation::get_number_subproblems_solved() const {
       return this->inequality_handling_method->number_subproblems_solved;
    }
+
+   // protected member functions
+
+   // v: l1 norm: ||c(x_k)||_1
+   double l1Relaxation::v(const Iterate& current_iterate) const {
+      return current_iterate.progress.infeasibility;
+   }
+
+   // l: constraint violation of linearized constraint: ||c(x_k) + \nabla c(x_k)^T d||_1
+   double l1Relaxation::l(const Direction& direction, const Iterate& current_iterate) const {
+      // TODO preallocate
+      Vector<double> constraint_directional_derivative(this->relaxed_problem->model.number_constraints);
+      const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space();
+      evaluation_space.compute_constraint_jacobian_vector_product(direction.primals, constraint_directional_derivative);
+      const auto linearized_constraint = current_iterate.evaluations.constraints + constraint_directional_derivative;
+      return this->relaxed_problem->model.constraint_violation(linearized_constraint, Norm::L1);
+   }
+
+   // delta_l: predicted (linear model) reduction of constraint violation
+   // ||c(x_k)||_1 - ||c(x_k) + \alpha \nabla c(x_k)^T d||_1
+   double l1Relaxation::delta_l(const Direction& direction, const Iterate& current_iterate) const {
+      return this->v(current_iterate) - this->l(direction, current_iterate);
+   }
+
+   // Ropt: measure that combines stationarity error and complementarity error for l1 relaxed problem
+   double l1Relaxation::Ropt(Iterate& current_iterate, double objective_multiplier, const Multipliers& multipliers) const {
+      /*
+      // stationarity error
+      this->relaxed_problem.evaluate_lagrangian_gradient(current_iterate.residuals.lagrangian_gradient, current_iterate, multipliers);
+      const auto scaled_lagrangian = objective_multiplier * current_iterate.residuals.lagrangian_gradient.objective_contribution +
+         current_iterate.residuals.lagrangian_gradient.constraints_contribution;
+      double error = norm_inf(scaled_lagrangian);
+
+      // complementarity error
+      error += this->relaxed_problem.complementarity_error(current_iterate.primals, current_iterate.evaluations.constraints, multipliers, 0., Norm::INF);
+      return error;
+      */
+      return 0.;
+   }
+
+   // Rinf: measure that combines stationarity error and complementarity error for feasibility problem
+   double l1Relaxation::Rinf(Iterate& current_iterate, const Multipliers& multipliers) const {
+      /*
+      // stationarity error
+      this->feasibility_problem.evaluate_lagrangian_gradient(current_iterate.feasibility_residuals.lagrangian_gradient, current_iterate, multipliers);
+      const auto scaled_lagrangian = current_iterate.feasibility_residuals.lagrangian_gradient.constraints_contribution;
+      double error = norm_inf(scaled_lagrangian);
+
+      // complementarity error
+      error += this->feasibility_problem.complementarity_error(current_iterate.primals, current_iterate.evaluations.constraints, multipliers, 0., Norm::INF);
+      return error;
+      */
+      return 0.;
+   }
+
+   // m: predicted reduction of the l1 merit function
+   double l1Relaxation::delta_m(const Direction& direction, const Iterate& current_iterate, double objective_multiplier) const {
+      return -objective_multiplier * dot(direction.primals, current_iterate.evaluations.objective_gradient) +
+         this->delta_l(direction, current_iterate);
+   }
+
+   double l1Relaxation::compute_w(const Direction& feasibility_direction, const Direction& optimality_direction,
+         const Iterate& current_iterate) {
+      const double delta_l_dbar = this->delta_l(feasibility_direction, current_iterate);
+      // TODO preallocate
+      Vector<double> trial_direction(optimality_direction.primals.size());
+      Vector<double> constraint_directional_derivative(this->relaxed_problem->model.number_constraints);
+      const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space();
+      double weight = 0.;
+      double lower_bound = 0.;
+      double upper_bound = 1.;
+      bool termination = false;
+      while (not termination) {
+         DEBUG2 << "Testing the interpolation weight " << weight << '\n';
+         trial_direction = weight * feasibility_direction.primals + (1 - weight) * optimality_direction.primals;
+         // update reduction in linearized feasibility model
+         evaluation_space.compute_constraint_jacobian_vector_product(trial_direction, constraint_directional_derivative);
+         const auto trial_linearized_constraints = current_iterate.evaluations.constraints + constraint_directional_derivative;
+         double delta_l_d = current_iterate.progress.infeasibility - this->relaxed_problem->model.constraint_violation(trial_linearized_constraints, Norm::L1);
+         DEBUG2 << "Trial predicted infeasibility reduction = " << delta_l_d << '\n';
+
+         // sufficient decrease condition
+         if (delta_l_d >= this->parameters.beta * delta_l_dbar) {
+            // test if bisection has succeeded
+            if (weight == 0. || upper_bound - lower_bound <= 1e-8) {
+               termination = true;
+               DEBUG << "Decrease condition satisfied, terminate with interpolation weight " << weight << '\n';
+               DEBUG2 << "Interpolated direction:  ";
+               print_vector(DEBUG2, trial_direction);
+            }
+            else {
+               // keep the first half
+               upper_bound = weight;
+               weight = (lower_bound + upper_bound) / 2.;
+               DEBUG2 << "Decrease condition satisfied, decreasing the interpolation weight\n";
+            }
+         }
+         else if (weight == 1.) {
+            termination = true;
+            DEBUG << "Terminate with interpolation weight " << weight << '\n';
+            DEBUG2 << "Interpolated direction:  ";
+            print_vector(DEBUG2, trial_direction);
+         }
+         else {
+            // sufficient decrease condition violated: keep the second half
+            lower_bound = weight;
+            weight = (lower_bound + upper_bound) / 2.;
+            DEBUG2 << "Decrease condition violated, increasing the interpolation weight\n";
+         }
+      }
+      return weight;
+   }
+
+   // zeta: upper bound on objective multiplier
+   double l1Relaxation::compute_zeta(const Direction& direction, const Iterate& current_iterate) const {
+      const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space();
+      const double numerator = (1. - this->parameters.epsilon) * this->delta_l(direction, current_iterate);
+      const double denominator = dot(direction.primals, current_iterate.evaluations.objective_gradient);
+         /*+ this->inequality_handling_method->compute_hessian_quadratic_product()
+         evaluation_space.compute_hessian_quadratic_product(subproblem, direction.primals) / 2.;*/
+      return numerator / denominator;
+   }
+
+   // for information purposes, check that l1 is an exact relaxation
+   void l1Relaxation::check_exact_relaxation(const Iterate& iterate) const {
+      const double norm_inf_multipliers = norm_inf(iterate.multipliers.constraints);
+      if (0. < norm_inf_multipliers && this->penalty_parameter <= 1. / norm_inf_multipliers) {
+         DEBUG << "The value of the penalty parameter is consistent with an exact relaxation\n";
+      }
+   }
 } // namespace
\ No newline at end of file
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
index 31dba79f1..725bac455 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
@@ -11,6 +11,16 @@
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
 
 namespace uno {
+   struct l1RelaxationParameters {
+      double beta;
+      double theta;
+      double kappa_rho;
+      double kappa_lambda;
+      double epsilon;
+      double omega;
+      double delta;
+   };
+
    class l1Relaxation : public ConstraintRelaxationStrategy {
    public:
       explicit l1Relaxation(const Options& options);
@@ -35,14 +45,28 @@ namespace uno {
       [[nodiscard]] std::string get_name() const override;
       [[nodiscard]] size_t get_number_subproblems_solved() const override;
 
-   private:
-      std::unique_ptr<const l1RelaxedProblem> relaxed_problem{};
+   protected:
+      std::unique_ptr<l1RelaxedProblem> relaxed_problem{};
       std::unique_ptr<const l1RelaxedProblem> feasibility_problem{};
       std::unique_ptr<InequalityHandlingMethod> inequality_handling_method;
+      std::unique_ptr<InequalityHandlingMethod> feasibility_inequality_handling_method;
       std::unique_ptr<HessianModel> hessian_model{};
       std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
-      const double constraint_violation_coefficient;
-      double penalty_parameter{1.}; // TODO
+      Multipliers feasibility_multipliers;
+      double penalty_parameter;
+      const double tolerance;
+      const l1RelaxationParameters parameters;
+
+      // auxiliary functions with the notations of the paper
+      [[nodiscard]] double v(const Iterate& current_iterate) const;
+      [[nodiscard]] double l(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double delta_l(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double Ropt(Iterate& current_iterate, double objective_multiplier, const Multipliers& multipliers) const;
+      [[nodiscard]] double Rinf(Iterate& current_iterate, const Multipliers& multipliers) const;
+      [[nodiscard]] double delta_m(const Direction& direction, const Iterate& current_iterate, double objective_multiplier) const;
+      [[nodiscard]] double compute_zeta(const Direction& direction, const Iterate& current_iterate) const;
+      [[nodiscard]] double compute_w(const Direction& feasibility_direction, const Direction& optimality_direction, const Iterate& current_iterate);
+      void check_exact_relaxation(const Iterate& iterate) const;
    };
 } // namespace
 
diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
index 6a78e0862..9132eaf69 100644
--- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
@@ -27,16 +27,20 @@ namespace uno {
          constraint_violation_coefficient(constraint_violation_coefficient) {
    }
 
+   void l1RelaxedProblem::set_objective_multiplier(double new_objective_multiplier) {
+      this->objective_multiplier = new_objective_multiplier;
+   }
+
    double l1RelaxedProblem::get_objective_multiplier() const {
       return this->objective_multiplier;
    }
 
-   void l1RelaxedProblem::set_proximal_coefficient(double proximal_coefficient) {
-      this->proximal_coefficient = proximal_coefficient;
+   void l1RelaxedProblem::set_proximal_coefficient(double new_proximal_coefficient) {
+      this->proximal_coefficient = new_proximal_coefficient;
    }
 
-   void l1RelaxedProblem::set_proximal_center(double* proximal_center) {
-      this->proximal_center = proximal_center;
+   void l1RelaxedProblem::set_proximal_center(double* new_proximal_center) {
+      this->proximal_center = new_proximal_center;
    }
 
    void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const {
diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
index 92ef5eeee..aceb57110 100644
--- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
@@ -16,9 +16,10 @@ namespace uno {
          bool relax_linear_constraints);
       ~l1RelaxedProblem() override = default;
 
+      void set_objective_multiplier(double new_objective_multiplier);
       [[nodiscard]] double get_objective_multiplier() const override;
-      void set_proximal_coefficient(double proximal_coefficient);
-      void set_proximal_center(double* proximal_center);
+      void set_proximal_coefficient(double new_proximal_coefficient);
+      void set_proximal_center(double* new_proximal_center);
 
       // constraint evaluations
       void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
@@ -56,7 +57,7 @@ namespace uno {
       [[nodiscard]] const Collection<size_t>& get_dual_regularization_constraints() const override;
 
       [[nodiscard]] SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
-         double dual_tolerance) const;
+         double dual_tolerance) const override;
 
       void set_elastic_variable_values(Iterate& iterate, const std::function<void(Iterate&, size_t, size_t, double)>& elastic_setting_function) const;
 
@@ -68,7 +69,7 @@ namespace uno {
    protected:
       ElasticVariables elastic_variables;
       const size_t number_elastic_variables;
-      const double objective_multiplier;
+      double objective_multiplier;
       const double constraint_violation_coefficient;
       double proximal_coefficient{INF<double>};
       double* proximal_center{};
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
index 3f7eff359..f1f366770 100644
--- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
@@ -49,7 +49,6 @@ namespace uno {
       virtual void evaluate_constraint_jacobian(Iterate& iterate) = 0;
       virtual void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const = 0;
       virtual void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const = 0;
-      [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const = 0;
 
       // progress measures
       [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp
index c281af211..06790f40b 100644
--- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp
+++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp
@@ -100,11 +100,6 @@ namespace uno {
       evaluation_space.compute_constraint_jacobian_transposed_vector_product(vector, result);
    }
 
-   double InequalityConstrainedMethod::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const {
-      const auto& evaluation_space = this->solver->get_evaluation_space();
-      return evaluation_space.compute_hessian_quadratic_product(subproblem, vector);
-   }
-
    // compute dual *displacements*
    // because of the way we form LPs/QPs, we get the new *multipliers* back from the solver. To get the dual displacements/direction,
    // we need to subtract the current multipliers
diff --git a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp
index 95af11afd..5de2d3a36 100644
--- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp
@@ -35,7 +35,6 @@ namespace uno {
       void evaluate_constraint_jacobian(Iterate& iterate) override;
       void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
       void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
-      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const override;
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp
index 63789f481..7dfb7ad74 100644
--- a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp
@@ -41,7 +41,6 @@ namespace uno {
       void evaluate_constraint_jacobian(Iterate& iterate) override;
       void compute_constraint_jacobian_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
       void compute_constraint_jacobian_transposed_vector_product(const Vector<double>& vector, Vector<double>& result) const override;
-      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const override;
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
@@ -259,12 +258,6 @@ namespace uno {
       evaluation_space.compute_constraint_jacobian_transposed_vector_product(vector, result);
    }
 
-   template <typename BarrierProblem>
-   double InteriorPointMethod<BarrierProblem>::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const {
-      const auto& evaluation_space = this->linear_solver->get_evaluation_space();
-      return evaluation_space.compute_hessian_quadratic_product(subproblem, vector);
-   }
-
    template <typename BarrierProblem>
    bool InteriorPointMethod<BarrierProblem>::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
          HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy, double trust_region_radius,
diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp
index d2d4de3dd..91294c5a6 100644
--- a/uno/optimization/OptimizationProblem.hpp
+++ b/uno/optimization/OptimizationProblem.hpp
@@ -79,7 +79,7 @@ namespace uno {
       [[nodiscard]] virtual double complementarity_error(const Vector<double>& primals, const Vector<double>& constraints,
          const Multipliers& multipliers, double shift_value, Norm residual_norm) const;
 
-      [[nodiscard]] SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
+      [[nodiscard]] virtual SolutionStatus check_first_order_convergence(const Iterate& current_iterate, double primal_tolerance,
          double dual_tolerance) const;
 
       // progress measures

From 148c8f6bbb170002fe4113b8f13db5914c7223fa Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Wed, 5 Nov 2025 12:34:57 +0100
Subject: [PATCH 03/10] Replaced not with ! for MSVC

---
 .../constraint_relaxation_strategies/l1Relaxation.cpp           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index 289c0ad6d..6e8dbb277 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -291,7 +291,7 @@ namespace uno {
       double lower_bound = 0.;
       double upper_bound = 1.;
       bool termination = false;
-      while (not termination) {
+      while (!termination) {
          DEBUG2 << "Testing the interpolation weight " << weight << '\n';
          trial_direction = weight * feasibility_direction.primals + (1 - weight) * optimality_direction.primals;
          // update reduction in linearized feasibility model

From 7e0e0f8bfc86c4a028fbaf8510c0d97c204d7e69 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Thu, 6 Nov 2025 10:56:44 +0100
Subject: [PATCH 04/10] Forgot initialization in constructor

---
 .../l1Relaxation.cpp                          | 23 +++++++++++++++----
 .../l1Relaxation.hpp                          |  6 +++--
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index 6e8dbb277..c84be6eda 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -30,6 +30,7 @@ namespace uno {
          ConstraintRelaxationStrategy(options),
          inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
          inertia_correction_strategy(InertiaCorrectionStrategyFactory::create(options)),
+         feasibility_inequality_handling_method(InequalityHandlingMethodFactory::create(options)),
          penalty_parameter(0.1), // TODO add option
          tolerance(options.get_double("primal_tolerance")),
          parameters({
@@ -41,6 +42,9 @@ namespace uno {
             1. - 1e-18, /* omega */
             0.5 /* delta */
          }) {
+      assert(this->inequality_handling_method != nullptr);
+      assert(this->inertia_correction_strategy != nullptr);
+      assert(this->feasibility_inequality_handling_method != nullptr);
    }
 
    void l1Relaxation::initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate,
@@ -50,14 +54,23 @@ namespace uno {
       this->feasibility_problem = std::make_unique<const l1RelaxedProblem>(model, 0., 1.);
       assert(this->feasibility_problem != nullptr);
 
-      // Hessian model
+      // Hessian models
       this->hessian_model = HessianModelFactory::create(model, options);
+      this->feasibility_hessian_model = HessianModelFactory::create(model, options);
+
+      std::cout << "step 1\n";
 
       // memory allocation
       this->inequality_handling_method->initialize(*this->relaxed_problem, initial_iterate, *this->hessian_model,
          *this->inertia_correction_strategy, trust_region_radius);
+      std::cout << "step 2\n";
+      this->feasibility_inequality_handling_method->initialize(*this->feasibility_problem, initial_iterate,
+         *this->feasibility_hessian_model, *this->inertia_correction_strategy, trust_region_radius);
+      std::cout << "step 3\n";
       direction = Direction(this->relaxed_problem->number_variables, this->relaxed_problem->number_constraints);
 
+      std::cout << "step 4\n";
+
       // statistics
       this->inertia_correction_strategy->initialize_statistics(statistics, options);
       this->inequality_handling_method->initialize_statistics(statistics, options);
@@ -80,8 +93,8 @@ namespace uno {
       initial_iterate.evaluate_objective_gradient(model);
       initial_iterate.evaluate_constraints(model);
       this->inequality_handling_method->evaluate_constraint_jacobian(initial_iterate);
-      this->relaxed_problem->evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, *this->inequality_handling_method,
-         initial_iterate);
+      this->relaxed_problem->evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient,
+         *this->inequality_handling_method, initial_iterate);
       this->compute_primal_dual_residuals(*this->relaxed_problem, initial_iterate);
    }
 
@@ -94,8 +107,8 @@ namespace uno {
       DEBUG << "Solving the l1 feasibility problem\n";
       Direction feasibility_direction(this->feasibility_problem->number_variables, this->feasibility_problem->number_constraints);
       std::swap(current_iterate.multipliers, this->feasibility_multipliers);
-      this->feasibility_inequality_handling_method->solve(statistics, current_iterate, feasibility_direction, *this->hessian_model,
-         *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
+      this->feasibility_inequality_handling_method->solve(statistics, current_iterate, feasibility_direction,
+         *this->feasibility_hessian_model, *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
       assert(direction.status == SubproblemStatus::OPTIMAL && "The feasibility subproblem was not solved to optimality");
       std::swap(current_iterate.multipliers, this->feasibility_multipliers);
       // assemble multipliers for feasibility problem
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
index 725bac455..2c1e0cf95 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
@@ -47,11 +47,13 @@ namespace uno {
 
    protected:
       std::unique_ptr<l1RelaxedProblem> relaxed_problem{};
-      std::unique_ptr<const l1RelaxedProblem> feasibility_problem{};
       std::unique_ptr<InequalityHandlingMethod> inequality_handling_method;
-      std::unique_ptr<InequalityHandlingMethod> feasibility_inequality_handling_method;
       std::unique_ptr<HessianModel> hessian_model{};
       std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
+      // feasibility problem
+      std::unique_ptr<const l1RelaxedProblem> feasibility_problem{};
+      std::unique_ptr<InequalityHandlingMethod> feasibility_inequality_handling_method;
+      std::unique_ptr<HessianModel> feasibility_hessian_model{};
       Multipliers feasibility_multipliers;
       double penalty_parameter;
       const double tolerance;

From 3527c30f0aa5d62a2a50e039d9bf9d0c9d609828 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Thu, 6 Nov 2025 11:22:42 +0100
Subject: [PATCH 05/10] Properly print the Hessian in
 PrimalInertiaCorrection.hpp

---
 .../inertia_correction_strategies/PrimalInertiaCorrection.hpp   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp b/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
index 8e38bab8c..72816dc68 100644
--- a/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
+++ b/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
@@ -87,7 +87,7 @@ namespace uno {
          DirectSymmetricIndefiniteLinearSolver<double>& linear_solver, double* primal_regularization_values) {
       assert(hessian_values != nullptr);
       const size_t number_hessian_nonzeros = subproblem.number_regularized_hessian_nonzeros();
-
+      
       this->regularization_factor = 0.;
       bool good_inertia = false;
       while (!good_inertia) {

From 6aa62974ae1821bf400075cf37cb7d6d8ae41c80 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Thu, 6 Nov 2025 11:22:59 +0100
Subject: [PATCH 06/10] Bug fix: do not modify WarmstartInformation object

---
 .../l1Relaxation.cpp                              | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index c84be6eda..15a23e8f9 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -58,19 +58,13 @@ namespace uno {
       this->hessian_model = HessianModelFactory::create(model, options);
       this->feasibility_hessian_model = HessianModelFactory::create(model, options);
 
-      std::cout << "step 1\n";
-
       // memory allocation
       this->inequality_handling_method->initialize(*this->relaxed_problem, initial_iterate, *this->hessian_model,
          *this->inertia_correction_strategy, trust_region_radius);
-      std::cout << "step 2\n";
       this->feasibility_inequality_handling_method->initialize(*this->feasibility_problem, initial_iterate,
          *this->feasibility_hessian_model, *this->inertia_correction_strategy, trust_region_radius);
-      std::cout << "step 3\n";
       direction = Direction(this->relaxed_problem->number_variables, this->relaxed_problem->number_constraints);
 
-      std::cout << "step 4\n";
-
       // statistics
       this->inertia_correction_strategy->initialize_statistics(statistics, options);
       this->inequality_handling_method->initialize_statistics(statistics, options);
@@ -128,8 +122,7 @@ namespace uno {
          DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
 
          // update current multipliers according to (3.14)
-         const double multipliers_distance = norm_2(
-               current_iterate.multipliers.constraints + (-1) * this->feasibility_multipliers.constraints);
+         const double multipliers_distance = norm_2(current_iterate.multipliers.constraints + (-1) * this->feasibility_multipliers.constraints);
          const double alpha = std::min(1., this->parameters.kappa_lambda * feasibility_dual_error * feasibility_dual_error / multipliers_distance);
          current_iterate.multipliers.constraints = alpha * current_iterate.multipliers.constraints +
             (1. - alpha) * this->feasibility_multipliers.constraints;
@@ -140,7 +133,7 @@ namespace uno {
       DEBUG << "\nSolving the regular l1 relaxed problem\n";
       this->relaxed_problem->set_objective_multiplier(this->penalty_parameter);
       //this->subproblem->set_initial_point(feasibility_direction.primals);
-      warmstart_information.only_objective_changed();
+      // warmstart_information.only_objective_changed();
       Direction optimality_direction(direction.number_variables, direction.number_constraints);
       this->inequality_handling_method->solve(statistics, current_iterate, optimality_direction, *this->hessian_model,
          *this->inertia_correction_strategy, trust_region_radius, warmstart_information);
@@ -176,12 +169,14 @@ namespace uno {
       const double delta_l = this->delta_l(direction, current_iterate);
       if (delta_m >= this->parameters.epsilon * delta_l && w >= this->parameters.omega) {
          this->penalty_parameter *= this->parameters.delta;
+         DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
       }
       else if (delta_m < this->parameters.epsilon * delta_l) {
          const double zeta = this->compute_zeta(direction, current_iterate);
          this->penalty_parameter = std::min(this->parameters.delta * this->penalty_parameter, zeta);
+         DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
       }
-      DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+      assert(this->penalty_parameter >= 0);
 
       // construct the final direction
       direction.multipliers = optimality_direction.multipliers;

From 9bf1813e6c516793c01a628dd817cdff9d187e27 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Thu, 6 Nov 2025 16:17:07 +0100
Subject: [PATCH 07/10] Implemented Rinf + added a few assertions

---
 .../l1Relaxation.cpp                          | 23 ++++++++++---------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index 15a23e8f9..506ba9344 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -163,6 +163,7 @@ namespace uno {
       if (this->penalty_parameter * multipliers_inf_norm > 1.) {
          this->penalty_parameter = std::min(this->parameters.delta * this->penalty_parameter, (1. - this->parameters.epsilon) / multipliers_inf_norm);
          DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+         assert(this->penalty_parameter >= 0);
       }
       // update the penalty parameter (3.19)
       const double delta_m = this->delta_m(direction, current_iterate, this->penalty_parameter);
@@ -170,13 +171,16 @@ namespace uno {
       if (delta_m >= this->parameters.epsilon * delta_l && w >= this->parameters.omega) {
          this->penalty_parameter *= this->parameters.delta;
          DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+         assert(this->penalty_parameter >= 0);
       }
       else if (delta_m < this->parameters.epsilon * delta_l) {
-         const double zeta = this->compute_zeta(direction, current_iterate);
+         double zeta = this->compute_zeta(direction, current_iterate);
+         // TODO remove this line:
+         zeta = std::max(zeta, this->penalty_parameter);
          this->penalty_parameter = std::min(this->parameters.delta * this->penalty_parameter, zeta);
          DEBUG << "Penalty parameter updated to " << this->penalty_parameter << '\n';
+         assert(this->penalty_parameter >= 0);
       }
-      assert(this->penalty_parameter >= 0);
 
       // construct the final direction
       direction.multipliers = optimality_direction.multipliers;
@@ -190,8 +194,7 @@ namespace uno {
    }
 
    void l1Relaxation::switch_to_feasibility_problem(Statistics& /*statistics*/, GlobalizationStrategy& /*globalization_strategy*/,
-         Iterate& /*current_iterate*/, double /*trust_region_radius*/,
-         WarmstartInformation& /*warmstart_information*/) {
+         Iterate& /*current_iterate*/, double /*trust_region_radius*/, WarmstartInformation& /*warmstart_information*/) {
       throw std::runtime_error("l1Relaxation::switch_to_feasibility_problem is not implemented\n");
    }
 
@@ -269,17 +272,15 @@ namespace uno {
 
    // Rinf: measure that combines stationarity error and complementarity error for feasibility problem
    double l1Relaxation::Rinf(Iterate& current_iterate, const Multipliers& multipliers) const {
-      /*
       // stationarity error
-      this->feasibility_problem.evaluate_lagrangian_gradient(current_iterate.feasibility_residuals.lagrangian_gradient, current_iterate, multipliers);
-      const auto scaled_lagrangian = current_iterate.feasibility_residuals.lagrangian_gradient.constraints_contribution;
-      double error = norm_inf(scaled_lagrangian);
+      this->relaxed_problem->evaluate_lagrangian_gradient(current_iterate.residuals.lagrangian_gradient,
+         *this->inequality_handling_method, current_iterate);
+      double error = OptimizationProblem::stationarity_error(current_iterate.residuals.lagrangian_gradient, 0., this->residual_norm);
 
       // complementarity error
-      error += this->feasibility_problem.complementarity_error(current_iterate.primals, current_iterate.evaluations.constraints, multipliers, 0., Norm::INF);
+      error += this->feasibility_problem->complementarity_error(current_iterate.primals, current_iterate.evaluations.constraints,
+         multipliers, 0., this->residual_norm);
       return error;
-      */
-      return 0.;
    }
 
    // m: predicted reduction of the l1 merit function

From bd9970ef7c06b78b7be7a00eca78bef8a88696ce Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Fri, 7 Nov 2025 10:59:02 +0100
Subject: [PATCH 08/10] Added squid preset

---
 uno/options/Presets.cpp | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/uno/options/Presets.cpp b/uno/options/Presets.cpp
index 72316741e..b648b3f4c 100644
--- a/uno/options/Presets.cpp
+++ b/uno/options/Presets.cpp
@@ -133,6 +133,19 @@ namespace uno {
          preset_options.set_bool("switch_to_optimality_requires_linearized_feasibility", true);
          preset_options.set_bool("protect_actual_reduction_against_roundoff", false);
       }
+      else if (preset_name == "squid") {
+         preset_options.set_string("constraint_relaxation_strategy", "l1_relaxation");
+         preset_options.set_string("inequality_handling_method", "inequality_constrained");
+         preset_options.set_string("hessian_model", "exact");
+         preset_options.set_string("inertia_correction_strategy", "primal");
+         preset_options.set_string("globalization_mechanism", "LS");
+         preset_options.set_string("globalization_strategy", "l1_merit");
+         preset_options.set_string("progress_norm", "L1");
+         preset_options.set_string("residual_norm", "INF");
+         preset_options.set_double("l1_constraint_violation_coefficient", 1.);
+         preset_options.set_double("primal_tolerance", 1e-6);
+         preset_options.set_double("dual_tolerance", 1e-6);
+      }
       else {
          throw std::runtime_error("The preset " + preset_name + " is not known");
       }

From acb98b947e18f36f051f477ed07effcbf9c2d424 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Fri, 7 Nov 2025 11:00:05 +0100
Subject: [PATCH 09/10] Start primal inertia correction with 0 coefficient

---
 .../constraint_relaxation_strategies/l1Relaxation.cpp           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index 506ba9344..d7f3dc47b 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -249,7 +249,7 @@ namespace uno {
    }
 
    // delta_l: predicted (linear model) reduction of constraint violation
-   // ||c(x_k)||_1 - ||c(x_k) + \alpha \nabla c(x_k)^T d||_1
+   // ||c(x_k)||_1 - ||c(x_k) + \nabla c(x_k)^T d||_1
    double l1Relaxation::delta_l(const Direction& direction, const Iterate& current_iterate) const {
       return this->v(current_iterate) - this->l(direction, current_iterate);
    }

From 25aa3a8087a1693f2518c6e2cd971b6ce44abcfc Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Fri, 7 Nov 2025 18:10:13 +0100
Subject: [PATCH 10/10] Fixes after rebasing

---
 .../constraint_relaxation_strategies/l1Relaxation.cpp       | 6 ++++--
 .../constraint_relaxation_strategies/l1Relaxation.hpp       | 2 +-
 .../PrimalInertiaCorrection.hpp                             | 2 +-
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
index d7f3dc47b..b74dfa040 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp
@@ -49,9 +49,11 @@ namespace uno {
 
    void l1Relaxation::initialize(Statistics& statistics, const Model& model, Iterate& initial_iterate,
          Direction& direction, double trust_region_radius, const Options& options) {
-      this->relaxed_problem = std::make_unique<l1RelaxedProblem>(model, this->penalty_parameter, 1.);
+      // relax the linear constraints in the l1 relaxed problem only if we are using a trust-region constraint
+      const bool relax_linear_constraints = (trust_region_radius < INF<double>);
+      this->relaxed_problem = std::make_unique<l1RelaxedProblem>(model, this->penalty_parameter, 1., relax_linear_constraints);
       assert(this->relaxed_problem != nullptr);
-      this->feasibility_problem = std::make_unique<const l1RelaxedProblem>(model, 0., 1.);
+      this->feasibility_problem = std::make_unique<const l1RelaxedProblem>(model, 0., 1., relax_linear_constraints);
       assert(this->feasibility_problem != nullptr);
 
       // Hessian models
diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
index 2c1e0cf95..b15c574e9 100644
--- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp
@@ -6,7 +6,7 @@
 
 #include <memory>
 #include "ConstraintRelaxationStrategy.hpp"
-#include "l1RelaxedProblem.hpp"
+#include "relaxed_problems/l1RelaxedProblem.hpp"
 #include "ingredients/hessian_models/HessianModel.hpp"
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
 
diff --git a/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp b/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
index 72816dc68..8e38bab8c 100644
--- a/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
+++ b/uno/ingredients/inertia_correction_strategies/PrimalInertiaCorrection.hpp
@@ -87,7 +87,7 @@ namespace uno {
          DirectSymmetricIndefiniteLinearSolver<double>& linear_solver, double* primal_regularization_values) {
       assert(hessian_values != nullptr);
       const size_t number_hessian_nonzeros = subproblem.number_regularized_hessian_nonzeros();
-      
+
       this->regularization_factor = 0.;
       bool good_inertia = false;
       while (!good_inertia) {