diff --git a/bluemira/optimisation/_nlopt/functions.py b/bluemira/optimisation/_nlopt/functions.py
index 4b99f4083f..f858b482bd 100644
--- a/bluemira/optimisation/_nlopt/functions.py
+++ b/bluemira/optimisation/_nlopt/functions.py
@@ -123,6 +123,7 @@ def __init__(
         self.constraint_type = constraint_type
         self.tolerance = tolerance
         self.df = df if df is not None else self._approx_derivative
+        self.history: list[tuple[np.ndarray, ...]] = []
 
     def call(self, result: np.ndarray, x: np.ndarray, grad: np.ndarray) -> None:
         """
@@ -136,3 +137,9 @@ def call(self, result: np.ndarray, x: np.ndarray, grad: np.ndarray) -> None:
         self.f0 = result
         if grad.size > 0:
             grad[:] = self.df(x)
+
+    def call_with_history(
+        self, result: np.ndarray, x: np.ndarray, grad: np.ndarray
+    ) -> None:
+        self.call(result, x, grad)
+        self.history.append((np.copy(x), np.copy(result)))
diff --git a/bluemira/optimisation/_nlopt/optimiser.py b/bluemira/optimisation/_nlopt/optimiser.py
index 97dc647aed..ad2c39fe2d 100644
--- a/bluemira/optimisation/_nlopt/optimiser.py
+++ b/bluemira/optimisation/_nlopt/optimiser.py
@@ -169,7 +169,10 @@ def add_eq_constraint(
             df_constraint,
             bounds=(self.lower_bounds, self.upper_bounds),
         )
-        self._opt.add_equality_mconstraint(constraint.call, constraint.tolerance)
+        self._opt.add_equality_mconstraint(
+            constraint.call_with_history if self._keep_history else constraint.call,
+            constraint.tolerance,
+        )
         self._eq_constraints.append(constraint)
 
     def add_ineq_constraint(
@@ -200,7 +203,10 @@ def add_ineq_constraint(
             df_constraint,
             bounds=(self.lower_bounds, self.upper_bounds),
         )
-        self._opt.add_inequality_mconstraint(constraint.call, constraint.tolerance)
+        self._opt.add_inequality_mconstraint(
+            constraint.call_with_history if self._keep_history else constraint.call,
+            constraint.tolerance,
+        )
         self._ineq_constraints.append(constraint)
 
     def optimise(self, x0: np.ndarray | None = None) -> OptimiserResult:
@@ -247,8 +253,15 @@ def optimise(self, x0: np.ndarray | None = None) -> OptimiserResult:
             x=x_star,
             n_evals=self._opt.get_numevals(),
             history=self._objective.history,
+            constraint_history=self._get_constraint_history(),
         )
 
+    def _get_constraint_history(self):
+        return [
+            constraint.history
+            for constraint in self._eq_constraints + self._ineq_constraints
+        ]
+
     def set_lower_bounds(self, bounds: npt.ArrayLike) -> None:
         """
         Set the lower bound for each optimisation parameter.
diff --git a/bluemira/optimisation/_optimiser.py b/bluemira/optimisation/_optimiser.py
index 86814e893f..a5349711b3 100644
--- a/bluemira/optimisation/_optimiser.py
+++ b/bluemira/optimisation/_optimiser.py
@@ -30,6 +30,10 @@ class OptimiserResult:
     The first element of each tuple is the parameterisation (x), the
     second is the evaluation of the objective function at x (f(x)).
     """
+    constraint_history: list[tuple[np.ndarray, ...]] = field(repr=False)
+    """
+    Constraint history
+    """
     constraints_satisfied: bool | None = None
     """
     Whether all constraints have been satisfied to within the required tolerance.