diff --git a/interfaces/AMPL/AMPLModel.cpp b/interfaces/AMPL/AMPLModel.cpp
index 7601a5cb6..ad77d9b3c 100644
--- a/interfaces/AMPL/AMPLModel.cpp
+++ b/interfaces/AMPL/AMPLModel.cpp
@@ -91,12 +91,13 @@ namespace uno {
 
    double AMPLModel::evaluate_objective(const Vector<double>& x) const {
       fint error_flag = 0;
-      const double result = this->optimization_sense * (*(this->asl)->p.Objval)(this->asl, 0, const_cast<double*>(x.data()), &error_flag);
+      double objective_value = (*(this->asl)->p.Objval)(this->asl, 0, const_cast<double*>(x.data()), &error_flag);
       if (0 < error_flag) {
          throw FunctionEvaluationError();
       }
+      objective_value *= this->optimization_sense;
       ++this->number_model_evaluations.objective;
-      return result;
+      return objective_value;
    }
 
    /*
diff --git a/uno/Uno.cpp b/uno/Uno.cpp
index d61b10660..f9363fb4f 100644
--- a/uno/Uno.cpp
+++ b/uno/Uno.cpp
@@ -2,7 +2,6 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "Uno.hpp"
-#include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp"
 #include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.hpp"
 #include "ingredients/globalization_mechanisms/GlobalizationMechanismFactory.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp"
diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp
index cf07320ad..3b2ed8d45 100644
--- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp
@@ -69,7 +69,14 @@ namespace uno {
       const auto& evaluation_space = this->optimality_inequality_handling_method->get_evaluation_space();
       this->optimality_inequality_handling_method->evaluate_constraint_jacobian(initial_iterate);
       this->optimality_problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient,
-         evaluation_space, initial_iterate);
+         evaluation_space, initial_iterate, this->scaling);
+
+      // optional scaling
+      if (options.get_bool("use_function_scaling")) {
+         this->scaling.emplace(initial_iterate, this->optimality_inequality_handling_method->get_evaluation_space(),
+            options.get_double("function_scaling_threshold"));
+      }
+      this->optimality_inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling);
       ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, initial_iterate);
       this->optimality_globalization_strategy->initialize(statistics, initial_iterate, options);
       this->feasibility_globalization_strategy.initialize(statistics, initial_iterate, options);
@@ -157,7 +164,8 @@ namespace uno {
          const OptimizationProblem& problem, Iterate& current_iterate, Direction& direction, double trust_region_radius,
          WarmstartInformation& warmstart_information) {
       direction.set_dimensions(problem.number_variables, problem.number_constraints);
-      inequality_handling_method.solve(statistics, current_iterate, direction, trust_region_radius, warmstart_information);
+      inequality_handling_method.solve(statistics, current_iterate, direction, trust_region_radius, this->scaling,
+         warmstart_information);
       direction.norm = norm_inf(view(direction.primals, 0, problem.get_number_original_variables()));
       DEBUG3 << direction << '\n';
    }
@@ -200,11 +208,13 @@ namespace uno {
       // determine acceptability, depending on the current phase
       if (this->current_phase == Phase::OPTIMALITY) {
          accept_iterate = this->optimality_inequality_handling_method->is_iterate_acceptable(statistics,
-            *this->optimality_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks);
+            *this->optimality_globalization_strategy, this->scaling, current_iterate, trial_iterate, direction,
+            step_length, user_callbacks);
       }
       else {
          accept_iterate = this->feasibility_inequality_handling_method->is_iterate_acceptable(statistics,
-            this->feasibility_globalization_strategy, current_iterate, trial_iterate, direction, step_length, user_callbacks);
+            this->feasibility_globalization_strategy, this->scaling, current_iterate, trial_iterate, direction,
+            step_length, user_callbacks);
       }
       trial_iterate.status = this->check_termination(model, trial_iterate);
 
@@ -227,13 +237,13 @@ namespace uno {
 
       if (this->current_phase == Phase::OPTIMALITY) {
          this->optimality_problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient,
-            this->optimality_inequality_handling_method->get_evaluation_space(), iterate);
+            this->optimality_inequality_handling_method->get_evaluation_space(), iterate, this->scaling);
          ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->optimality_problem, iterate);
          return ConstraintRelaxationStrategy::check_termination(this->optimality_problem, iterate);
       }
       else {
          this->feasibility_problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient,
-            this->feasibility_inequality_handling_method->get_evaluation_space(), iterate);
+            this->feasibility_inequality_handling_method->get_evaluation_space(), iterate, this->scaling);
          ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->feasibility_problem, iterate);
          return ConstraintRelaxationStrategy::check_termination(this->feasibility_problem, iterate);
       }
diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp
index fa420dca3..59171083f 100644
--- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp
@@ -5,6 +5,7 @@
 #define UNO_FEASIBILITYRESTORATION_H
 
 #include <memory>
+#include <optional>
 #include "ConstraintRelaxationStrategy.hpp"
 #include "relaxed_problems/l1RelaxedProblem.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
@@ -12,6 +13,7 @@
 #include "ingredients/globalization_strategies/ProgressMeasures.hpp"
 #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategy.hpp"
 #include "linear_algebra/Vector.hpp"
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    // forward declaration
@@ -56,6 +58,7 @@ namespace uno {
       std::unique_ptr<InequalityHandlingMethod> feasibility_inequality_handling_method;
       std::unique_ptr<GlobalizationStrategy> optimality_globalization_strategy;
       l1MeritFunction feasibility_globalization_strategy;
+      std::optional<Scaling> scaling{};
       // the class maintains multipliers for the other phase (feasibility multipliers if we are in the optimality phase,
       // and vice versa). These multipliers and those of the iterate are swapped whenever we switch phases.
       Multipliers other_phase_multipliers;
diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp
index 85adcc7ae..33afde21a 100644
--- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp
@@ -9,6 +9,7 @@
 #include "optimization/Iterate.hpp"
 #include "optimization/OptimizationProblem.hpp"
 #include "optimization/WarmstartInformation.hpp"
+#include "options/Options.hpp"
 #include "symbolic/VectorView.hpp"
 #include "tools/Logger.hpp"
 
@@ -36,10 +37,19 @@ namespace uno {
       // 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);
       const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space();
-      this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space, initial_iterate);
+      this->problem.evaluate_lagrangian_gradient(initial_iterate.residuals.lagrangian_gradient, evaluation_space,
+         initial_iterate, this->scaling);
+      this->compute_primal_dual_residuals(this->problem, initial_iterate);
+      this->globalization_strategy.initialize(statistics, initial_iterate, options);
+
+      // optional scaling
+      if (options.get_bool("use_function_scaling")) {
+         this->scaling.emplace(initial_iterate, this->inequality_handling_method->get_evaluation_space(),
+            options.get_double("function_scaling_threshold"));
+      }
+      this->inequality_handling_method->evaluate_progress_measures(initial_iterate, this->scaling);
       this->compute_primal_dual_residuals(this->problem, initial_iterate);
       this->globalization_strategy.initialize(statistics, initial_iterate, options);
    }
@@ -49,7 +59,8 @@ namespace uno {
       direction.reset();
       DEBUG << "Solving the subproblem\n";
       direction.set_dimensions(this->problem.number_variables, this->problem.number_constraints);
-      this->inequality_handling_method->solve(statistics, current_iterate, direction, trust_region_radius, warmstart_information);
+      this->inequality_handling_method->solve(statistics, current_iterate, direction, trust_region_radius, this->scaling,
+         warmstart_information);
       direction.norm = norm_inf(view(direction.primals, 0, this->problem.get_number_original_variables()));
       DEBUG3 << direction << '\n';
       warmstart_information.no_changes();
@@ -68,7 +79,7 @@ namespace uno {
          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, this->globalization_strategy,
-         current_iterate, trial_iterate, direction, step_length, user_callbacks);
+         this->scaling, 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;
@@ -79,7 +90,7 @@ namespace uno {
       iterate.evaluate_constraints(model);
 
       const auto& evaluation_space = this->inequality_handling_method->get_evaluation_space();
-      this->problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, evaluation_space, iterate);
+      this->problem.evaluate_lagrangian_gradient(iterate.residuals.lagrangian_gradient, evaluation_space, iterate, this->scaling);
       ConstraintRelaxationStrategy::compute_primal_dual_residuals(this->problem, iterate);
       return ConstraintRelaxationStrategy::check_termination(this->problem, iterate);
    }
diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp
index a37f18b53..6edfc9db5 100644
--- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.hpp
@@ -5,11 +5,13 @@
 #define UNO_NORELAXATION_H
 
 #include <memory>
+#include <optional>
 #include "ConstraintRelaxationStrategy.hpp"
 #include "ingredients/globalization_strategies/l1MeritFunction.hpp"
 #include "ingredients/hessian_models/HessianModel.hpp"
 #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
 #include "optimization/OptimizationProblem.hpp"
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    class NoRelaxation : public ConstraintRelaxationStrategy {
@@ -42,6 +44,7 @@ namespace uno {
       std::unique_ptr<HessianModel> hessian_model;
       std::unique_ptr<InertiaCorrectionStrategy<double>> inertia_correction_strategy;
       l1MeritFunction globalization_strategy;
+      std::optional<Scaling> scaling{};
    };
 } // namespace
 
diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
index af60149a6..c8a11b70f 100644
--- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
+++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp
@@ -53,13 +53,14 @@ namespace uno {
       }
    }
 
-   void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const {
+   void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional<Scaling>& scaling) const {
       // scale nabla f(x) by rho
       if (this->objective_multiplier != 0.) {
          iterate.evaluate_objective_gradient(this->model);
          // TODO change this
+         const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.;
          for (size_t index: Range(this->model.number_variables)) {
-            objective_gradient[index] = this->objective_multiplier * iterate.evaluations.objective_gradient[index];
+            objective_gradient[index] = this->objective_multiplier * objective_scaling * iterate.evaluations.objective_gradient[index];
          }
       }
       else {
@@ -151,7 +152,7 @@ namespace uno {
 
    // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
    void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const {
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& /*scaling*/) const {
       lagrangian_gradient.objective_contribution.fill(0.);
       lagrangian_gradient.constraints_contribution.fill(0.);
 
@@ -191,9 +192,9 @@ namespace uno {
    }
 
    void l1RelaxedProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const {
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& scaling) const {
       hessian_model.evaluate_hessian(statistics, primal_variables, this->get_objective_multiplier(),
-         multipliers.constraints, hessian_values);
+         multipliers.constraints, hessian_values, scaling);
 
       // proximal contribution
       size_t nonzero_index = hessian_model.number_nonzeros();
@@ -208,8 +209,9 @@ namespace uno {
    }
 
    void l1RelaxedProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector,
-         const Multipliers& multipliers, double* result) const {
-      hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result);
+         const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const {
+      hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints,
+         result, scaling);
 
       // proximal contribution
       if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) {
diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
index b760826fb..4e79df31e 100644
--- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
+++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp
@@ -24,7 +24,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
 
       // dense objective gradient
-      void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const override;
+      void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional<Scaling>& scaling) const override;
 
       // sparsity patterns of Jacobian and Hessian
       [[nodiscard]] size_t number_jacobian_nonzeros() const override;
