diff --git a/CMakeLists.txt b/CMakeLists.txt
index e3e10d42e..02e4f3b62 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -90,6 +90,7 @@ file(GLOB UNO_SOURCE_FILES CONFIGURE_DEPENDS
    uno/ingredients/inertia_correction_strategies/*.cpp
    uno/ingredients/subproblem/*.cpp
    uno/ingredients/subproblem_solvers/*.cpp
+   uno/ingredients/subproblem_solvers/TRON/*.cpp
    uno/model/*.cpp
    uno/optimization/*.cpp
    uno/options/*.cpp
diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp
index c4532d583..a55ade23c 100644
--- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp
+++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp
@@ -28,6 +28,13 @@ namespace uno {
          }
       }
       // from now on, the problem has inequalities
+      /*
+      // bound-constrained problem: do not reformulate
+      if (problem.number_constraints == 0 && problem.has_bound_constraints()) {
+         return std::make_unique<NoInequalityReformulation>("bound-constrained method");
+      }
+      */
+
       const std::string inequality_handling_method = options.get_string("inequality_handling_method");
       // inequality-constrained methods
       if (inequality_handling_method == "inequality_constrained") {
diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp
index 03aeed312..ae0d31494 100644
--- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp
+++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp
@@ -11,6 +11,7 @@
 #include "InverseNewtonSolver.hpp"
 #include "QPSolverFactory.hpp"
 #include "WoodburyEQPSolver.hpp"
+#include "TRON/TRONSolver.hpp"
 #include "ingredients/subproblem/Subproblem.hpp"
 #include "options/Options.hpp"
 #include "tools/Logger.hpp"
@@ -52,6 +53,14 @@ namespace uno {
             return subproblem_solver;
          }
       }
+      /*
+      // if only bound constraints, allocate bound-constrained solver
+      else if (subproblem.number_constraints == 0 && subproblem.has_bound_constraints()) {
+         auto subproblem_solver = std::make_unique<TRONSolver>();
+         subproblem_solver->initialize_memory(subproblem);
+         return subproblem_solver;
+      }
+      */
       // if no inequality constraint and no trust region, allocate EQP solver
       else if (!subproblem.has_inequality_constraints() && !subproblem.has_bound_constraints() && !uses_trust_region) {
          if constexpr (std::is_same_v<HessianType, InverseLBFGSHessian>) { // unconstrained
diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp
new file mode 100644
index 000000000..8ac9c41b6
--- /dev/null
+++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp
@@ -0,0 +1,368 @@
+// Copyright (c) 2026 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <optional>
+#include <stdexcept>
+#include "TRONSolver.hpp"
+#include "ingredients/subproblem/Subproblem.hpp"
+#include "linear_algebra/BLAS.hpp"
+#include "linear_algebra/Vector.hpp"
+#include "symbolic/UnaryNegation.hpp"
+
+namespace uno {
+   void TRONSolver::initialize_memory(const Subproblem& subproblem) {
+      this->lower_bounds.resize(subproblem.number_variables);
+      this->upper_bounds.resize(subproblem.number_variables);
+      this->workspace.objective_gradient.resize(subproblem.number_variables);
+   }
+
+   void TRONSolver::solve(Statistics& /*statistics*/, const Subproblem& subproblem, double trust_region_radius,
+         const Vector<double>& /*initial_point*/, Direction& direction, Evaluations& current_evaluations,
+         const WarmstartInformation& /*warmstart_information*/) {
+      if (0 < subproblem.number_constraints) {
+         throw std::runtime_error("TRONSolver cannot solve problems with general constraints");
+      }
+      std::cout << "TRON solver created\n";
+      throw std::runtime_error("TRONSolver created");
+   }
+
+   SolverWorkspace& TRONSolver::get_workspace() {
+      return this->workspace;
+   }
+
+   void TRONSolver::solve(const ObjectiveOperator& objective_operator, const GradientOperator& gradient_operator,
+         const MatrixOperator& hessian_operator, const Vector<double>& initial_point) {
+      SolveStatus status = SolveStatus::Unknown;
+
+      // initialize x within bounds
+      this->x = initial_point;
+      project_onto_bounds(x);
+      gradient_operator(this->x, gx);
+
+      // projected-gradient norm at x0
+      project_step(x, -1., gx, this->gpx); // projected_gradient = P(x - g) - x
+      double pix = norm_2(this->gpx);
+      double epsilon = atol + rtol * pix;
+
+      // trust-region radius
+      double radius = std::min(std::max(1., pix / 10.), this->max_radius);
+
+      if (pix <= epsilon) {
+         status = SolveStatus::Optimal;
+         return;
+      }
+
+      // TODO change point at which Hessian operator is evaluated
+
+      // compute Cauchy step and store it in s
+      double alpha_c = 1.;
+      const bool cauchy_success = compute_cauchy_step(hessian_operator, this->gx, alpha_c, radius);
+      if (!cauchy_success) {
+         status = SolveStatus::Error;
+         return;
+      }
+
+      // projected Newton refinement (CG)
+      std::string cg_info = projected_newton(hessian_operator, gx, radius);
+   }
+
+   // protected member functions
+
+   // project v component-wise onto [ℓ, u]
+   void TRONSolver::project_onto_bounds(Vector<double>& v) const {
+      for (size_t i = 0; i < v.size(); ++i)
+         v[i] = std::max(this->lower_bounds[i], std::min(v[i], this->upper_bounds[i]));
+   }
+
+   // active-set indicator: active[i] = true if x[i] at a bound
+   void TRONSolver::compute_active_set(std::vector<bool>& active, const Vector<double>& x) const {
+      for (size_t i = 0; i < x.size(); ++i) {
+         const double delta = (-INF<double> < this->lower_bounds[i] && this->lower_bounds[i] < this->upper_bounds[i] &&
+            this->upper_bounds[i] < INF<double>) ? std::min(atol, rtol * (this->upper_bounds[i] - this->lower_bounds[i])) : atol;
+         if (x[i] == this->lower_bounds[i] && x[i] == this->upper_bounds[i]) {
+            active[i] = true;
+         }
+         else if (x[i] <= this->lower_bounds[i] + delta) {
+            active[i] = true;
+         }
+         else if (x[i] >= this->upper_bounds[i] - delta) {
+            active[i] = true;
+         }
+         else {
+            active[i] = false;
+         }
+      }
+   }
+
+   // s = P(x + alpha*d) - x
+   void TRONSolver::project_step(const Vector<double>& x, double alpha, const Vector<double>& d, Vector<double>& s) const {
+      for (size_t i = 0; i < x.size(); ++i) {
+         if (x[i] + alpha * d[i] < this->lower_bounds[i]) {
+            s[i] = this->lower_bounds[i] - x[i];
+         }
+         else if (this->upper_bounds[i] < x[i] + alpha * d[i]) {
+            s[i] = this->upper_bounds[i] - x[i];
+         }
+         else {
+            s[i] = alpha * d[i];
+         }
+      }
+   }
+
+   // Hs = H*s, slope = gᵀs, qs = ½sᵀHs + gᵀs
+   std::pair<double, double> TRONSolver::compute_Hs_slope_qs(const MatrixOperator& hessian_operator, const Vector<double>& s,
+         const Vector<double>& g) {
+      hessian_operator(s, Hs); // at xc_
+      double slope = dot(g, s);
+      double qs = dot(s, Hs) / 2. + slope;
+      return {slope, qs};
+   }
+
+   // compute the minimal and maximal break-points of the projection of x + alpha*d on the n-dimensional interval [xl,xu].
+   BreakPoints TRONSolver::compute_break_points(const Vector<double>& x, const Vector<double>& d) const {
+      BreakPoints break_points{INF<double>, 0.};
+
+      for (size_t i = 0; i < x.size(); ++i) {
+         std::optional<double> break_point = std::nullopt;
+         if (x[i] < this->upper_bounds[i] && d[i] > 0.) {
+            break_point = (this->upper_bounds[i] - x[i]) / d[i];
+         }
+         else if (x[i] > this->lower_bounds[i] && d[i] < 0.) {
+            break_point = (this->lower_bounds[i] - x[i]) / d[i];
+         }
+         if (break_point.has_value()) {
+            break_points.min = std::min(*break_point, break_points.min);
+            break_points.max = std::max(*break_point, break_points.max);
+         }
+      }
+      return break_points;
+   }
+
+   // backtracking projected line search: find smallest t = 2^{-k} s.t.  q(s) ≤ μ₀ gᵀs, where s = P(x + t*d) - x.
+   void TRONSolver::projected_line_search(const MatrixOperator& hessian_operator, Vector<double>& x, const Vector<double>& d,
+         const Vector<double>& g, Vector<double>& s) {
+      double alpha = 1.;
+      const BreakPoints break_points = compute_break_points(x, d);
+
+      bool search = true;
+      while (search && alpha > break_points.min) {
+         project_step(x, alpha, d, s);
+         const auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g);
+         if (qs <= mu0 * slope) {
+            search = false;
+         }
+         else {
+            alpha /= 2.;
+         }
+      }
+      if (alpha < std::min(1., break_points.min)) {
+         alpha = break_points.min;
+         project_step(x, alpha, d, s);
+         hessian_operator(s, this->Hs);
+      }
+      project_step(x, alpha, d, s);
+      x += s;
+   }
+
+   // compute Cauchy step s = P(x - α g) - x satisfying sufficient decrease.
+   // returns true upon success, false upon failure
+   bool TRONSolver::compute_cauchy_step(const MatrixOperator& hessian_operator, const Vector<double>& g, double& alpha,
+         double radius) {
+      // compute breakpoints along negative gradient direction
+      s = -g;
+      const BreakPoints break_points = compute_break_points(x, s);
+
+      s.fill(0.);
+      Hs.fill(0.);
+
+      project_step(x, -alpha, g, s);
+
+      // interpolate or extrapolate
+      bool interpolate = true;
+      if (norm_2(s) <= mu1 * radius) {
+         auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g);
+         interpolate = (qs >= mu0 * slope);
+      }
+
+      if (interpolate) {
+         bool search = true;
+         while (search) {
+            alpha /= sigma;
+            project_step(x, -alpha, g, s);
+            if (norm_2(s) <= mu1 * radius) {
+               auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g);
+               search = (qs >= mu0 * slope);
+            }
+            // TODO: correctly assess why this fails
+            if (alpha < std::sqrt(std::nextafter(0., 1.))) {
+               return false;
+            }
+         }
+      }
+      else { // extrapolation
+         double alpha_success = alpha;
+         bool search = true;
+         while (search && alpha <= break_points.max) {
+            alpha *= sigma;
+            project_step(x, -alpha, g, s);
+            if (norm_2(s) <= mu1 * radius) {
+               auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g);
+               if (qs <= mu0 * slope) {
+                  alpha_success = alpha;
+               }
+            }
+            else {
+               search = false;
+            }
+         }
+         // recover the last successful step
+         alpha = alpha_success;
+         project_step(x, -alpha, g, s);
+      }
+      return true;
+   }
+
+   // find t ≥ 0 so that ‖d + t*p‖ = Delta (solve quadratic for t ≥ 0)
+   double TRONSolver::compute_distance_to_trust_region(const Vector<double>& d, const Vector<double>& p, double radius) {
+      const double dTd = dot(d, d);
+      const double dTp = dot(d, p);
+      const double pTp = dot(p, p);
+      double discriminant = dTp * dTp - pTp * (dTd - radius * radius);
+      if (discriminant < 0.0) {
+         discriminant = 0.0;
+      }
+      return (-dTp + std::sqrt(discriminant)) / pTp;
+   }
+
+   // Conjugate Gradient (Steihaug-Toint style)
+   CGStatus TRONSolver::CG(Vector<double>& d, const MatrixOperator& matrix_operator, const Vector<double>& rhs, double radius,
+         double gfnorm_sqrt) {
+      const size_t n = d.size();
+      Vector<double> r(n), p(n), Hp(n);
+      d.fill(0.);
+      r = rhs;
+      p = r;
+      double norm_r = norm_2(r);
+
+      for (size_t cg_it = 0; cg_it < max_cgiter; ++cg_it) {
+         matrix_operator(p, Hp);
+         double pHp = dot(p, Hp);
+         if (pHp <= 0.0) {
+            // Negative curvature: go to boundary
+            const double t = compute_distance_to_trust_region(d, p, radius);
+            d += t*p;
+            return CGStatus::ON_TR_BOUNDARY;
+         }
+         double alpha_cg = norm_r / pHp;
+
+         // trial step
+         for (size_t i = 0; i < n; ++i) {
+            w[i] = d[i] + alpha_cg * p[i];
+         }
+         if (norm_2(w) >= radius) {
+            const double t = compute_distance_to_trust_region(d, p, radius);
+            d += t*p;
+            return CGStatus::ON_TR_BOUNDARY;
+         }
+         d = w;
+
+         // Update residual r = r - alpha * Hp
+         // TODO
+         // r -= alpha_cg*Hp;
+         double rr_new = dot(r, r);
+
+         if (std::sqrt(rr_new) <= cgtol * gfnorm_sqrt) {
+            return CGStatus::SUCCESS;
+         }
+         double beta = rr_new / norm_r;
+         // p = r + beta*p
+         for (size_t i = 0; i < n; ++i) {
+            p[i] = r[i] + beta * p[i];
+         }
+         norm_r = rr_new;
+      }
+      return CGStatus::MAX_ITERATIONS;
+   }
+
+   std::string TRONSolver::projected_newton(const MatrixOperator& hessian_operator, const Vector<double>& g, double radius) {
+      const size_t n = this->x.size();
+
+      // Update Hs = H * s
+      hessian_operator(s, Hs);
+
+      // projected Newton step
+      bool exit_optimal = false, exit_pcg = false, exit_itmax = false;
+      std::string exit_status = "maximum number of iterations";
+      size_t iters = 0;
+      x += s;
+      project_onto_bounds(x);
+      std::vector<bool> active(n, false);
+      while (!exit_optimal && !exit_pcg && !exit_itmax) {
+         compute_active_set(active, x);
+         // stop if all bounds are active
+         if (std::all_of(active.begin(), active.end(), [](bool active) { return active; })) {
+            exit_optimal = true;
+            continue;
+         }
+
+         // build RHS = -(g + Hs) for free variables
+         double gfnorm = 0.0;
+         for (size_t i = 0; i < n; ++i) {
+            this->quadratic_gradient[i] = active[i] ? 0. : -g[i];
+            gfnorm += this->quadratic_gradient[i] * this->quadratic_gradient[i];
+            this->quadratic_gradient[i] -= active[i] ? 0. : Hs[i];
+         }
+         double gfnorm_sqrt = std::sqrt(gfnorm);
+
+         // define ZHZ, the Hessian-vector product with free-variable masking
+         // (Hp)_i = ifix[i] ? 0 : (H * (ifix-masked d))_i
+         const auto ZHZ = [&](const Vector<double>& d, Vector<double>& result) {
+            // zero out fixed components of d, then apply H
+            for (size_t i = 0; i < n; ++i) {
+               this->d_masked[i] = active[i] ? 0. : d[i];
+            }
+            hessian_operator(this->d_masked, result);
+            // zero out fixed components of the result
+            for (size_t i = 0; i < n; ++i) {
+               if (active[i]) {
+                  result[i] = 0.;
+               }
+            }
+         };
+
+         const CGStatus cg_status = CG(this->d_cg, ZHZ, this->quadratic_gradient, radius, gfnorm_sqrt);
+         iters++;
+
+         // projected line search
+         this->quadratic_gradient.scale(-1.);
+         projected_line_search(ZHZ, x, d_cg, this->quadratic_gradient, w);
+         s += w;
+         hessian_operator(s, Hs);
+
+         // Check optimality: ‖(g + Hs) restricted‖ ≤ cgtol * gfnorm_sqrt
+         double new_norm = 0.0;
+         for (size_t i = 0; i < n; ++i) {
+            const double ri = active[i] ? 0. : Hs[i] + g[i];
+            new_norm += ri * ri;
+         }
+         if (std::sqrt(new_norm) <= cgtol * gfnorm_sqrt) {
+            exit_optimal = true;
+         }
+         else if (cg_status == CGStatus::ON_TR_BOUNDARY) {
+            exit_pcg = true;
+         }
+         else if (iters >= max_cgiter) {
+            exit_itmax = true;
+         }
+      }
+
+      if (exit_optimal) return "stationary point found";
+      if (exit_pcg) return "on trust-region boundary";
+      if (exit_itmax) return "maximum number of iterations";
+      return exit_status;
+   }
+} // namespace
\ No newline at end of file
diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp
new file mode 100644
index 000000000..408746848
--- /dev/null
+++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp
@@ -0,0 +1,125 @@
+// Copyright (c) 2026 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_TRONSOLVER_H
+#define UNO_TRONSOLVER_H
+
+#include <cmath>
+#include <functional>
+#include <vector>
+#include "../SubproblemSolver.hpp"
+#include "../SolverWorkspace.hpp"
+#include "linear_algebra/Vector.hpp"
+
+namespace uno {
+   class TRONSolverWorkspace: public SolverWorkspace {
+   public:
+      TRONSolverWorkspace() = default;
+
+      [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, const Vector<double>& /*vector*/) const override {
+         return 0.;
+      }
+
+      Vector<double> objective_gradient;
+   };
+
+   enum class SolveStatus {
+      Unknown,
+      Optimal,
+      MaxIter,
+      SmallStep,
+      NegPred,
+      Error
+   };
+
+   enum class CGStatus {
+      SUCCESS,
+      ON_TR_BOUNDARY,
+      MAX_ITERATIONS
+   };
+
+   struct TRONStats {
+      SolveStatus status = SolveStatus::Unknown;
+      int iter = 0;
+      double objective = 0.0;
+      double dual_residual = 0.0; ///< ‖P(x - g) - x‖ (projected gradient)
+      std::vector<double> solution;
+   };
+
+   struct BreakPoints {
+      double min; // minimal break-point
+      double max; // maximal break-point
+   };
+
+   class TRONSolver: public SubproblemSolver {
+   public:
+      using ObjectiveOperator = std::function<double(const Vector<double>& x)>;
+      using GradientOperator = std::function<void(const Vector<double>& x, Vector<double>& gradient)>;
+      using MatrixOperator = std::function<void(const Vector<double>& x, Vector<double>& Hv)>;
+
+      TRONSolver() = default;
+      ~TRONSolver() override = default;
+
+      void initialize_memory(const Subproblem& subproblem) override;
+
+      void solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, const Vector<double>& initial_point,
+         Direction& direction, Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) override;
+      /**
+       * @param x   Initial guess (length n); need not be feasible.
+       * @param opts Solver options.
+       * @return     Execution statistics including the solution vector.
+       */
+      void solve(const ObjectiveOperator& objective_operator, const GradientOperator& gradient_operator,
+         const MatrixOperator& hessian_operator, const Vector<double>& initial_point);
+
+      [[nodiscard]] SolverWorkspace& get_workspace() override;
+
+   protected:
+      std::vector<double> lower_bounds;
+      std::vector<double> upper_bounds;
+      TRONSolverWorkspace workspace{};
+
+      // Workspace vectors (all length n_)
+      Vector<double> x, xc, x_copy, gpx, gx, s, Hs, temp_, w;
+      Vector<double> d_masked;
+      Vector<double> quadratic_gradient;
+      Vector<double> d_cg; // CG point
+      double fc;
+
+      // parameters
+      double mu0 = 1.0 / 100.0; ///< sufficient-decrease parameter  ∈ (0, 0.5)
+      double mu1 = 1.0; ///< trust-region scaling            ∈ (0, ∞)
+      double sigma = 10.0; ///< step-size update factor         ∈ (1, ∞)
+      // trust-region constants
+      double eta1_ = 0.1; // acceptance threshold (ratio)
+      double eta2_ = 0.75; // "very successful" threshold
+      double min_radius_ = 1e-10;
+      double max_radius = std::min(1.0 / std::sqrt(2.0 * std::numeric_limits<double>::epsilon()), 100.0);
+      // options
+      size_t max_iter = 100000;
+      size_t max_cgiter = 50;
+      double cgtol = 0.1; ///< CG sub-problem tolerance
+      const double atol = std::sqrt(std::numeric_limits<double>::epsilon());
+      const double rtol = std::sqrt(std::numeric_limits<double>::epsilon());
+
+      std::function<std::pair<double, Vector<double>> (const Vector<double>& x)> fg_;
+
+      void project_onto_bounds(Vector<double>& v) const;
+      void compute_active_set(std::vector<bool>& active, const Vector<double>& x) const;
+      void project_step(const Vector<double>& x, double alpha, const Vector<double>& d, Vector<double>& s) const;
+      [[nodiscard]] std::pair<double, double> compute_Hs_slope_qs(const MatrixOperator& hessian_operator,
+         const Vector<double>& s, const Vector<double>& g);
+      [[nodiscard]] BreakPoints compute_break_points(const Vector<double>& x, const Vector<double>& d) const;
+      void projected_line_search(const MatrixOperator& hessian_operator, Vector<double>& x, const Vector<double>& d,
+         const Vector<double>& g, Vector<double>& s);
+      [[nodiscard]] bool compute_cauchy_step(const MatrixOperator& hessian_operator, const Vector<double>& g, double& alpha,
+         double radius);
+      [[nodiscard]] static double compute_distance_to_trust_region(const Vector<double>& d, const Vector<double>& p,
+         double radius);
+      [[nodiscard]] CGStatus CG(Vector<double>& d, const MatrixOperator& matrix_operator, const Vector<double>& rhs,
+         double radius, double gfnorm_sqrt);
+      std::string projected_newton(const MatrixOperator& hessian_operator, const Vector<double>& g, double radius);
+   };
+} // namespace
+
+#endif // UNO_TRONSOLVER_H
\ No newline at end of file
diff --git a/uno/symbolic/UnaryNegation.hpp b/uno/symbolic/UnaryNegation.hpp
index 09994fde6..7f6baba05 100644
--- a/uno/symbolic/UnaryNegation.hpp
+++ b/uno/symbolic/UnaryNegation.hpp
@@ -26,6 +26,8 @@ namespace uno {
          return -this->expression[index];
       }
 
+		UNO_FORWARD_ACCESSOR(get_expression, this->expression)
+
    protected:
       storage_t<Expression> expression;
    };