diff --git a/docs/conf.py b/docs/conf.py
index ca479a2..d11eca2 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -14,7 +14,7 @@
 project = "pygradflow"
 copyright = "2023, Christoph Hansknecht"
 author = "Christoph Hansknecht"
-release = "0.5.18"
+release = "0.5.19"
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
diff --git a/pygradflow/implicit_func.py b/pygradflow/implicit_func.py
index 912c3a4..97e3454 100644
--- a/pygradflow/implicit_func.py
+++ b/pygradflow/implicit_func.py
@@ -100,6 +100,21 @@ def apply_project_deriv(self, mat: sp.sparse.spmatrix, active_set: np.ndarray):
 
 
 class ImplicitFunc(StepFunc):
+    """
+    Standard residual function for implicit Euler steps
+    with the value
+
+    .. math::
+        F(x, y; \\hat{x}, \\hat{y}) =
+        \\begin{pmatrix}
+            x - P_{C} (\\hat{x} - \\Delta t \\nabla_{x} \\mathcal{L}^{\\rho}(x, y)) \\\\
+            y - (\\hat{y} + \\Delta t \\nabla_{y} \\mathcal{L}^{\\rho}(x, y)
+        \\end{pmatrix}
+
+    The values :math:`\\hat{x}` and :math:`\\hat{y}` are obtained
+    from the iteratre provided in the constructor
+    """
+
     def __init__(self, problem: Problem, iterate: Iterate, dt: float) -> None:
         super().__init__(problem, iterate, dt)
 