@@ -38,11 +38,11 @@ namespace uno {
       // numerical evaluations of Jacobian and Hessian
       void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const override;
       void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const override;
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const override;
       void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const override;
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& scaling) const override;
       void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector,
-         const Multipliers& multipliers, double* result) const override;
+         const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp
index 68f1646b2..84f0b15d8 100644
--- a/uno/ingredients/hessian_models/ExactHessian.cpp
+++ b/uno/ingredients/hessian_models/ExactHessian.cpp
@@ -34,13 +34,20 @@ namespace uno {
    }
 
    void ExactHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& primal_variables,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) {
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) {
+      if (scaling.has_value()) {
+         objective_multiplier *= scaling->get_objective_scaling();
+      }
       this->model.evaluate_lagrangian_hessian(primal_variables, objective_multiplier, constraint_multipliers, hessian_values);
       ++this->evaluation_count;
    }
 
-   void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result) {
+   void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& scaling) {
+      if (scaling.has_value()) {
+         objective_multiplier *= scaling->get_objective_scaling();
+      }
       this->model.compute_hessian_vector_product(x, vector, objective_multiplier, constraint_multipliers, result);
       ++this->evaluation_count;
    }
diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp
index 9022ce0b9..9128668c3 100644
--- a/uno/ingredients/hessian_models/ExactHessian.hpp
+++ b/uno/ingredients/hessian_models/ExactHessian.hpp
@@ -22,9 +22,9 @@ namespace uno {
       [[nodiscard]] bool is_positive_definite() const override;
 
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* hessian_values) override;
+         const Vector<double>& constraint_multipliers, double* hessian_values, const std::optional<Scaling>& scaling) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* result) override;
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& scaling) override;
 
    protected:
       const Model& model;
diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp
index 202a2d6b7..a1f33e0cb 100644
--- a/uno/ingredients/hessian_models/HessianModel.hpp
+++ b/uno/ingredients/hessian_models/HessianModel.hpp
@@ -5,9 +5,11 @@
 #define UNO_HESSIANMODEL_H
 
 #include <cstddef>
+#include <optional>
 #include <string>
 #include <string_view>
 #include "../interfaces/C/uno_int.h"
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    // forward declarations
@@ -32,9 +34,11 @@ namespace uno {
       [[nodiscard]] virtual bool is_positive_definite() const = 0;
 
       virtual void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) = 0;
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) = 0;
       virtual void compute_hessian_vector_product(const double* x, const double* vector,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result) = 0;
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* result,
+         const std::optional<Scaling>& scaling) = 0;
    };
 } // namespace
 
diff --git a/uno/ingredients/hessian_models/IdentityHessian.cpp b/uno/ingredients/hessian_models/IdentityHessian.cpp
index b7d66228d..22af946a7 100644
--- a/uno/ingredients/hessian_models/IdentityHessian.cpp
+++ b/uno/ingredients/hessian_models/IdentityHessian.cpp
@@ -38,7 +38,8 @@ namespace uno {
    }
 
    void IdentityHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& /*primal_variables*/,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* hessian_values) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* hessian_values,
+         const std::optional<Scaling>& /*scaling*/) {
       DEBUG << "Setting identity Hessian\n";
       for (size_t variable_index: Range(this->number_variables)) {
          hessian_values[variable_index] = 1.;
@@ -46,7 +47,8 @@ namespace uno {
    }
 
    void IdentityHessian::compute_hessian_vector_product(const double* /*x*/, const double* vector,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result,
+         const std::optional<Scaling>& /*scaling*/) {
       for (size_t variable_index: Range(this->number_variables)) {
          result[variable_index] = vector[variable_index];
       }
diff --git a/uno/ingredients/hessian_models/IdentityHessian.hpp b/uno/ingredients/hessian_models/IdentityHessian.hpp
index 5fde92d9b..7332e28e5 100644
--- a/uno/ingredients/hessian_models/IdentityHessian.hpp
+++ b/uno/ingredients/hessian_models/IdentityHessian.hpp
@@ -19,9 +19,10 @@ namespace uno {
       [[nodiscard]] bool is_positive_definite() const override;
 
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables,
-         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values) override;
+         double objective_multiplier, const Vector<double>& constraint_multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* result) override;
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& scaling) override;
 
    protected:
       const size_t number_variables;
diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp
index ecccc6744..299aa0471 100644
--- a/uno/ingredients/hessian_models/ZeroHessian.cpp
+++ b/uno/ingredients/hessian_models/ZeroHessian.cpp
@@ -33,12 +33,14 @@ namespace uno {
    }
 
    void ZeroHessian::evaluate_hessian(Statistics& /*statistics*/, const Vector<double>& /*primal_variables*/,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* /*hessian_values*/) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* /*hessian_values*/,
+         const std::optional<Scaling>& /*scaling*/) {
       // do nothing
    }
 
    void ZeroHessian::compute_hessian_vector_product(const double* /*x*/, const double* /*vector*/,
-         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result) {
+         double /*objective_multiplier*/, const Vector<double>& /*constraint_multipliers*/, double* result,
+         const std::optional<Scaling>& /*scaling*/) {
       for (size_t variable_index: Range(this->number_variables)) {
          result[variable_index] = 0.;
       }
diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp
index e8c11858d..7ce996973 100644
--- a/uno/ingredients/hessian_models/ZeroHessian.hpp
+++ b/uno/ingredients/hessian_models/ZeroHessian.hpp
@@ -19,9 +19,9 @@ namespace uno {
       [[nodiscard]] bool is_positive_definite() const override;
 
       void evaluate_hessian(Statistics& statistics, const Vector<double>& primal_variables, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* hessian_values) override;
+         const Vector<double>& constraint_multipliers, double* hessian_values, const std::optional<Scaling>& scaling) override;
       void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier,
-         const Vector<double>& constraint_multipliers, double* result) override;
+         const Vector<double>& constraint_multipliers, double* result, const std::optional<Scaling>& scaling) override;
 
    protected:
       const size_t number_variables;
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp
index f2cca245d..dbf30aef0 100644
--- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp
+++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.cpp
@@ -17,15 +17,17 @@ namespace uno {
       progress_norm(norm_from_string(options.get_string("progress_norm"))) {
    }
 
-   void InequalityHandlingMethod::evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate) const {
+   void InequalityHandlingMethod::evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate,
+         const std::optional<Scaling>& scaling) const {
       problem.set_infeasibility_measure(iterate, this->progress_norm);
-      problem.set_objective_measure(iterate);
+      problem.set_objective_measure(iterate, scaling);
       problem.set_auxiliary_measure(iterate);
    }
 
    bool InequalityHandlingMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         const Subproblem& subproblem, const EvaluationSpace& evaluation_space, Iterate& current_iterate, Iterate& trial_iterate,
-         const Direction& direction, double step_length, UserCallbacks& user_callbacks) {
+         const Subproblem& subproblem, const std::optional<Scaling>& scaling, const EvaluationSpace& evaluation_space,
+         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
+         UserCallbacks& user_callbacks) {
       this->postprocess_iterate(trial_iterate);
       const double objective_multiplier = subproblem.problem.get_objective_multiplier();
 
@@ -37,7 +39,7 @@ namespace uno {
          subproblem.problem.set_auxiliary_measure(current_iterate);
          this->subproblem_definition_changed = false;
       }
-      this->evaluate_progress_measures(subproblem.problem, trial_iterate);
+      this->evaluate_progress_measures(subproblem.problem, trial_iterate, scaling);
 
       bool accept_iterate = false;
       if (direction.norm == 0.) {
@@ -48,7 +50,7 @@ namespace uno {
       }
       else {
          const ProgressMeasures predicted_reductions = subproblem.compute_predicted_reductions(evaluation_space, direction,
-            step_length, this->progress_norm);
+            step_length, this->progress_norm, scaling);
          accept_iterate = globalization_strategy.is_iterate_acceptable(statistics, current_iterate.progress, trial_iterate.progress,
             predicted_reductions, objective_multiplier);
       }
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
index ac66a1909..4b7d739bd 100644
--- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp
@@ -4,8 +4,10 @@
 #ifndef UNO_INEQUALITYHANDLINGMETHOD_H
 #define UNO_INEQUALITYHANDLINGMETHOD_H
 
+#include <optional>
 #include <string>
 #include "linear_algebra/Norm.hpp"
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    // forward declarations
@@ -33,11 +35,13 @@ namespace uno {
       virtual ~InequalityHandlingMethod() = default;
 
       virtual void initialize(const OptimizationProblem& problem, Iterate& current_iterate,
-         HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy, double trust_region_radius) = 0;
+         HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy,
+         double trust_region_radius) = 0;
       virtual void initialize_statistics(Statistics& statistics, const Options& options) = 0;
       virtual void generate_initial_iterate(Iterate& initial_iterate) = 0;
+      virtual void evaluate_progress_measures(Iterate& iterate, const std::optional<Scaling>& scaling) const = 0;
       virtual void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius,
-         WarmstartInformation& warmstart_information) = 0;
+         const std::optional<Scaling>& scaling, WarmstartInformation& warmstart_information) = 0;
 
       virtual void initialize_feasibility_problem(Iterate& current_iterate) = 0;
       virtual void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) = 0;
@@ -49,8 +53,8 @@ namespace uno {
 
       // progress measures
       [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
-         UserCallbacks& user_callbacks) = 0;
+         const std::optional<Scaling>& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, UserCallbacks& user_callbacks) = 0;
 
       virtual void postprocess_iterate(Iterate& iterate) = 0;
 
@@ -65,10 +69,11 @@ namespace uno {
       // when the parameterization of the subproblem (e.g. penalty or barrier parameter) is updated, signal it
       bool subproblem_definition_changed{false};
 
-      void evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate) const;
+      void evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate, const std::optional<Scaling>& scaling) const;
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         const Subproblem& subproblem, const EvaluationSpace& evaluation_space, Iterate& current_iterate, Iterate& trial_iterate,
-         const Direction& direction, double step_length, UserCallbacks& user_callbacks);
+         const Subproblem& subproblem, const std::optional<Scaling>& scaling, const EvaluationSpace& evaluation_space,
+         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
+         UserCallbacks& user_callbacks);
    };
 } // namespace
 
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 2f79a8a89..fa0eb2203 100644
--- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp
+++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.cpp
@@ -5,7 +5,6 @@
 #include "InequalityConstrainedMethod.hpp"
 #include "optimization/Iterate.hpp"
 #include "ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp"
-#include "ingredients/hessian_models/HessianModel.hpp"
 #include "ingredients/subproblem/Subproblem.hpp"
 #include "ingredients/subproblem_solvers/BoxLPSolverFactory.hpp"
 #include "ingredients/subproblem_solvers/LPSolverFactory.hpp"
@@ -50,14 +49,19 @@ namespace uno {
       // do nothing
    }
 
-   void InequalityConstrainedMethod::generate_initial_iterate(Iterate& initial_iterate) {
-      this->evaluate_progress_measures(*this->problem, initial_iterate);
+   void InequalityConstrainedMethod::generate_initial_iterate(Iterate& /*iterate*/) {
+      // do nothing
+   }
+
+   void InequalityConstrainedMethod::evaluate_progress_measures(Iterate& iterate, const std::optional<Scaling>& scaling) const {
+      InequalityHandlingMethod::evaluate_progress_measures(*this->problem, iterate, scaling);
    }
 
    void InequalityConstrainedMethod::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction,
-        double trust_region_radius, WarmstartInformation& warmstart_information) {
+        double trust_region_radius, const std::optional<Scaling>& scaling, WarmstartInformation& warmstart_information) {
       // solve the subproblem
-      this->solver->solve(statistics, *this->subproblem, trust_region_radius, this->initial_point, direction, warmstart_information);
+      this->solver->solve(statistics, *this->subproblem, trust_region_radius, scaling, this->initial_point, direction,
+         warmstart_information);
       InequalityConstrainedMethod::compute_dual_displacements(current_iterate.multipliers, direction.multipliers);
       ++this->number_subproblems_solved;
       // reset the initial point
@@ -99,10 +103,10 @@ namespace uno {
    }
 
    bool InequalityConstrainedMethod::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
-         UserCallbacks& user_callbacks) {
+         const std::optional<Scaling>& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, UserCallbacks& user_callbacks) {
       return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem,
-         this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks);
+         scaling, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks);
    }
 
    void InequalityConstrainedMethod::postprocess_iterate(Iterate& /*iterate*/) {
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 f0c4360cc..668109009 100644
--- a/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/inequality_constrained_methods/InequalityConstrainedMethod.hpp
@@ -19,11 +19,13 @@ namespace uno {
       ~InequalityConstrainedMethod() override = default;
 
       void initialize(const OptimizationProblem& problem, Iterate& current_iterate,
-         HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy, double trust_region_radius) override;
+         HessianModel& hessian_model, InertiaCorrectionStrategy<double>& inertia_correction_strategy,
+         double trust_region_radius) override;
       void initialize_statistics(Statistics& statistics, const Options& options) override;
       void generate_initial_iterate(Iterate& initial_iterate) override;
+      void evaluate_progress_measures(Iterate& iterate, const std::optional<Scaling>& scaling) const override;
       void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius,
-         WarmstartInformation& warmstart_information) override;
+         const std::optional<Scaling>& scaling, WarmstartInformation& warmstart_information) override;
 
       void initialize_feasibility_problem(Iterate& current_iterate) override;
       void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& current_iterate) override;
@@ -35,8 +37,8 @@ namespace uno {
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
-         UserCallbacks& user_callbacks) override;
+         const std::optional<Scaling>& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, UserCallbacks& user_callbacks) override;
 
       void postprocess_iterate(Iterate& iterate) override;
 
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 45eb8547c..0400e0ba3 100644
--- a/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp
+++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/InteriorPointMethod.hpp
@@ -28,8 +28,9 @@ namespace uno {
          InertiaCorrectionStrategy<double>& inertia_correction_strategy, double trust_region_radius) override;
       void initialize_statistics(Statistics& statistics, const Options& options) override;
       void generate_initial_iterate(Iterate& initial_iterate) override;
+      void evaluate_progress_measures(Iterate& iterate, const std::optional<Scaling>& scaling) const override;
       void solve(Statistics& statistics, Iterate& current_iterate, Direction& direction, double trust_region_radius,
-         WarmstartInformation& warmstart_information) override;
+         const std::optional<Scaling>& scaling, WarmstartInformation& warmstart_information) override;
 
       void initialize_feasibility_problem(Iterate& current_iterate) override;
       void set_elastic_variable_values(const l1RelaxedProblem& problem, Iterate& iterate) override;
@@ -41,8 +42,8 @@ namespace uno {
 
       // acceptance
       [[nodiscard]] bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
-         UserCallbacks& user_callbacks) override;
+         const std::optional<Scaling>& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, UserCallbacks& user_callbacks) override;
 
       void postprocess_iterate(Iterate& iterate) override;
 