@@ -185,6 +200,14 @@ def deriv_at(
 
 
 class ScaledImplicitFunc(StepFunc):
+    """
+    The *scaled* residual function for implicit Euler steps
+    is defined as the normal residual function
+    given by :py:class:`pygradflow.implicit_func.ImplicitFunc`
+    scaled by a factor
+    of :math:`\\lambda = 1 / \\Delta t`.
+    """
+
     def __init__(self, problem: Problem, iterate: Iterate, dt: float) -> None:
         super().__init__(problem, iterate, dt)
         self.lamb = 1.0 / dt
diff --git a/pygradflow/newton.py b/pygradflow/newton.py
index 7b6ff9f..12d5009 100644
--- a/pygradflow/newton.py
+++ b/pygradflow/newton.py
@@ -200,11 +200,17 @@ def __init__(
 
         self.step_solver = step_solver
         self.step_solver.update_derivs(orig_iterate)
+        self._curr_active_set = None
 
     def step(self, iterate):
         active_set = self.func.compute_active_set(iterate, self.rho, self.tau)
 
-        self.step_solver.update_active_set(active_set)
+        if self._curr_active_set is None:
+            self.step_solver.update_active_set(active_set)
+        elif (self._curr_active_set != active_set).any():
+            self.step_solver.update_active_set(active_set)
+
+        self._curr_active_set = active_set
 
         return self.step_solver.solve(iterate)
 
@@ -247,6 +253,9 @@ def step(self, iterate):
         func_value = self.func.value_at(iterate, self.rho)
         res_value = 0.5 * np.dot(func_value, func_value)
 
+        if res_value <= params.newton_tol:
+            return step_result
+
         # TODO: Create a method in the step function
         # to compute a forward product instead to make
         # everything more efficient
@@ -271,7 +280,7 @@ def step(self, iterate):
             next_func_value = self.func.value_at(next_iterate, self.rho)
             next_res_value = 0.5 * np.dot(next_func_value, next_func_value)
 
-            if np.isclose(next_res_value, 0.0):
+            if next_res_value <= params.newton_tol:
                 break
 
             if next_res_value <= res_value + (1e-4 * alpha * inner_product):
diff --git a/pygradflow/result.py b/pygradflow/result.py
index e3aca51..07666e2 100644
--- a/pygradflow/result.py
+++ b/pygradflow/result.py
@@ -23,7 +23,8 @@ def __init__(
         dist_factor: float,
         **attrs
     ):
-        self.problem = problem
+        self.num_vars = problem.num_vars
+        self.num_cons = problem.num_cons
         self._attrs = attrs
 
         self._x = x
@@ -39,8 +40,8 @@ def _set_path(self, path, model_times):
         self._attrs["path"] = path
         self._attrs["model_times"] = model_times
 
-        num_vars = self.problem.num_vars
-        num_cons = self.problem.num_cons
+        num_vars = self.num_vars
+        num_cons = self.num_cons
 
         assert model_times.ndim == 1
         assert path.shape == (num_vars + num_cons, len(model_times))
diff --git a/pygradflow/runners/instance.py b/pygradflow/runners/instance.py
index 12e4b32..4bb138c 100644
--- a/pygradflow/runners/instance.py
+++ b/pygradflow/runners/instance.py
@@ -9,6 +9,9 @@ def __init__(self, name, num_vars, num_cons):
         self.num_vars = num_vars
         self.num_cons = num_cons
 
+    def __repr__(self):
+        return f"{self.__class__.__name__}({self.name})"
+
     @property
     def size(self):
         return self.num_vars + self.num_cons
diff --git a/pygradflow/step/box_control.py b/pygradflow/step/box_control.py
index 857c521..8cffe80 100644
--- a/pygradflow/step/box_control.py
+++ b/pygradflow/step/box_control.py
@@ -216,9 +216,10 @@ def solve_step_ipopt(self, iterate, rho, dt, timer):
         return reduced_problem.solve(timer=timer)
 
     def solve_step_box(self, iterate, rho, dt, timer):
-        from .box_solver import BoxSolverError, solve_box_constrained
+        from .box_solver import BoxSolverError, BoxSolverStatus, solve_box_constrained
 
         problem = self.problem
+        params = self.params
         lamb = 1.0 / dt
 
         def objective(x):
@@ -231,9 +232,28 @@ def hessian(x):
             return self.hessian(iterate, x, lamb, rho)
 
         try:
-            return solve_box_constrained(
-                iterate.x, objective, gradient, hessian, problem.var_lb, problem.var_ub
+            result = solve_box_constrained(
+                iterate.x,
+                objective,
+                gradient,
+                hessian,
+                problem.var_lb,
+                problem.var_ub,
+                obj_lower=params.obj_lower_limit,
             )
+
+            if result.success:
+                return result.x
+
+            elif result.status == BoxSolverStatus.Unbounded:
+                return result.x
+            else:
+                raise StepSolverError(
+                    "Box-constrained solver failed to converge (status = {})".format(
+                        result.status
+                    )
+                )
+
         except BoxSolverError as e:
             raise StepSolverError("Box-constrained solver failed to converge") from e
 
diff --git a/pygradflow/step/box_solver.py b/pygradflow/step/box_solver.py
index 42484bc..84d2752 100644
--- a/pygradflow/step/box_solver.py
+++ b/pygradflow/step/box_solver.py
@@ -1,3 +1,4 @@
+from enum import Enum, auto
 from typing import Callable
 
 import numpy as np
@@ -10,6 +11,22 @@ class BoxSolverError(Exception):
     pass
 
 
+class BoxSolverStatus(Enum):
+    Optimal = auto()
+    Unbounded = auto()
+    IterationLimit = auto()
+
+
+class BoxSolverResult:
+    def __init__(self, x, status, iterations):
+        self.x = x
+        self.status = status
+
+    @property
+    def success(self):
+        return self.status == BoxSolverStatus.Optimal
+
+
 # Based on "Projected Newton Methods for Optimization Problems with Simple Constraints"
 def solve_box_constrained(
     x0: np.ndarray,
@@ -18,10 +35,11 @@ def solve_box_constrained(
     hess: Callable[[np.ndarray], sp.sparse.spmatrix],
     lb: np.ndarray,
     ub: np.ndarray,
+    obj_lower: float,
     max_it=1000,
-    atol: float = 1e-8,
-    rtol: float = 1e-8,
-):
+    atol: float = 1e-6,
+    rtol: float = 1e-6,
+) -> BoxSolverResult:
 
     (n,) = x0.shape
     assert lb.shape == (n,)
@@ -32,10 +50,16 @@ def solve_box_constrained(
     beta = 0.5
     sigma = 1e-3
 
+    status = BoxSolverStatus.IterationLimit
+
     for iteration in range(max_it):
         curr_func = func(curr_x)
         curr_grad = grad(curr_x)
 
+        if curr_func <= obj_lower:
+            status = BoxSolverStatus.Unbounded
+            break
+
         assert curr_grad.shape == (n,)
 
         at_lower = np.isclose(curr_x, lb)
@@ -52,9 +76,11 @@ def solve_box_constrained(
         grad_norm = np.linalg.norm(curr_grad, ord=np.inf)
 
         if grad_norm < atol:
+            status = BoxSolverStatus.Optimal
             break
 
         if (residuum < atol) or (residuum / grad_norm) < rtol:
+            status = BoxSolverStatus.Optimal
             break
 
         active = np.logical_or(active_lower, active_upper)
@@ -105,4 +131,4 @@ def solve_box_constrained(
     else:
         raise BoxSolverError(f"Did not converge after {max_it} iterations")
 
-    return curr_x
+    return BoxSolverResult(curr_x, status, iteration)
diff --git a/pyproject.toml b/pyproject.toml
index cf78fbb..8b52a91 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pygradflow"
-version = "0.5.18"
+version = "0.5.19"
 description = "PyGradFlow is a simple implementation of the sequential homotopy method to be used to solve general nonlinear programs."
 authors = ["Christoph Hansknecht <christoph.hansknecht@tu-clausthal.de>"]
 readme = "README.md"
diff --git a/tests/pygradflow/qp.py b/tests/pygradflow/qp.py
new file mode 100644
index 0000000..3fb062a
--- /dev/null
+++ b/tests/pygradflow/qp.py
@@ -0,0 +1,30 @@
+from pygradflow.problem import Problem
+
+
+class QP(Problem):
+    def __init__(self, hess, jac, grad, res, lb, ub):
+        self.m, self.n = jac.shape
+        assert hess.shape[0] == self.n
+        assert hess.shape[1] == self.n
+        assert grad.shape[0] == self.n
+        assert res.shape[0] == self.m
+        super().__init__(lb, ub, num_cons=self.m)
+        self.A = hess
+        self.B = jac
+        self.a = grad
+        self.b = res
+
+    def obj(self, x):
+        return 0.5 * x.T @ self.A @ x + self.a.T @ x
+
+    def obj_grad(self, x):
+        return self.A @ x + self.a
+
+    def cons(self, x):
+        return self.B @ x + self.b
+
+    def cons_jac(self, x):
+        return self.B
+
+    def lag_hess(self, x, _):
+        return self.A
diff --git a/tests/pygradflow/test_qp.py b/tests/pygradflow/test_qp.py
new file mode 100644
index 0000000..fc22f12
--- /dev/null
+++ b/tests/pygradflow/test_qp.py
@@ -0,0 +1,129 @@
+import numpy as np
+import pytest
+import scipy as sp
+
+from pygradflow.params import NewtonType, Params, StepControlType
+from pygradflow.solver import Solver
+from pygradflow.status import SolverStatus
+
+from .qp import QP
+
+
+@pytest.fixture
+def unbounded_qp():
+    n = 199
+    h = 1 / n
+    e = np.ones(n)
+    H = (1 / h**2) * sp.sparse.spdiags([e, -2 * e, e], [-1, 0, 1], n, n, format="coo")
+    A = sp.sparse.coo_matrix((0, n))
+    g = -e
+    b = np.array([])
+    lb = -np.inf * e
+    lb[n // 4] = 0.0
+    lb[3 * n // 4] = 0.0
+    lb[n // 2] = 0.0
+    ub = np.inf * e
+    return QP(H, sp.sparse.coo_matrix((0, n)), g, np.array([]), lb=lb, ub=ub)
+
+
+@pytest.fixture
+def boxed_qp():
+    n = 49
+    h = 1 / n
+    e = np.ones(n)
+    H = (1 / h**2) * sp.sparse.spdiags([-e, 2 * e, -e], [-1, 0, 1], n, n, format="coo")
+    g = e
+
+    lb = np.linspace(0, -0.01, n + 2)[1:-1]
+    lb[n // 4] = 0.0
+    lb[3 * n // 4] = 0.0
+    lb[n // 2] = 0.0
+    ub = np.inf * e
+    return QP(H, sp.sparse.coo_matrix((0, n)), g, np.array([]), lb=lb, ub=ub)
+
+
+@pytest.fixture
+def unbounded_x0():
+    return 0.0
+
+
+@pytest.mark.parametrize(
+    "step_control_type",
+    [
+        StepControlType.Exact,
+        StepControlType.Optimizing,
+        StepControlType.BoxReduced,
+        StepControlType.ResiduumRatio,
+        # StepControlType.DistanceRatio,
+    ],
+)
+def test_unbounded(unbounded_qp, unbounded_x0, step_control_type):
+
+    problem = unbounded_qp
+
+    params = Params(
+        step_control_type=step_control_type, iteration_limit=1000, display_interval=0.0
+    )
+
+    solver = Solver(problem, params)
+
+    result = solver.solve(x0=np.amax(problem.var_lb, 0), y0=np.array([]))
+
+    assert result.status == SolverStatus.Unbounded
+
+
+def test_boxed(boxed_qp, unbounded_x0):
+    problem = boxed_qp
+
+    params = Params(
+        iteration_limit=1000,
+        display_interval=0.0,
+        # newton_type=NewtonType.Full,
+        # lamb_init=1e-12,
+    )
+
+    solver = Solver(problem, params)
+
+    result = solver.solve(x0=0.0, y0=np.array([]))
+
+    assert result.status == SolverStatus.Optimal
+
+
+@pytest.mark.parametrize(
+    "newton_type",
+    [
+        NewtonType.ActiveSet,
+        NewtonType.Full,
+        # NewtonType.Globalized,
+        NewtonType.Simplified,
+    ],
+)
+def test_newton_types(newton_type):
+    n = 49
+    h = 1 / n
+    e = np.ones(n)
+    H = (1 / h**2) * sp.sparse.spdiags([-e, 2 * e, -e], [-1, 0, 1], n, n, format="coo")
+    A = sp.sparse.coo_matrix((0, n))
+    g = e
+    b = np.array([])
+    lb = -np.inf * e
+    ub = np.inf * e
+    qp = QP(H, A, g, b, lb, ub)
+
+    lb = np.linspace(0, -0.01, n + 2)[1:-1]
+    lb[n // 4] = 0.0
+    lb[3 * n // 4] = 0.0
+    lb[n // 2] = 0.0
+    ub = np.inf * e
+    qp = QP(H, sp.sparse.coo_matrix((0, n)), g, np.array([]), lb=lb, ub=ub)
+
+    params = Params()
+    params.display_interval = 1e-16
+    params.lamb_init = 1e-12
+    params.iteration_limit = 1000
+    params.newton_type = newton_type
+
+    solver = Solver(problem=qp, params=params)
+    result = solver.solve(x0=np.amax(lb, 0), y0=np.array([]))
+
+    assert result.success