@@ -123,12 +124,16 @@ namespace uno {
       // resize the initial iterate
       initial_iterate.set_number_variables(this->barrier_problem->number_variables);
       this->barrier_problem->generate_initial_iterate(initial_iterate);
-      this->evaluate_progress_measures(*this->barrier_problem, initial_iterate);
+   }
+
+   template <typename BarrierProblem>
+   void InteriorPointMethod<BarrierProblem>::evaluate_progress_measures(Iterate& iterate, const std::optional<Scaling>& scaling) const {
+      InequalityHandlingMethod::evaluate_progress_measures(*this->barrier_problem, iterate, scaling);
    }
 
    template <typename BarrierProblem>
    void InteriorPointMethod<BarrierProblem>::solve(Statistics& statistics, Iterate& current_iterate, Direction& direction,
-         double trust_region_radius, WarmstartInformation& warmstart_information) {
+         double trust_region_radius, const std::optional<Scaling>& scaling, WarmstartInformation& warmstart_information) {
       if (is_finite(trust_region_radius)) {
          throw std::runtime_error("The interior-point subproblem has a trust region. This is not implemented yet");
       }
@@ -143,7 +148,7 @@ namespace uno {
       statistics.set("barrier", this->barrier_parameter());
 
       // compute the primal-dual solution
-      this->linear_solver->solve_indefinite_system(statistics, *this->subproblem, direction, warmstart_information);
+      this->linear_solver->solve_indefinite_system(statistics, *this->subproblem, scaling, direction, warmstart_information);
       ++this->number_subproblems_solved;
 
       // check whether the augmented matrix was singular, in which case the subproblem is infeasible
@@ -208,7 +213,8 @@ namespace uno {
       // where jacobian_coefficient = -1 for p, +1 for n
       // Note: IPOPT uses a '+' sign because they define the Lagrangian as f(x) + \lambda^T c(x)
       const double mu = this->barrier_parameter();
-      const auto elastic_setting_function = [&](Iterate& iterate, size_t /*constraint_index*/, size_t elastic_index, double jacobian_coefficient) {
+      const auto elastic_setting_function = [&](Iterate& iterate, size_t /*constraint_index*/, size_t elastic_index,
+            double jacobian_coefficient) {
          // precomputations
          const double constraint_j = 0.; // TODO the constraints are not stored here any more... this->constraints[constraint_index];
          const double rho = this->l1_constraint_violation_coefficient;
@@ -244,10 +250,10 @@ namespace uno {
 
    template <typename BarrierProblem>
    bool InteriorPointMethod<BarrierProblem>::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
-         Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction, double step_length,
-         UserCallbacks& user_callbacks) {
+         const std::optional<Scaling>& scaling, Iterate& current_iterate, Iterate& trial_iterate, const Direction& direction,
+         double step_length, UserCallbacks& user_callbacks) {
       return InequalityHandlingMethod::is_iterate_acceptable(statistics, globalization_strategy, *this->subproblem,
-         this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks);
+         scaling, this->get_evaluation_space(), current_iterate, trial_iterate, direction, step_length, user_callbacks);
    }
 
    template <typename BarrierProblem>
diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp
index 5c71ce0c3..d24f8414f 100644
--- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp
+++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp
@@ -69,8 +69,9 @@ namespace uno {
       this->first_reformulation.evaluate_constraints(iterate, constraints);
    }
 
-   void PrimalDualInteriorPointProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const {
-      this->first_reformulation.evaluate_objective_gradient(iterate, objective_gradient);
+   void PrimalDualInteriorPointProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient,
+         const std::optional<Scaling>& scaling) const {
+      this->first_reformulation.evaluate_objective_gradient(iterate, objective_gradient, scaling);
 
       // barrier terms
       for (size_t variable_index: Range(this->first_reformulation.number_variables)) {
@@ -153,8 +154,8 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const {
-      this->first_reformulation.evaluate_lagrangian_gradient(lagrangian_gradient, evaluation_space, iterate);
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const {
+      this->first_reformulation.evaluate_lagrangian_gradient(lagrangian_gradient, evaluation_space, iterate, scaling);
 
       // barrier terms
       for (size_t variable_index: Range(this->first_reformulation.number_variables)) {
@@ -179,9 +180,10 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const {
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& scaling) const {
       // original Lagrangian Hessian
-      this->first_reformulation.evaluate_lagrangian_hessian(statistics, hessian_model, primal_variables, multipliers, hessian_values);
+      this->first_reformulation.evaluate_lagrangian_hessian(statistics, hessian_model, primal_variables, multipliers,
+         hessian_values, scaling);
 
       // barrier terms
       size_t nonzero_index = this->first_reformulation.number_hessian_nonzeros(hessian_model);
@@ -205,9 +207,9 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x,
-         const double* vector, const Multipliers& multipliers, double* result) const {
+         const double* vector, const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const {
       // original Lagrangian Hessian
-      this->first_reformulation.compute_hessian_vector_product(hessian_model, x, vector, multipliers, result);
+      this->first_reformulation.compute_hessian_vector_product(hessian_model, x, vector, multipliers, result, scaling);
 
       // barrier terms
       for (size_t variable_index: Range(this->first_reformulation.number_variables)) {
diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp
index 1da605d4b..5f3add3de 100644
--- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp
+++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp
@@ -23,7 +23,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const override;
 
       // dense objective gradient
-      void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const override;
+      void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional<Scaling>& scaling) const override;
 
       // sparsity patterns of Jacobian and Hessian
       void compute_constraint_jacobian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing,
@@ -37,11 +37,11 @@ namespace uno {
       [[nodiscard]] size_t number_hessian_nonzeros(const HessianModel& hessian_model) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const override;
       void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const override;
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const override;
       void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const override;
+         const Multipliers& multipliers, double* hessian_values, const std::optional<Scaling>& scaling) const override;
       void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector,
-         const Multipliers& multipliers, double* result) const override;
+         const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp
index cf74d286a..bdb0b79b4 100644
--- a/uno/ingredients/subproblem/Subproblem.cpp
+++ b/uno/ingredients/subproblem/Subproblem.cpp
@@ -74,10 +74,10 @@ namespace uno {
       }
    }
 
-   void Subproblem::evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values) const {
+   void Subproblem::evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values, const std::optional<Scaling>& scaling) const {
       // evaluate the Lagrangian Hessian of the problem at the current primal-dual point
       this->problem.evaluate_lagrangian_hessian(statistics, this->hessian_model, this->current_iterate.primals,
-         this->current_iterate.multipliers, hessian_values);
+         this->current_iterate.multipliers, hessian_values, scaling);
    }
 
    void Subproblem::regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const {
@@ -92,9 +92,11 @@ namespace uno {
       }
    }
 
-   void Subproblem::compute_hessian_vector_product(const double* x, const double* vector, double* result) const {
+   void Subproblem::compute_hessian_vector_product(const double* x, const double* vector, double* result,
+         const std::optional<Scaling>& scaling) const {
       // unregularized Hessian-vector product
-      this->problem.compute_hessian_vector_product(this->hessian_model, x, vector, this->current_iterate.multipliers, result);
+      this->problem.compute_hessian_vector_product(this->hessian_model, x, vector, this->current_iterate.multipliers,
+         result, scaling);
 
       // contribution of the regularization strategy
       const double regularization_factor = this->inertia_correction_strategy.get_primal_regularization_factor();
@@ -105,10 +107,10 @@ namespace uno {
       }
    }
 
-   void Subproblem::assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values) const {
+   void Subproblem::assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, const std::optional<Scaling>& scaling) const {
       // evaluate the Lagrangian Hessian of the problem at the current primal-dual point
       this->problem.evaluate_lagrangian_hessian(statistics, this->hessian_model, this->current_iterate.primals,
-         this->current_iterate.multipliers, augmented_matrix_values);
+         this->current_iterate.multipliers, augmented_matrix_values, scaling);
 
       // Jacobian of general constraints
       this->problem.evaluate_constraint_jacobian(this->current_iterate, augmented_matrix_values + this->number_hessian_nonzeros());
@@ -241,25 +243,27 @@ namespace uno {
    }
 
    std::function<double(double)> Subproblem::compute_predicted_objective_reduction(const EvaluationSpace& evaluation_space,
-         const Vector<double>& primal_direction, double step_length) const {
+         const Vector<double>& primal_direction, double step_length, const std::optional<Scaling>& scaling) const {
       // predicted objective reduction: "-∇f(x)^T (αd) - α^2/2 d^T H d"
       const double directional_derivative = dot(primal_direction, this->current_iterate.evaluations.objective_gradient);
       // if the regularized Hessian is positive definite (as it usually is in line-search methods), we can compute the
       // predicted reduction with only first-order information (the directional derivative)
       const bool is_regularized_hessian_positive_definite = this->hessian_model.is_positive_definite() && this->performs_primal_regularization();
       const double quadratic_term = is_regularized_hessian_positive_definite ? 0. :
-         evaluation_space.compute_hessian_quadratic_product(*this, primal_direction);
+         evaluation_space.compute_hessian_quadratic_product(*this, scaling, primal_direction);
+      const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.;
       return [=](double objective_multiplier) {
-         return step_length * (-objective_multiplier*directional_derivative) - step_length*step_length/2. * quadratic_term;
+         return step_length * (-objective_multiplier * objective_scaling * directional_derivative) -
+            step_length*step_length/2. * quadratic_term;
       };
    }
 
    ProgressMeasures Subproblem::compute_predicted_reductions(const EvaluationSpace& evaluation_space, const Direction& direction,
-         double step_length, Norm norm) const {
+         double step_length, Norm norm, const std::optional<Scaling>& scaling) const {
       return {
          this->compute_predicted_infeasibility_reduction(evaluation_space, this->problem.model, direction.primals,
             step_length, norm),
-         this->compute_predicted_objective_reduction(evaluation_space, direction.primals, step_length),
+         this->compute_predicted_objective_reduction(evaluation_space, direction.primals, step_length, scaling),
          this->problem.compute_predicted_auxiliary_reduction(this->current_iterate, direction.primals, step_length)
       };
    }
diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp
index 7056096ad..3a1d8ca2d 100644
--- a/uno/ingredients/subproblem/Subproblem.hpp
+++ b/uno/ingredients/subproblem/Subproblem.hpp
@@ -5,11 +5,13 @@
 #define UNO_SUBPROBLEM_H
 
 #include <functional>
+#include <optional>
 #include "optimization/Iterate.hpp"
 #include "linear_algebra/Matrix.hpp"
 #include "linear_algebra/MatrixOrder.hpp"
 #include "optimization/Multipliers.hpp"
 #include "optimization/OptimizationProblem.hpp"
+#include "optimization/Scaling.hpp"
 #include "symbolic/Range.hpp"
 #include "symbolic/UnaryNegation.hpp"
 #include "symbolic/VectorView.hpp"
@@ -42,12 +44,12 @@ namespace uno {
          const uno_int* jacobian_row_indices, const uno_int* jacobian_column_indices, uno_int solver_indexing) const;
 
       // regularized Hessian
-      void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values) const;
+      void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values, const std::optional<Scaling>& scaling) const;
       void regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const;
-      void compute_hessian_vector_product(const double* x, const double* vector, double* result) const;
+      void compute_hessian_vector_product(const double* x, const double* vector, double* result, const std::optional<Scaling>& scaling) const;
 
       // augmented system
-      void assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values) const;
+      void assemble_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, const std::optional<Scaling>& scaling) const;
       void regularize_augmented_matrix(Statistics& statistics, double* augmented_matrix_values,
          double dual_regularization_parameter, DirectSymmetricIndefiniteLinearSolver<double>& linear_solver) const;
       template <typename IndexType>
@@ -86,9 +88,9 @@ namespace uno {
       [[nodiscard]] double compute_predicted_infeasibility_reduction(const EvaluationSpace& evaluation_space, const Model& model,
          const Vector<double>& primal_direction, double step_length, Norm norm) const;
       [[nodiscard]] std::function<double(double)> compute_predicted_objective_reduction(const EvaluationSpace& evaluation_space,
-         const Vector<double>& primal_direction, double step_length) const;
+         const Vector<double>& primal_direction, double step_length, const std::optional<Scaling>& scaling) const;
       [[nodiscard]] ProgressMeasures compute_predicted_reductions(const EvaluationSpace& evaluation_space,
-         const Direction& direction, double step_length, Norm norm) const;
+         const Direction& direction, double step_length, Norm norm, const std::optional<Scaling>& scaling) const;
 
       const OptimizationProblem& problem;
       Iterate& current_iterate;
diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp
index d78898a81..ed4d4b691 100644
--- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp
+++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.cpp
@@ -76,12 +76,23 @@ namespace uno {
       }
    }
 
-   double BQPDEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& subproblem, const Vector<double>& vector) const {
+   void BQPDEvaluationSpace::compute_constraint_jacobian_norms(Vector<double>& row_norms) const {
+      row_norms.fill(0.);
+      const size_t number_constraint_jacobian_nonzeros = this->jacobian_row_indices.size();
+      for (size_t nonzero_index: Range(number_constraint_jacobian_nonzeros)) {
+         const size_t constraint_index = static_cast<size_t>(this->jacobian_row_indices[nonzero_index]);
+         const double derivative = this->jacobian_values[nonzero_index];
+         row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative));
+      }
+   }
+
+   double BQPDEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const Vector<double>& vector) const {
       if (subproblem.has_hessian_operator()) { // linear operator
          // TODO preallocate
          Vector<double> result(subproblem.number_variables);
          // compute Hv
-         subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), vector.data(), result.data());
+         subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), vector.data(), result.data(), scaling);
          // compute the dot product <v, Hv>
          return dot(vector, result);
       }
@@ -105,16 +116,16 @@ namespace uno {
       }
    }
 
-   void BQPDEvaluationSpace::evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate,
-         const WarmstartInformation& warmstart_information) {
+   void BQPDEvaluationSpace::evaluate_functions(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Iterate& current_iterate, const WarmstartInformation& warmstart_information) {
       // evaluate the functions based on warmstart information
       // gradients is a concatenation of the dense objective gradient and the sparse Jacobian
       if (warmstart_information.objective_changed) {
-         problem.evaluate_objective_gradient(current_iterate, this->gradients.data());
+         subproblem.problem.evaluate_objective_gradient(current_iterate, this->gradients.data(), scaling);
       }
       if (warmstart_information.constraints_changed) {
-         problem.evaluate_constraints(current_iterate, this->constraints);
-         this->evaluate_constraint_jacobian(problem, current_iterate);
+         subproblem.problem.evaluate_constraints(current_iterate, this->constraints);
+         this->evaluate_constraint_jacobian(subproblem.problem, current_iterate);
       }
       if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
          this->evaluate_hessian = true;
diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp
index 699aaf39a..84f91f0c8 100644
--- a/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp
+++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDEvaluationSpace.hpp
@@ -5,8 +5,10 @@
 #define UNO_BQPDEVALUATIONSPACE_H
 
 #include <cstddef>
+#include <optional>
 #include "linear_algebra/Vector.hpp"
 #include "optimization/EvaluationSpace.hpp"
+#include "optimization/Scaling.hpp"
 #include "../interfaces/C/uno_int.h"
 
 namespace uno {
@@ -25,9 +27,12 @@ namespace uno {
       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;
+      void compute_constraint_jacobian_norms(Vector<double>& row_norms) const override;
+      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const Vector<double>& vector) const override;
 
-      void evaluate_functions(const OptimizationProblem& problem, Iterate& current_iterate, const WarmstartInformation& warmstart_information);
+      void evaluate_functions(const Subproblem& subproblem, const std::optional<Scaling>& scaling, Iterate& current_iterate,
+         const WarmstartInformation& warmstart_information);
 
       Vector<double> constraints{};
       Vector<double> gradients{};
diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp
index 2fa610282..1b9559944 100644
--- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp
+++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp
@@ -98,17 +98,17 @@ namespace uno {
       this->mxws = static_cast<size_t>(this->kmax * (this->kmax + 9) / 2) + 2 * subproblem.number_variables +
          subproblem.number_constraints /* (required by bqpd.f) */ + 5 * subproblem.number_variables + this->nprof /* (required
          by sparseL.f) */;
-      // 6 pointers hidden in lws
-      constexpr size_t hidden_pointers_size = 6*sizeof(intptr_t);
+      // 7 pointers hidden in lws
+      constexpr size_t hidden_pointers_size = 7*sizeof(intptr_t);
       this->mxlws = hidden_pointers_size + static_cast<size_t>(this->kmax) /* (required by bqpd.f) */ +
          9 * subproblem.number_variables + subproblem.number_constraints /* (required by sparseL.f) */;
       this->ws.resize(this->mxws);
       this->lws.resize(this->mxlws);
    }
 
-   void BQPDSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius,
+   void BQPDSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
          const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) {
-      this->set_up_subproblem(statistics, subproblem, trust_region_radius, warmstart_information);
+      this->set_up_subproblem(statistics, subproblem, trust_region_radius, scaling, warmstart_information);
       if (this->print_subproblem) {
          this->display_subproblem(subproblem, initial_point);
       }
@@ -122,14 +122,14 @@ namespace uno {
    // protected member functions
 
    void BQPDSolver::set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
-         const WarmstartInformation& warmstart_information) {
+         const std::optional<Scaling>& scaling, const WarmstartInformation& warmstart_information) {
       // initialize wsc_ common block (Hessian & workspace for BQPD)
       // setting the common block here ensures that several instances of BQPD can run simultaneously
       WSC.mxws = static_cast<int>(this->mxws);
       WSC.mxlws = static_cast<int>(this->mxlws);
 
       // evaluate the functions and derivatives
-      this->evaluation_space.evaluate_functions(subproblem.problem, subproblem.current_iterate, warmstart_information);
+      this->evaluation_space.evaluate_functions(subproblem, scaling, subproblem.current_iterate, warmstart_information);
 
       // variable bounds
       if (warmstart_information.variable_bounds_changed) {
@@ -149,7 +149,7 @@ namespace uno {
          this->upper_bounds[variable_index] = std::min(BIG, this->upper_bounds[variable_index]);
       }
 
-      this->hide_pointers_in_workspace(statistics, subproblem);
+      this->hide_pointers_in_workspace(statistics, subproblem, scaling);
    }
 
    void BQPDSolver::display_subproblem(const Subproblem& subproblem, const Vector<double>& initial_point) const {
@@ -242,7 +242,8 @@ namespace uno {
    }
 
    // hide pointers to arbitrary objects into this->workspace_sparsity (BQPD's lws)
-   void BQPDSolver::hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem) {
+   void BQPDSolver::hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem,
+         const std::optional<Scaling>& scaling) {
       WSC.kk = 0; // length of ws that is used by gdotx
       WSC.ll = 0; // length of lws that is used by gdotx
 
@@ -253,11 +254,13 @@ namespace uno {
       WSC.ll += sizeof(intptr_t);
       hide_pointer(2, this->lws.data(), subproblem);
       WSC.ll += sizeof(intptr_t);
-      hide_pointer(3, this->lws.data(), this->evaluation_space.hessian_row_indices);
+      hide_pointer(3, this->lws.data(), scaling);
+      WSC.ll += sizeof(intptr_t);
+      hide_pointer(4, this->lws.data(), this->evaluation_space.hessian_row_indices);
       WSC.ll += sizeof(intptr_t);
-      hide_pointer(4, this->lws.data(), this->evaluation_space.hessian_column_indices);
+      hide_pointer(5, this->lws.data(), this->evaluation_space.hessian_column_indices);
       WSC.ll += sizeof(intptr_t);
-      hide_pointer(5, this->lws.data(), this->evaluation_space.hessian_values);
+      hide_pointer(6, this->lws.data(), this->evaluation_space.hessian_values);
       WSC.ll += sizeof(intptr_t);
    }
 
@@ -343,12 +346,14 @@ void hessian_vector_product(int* dimension, const double vector[], const double
    bool* evaluate_hessian = uno::retrieve_pointer<bool>(0, lws);
    uno::Statistics* statistics = uno::retrieve_pointer<uno::Statistics>(1, lws);
    uno::Subproblem* subproblem = uno::retrieve_pointer<uno::Subproblem>(2, lws);
-   uno::Vector<int>* hessian_row_indices = uno::retrieve_pointer<uno::Vector<int>>(3, lws);
-   uno::Vector<int>* hessian_column_indices = uno::retrieve_pointer<uno::Vector<int>>(4, lws);
-   uno::Vector<double>* hessian_values = uno::retrieve_pointer<uno::Vector<double>>(5, lws);
+   const std::optional<uno::Scaling>* scaling = uno::retrieve_pointer<const std::optional<uno::Scaling>>(3, lws);
+   uno::Vector<int>* hessian_row_indices = uno::retrieve_pointer<uno::Vector<int>>(4, lws);
+   uno::Vector<int>* hessian_column_indices = uno::retrieve_pointer<uno::Vector<int>>(5, lws);
+   uno::Vector<double>* hessian_values = uno::retrieve_pointer<uno::Vector<double>>(6, lws);
    assert(evaluate_hessian != nullptr);
    assert(statistics != nullptr);
    assert(subproblem != nullptr);
+   assert(scaling != nullptr);
    assert(hessian_row_indices != nullptr);
    assert(hessian_column_indices != nullptr);
    assert(hessian_values != nullptr);
@@ -360,7 +365,7 @@ void hessian_vector_product(int* dimension, const double vector[], const double
       if (subproblem->has_hessian_matrix()) {
          // if the Hessian has not been evaluated at the current point, evaluate it
          if (*evaluate_hessian) {
-            subproblem->evaluate_lagrangian_hessian(*statistics, hessian_values->data());
+            subproblem->evaluate_lagrangian_hessian(*statistics, hessian_values->data(), *scaling);
             subproblem->regularize_lagrangian_hessian(*statistics, hessian_values->data());
             *evaluate_hessian = false;
          }
@@ -381,6 +386,6 @@ void hessian_vector_product(int* dimension, const double vector[], const double
    }
    // otherwise, try to perform a Hessian-vector product if possible
    else if (subproblem->has_hessian_operator()) {
-      subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.data(), vector, result);
+      subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.data(), vector, result, *scaling);
    }
 }
\ No newline at end of file
diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp
index d52e95d6d..f6ece1d6f 100644
--- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.hpp
@@ -45,8 +45,8 @@ namespace uno {
 
       void initialize_memory(const Subproblem& subproblem) override;
 
-      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
-         Direction& direction, const WarmstartInformation& warmstart_information) override;
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
+         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] EvaluationSpace& get_evaluation_space() override;
 
@@ -74,12 +74,12 @@ namespace uno {
       const bool print_subproblem;
 
       void set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
-         const WarmstartInformation& warmstart_information);
+         const std::optional<Scaling>& scaling, const WarmstartInformation& warmstart_information);
       void display_subproblem(const Subproblem& subproblem, const Vector<double>& initial_point) const;
       void solve_subproblem(const Subproblem& subproblem, const Vector<double>& initial_point, Direction& direction,
          const WarmstartInformation& warmstart_information);
       [[nodiscard]] static BQPDMode determine_mode(const WarmstartInformation& warmstart_information);
-      void hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem);
+      void hide_pointers_in_workspace(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling);
       void set_multipliers(size_t number_variables, Multipliers& direction_multipliers) const;
       [[nodiscard]] static BQPDStatus bqpd_status_from_int(int ifail);
       [[nodiscard]] bool check_sufficient_workspace_size(BQPDStatus bqpd_status);
diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp
index 8259c0566..fe030801f 100644
--- a/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp
+++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.cpp
@@ -14,12 +14,14 @@ namespace uno {
    }
 
    void BoxLPSolver::solve(Statistics& /*statistics*/, Subproblem& subproblem, double trust_region_radius,
-         const Vector<double>& /*initial_point*/, Direction& direction, const WarmstartInformation& /*warmstart_information*/) {
+         const std::optional<Scaling>& scaling, const Vector<double>& /*initial_point*/, Direction& direction,
+         const WarmstartInformation& /*warmstart_information*/) {
       if (0 < subproblem.number_constraints) {
          throw std::runtime_error("BoxLPSolver cannot solve problems with general constraints");
       }
       // compute the objective gradient
-      subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->evaluation_space.objective_gradient.data());
+      subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->evaluation_space.objective_gradient.data(),
+         scaling);
 
       // compute the variables bounds
       subproblem.set_variables_bounds(this->variable_lower_bounds, this->variable_upper_bounds, trust_region_radius);
diff --git a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp
index bb4a4c866..4f7fae960 100644
--- a/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/BoxLPSolver.hpp
@@ -18,7 +18,11 @@ namespace uno {
       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 {
+      void compute_constraint_jacobian_norms(Vector<double>& row_norms) const override {
+         row_norms.fill(0.);
+      }
+      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional<Scaling>& /*scaling*/,
+            const Vector<double>& /*vector*/) const override {
          return 0.;
       }
 
@@ -32,8 +36,8 @@ namespace uno {
 
       void initialize_memory(const Subproblem& subproblem) override;
 
-      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
-         Direction& direction, const WarmstartInformation& warmstart_information) override;
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
+         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] EvaluationSpace& get_evaluation_space() override;
 
diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp
index e7441d192..12f7e7d9c 100644
--- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp
+++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.cpp
@@ -95,15 +95,27 @@ namespace uno {
       }
    }
 
-   double COOEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector<double>& /*vector*/) const {
+   // compute the inf norm of each row
+   void COOEvaluationSpace::compute_constraint_jacobian_norms(Vector<double>& row_norms) const {
+      row_norms.fill(0.);
+      const size_t offset = this->number_hessian_nonzeros;
+      for (size_t nonzero_index: Range(this->number_jacobian_nonzeros)) {
+         const size_t constraint_index = static_cast<size_t>(this->jacobian_row_indices[nonzero_index]);
+         const double derivative = this->matrix_values[offset + nonzero_index];
+         row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative));
+      }
+   }
+
+   double COOEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional<Scaling>& /*scaling*/,
+         const Vector<double>& /*vector*/) const {
       return 0.;
    }
 
-   void COOEvaluationSpace::set_up_linear_system(Statistics& statistics, const Subproblem& subproblem,
+   void COOEvaluationSpace::set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
          DirectSymmetricIndefiniteLinearSolver<double>& linear_solver, const WarmstartInformation& warmstart_information) {
       // evaluate the functions at the current iterate
       if (warmstart_information.objective_changed) {
-         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data());
+         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->objective_gradient.data(), scaling);
       }
       if (warmstart_information.constraints_changed) {
          subproblem.problem.evaluate_constraints(subproblem.current_iterate, this->constraints);
@@ -117,7 +129,7 @@ namespace uno {
             this->analysis_performed = true;
          }
          // assemble the augmented matrix
-         subproblem.assemble_augmented_matrix(statistics, this->matrix_values.data());
+         subproblem.assemble_augmented_matrix(statistics, this->matrix_values.data(), scaling);
          // regularize the augmented matrix (this calls the analysis and the factorization)
          subproblem.regularize_augmented_matrix(statistics, this->matrix_values.data(),
             subproblem.dual_regularization_factor(), linear_solver);
diff --git a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp
index cb471526d..cf031de9a 100644
--- a/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp
+++ b/uno/ingredients/subproblem_solvers/COOEvaluationSpace.hpp
@@ -5,9 +5,12 @@
 #define UNO_COOEVALUATIONSPACE_H
 
 #include <cstddef>
+#include <optional>
 #include <vector>
+#include "linear_algebra/Norm.hpp"
 #include "linear_algebra/Vector.hpp"
 #include "optimization/EvaluationSpace.hpp"
+#include "optimization/Scaling.hpp"
 #include "../interfaces/C/uno_int.h"
 
 namespace uno {
@@ -23,10 +26,12 @@ namespace uno {
       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;
+      void compute_constraint_jacobian_norms(Vector<double>& row_norms) const override;
+      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const Vector<double>& vector) const override;
 
-      void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, DirectSymmetricIndefiniteLinearSolver<double>& linear_solver,
-         const WarmstartInformation& warmstart_information);
+      void set_up_linear_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         DirectSymmetricIndefiniteLinearSolver<double>& linear_solver, const WarmstartInformation& warmstart_information);
 
       Vector<double> objective_gradient{}; /*!< Sparse Jacobian of the objective */
       Vector<double> constraints{}; /*!< Constraint values (size \f$m)\f$ */
diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp
index d2dbc723d..cd0e70e4a 100644
--- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp
+++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.cpp
@@ -75,7 +75,19 @@ namespace uno {
       }
    }
 
-   double HiGHSEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const Vector<double>& vector) const {
+   void HiGHSEvaluationSpace::compute_constraint_jacobian_norms(Vector<double>& row_norms) const {
+      row_norms.fill(0.);
+      const size_t number_constraint_jacobian_nonzeros = this->jacobian_row_indices.size();
+      for (size_t nonzero_index: Range(number_constraint_jacobian_nonzeros)) {
+         const size_t permuted_nonzero_index = this->jacobian_permutation_vector[nonzero_index];
+         const size_t constraint_index = static_cast<size_t>(this->jacobian_row_indices[permuted_nonzero_index]);
+         const double derivative = this->model.lp_.a_matrix_.value_[nonzero_index];
+         row_norms[constraint_index] = std::max(row_norms[constraint_index], std::abs(derivative));
+      }
+   }
+
+   double HiGHSEvaluationSpace::compute_hessian_quadratic_product(const Subproblem& /*subproblem*/, const std::optional<Scaling>& /*scaling*/,
+         const Vector<double>& vector) const {
       double quadratic_product = 0.;
       const size_t number_hessian_nonzeros = this->hessian_values.size();
       for (size_t nonzero_index: Range(number_hessian_nonzeros)) {
@@ -92,10 +104,10 @@ namespace uno {
    }
 
    void HiGHSEvaluationSpace::evaluate_functions(Statistics& statistics, const Subproblem& subproblem,
-         const WarmstartInformation& warmstart_information) {
+         const std::optional<Scaling>& scaling, const WarmstartInformation& warmstart_information) {
       // evaluate the functions based on warmstart information
       if (warmstart_information.objective_changed) {
-         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->model.lp_.col_cost_.data());
+         subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->model.lp_.col_cost_.data(), scaling);
       }
       if (warmstart_information.constraints_changed) {
          subproblem.problem.evaluate_constraints(subproblem.current_iterate, this->constraints);
@@ -103,7 +115,7 @@ namespace uno {
       }
       // evaluate the Hessian and regularize it
       if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
-         subproblem.evaluate_lagrangian_hessian(statistics, this->hessian_values.data());
+         subproblem.evaluate_lagrangian_hessian(statistics, this->hessian_values.data(), scaling);
          // copy the Hessian with permutation into this->model.hessian_.value_
          for (size_t nonzero_index: Range(subproblem.number_regularized_hessian_nonzeros())) {
             const size_t permuted_nonzero_index = this->hessian_permutation_vector[nonzero_index];
diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp
index e2ee04fd0..d599e41a7 100644
--- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp
+++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSEvaluationSpace.hpp
@@ -5,9 +5,11 @@
 #define UNO_HIGHSEVALUATIONSPACE_H
 
 #include <cstddef>
+#include <optional>
 #include "optimization/EvaluationSpace.hpp"
 #include "Highs.h"
 #include "linear_algebra/Vector.hpp"
+#include "optimization/Scaling.hpp"
 #include "../interfaces/C/uno_int.h"
 
 namespace uno {
@@ -27,9 +29,12 @@ namespace uno {
       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;
+      void compute_constraint_jacobian_norms(Vector<double>& row_norms) const override;
+      [[nodiscard]] double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const Vector<double>& vector) const override;
 
-      void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const WarmstartInformation& warmstart_information);
+      void evaluate_functions(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const WarmstartInformation& warmstart_information);
 
       HighsModel model;
       Vector<double> constraints{};
diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp
index 1cabd11d3..0783137fd 100644
--- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp
+++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp
@@ -20,8 +20,9 @@ namespace uno {
    }
 
    void HiGHSSolver::solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius,
-         const Vector<double>& /*initial_point*/, Direction& direction, const WarmstartInformation& warmstart_information) {
-      this->set_up_subproblem(statistics, subproblem, trust_region_radius, warmstart_information);
+         const std::optional<Scaling>& scaling, const Vector<double>& /*initial_point*/, Direction& direction,
+         const WarmstartInformation& warmstart_information) {
+      this->set_up_subproblem(statistics, subproblem, trust_region_radius, scaling, warmstart_information);
       this->solve_subproblem(subproblem, direction);
    }
 
@@ -32,9 +33,9 @@ namespace uno {
    // protected member functions
 
    void HiGHSSolver::set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
-         const WarmstartInformation& warmstart_information) {
+         const std::optional<Scaling>& scaling, const WarmstartInformation& warmstart_information) {
       // evaluate the functions and derivatives
-      this->evaluation_space.evaluate_functions(statistics, subproblem, warmstart_information);
+      this->evaluation_space.evaluate_functions(statistics, subproblem, scaling, warmstart_information);
 
       // variable bounds
       if (warmstart_information.variable_bounds_changed) {
diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp
index 5c67a5a35..48109213a 100644
--- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.hpp
@@ -18,8 +18,8 @@ namespace uno {
 
       void initialize_memory(const Subproblem& subproblem) override;
 
-      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
-         Direction& direction, const WarmstartInformation& warmstart_information) override;
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
+         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] EvaluationSpace& get_evaluation_space() override;
 
@@ -30,7 +30,7 @@ namespace uno {
       const bool print_subproblem;
 
       void set_up_subproblem(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius,
-         const WarmstartInformation& warmstart_information);
+         const std::optional<Scaling>& scaling, const WarmstartInformation& warmstart_information);
       void solve_subproblem(const Subproblem& subproblem, Direction& direction);
    };
 } // namespace
diff --git a/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp b/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp
index 5712efb7f..6d6807a11 100644
--- a/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/InequalityConstrainedSolver.hpp
@@ -4,6 +4,9 @@
 #ifndef UNO_INEQUALITYCONSTRAINEDSOLVER_H
 #define UNO_INEQUALITYCONSTRAINEDSOLVER_H
 
+#include <optional>
+#include "optimization/Scaling.hpp"
+
 namespace uno {
    // forward declarations
    class Direction;
@@ -22,7 +25,8 @@ namespace uno {
       virtual void initialize_memory(const Subproblem& subproblem) = 0;
 
       virtual void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius,
-         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) = 0;
+         const std::optional<Scaling>& scaling, const Vector<double>& initial_point, Direction& direction,
+         const WarmstartInformation& warmstart_information) = 0;
 
       [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0;
    };
diff --git a/uno/ingredients/subproblem_solvers/LPSolver.hpp b/uno/ingredients/subproblem_solvers/LPSolver.hpp
index b2664f47e..80ec85216 100644
--- a/uno/ingredients/subproblem_solvers/LPSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/LPSolver.hpp
@@ -14,8 +14,8 @@ namespace uno {
 
       void initialize_memory(const Subproblem& subproblem) override = 0;
 
-      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
-         Direction& direction, const WarmstartInformation& warmstart_information) override = 0;
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
+         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override = 0;
 
       [[nodiscard]] EvaluationSpace& get_evaluation_space() override = 0;
    };
diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp
index 6cb3f449a..1774e520e 100644
--- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp
+++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp
@@ -231,10 +231,10 @@ namespace uno {
       }
    }
 
-   void MA27Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) {
+   void MA27Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) {
       // set up the linear system by evaluating the functions at the current iterate
-      this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information);
+      this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information);
       // solve the linear system
       this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution);
       // assemble the full primal-dual direction
diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp
index 3a58f1b91..93cfd182a 100644
--- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp
+++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp
@@ -45,8 +45,8 @@ namespace uno {
       void do_symbolic_analysis() override;
       void do_numerical_factorization(const double* matrix_values) override;
       void solve_indefinite_system(const Vector<double>& matrix_values, const Vector<double>& rhs, Vector<double>& result) override;
-      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) override;
+      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] Inertia get_inertia() const override;
       [[nodiscard]] size_t number_negative_eigenvalues() const override;
diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp
index a5a81ab3e..e74e2cf5b 100644
--- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp
+++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp
@@ -203,10 +203,10 @@ namespace uno {
       }
    }
 
-   void MA57Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) {
+   void MA57Solver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) {
       // set up the linear system by evaluating the functions at the current iterate
-      this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information);
+      this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information);
       // solve the linear system
       this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution);
       // assemble the full primal-dual direction
diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp
index 13020acc5..921326345 100644
--- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp
+++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp
@@ -52,8 +52,8 @@ namespace uno {
       void do_symbolic_analysis() override;
       void do_numerical_factorization(const double* matrix_values) override;
       void solve_indefinite_system(const Vector<double>& matrix_values, const Vector<double>& rhs, Vector<double>& result) override;
-      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) override;
+      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] Inertia get_inertia() const override;
       [[nodiscard]] size_t number_negative_eigenvalues() const override;
diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp
index dc3c994f9..4341e31b2 100644
--- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp
+++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp
@@ -94,10 +94,10 @@ namespace uno {
       dmumps_c(&this->workspace);
    }
 
-   void MUMPSSolver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) {
+   void MUMPSSolver::solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) {
       // set up the linear system by evaluating the functions at the current iterate
-      this->evaluation_space.set_up_linear_system(statistics, subproblem, *this, warmstart_information);
+      this->evaluation_space.set_up_linear_system(statistics, subproblem, scaling, *this, warmstart_information);
       // solve the linear system
       this->solve_indefinite_system(this->evaluation_space.matrix_values, this->evaluation_space.rhs, this->evaluation_space.solution);
       // assemble the full primal-dual direction
diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp
index 06bc93783..ed8ff18a0 100644
--- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp
@@ -21,8 +21,8 @@ namespace uno {
       void do_symbolic_analysis() override;
       void do_numerical_factorization(const double* matrix_values) override;
       void solve_indefinite_system(const Vector<double>& matrix_values, const Vector<double>& rhs, Vector<double>& result) override;
-      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) override;
+      void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) override;
 
       [[nodiscard]] Inertia get_inertia() const override;
       [[nodiscard]] size_t number_negative_eigenvalues() const override;
diff --git a/uno/ingredients/subproblem_solvers/QPSolver.hpp b/uno/ingredients/subproblem_solvers/QPSolver.hpp
index 6e5faf874..6ed02e522 100644
--- a/uno/ingredients/subproblem_solvers/QPSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/QPSolver.hpp
@@ -14,8 +14,8 @@ namespace uno {
 
       void initialize_memory(const Subproblem& subproblem) override = 0;
 
-      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
-         Direction& direction, const WarmstartInformation& warmstart_information) override = 0;
+      void solve(Statistics& statistics, Subproblem& subproblem, double trust_region_radius, const std::optional<Scaling>& scaling,
+         const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) override = 0;
    };
 } // namespace
 
diff --git a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp
index 20390d524..101dd0216 100644
--- a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp
+++ b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp
@@ -5,6 +5,8 @@
 #define UNO_SYMMETRICINDEFINITELINEARSOLVER_H
 
 #include <cstddef>
+#include <optional>
+#include "optimization/Scaling.hpp"
 
 namespace uno {
    // forward declarations
@@ -27,8 +29,8 @@ namespace uno {
 
       virtual void solve_indefinite_system(const Vector<double>& matrix_values, const Vector<ElementType>& rhs,
          Vector<ElementType>& result) = 0;
-      virtual void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, Direction& direction,
-         const WarmstartInformation& warmstart_information) = 0;
+      virtual void solve_indefinite_system(Statistics& statistics, const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         Direction& direction, const WarmstartInformation& warmstart_information) = 0;
 
       [[nodiscard]] virtual EvaluationSpace& get_evaluation_space() = 0;
    };
diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp
index 83d562175..f1e240cc3 100644
--- a/uno/model/Model.hpp
+++ b/uno/model/Model.hpp
@@ -4,11 +4,13 @@
 #ifndef UNO_MODEL_H
 #define UNO_MODEL_H
 
+#include <optional>
 #include <string>
 #include <vector>
 #include "linear_algebra/MatrixOrder.hpp"
 #include "linear_algebra/Norm.hpp"
 #include "optimization/ProblemType.hpp"
+#include "optimization/Scaling.hpp"
 #include "symbolic/VectorExpression.hpp"
 #include "../interfaces/C/uno_int.h"
 
diff --git a/uno/optimization/EvaluationSpace.hpp b/uno/optimization/EvaluationSpace.hpp
index d42a867b7..4afa88198 100644
--- a/uno/optimization/EvaluationSpace.hpp
+++ b/uno/optimization/EvaluationSpace.hpp
@@ -4,6 +4,9 @@
 #ifndef UNO_EVALUATIONSPACE_H
 #define UNO_EVALUATIONSPACE_H
 
+#include <optional>
+#include "optimization/Scaling.hpp"
+
 namespace uno {
    // forward declaration
    template <typename ElementType>
@@ -25,7 +28,9 @@ namespace uno {
       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;
+      virtual void compute_constraint_jacobian_norms(Vector<double>& row_norms) const = 0;
+      [[nodiscard]] virtual double compute_hessian_quadratic_product(const Subproblem& subproblem, const std::optional<Scaling>& scaling,
+         const Vector<double>& vector) const = 0;
    };
 } // namespace
 
diff --git a/uno/optimization/Iterate.hpp b/uno/optimization/Iterate.hpp
index c0906b5be..9996ee630 100644
--- a/uno/optimization/Iterate.hpp
+++ b/uno/optimization/Iterate.hpp
@@ -7,8 +7,8 @@
 #include "Evaluations.hpp"
 #include "SolutionStatus.hpp"
 #include "ingredients/globalization_strategies/ProgressMeasures.hpp"
-#include "optimization/Multipliers.hpp"
 #include "optimization/DualResiduals.hpp"
+#include "optimization/Multipliers.hpp"
 
 namespace uno {
    // forward declaration
diff --git a/uno/optimization/OptimizationProblem.cpp b/uno/optimization/OptimizationProblem.cpp
index 9500e4a50..4be1c23ce 100644
--- a/uno/optimization/OptimizationProblem.cpp
+++ b/uno/optimization/OptimizationProblem.cpp
@@ -30,10 +30,11 @@ namespace uno {
       constraints = iterate.evaluations.constraints;
    }
 
-   void OptimizationProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const {
+   void OptimizationProblem::evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional<Scaling>& scaling) const {
       iterate.evaluate_objective_gradient(this->model);
+      const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.;
       for (size_t index: Range(this->number_variables)) {
-         objective_gradient[index] = iterate.evaluations.objective_gradient[index];
+         objective_gradient[index] = objective_scaling * iterate.evaluations.objective_gradient[index];
       }
    }
 
@@ -50,12 +51,12 @@ namespace uno {
    }
 
    void OptimizationProblem::compute_constraint_jacobian_sparsity(uno_int *row_indices, uno_int *column_indices, uno_int solver_indexing,
-                                                                  MatrixOrder matrix_order) const {
+         MatrixOrder matrix_order) const {
       this->model.compute_constraint_jacobian_sparsity(row_indices, column_indices, solver_indexing, matrix_order);
    }
 
    void OptimizationProblem::compute_hessian_sparsity(const HessianModel& hessian_model, uno_int *row_indices,
-                                                      uno_int *column_indices, uno_int solver_indexing) const {
+         uno_int *column_indices, uno_int solver_indexing) const {
       hessian_model.compute_sparsity(row_indices, column_indices, solver_indexing);
    }
 
@@ -66,12 +67,12 @@ namespace uno {
    // Lagrangian gradient ∇f(x_k) - ∇c(x_k) y_k - z_k
    // split in two parts: objective contribution and constraints' contribution
    void OptimizationProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const {
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const {
       lagrangian_gradient.objective_contribution.fill(0.);
       lagrangian_gradient.constraints_contribution.fill(0.);
 
       // ∇f(x_k)
-      this->evaluate_objective_gradient(iterate, lagrangian_gradient.objective_contribution.data());
+      this->evaluate_objective_gradient(iterate, lagrangian_gradient.objective_contribution.data(), scaling);
 
       // ∇c(x_k) λ_k
       evaluation_space.compute_constraint_jacobian_transposed_vector_product(iterate.multipliers.constraints,
@@ -86,14 +87,16 @@ namespace uno {
    }
 
    void OptimizationProblem::evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model,
-         const Vector<double>& primal_variables, const Multipliers& multipliers, double* hessian_values) const {
+         const Vector<double>& primal_variables, const Multipliers& multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) const {
       hessian_model.evaluate_hessian(statistics, primal_variables, this->get_objective_multiplier(),
-         multipliers.constraints, hessian_values);
+         multipliers.constraints, hessian_values, scaling);
    }
 
    void OptimizationProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector,
-         const Multipliers& multipliers, double* result) const {
-      hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result);
+         const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const {
+      hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints,
+         result, scaling);
    }
 
    size_t OptimizationProblem::get_number_original_variables() const {
@@ -213,11 +216,12 @@ namespace uno {
    }
 
    // objective measure: scaled objective
-   void OptimizationProblem::set_objective_measure(Iterate& iterate) const {
+   void OptimizationProblem::set_objective_measure(Iterate& iterate, const std::optional<Scaling>& scaling) const {
       iterate.evaluate_objective(this->model);
-      const double objective = iterate.evaluations.objective;
+      const double objective_value = iterate.evaluations.objective;
+      const double objective_scaling = scaling.has_value() ? scaling->get_objective_scaling() : 1.;
       iterate.progress.objective = [=](double objective_multiplier) {
-         return objective_multiplier * objective;
+         return objective_multiplier * objective_scaling * objective_value;
       };
    }
 
diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp
index 3ffaf62ec..89152d866 100644
--- a/uno/optimization/OptimizationProblem.hpp
+++ b/uno/optimization/OptimizationProblem.hpp
@@ -37,7 +37,7 @@ namespace uno {
       virtual void evaluate_constraints(Iterate& iterate, Vector<double>& constraints) const;
 
       // dense objective gradient
-      virtual void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient) const;
+      virtual void evaluate_objective_gradient(Iterate& iterate, double* objective_gradient, const std::optional<Scaling>& scaling) const;
 
       // sparsity patterns of Jacobian and Hessian
       [[nodiscard]] virtual size_t number_jacobian_nonzeros() const;
@@ -51,11 +51,12 @@ namespace uno {
       // numerical evaluations of Jacobian and Hessian
       virtual void evaluate_constraint_jacobian(Iterate& iterate, double* jacobian_values) const;
       virtual void evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient,
-         const EvaluationSpace& evaluation_space, Iterate& iterate) const;
-      virtual void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model, const Vector<double>& primal_variables,
-         const Multipliers& multipliers, double* hessian_values) const;
+         const EvaluationSpace& evaluation_space, Iterate& iterate, const std::optional<Scaling>& scaling) const;
+      virtual void evaluate_lagrangian_hessian(Statistics& statistics, HessianModel& hessian_model,
+         const Vector<double>& primal_variables, const Multipliers& multipliers, double* hessian_values,
+         const std::optional<Scaling>& scaling) const;
       virtual void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector,
-         const Multipliers& multipliers, double* result) const;
+         const Multipliers& multipliers, double* result, const std::optional<Scaling>& scaling) const;
 
       [[nodiscard]] size_t get_number_original_variables() const;
       [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const;
@@ -69,11 +70,12 @@ namespace uno {
       [[nodiscard]] virtual const Collection<size_t>& get_inequality_constraints() const;
       [[nodiscard]] virtual const Collection<size_t>& get_dual_regularization_constraints() const;
 
-      virtual void assemble_primal_dual_direction(const Iterate& current_iterate, const Vector<double>& solution, Direction& direction) const;
+      virtual void assemble_primal_dual_direction(const Iterate& current_iterate, const Vector<double>& solution,
+         Direction& direction) const;
       [[nodiscard]] virtual double dual_regularization_factor() const;
 
-      [[nodiscard]] static double stationarity_error(const LagrangianGradient<double>& lagrangian_gradient, double objective_multiplier,
-         Norm residual_norm);
+      [[nodiscard]] static double stationarity_error(const LagrangianGradient<double>& lagrangian_gradient,
+         double objective_multiplier, Norm residual_norm);
       [[nodiscard]] virtual double complementarity_error(const Vector<double>& primals, const Vector<double>& constraints,
          const Multipliers& multipliers, double shift_value, Norm residual_norm) const;
 
@@ -82,7 +84,7 @@ namespace uno {
 
       // progress measures
       virtual void set_infeasibility_measure(Iterate& iterate, Norm norm) const;
-      virtual void set_objective_measure(Iterate& iterate) const;
+      virtual void set_objective_measure(Iterate& iterate, const std::optional<Scaling>& scaling) const;
       virtual void set_auxiliary_measure(Iterate& iterate) const;
       [[nodiscard]] virtual double compute_predicted_auxiliary_reduction(const Iterate& current_iterate,
          const Vector<double>& primal_direction, double step_length) const;
diff --git a/uno/optimization/Scaling.cpp b/uno/optimization/Scaling.cpp
new file mode 100644
index 000000000..ab8f5b058
--- /dev/null
+++ b/uno/optimization/Scaling.cpp
@@ -0,0 +1,38 @@
+// Copyright (c) 2018-2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include "Scaling.hpp"
+#include "Iterate.hpp"
+#include "linear_algebra/Norm.hpp"
+#include "linear_algebra/SparseVector.hpp"
+#include "linear_algebra/Vector.hpp"
+#include "optimization/EvaluationSpace.hpp"
+#include "tools/Logger.hpp"
+
+namespace uno {
+   Scaling::Scaling(const Iterate& initial_iterate, const EvaluationSpace& evaluation_space, double gradient_threshold):
+         gradient_threshold(gradient_threshold), objective_scaling(1.),
+         constraint_scaling(initial_iterate.evaluations.constraints.size(), 1.) {
+      // objective
+      this->objective_scaling = std::min(1., this->gradient_threshold / norm_inf(initial_iterate.evaluations.objective_gradient));
+
+      // constraints
+      Vector<double> row_norms(initial_iterate.evaluations.constraints.size());
+      evaluation_space.compute_constraint_jacobian_norms(row_norms);
+      for (size_t constraint_index: Range(this->constraint_scaling.size())) {
+         this->constraint_scaling[constraint_index] = std::min(1., this->gradient_threshold / row_norms[constraint_index]);
+      }
+
+      DEBUG << "Objective scaling: " << this->objective_scaling << '\n';
+      DEBUG << "Constraint scaling: "; print_vector(DEBUG, this->constraint_scaling);
+   }
+
+   double Scaling::get_objective_scaling() const {
+      return this->objective_scaling;
+   }
+
+   double Scaling::get_constraint_scaling(size_t constraint_index) const {
+      assert(constraint_index < this->constraint_scaling.size() && "The constraint index is not valid.");
+      return this->constraint_scaling[constraint_index];
+   }
+} // namespace
\ No newline at end of file
diff --git a/uno/optimization/Scaling.hpp b/uno/optimization/Scaling.hpp
new file mode 100644
index 000000000..a6853bdfc
--- /dev/null
+++ b/uno/optimization/Scaling.hpp
@@ -0,0 +1,30 @@
+// Copyright (c) 2018-2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_SCALING_H
+#define UNO_SCALING_H
+
+#include <cstddef>
+#include <vector>
+
+namespace uno {
+   // forward declarations
+   class EvaluationSpace;
+   class Iterate;
+
+   class Scaling {
+   public:
+      Scaling(const Iterate& initial_iterate, const EvaluationSpace& evaluation_space, double gradient_threshold);
+      ~Scaling() = default;
+
+      [[nodiscard]] double get_objective_scaling() const;
+      [[nodiscard]] double get_constraint_scaling(size_t constraint_index) const;
+
+   protected:
+      const double gradient_threshold;
+      double objective_scaling;
+      std::vector<double> constraint_scaling;
+   };
+} // namespace
+
+#endif // UNO_SCALING_H
\ No newline at end of file
diff --git a/uno/options/DefaultOptions.cpp b/uno/options/DefaultOptions.cpp
index f7ddf1d5a..6a45b264b 100644
--- a/uno/options/DefaultOptions.cpp
+++ b/uno/options/DefaultOptions.cpp
@@ -58,8 +58,8 @@ namespace uno {
       // Hessian model (exact|zero)
       options.set_string("hessian_model", "exact");
       options.set_string("inertia_correction_strategy", "primal");
-      // scale the functions (yes|no)
-      options.set_bool("scale_functions", false);
+      // use function scaling based on derivatives at the initial point (yes|no)
+      options.set_bool("use_function_scaling", false);
       options.set_double("function_scaling_threshold", 100.);
       // factor scaling
       options.set_double("function_scaling_factor", 100.);
diff --git a/uno/options/Options.hpp b/uno/options/Options.hpp
index 66d62e540..e5bd7d4b6 100644
--- a/uno/options/Options.hpp
+++ b/uno/options/Options.hpp
@@ -82,7 +82,7 @@ namespace uno {
          {"globalization_strategy", OptionType::STRING},
          {"hessian_model", OptionType::STRING},
          {"inertia_correction_strategy", OptionType::STRING},
-         {"scale_functions", OptionType::BOOL},
+         {"use_function_scaling", OptionType::BOOL},
          {"function_scaling_threshold", OptionType::DOUBLE},
          {"function_scaling_factor", OptionType::DOUBLE},
          {"scale_residuals", OptionType::BOOL},
diff --git a/uno/options/Presets.cpp b/uno/options/Presets.cpp
index 72316741e..4a2b8def0 100644
--- a/uno/options/Presets.cpp
+++ b/uno/options/Presets.cpp
@@ -58,7 +58,7 @@ namespace uno {
          preset_options.set_double("l1_constraint_violation_coefficient", 1000.);
          preset_options.set_string("progress_norm", "L1");
          preset_options.set_string("residual_norm", "INF");
-         preset_options.set_bool("scale_functions", true);
+         preset_options.set_bool("use_function_scaling", true);
          preset_options.set_double("primal_tolerance", 1e-8);
          preset_options.set_double("dual_tolerance", 1e-8);
          preset_options.set_double("loose_primal_tolerance", 1e-6);