diff --git a/bofire/data_models/constraints/interpoint.py b/bofire/data_models/constraints/interpoint.py
index ec03e0b8f..f7b1f09f3 100644
--- a/bofire/data_models/constraints/interpoint.py
+++ b/bofire/data_models/constraints/interpoint.py
@@ -1,5 +1,5 @@
 import math
-from typing import Annotated, Literal, Optional
+from typing import Annotated, Dict, Literal, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -83,3 +83,6 @@ def __call__(self, experiments: pd.DataFrame) -> pd.Series:
 
     def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
         raise NotImplementedError("Method `jacobian` currently not implemented.")
+
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        raise NotImplementedError("Method `hessian` currently not implemented.")
diff --git a/bofire/data_models/constraints/linear.py b/bofire/data_models/constraints/linear.py
index 264c9fe3d..9c94fad35 100644
--- a/bofire/data_models/constraints/linear.py
+++ b/bofire/data_models/constraints/linear.py
@@ -1,4 +1,4 @@
-from typing import Annotated, List, Literal, Tuple
+from typing import Annotated, Dict, List, Literal, Tuple, Union
 
 import numpy as np
 import pandas as pd
@@ -64,6 +64,9 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             columns=[f"dg/d{name}" for name in self.features],
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[int, str], float]:
+        return {idx: 0.0 for idx in range(experiments.shape[0])}
+
 
 class LinearEqualityConstraint(LinearConstraint, EqualityConstraint):
     """Linear equality constraint of the form `coefficients * x = rhs`.
diff --git a/bofire/data_models/constraints/nchoosek.py b/bofire/data_models/constraints/nchoosek.py
index d49caaf13..d6f2d85ca 100644
--- a/bofire/data_models/constraints/nchoosek.py
+++ b/bofire/data_models/constraints/nchoosek.py
@@ -1,4 +1,4 @@
-from typing import Literal
+from typing import Dict, Literal, Union
 
 import numpy as np
 import pandas as pd
@@ -119,3 +119,6 @@ def is_fulfilled(self, experiments: pd.DataFrame, tol: float = 1e-6) -> pd.Serie
 
     def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
         raise NotImplementedError("Jacobian not implemented for NChooseK constraints.")
+
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        raise NotImplementedError("Hessian not implemented for NChooseK constraints.")
diff --git a/bofire/data_models/constraints/nonlinear.py b/bofire/data_models/constraints/nonlinear.py
index a884bbed5..f868d650a 100644
--- a/bofire/data_models/constraints/nonlinear.py
+++ b/bofire/data_models/constraints/nonlinear.py
@@ -1,6 +1,6 @@
 import inspect
 import warnings
-from typing import Callable, Literal, Optional, Union
+from typing import Callable, Dict, Literal, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -9,6 +9,7 @@
 
 try:
     import torch
+    from torch.autograd.functional import hessian as torch_hessian
     from torch.autograd.functional import jacobian as torch_jacobian
 
     torch_tensor = torch.tensor
@@ -21,6 +22,7 @@ def error_func(*args, **kwargs):
     torch_jacobian = error_func
     torch_tensor = error_func
     torch_diag = error_func
+    torch_hessian = error_func
 
 from bofire.data_models.constraints.constraint import (
     EqualityConstraint,
@@ -47,6 +49,9 @@ class NonlinearConstraint(IntrapointConstraint):
     jacobian_expression: Optional[Union[str, Callable]] = Field(
         default=None, validate_default=True
     )
+    hessian_expression: Optional[Union[str, Callable]] = Field(
+        default=None, validate_default=True
+    )
 
     def validate_inputs(self, inputs: Inputs):
         if self.features is not None:
@@ -100,12 +105,53 @@ def set_jacobian_expression(cls, jacobian_expression, info) -> Union[str, Callab
 
         return jacobian_expression
 
+    @field_validator("hessian_expression")
+    @classmethod
+    def set_hessian_expression(cls, hessian_expression, info) -> Union[str, Callable]:
+        if (
+            hessian_expression is None
+            and "features" in info.data.keys()
+            and "expression" in info.data.keys()
+        ):
+            try:
+                import sympy
+            except ImportError as e:
+                warnings.warn(e.msg)
+                warnings.warn("please run `pip install sympy` for this functionality.")
+                return hessian_expression
+            if info.data["features"] is not None:
+                if isinstance(info.data["expression"], str):
+                    return (
+                        "["
+                        + ", ".join(
+                            [
+                                "["
+                                + ", ".join(
+                                    [
+                                        str(
+                                            sympy.S(info.data["expression"])
+                                            .diff(key1)
+                                            .diff(key2)
+                                        )
+                                        for key1 in info.data["features"]
+                                    ]
+                                )
+                                + "]"
+                                for key2 in info.data["features"]
+                            ]
+                        )
+                        + "]"
+                    )
+
+        return hessian_expression
+
     def __call__(self, experiments: pd.DataFrame) -> pd.Series:
         if isinstance(self.expression, str):
             return experiments.eval(self.expression)
         elif isinstance(self.expression, Callable):
             func_input = {
-                col: torch_tensor(experiments[col]) for col in experiments.columns
+                col: torch_tensor(experiments[col], requires_grad=False)
+                for col in experiments.columns
             }
             return pd.Series(self.expression(**func_input).cpu().numpy())
 
@@ -129,7 +175,10 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             elif isinstance(self.jacobian_expression, Callable):
                 args = inspect.getfullargspec(self.jacobian_expression).args
 
-                func_input = {arg: torch_tensor(experiments[arg]) for arg in args}
+                func_input = {
+                    arg: torch_tensor(experiments[arg], requires_grad=False)
+                    for arg in args
+                }
                 result = self.jacobian_expression(**func_input)
 
                 return pd.DataFrame(
@@ -141,12 +190,14 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
                             for col in experiments.columns
                         ]
                     ),
-                    index=["dg/d" + name for name in args],
+                    index=["dg/d" + name for name in experiments.columns],
                 ).transpose()
         elif isinstance(self.expression, Callable):
             args = inspect.getfullargspec(self.expression).args
 
-            func_input = tuple([torch_tensor(experiments[arg]) for arg in args])
+            func_input = tuple(
+                [torch_tensor(experiments[arg], requires_grad=False) for arg in args]
+            )
 
             result = torch_jacobian(self.expression, func_input)
             result = [torch_diag(result[i]).cpu().numpy() for i in range(len(args))]
@@ -160,6 +211,77 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             "The jacobian of a nonlinear constraint cannot be evaluated if jacobian_expression is None and expression is not Callable.",
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        """
+        Computes a dict of dataframes where the key dimension is the index of the experiments dataframe
+        and the value is the hessian matrix of the constraint evaluated at the corresponding experiment.
+
+        Args:
+            experiments (pd.DataFrame): Dataframe to evaluate the constraint Hessian on.
+
+        Returns:
+            Dict[pd.DataFrame]: Dictionary of dataframes where the key is the index of the experiments dataframe
+            and the value is the Hessian matrix of the constraint evaluated at the corresponding experiment.
+        """
+        if self.hessian_expression is not None:
+            if isinstance(self.hessian_expression, str):
+                res = experiments.eval(self.hessian_expression)
+            else:
+                if not isinstance(self.hessian_expression, Callable):
+                    raise ValueError(
+                        "The hessian_expression must be a string or a callable.",
+                    )
+                args = inspect.getfullargspec(self.hessian_expression).args
+
+                func_input = {
+                    arg: torch_tensor(experiments[arg], requires_grad=False)
+                    for arg in args
+                }
+                res = self.hessian_expression(**func_input)
+            for i, _ in enumerate(res):
+                for j, entry in enumerate(res[i]):
+                    if not hasattr(entry, "__iter__"):
+                        res[i][j] = pd.Series(np.repeat(entry, experiments.shape[0]))
+            res = np.array(res)
+            names = self.features or [f"x{i}" for i in range(experiments.shape[1])]
+
+            return {
+                idx: pd.DataFrame(
+                    res[..., i],
+                    columns=[f"d/d{name}" for name in names],
+                    index=[f"d/d{name}" for name in names],
+                )
+                for i, idx in enumerate(experiments.index)
+            }
+
+        elif isinstance(self.expression, Callable):
+            args = inspect.getfullargspec(self.expression).args
+
+            func_input = {
+                idx: tuple(
+                    [
+                        torch_tensor(experiments[arg][idx], requires_grad=False)
+                        for arg in args
+                    ]
+                )
+                for idx in experiments.index
+            }
+
+            names = self.features or [f"x{i}" for i in range(experiments.shape[1])]
+            res = {
+                idx: pd.DataFrame(
+                    np.array(torch_hessian(self.expression, func_input[idx])),
+                    columns=[f"d/d{name}" for name in names],
+                    index=[f"d/d{name}" for name in names],
+                )
+                for idx in experiments.index
+            }
+            return res
+
+        raise ValueError(
+            "The hessian of a nonlinear constraint cannot be evaluated if hessian_expression is None and expression is not Callable.",
+        )
+
 
 class NonlinearEqualityConstraint(NonlinearConstraint, EqualityConstraint):
     """Nonlinear equality constraint of the form 'expression == 0'.
diff --git a/bofire/data_models/constraints/product.py b/bofire/data_models/constraints/product.py
index 85b71d60a..f274f9e4a 100644
--- a/bofire/data_models/constraints/product.py
+++ b/bofire/data_models/constraints/product.py
@@ -83,6 +83,11 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             "Jacobian for product constraints is not yet implemented.",
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> List[pd.DataFrame]:
+        raise NotImplementedError(
+            "Hessian for product constraints is not yet implemented.",
+        )
+
 
 class ProductEqualityConstraint(ProductConstraint, EqualityConstraint):
     """Represents a product constraint of the form `sign * x1**e1 * x2**e2 * ... * xn**en == rhs`.
diff --git a/bofire/data_models/strategies/doe.py b/bofire/data_models/strategies/doe.py
index bcfbac797..a627edba2 100644
--- a/bofire/data_models/strategies/doe.py
+++ b/bofire/data_models/strategies/doe.py
@@ -112,6 +112,8 @@ class DoEStrategy(Strategy):
 
     verbose: bool = False  # get rid of this at a later stage
     ipopt_options: Optional[Dict] = None
+    use_hessian: bool = False
+    use_cyipopt: Optional[bool] = None
     sampling: Optional[List[List]] = None
 
     def is_constraint_implemented(self, my_type: Type[Constraint]) -> bool:
diff --git a/bofire/strategies/doe/branch_and_bound.py b/bofire/strategies/doe/branch_and_bound.py
index 7e3556114..edd7cc8cb 100644
--- a/bofire/strategies/doe/branch_and_bound.py
+++ b/bofire/strategies/doe/branch_and_bound.py
@@ -235,6 +235,8 @@ def find_local_max_ipopt_BaB(
     categorical_groups: Optional[List[List[ContinuousInput]]] = None,
     discrete_variables: Optional[Dict[str, Tuple[ContinuousInput, List[float]]]] = None,
     verbose: bool = False,
+    use_hessian: bool = False,
+    use_cyipopt: Optional[bool] = None,
 ) -> pd.DataFrame:
     """Function computing a d-optimal design" for a given domain and model.
     It allows for the problem to have categorical values which is solved by Branch-and-Bound
@@ -261,6 +263,9 @@ def find_local_max_ipopt_BaB(
             verbose (bool): if true, print information during the optimization process
             transform_range (Optional[Bounds]): range to which the input variables are transformed.
                 If None is provided, the features will not be scaled. Defaults to None.
+            use_hessian (bool): if true, the hessian will be used in the optimization process. Defaults to False.
+            use_cyipopt (Optional[bool]): if true, the cyipopt solver will be used, otherwise scipy.minimize(). Defaults to None.
+                If None is provided, the cyipopt solver will be used if available.
         Returns:
             A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
             local optimum.
@@ -324,6 +329,8 @@ def find_local_max_ipopt_BaB(
         sampling,
         None,
         partially_fixed_experiments=initial_branch,
+        use_hessian=use_hessian,
+        use_cyipopt=use_cyipopt,
     )
     initial_value = objective_function.evaluate(
         initial_design.to_numpy().flatten(),
@@ -367,6 +374,8 @@ def find_local_max_ipopt_exhaustive(
     categorical_groups: Optional[List[List[ContinuousInput]]] = None,
     discrete_variables: Optional[Dict[str, Tuple[ContinuousInput, List[float]]]] = None,
     verbose: bool = False,
+    use_hessian: bool = False,
+    use_cyipopt: Optional[bool] = None,
 ) -> pd.DataFrame:
     """Function computing a d-optimal design" for a given domain and model.
     It allows for the problem to have categorical values which is solved by exhaustive search
@@ -392,6 +401,9 @@ def find_local_max_ipopt_exhaustive(
                 with key:(relaxed variable, valid values). Defaults to None
             verbose (bool): if true, print information during the optimization process
             transform_range (Optional[Bounds]): range to which the input variables are transformed.
+            use_hessian (bool): if true, the hessian will be used in the optimization process. Defaults to False.
+            use_cyipopt (Optional[bool]): if true, the cyipopt solver will be used, otherwise scipy.minimize(). Defaults to None.
+                If None is provided, the cyipopt solver will be used if available.
         Returns:
             A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
             local optimum.
@@ -503,6 +515,8 @@ def find_local_max_ipopt_exhaustive(
                 sampling,
                 None,
                 one_set_of_experiments,
+                use_hessian=use_hessian,
+                use_cyipopt=use_cyipopt,
             )
             domain.validate_candidates(
                 candidates=current_design.apply(lambda x: np.round(x, 8)),
diff --git a/bofire/strategies/doe/design.py b/bofire/strategies/doe/design.py
index 5bdbdfea4..b55a3af06 100644
--- a/bofire/strategies/doe/design.py
+++ b/bofire/strategies/doe/design.py
@@ -4,7 +4,6 @@
 import numpy as np
 import pandas as pd
 from formulaic import Formula
-from scipy.optimize._minimize import standardize_constraints
 
 from bofire.data_models.constraints.api import (
     ConstraintNotFulfilledError,
@@ -17,6 +16,7 @@
 from bofire.data_models.strategies.doe import AnyOptimalityCriterion
 from bofire.strategies.doe.objective import get_objective_function
 from bofire.strategies.doe.utils import (
+    _minimize,
     constraints_as_scipy_constraints,
     nchoosek_constraints_as_bounds,
 )
@@ -31,6 +31,8 @@ def find_local_max_ipopt(
     sampling: Optional[pd.DataFrame] = None,
     fixed_experiments: Optional[pd.DataFrame] = None,
     partially_fixed_experiments: Optional[pd.DataFrame] = None,
+    use_hessian: bool = False,
+    use_cyipopt: Optional[bool] = None,
 ) -> pd.DataFrame:
     """Function computing an optimal design for a given domain and model.
 
@@ -48,6 +50,9 @@ def find_local_max_ipopt(
             Variables can be fixed to one value or can be set to a range by setting a tuple with lower and upper bound
             Non-fixed variables have to be set to None or nan.
         criterion (OptimalityCriterion): OptimalityCriterion object indicating which criterion function to use.
+        use_hessian (bool): If True, the hessian of the objective function is used. Default is False.
+        use_cyipopt (bool, optional): If True, cyipopt is used, otherwise scipy.minimize(). Default is None.
+            If None, cyipopt is used if available.
 
     Returns:
         A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
@@ -58,14 +63,6 @@ def find_local_max_ipopt(
     # Checks and preparation steps
     #
 
-    # warn user if IPOPT scipy interface is not available
-    try:
-        from cyipopt import minimize_ipopt  # type: ignore
-    except ImportError:
-        raise ImportError(
-            "cyipopt is not installed. Install it via `conda install -c conda-forge cyipopt`"
-        )
-
     objective_function = get_objective_function(
         criterion, domain=domain, n_experiments=n_experiments
     )
@@ -166,28 +163,27 @@ def find_local_max_ipopt(
     # set ipopt options
     if ipopt_options is None:
         ipopt_options = {}
-    _ipopt_options = {"maxiter": 500, "disp": 0}
+    _ipopt_options = {"max_iter": 500, "print_level": 0}
     for key in ipopt_options.keys():
         _ipopt_options[key] = ipopt_options[key]
-    if _ipopt_options["disp"] > 12:
-        _ipopt_options["disp"] = 0
+    if _ipopt_options["print_level"] > 12:
+        _ipopt_options["print_level"] = 0
 
     #
     # Do the optimization
     #
-
-    result = minimize_ipopt(
-        objective_function.evaluate,
+    x = _minimize(
+        objective_function=objective_function,
         x0=x0,
         bounds=bounds,
-        # "SLSQP" has no deeper meaning here and just ensures correct constraint standardization
-        constraints=standardize_constraints(constraints, x0, "SLSQP"),
-        options=_ipopt_options,
-        jac=objective_function.evaluate_jacobian,
+        constraints=constraints,
+        use_hessian=use_hessian,
+        ipopt_options=_ipopt_options,
+        use_cyipopt=use_cyipopt,
     )
 
     design = pd.DataFrame(
-        result["x"].reshape(n_experiments, len(domain.inputs)),
+        x.reshape(n_experiments, len(domain.inputs)),
         columns=domain.inputs.get_keys(),
         index=[f"exp{i}" for i in range(n_experiments)],
     )
diff --git a/bofire/strategies/doe/doe_problem.py b/bofire/strategies/doe/doe_problem.py
new file mode 100644
index 000000000..d69ad2413
--- /dev/null
+++ b/bofire/strategies/doe/doe_problem.py
@@ -0,0 +1,142 @@
+# warn user if IPOPT scipy interface is not available
+try:
+    from cyipopt import Problem  # type: ignore
+except ImportError:
+
+    class Problem:
+        def __init__(self, *args, **kwargs):
+            raise ImportError(
+                "cyipopt is not installed. Install it via `conda install -c conda-forge cyipopt`"
+            )
+
+
+from typing import List, Optional, Tuple, Union
+
+import numpy as np
+from scipy import sparse
+from scipy.optimize._minimize import LinearConstraint, NonlinearConstraint
+
+from bofire.strategies.doe.objective_base import Objective
+
+
+class FirstOrderDoEProblem(Problem):  # type: ignore
+    def __init__(
+        self,
+        doe_objective: Objective,
+        bounds: List[Tuple[float, float]],
+        constraints: Optional[
+            List[Union[NonlinearConstraint, LinearConstraint]]
+        ] = None,
+    ):
+        self.doe_objective = doe_objective
+        n_vars = doe_objective.n_vars * doe_objective.n_experiments
+
+        # assemble linear constraints to one big linear constraint
+        if constraints is None:
+            constraints = []
+        self.linear_constraints = [
+            c for c in constraints if isinstance(c, LinearConstraint)
+        ]
+        self.A = (
+            sparse.coo_array(np.vstack([c.A for c in self.linear_constraints]))
+            if len(self.linear_constraints) > 0
+            else sparse.coo_matrix([]).T
+        )
+
+        # assemble nonlinear constraints to one big nonlinear constraint
+        self.nonlinear_constraints = [
+            c for c in constraints if isinstance(c, NonlinearConstraint)
+        ]
+
+        # generate jacobian structure
+        start_row = self.A.shape[0]
+        nonlinear_constraint_row = []
+        nonlinear_constraint_col = []
+        for _ in self.nonlinear_constraints:
+            nonlinear_constraint_row.append(
+                start_row
+                + np.repeat(
+                    np.arange(doe_objective.n_experiments), doe_objective.n_vars
+                )
+            )
+            nonlinear_constraint_col.append(np.arange(n_vars))
+            start_row += doe_objective.n_experiments
+        nonlinear_constraint_row = (
+            np.concatenate(nonlinear_constraint_row)
+            if len(nonlinear_constraint_row) > 0
+            else np.array([])
+        )
+        nonlinear_constraint_col = (
+            np.concatenate(nonlinear_constraint_col)
+            if len(nonlinear_constraint_col) > 0
+            else np.array([])
+        )
+
+        self.row = np.concatenate((self.A.row, nonlinear_constraint_row))
+        self.col = np.concatenate((self.A.col, nonlinear_constraint_col))
+
+        self.hessian_structure = np.nonzero(np.ones((4, 4)))
+
+        super().__init__(
+            n=int(n_vars),
+            m=sum([len(c.lb) for c in constraints]),
+            lb=np.array([b[0] for b in bounds]),
+            ub=np.array([b[1] for b in bounds]),
+            cl=np.array(
+                [c.lb for c in self.linear_constraints]
+                + [c.lb for c in self.nonlinear_constraints]
+            ).flatten(),
+            cu=np.array(
+                [c.ub for c in self.linear_constraints]
+                + [c.ub for c in self.nonlinear_constraints]
+            ).flatten(),
+        )
+
+    def objective(self, x: np.ndarray) -> float:
+        return self.doe_objective.evaluate(x)
+
+    def gradient(self, x: np.ndarray) -> np.ndarray:
+        return self.doe_objective.evaluate_jacobian(x)
+
+    def constraints(self, x: np.ndarray) -> np.ndarray:
+        linear = (
+            np.array(self.A @ x) if len(self.linear_constraints) > 0 else np.array([])
+        )
+        if len(linear.shape) == 0:
+            linear = np.array([linear])
+        nonlinear = (
+            np.array([c.fun(x) for c in self.nonlinear_constraints]).flatten()
+            if len(self.nonlinear_constraints) > 0
+            else np.array([])
+        )
+        return np.concatenate((linear, nonlinear))
+
+    def jacobian(self, x: np.ndarray):
+        linear = self.A.data if len(self.linear_constraints) > 0 else np.array([])
+        nonlinear = (
+            np.concatenate([c.jac(x, sparse=True) for c in self.nonlinear_constraints])
+            if len(self.nonlinear_constraints) > 0
+            else np.array([])
+        )
+        return np.concatenate((linear, nonlinear))
+
+    def jacobianstructure(self) -> Tuple[np.ndarray, np.ndarray]:
+        return self.row, self.col
+
+
+class SecondOrderDoEProblem(FirstOrderDoEProblem):
+    def hessian(
+        self, x: np.ndarray, lagrange: np.ndarray, obj_factor: float
+    ) -> np.ndarray:
+        H = obj_factor * np.array(self.doe_objective.evaluate_hessian(x))
+
+        # linear constraints have a vanishing hessian
+        for i, c in enumerate(self.nonlinear_constraints):
+            H += lagrange[i] * c.hess(x)
+
+        row, col = self.hessianstructure()
+
+        return H[row, col]
+
+    def hessianstructure(self) -> Tuple[np.ndarray, np.ndarray]:
+        return self.hessian_structure
diff --git a/bofire/strategies/doe/objective.py b/bofire/strategies/doe/objective.py
index 2ebb0cae6..53390fb43 100644
--- a/bofire/strategies/doe/objective.py
+++ b/bofire/strategies/doe/objective.py
@@ -8,8 +8,8 @@
 import pandas as pd
 import torch
 from formulaic import Formula
-from scipy.optimize._minimize import standardize_constraints
 from torch import Tensor
+from torch.autograd.functional import hessian, jacobian
 
 from bofire.data_models.constraints.api import (
     ConstraintNotFulfilledError,
@@ -29,93 +29,32 @@
     SpaceFillingCriterion,
 )
 from bofire.data_models.types import Bounds
-from bofire.strategies.doe.transform import IndentityTransform, MinMaxTransform
+from bofire.strategies.doe.objective_base import Objective
 from bofire.strategies.doe.utils import (
+    _minimize,
     constraints_as_scipy_constraints,
+    convert_formula_to_string,
     get_formula_from_string,
     nchoosek_constraints_as_bounds,
 )
 from bofire.utils.torch_tools import tkwargs
 
 
-class Objective:
-    def __init__(
-        self,
-        domain: Domain,
-        n_experiments: int,
-        delta: float = 1e-6,
-        transform_range: Optional[Bounds] = None,
-    ) -> None:
-        """Args:
-        domain (Domain): A domain defining the DoE domain together with model_type.
-        model_type (str or Formula): A formula containing all model terms.
-        n_experiments (int): Number of experiments
-        delta (float): A regularization parameter for the information matrix. Default value is 1e-3.
-        transform_range (Bounds, optional): range to which the input variables are transformed before applying the objective function. Default is None.
-
-        """
-        self.domain = deepcopy(domain)
-
-        if transform_range is None:
-            self.transform = IndentityTransform()
-        else:
-            self.transform = MinMaxTransform(
-                inputs=self.domain.inputs,
-                feature_range=tuple(transform_range),  # type: ignore
-            )
-
-        self.n_experiments = n_experiments
-        self.delta = delta
-
-        self.vars = self.domain.inputs.get_keys()
-        self.n_vars = len(self.domain.inputs)
-
-    def __call__(self, x: np.ndarray) -> float:
-        return self.evaluate(x)
-
-    def evaluate(self, x: np.ndarray) -> float:
-        return self._evaluate(self.transform(x=x))
-
-    @abstractmethod
-    def _evaluate(self, x: np.ndarray) -> float:
-        pass
-
-    def evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        return self._evaluate_jacobian(self.transform(x)) * self.transform.jacobian(x=x)
-
-    @abstractmethod
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        pass
-
-    @abstractmethod
-    def _convert_input_to_model_tensor(
-        self,
-        x: np.ndarray,
-        requires_grad: bool = True,
-    ) -> Tensor:
-        """Args:
-        x: x (np.ndarray): values of design variables a 1d array.
-        """
-        assert x.ndim == 1, "values of design should be 1d array"
-        pass
-
-
 class ModelBasedObjective(Objective):
     def __init__(
         self,
         domain: Domain,
-        model: Formula,
+        formula: Formula,
         n_experiments: int,
-        delta: float = 1e-6,
+        delta: float = 1e-7,
         transform_range: Optional[Bounds] = None,
     ) -> None:
         """Args:
-        domain (Domain): A domain defining the DoE domain together with model_type.
-        model_type (str or Formula): A formula containing all model terms.
+        domain (Domain): A domain defining the DoE domain together with formula.
+        formula (str or Formula): A formula containing all model terms.
         n_experiments (int): Number of experiments
-        delta (float): A regularization parameter for the information matrix. Default value is 1e-3.
+        delta (float): A regularization parameter for the information matrix. Default value is 1e-7.
         transform_range (Bounds, optional): range to which the input variables are transformed before applying the objective function. Default is None.
-
         """
         super().__init__(
             domain=domain,
@@ -124,56 +63,63 @@ def __init__(
             transform_range=transform_range,
         )
 
-        self.model = deepcopy(model)
+        self.formula = deepcopy(formula)
+        self.model_terms_string_expression = convert_formula_to_string(
+            domain=domain, formula=formula
+        )
+        self.n_model_terms = len(list(np.array(formula, dtype=str)))
 
-        self.model_terms = list(np.array(model, dtype=str))
-        self.n_model_terms = len(self.model_terms)
+    def _evaluate(self, x: np.ndarray) -> float:
+        D = torch.tensor(
+            x.reshape(len(x.flatten()) // self.n_vars, self.n_vars),
+            **tkwargs,
+            requires_grad=False,
+        )
+        return float(self._evaluate_tensor(D))
 
-        # terms for model jacobian
-        self.terms_jacobian_t = []
-        for var in self.vars:
-            _terms = [
-                str(term).replace(":", "*") + f" + 0 * {self.vars[0]}"
-                for term in model.differentiate(var, use_sympy=True)
-            ]  # 0*vars[0] added to make sure terms are evaluated as series, not as number
-            terms = "["
-            for t in _terms:
-                terms += t + ", "
-            terms = terms[:-1] + "]"
+    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
+        D = torch.tensor(
+            x.reshape(len(x.flatten()) // self.n_vars, self.n_vars),
+            **tkwargs,
+        )
 
-            self.terms_jacobian_t.append(terms)
+        return (
+            torch.tensor(jacobian(self._evaluate_tensor, D)).detach().numpy().flatten()
+        )
 
-    def _convert_input_to_model_tensor(
-        self,
-        x: np.ndarray,
-        requires_grad: bool = True,
-    ) -> Tensor:
-        """Args:
-        x: x (np.ndarray): values of design variables a 1d array.
-        """
-        assert x.ndim == 1, "values of design should be 1d array"
-        X = pd.DataFrame(
-            x.reshape(len(x.flatten()) // self.n_vars, self.n_vars),
-            columns=self.vars,
+    # FIXME: currently not returning the hessian in a way that is compatible with ipopt
+    # also, the hessians of the constraints are missing
+    def _evaluate_hessian(self, x: np.ndarray) -> np.ndarray:
+        def _evaluate_from_flat_tensor(x: Tensor) -> Tensor:
+            D = x.reshape(len(x.flatten()) // self.n_vars, self.n_vars)
+            return self._evaluate_tensor(D)
+
+        return (
+            torch.tensor(
+                hessian(_evaluate_from_flat_tensor, torch.tensor(x, **tkwargs))
+            )
+            .detach()
+            .numpy()
         )
-        # scale to [0, 1]
-        # lower, upper = self.domain.inputs.get_bounds(specs={}, experiments=X)
-        # lower = np.array(lower)
-        # upper = np.array(upper)
-        # X = (X - lower) / (upper - lower)
-        # X = X * 2 - 1
-        # get model matrix
-        X = self.model.get_model_matrix(X)
-        return torch.tensor(X.values, requires_grad=requires_grad, **tkwargs)
 
-    def _model_jacobian_t(self, x: np.ndarray) -> np.ndarray:
-        """Computes the transpose of the model jacobian for each experiment in input x."""
-        X = pd.DataFrame(x.reshape(self.n_experiments, self.n_vars), columns=self.vars)
-        jacobians = np.swapaxes(X.eval(self.terms_jacobian_t), 0, 2)  # type: ignore
-        return np.swapaxes(jacobians, 1, 2)
+    def _evaluate_tensor(self, D: Tensor) -> Tensor:
+        """Evaluate the objective function on the design matrix as a tensor."""
+        var_dict = {var: D[:, i] for i, var in enumerate(self.vars)}
+        var_dict["torch"] = torch  # type: ignore
+        X = eval(str(self.model_terms_string_expression), {}, var_dict)
+        return self._criterion(X)
 
-    def get_model_matrix(self, design: pd.DataFrame) -> pd.DataFrame:
-        return self.model.get_model_matrix(design)
+    @abstractmethod
+    def _criterion(self, X: Tensor) -> Tensor:
+        """Function implementing the criterion acting on the design matrix X.
+
+        Args:
+            X (Tensor): Design matrix.
+
+        Returns:
+            A tensor containing a single number which is the value of the criterion.
+        """
+        pass
 
 
 class IOptimality(ModelBasedObjective):
@@ -184,9 +130,9 @@ class IOptimality(ModelBasedObjective):
     def __init__(
         self,
         domain: Domain,
-        model: Formula,
+        formula: Formula,
         n_experiments: int,
-        delta: float = 1e-6,
+        delta: float = 1e-7,
         transform_range: Optional[Bounds] = None,
         n_space_filling_points: Optional[int] = None,
         ipopt_options: Optional[dict] = None,
@@ -194,22 +140,15 @@ def __init__(
         """
         Args:
             domain (Domain): A domain defining the DoE domain together with model_type.
-            model_type (str or Formula): A formula containing all model terms.
+            formula (str or Formula): A formula containing all model terms.
             n_experiments (int): Number of experiments
             delta (float): A regularization parameter for the information matrix. Default value is 1e-3.
             transform_range (Bounds, optional): range to which the input variables are transformed before applying the objective function. Default is None.
             n_space_filling_points (int, optional): Number of space filling points. Only relevant if SpaceFilling is used
             ipopt_options (dict, optional): Options for the Ipopt solver to generate space filling point.
-                If None is provided, the default options (maxiter = 500) are used.
+                If None is provided, the default options (max_iter = 500) are used.
         """
 
-        try:
-            from cyipopt import minimize_ipopt  # type: ignore
-        except ImportError:
-            raise ImportError(
-                "cyipopt is not installed. Install it via `conda install -c conda-forge cyipopt`"
-            )
-
         if transform_range is not None:
             raise ValueError(
                 "IOptimality does not support transformations of the input variables."
@@ -217,7 +156,7 @@ def __init__(
 
         super().__init__(
             domain=domain,
-            model=model,
+            formula=formula,
             n_experiments=n_experiments,
             delta=delta,
             transform_range=transform_range,
@@ -239,7 +178,7 @@ def __init__(
                 .to_numpy()
                 .flatten()
             )
-            objective = SpaceFilling(
+            objective_function = SpaceFilling(
                 domain=domain,
                 n_experiments=n_space_filling_points,
                 delta=delta,
@@ -250,19 +189,17 @@ def __init__(
             )
             bounds = nchoosek_constraints_as_bounds(domain, n_space_filling_points)
 
-            result = minimize_ipopt(
-                objective.evaluate,
+            x = _minimize(
+                objective_function=objective_function,
                 x0=x0,
                 bounds=bounds,
-                constraints=standardize_constraints(constraints, x0, "SLSQP"),
-                options=ipopt_options
-                if ipopt_options is not None
-                else {"maxiter": 500, "disp": 0},
-                jac=objective.evaluate_jacobian,
+                constraints=constraints,
+                ipopt_options={"print_level": 0},
+                use_hessian=False,
             )
 
             self.Y = pd.DataFrame(
-                result["x"].reshape(n_space_filling_points, len(domain.inputs)),
+                x.reshape(n_space_filling_points, len(domain.inputs)),
                 columns=domain.inputs.get_keys(),
                 index=[f"gridpoint{i}" for i in range(n_space_filling_points)],
             )
@@ -308,162 +245,24 @@ def __init__(
                 UserWarning,
             )
 
-        X = model.get_model_matrix(self.Y).to_numpy()
+        X = formula.get_model_matrix(self.Y).to_numpy()
         self.YtY = torch.from_numpy(X.T @ X) / n_space_filling_points
         self.YtY.requires_grad = False
 
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes trace((Y.T@Y) / nY @ inv(X.T@X + delta)).
-        Where X is the model matrix corresponding to x, Y is the model matrix of points
-        uniformly filling up the feasible space and nY is the number of such points.
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-        Returns:
-            trace((Y.T@Y) / nY @ inv(X.T@X + delta))
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        return float(
-            torch.trace(
-                self.YtY.detach()
-                @ torch.linalg.inv(
-                    X.detach().T @ X.detach()
-                    + self.delta * torch.eye(self.n_model_terms)
-                )
-            )
-        )
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of trace((Y.T@Y) / nY @ inv(X.T@X + delta)).
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-        Returns:
-            The jacobian of trace((Y.T@Y) / nY @ inv(X.T@X + delta))
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.trace(
+    def _criterion(self, X: Tensor) -> Tensor:
+        return torch.trace(
             self.YtY.detach()
             @ torch.linalg.inv(X.T @ X + self.delta * torch.eye(self.n_model_terms))
-        ).backward()
-        J1 = X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments, self.n_vars, self.n_model_terms
         )
 
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
-
 
 class DOptimality(ModelBasedObjective):
-    """A class implementing the evaluation of logdet(X.T@X + delta) and its jacobian w.r.t. the inputs.
-    The Jacobian can be divided into two parts, one for logdet(X.T@ + delta) w.r.t. X (there is a simple
-    closed expression for this one) and one model dependent part for the jacobian of X.T@X
-    w.r.t. the inputs. Because each row of X only depends on the inputs of one experiment
-    the second part can be formulated in a simplified way. It is built up with n_experiment
-    blocks of the same structure which is represented by the attribute jacobian_building_block.
-
-    A nice derivation for the "first part" of the jacobian can be found [here](https://angms.science/doc/LA/logdet.pdf).
-    The second part consists of the partial derivatives of the model terms with
-    respect to the inputs. We denote the value of the i-th model term from the j-th experiment
-    with y_ij and the i-th input value of the j-th experiment with x_ij. N stands for the number
-    of model terms, n for the number of input terms and M for the number of experiments.
-    Here, we only consider models up to second order, but the computation can easily be extended
-    for higher-ordermodels.
-
-    To do the computation in the most basic way, we could compute the partial derivative of every
-    single model term and experiment with respect to every single input and experiment. We could write
-    this in one large matrix and multiply the first part of the gradient as a long vector from the right
-    side.
-    But because of the structure of the domain we can do the same computation with much smaller
-    matrices:
-    First, we write the first part of the jacobian as the matrix (df/dy_ij)_ij where i goes from 1 to N
-    and j goes from 1 to M.
-
-    Second, we compute a rank 3 tensor (K_kij)_kij. k goes from 1 to M, i from 1 to n and j from 1 to N.
-    For each k (K_kij)_ij contains the partial derivatives (dy_jk/dx_ik)_ij. Note that the values of the
-    entries of (dy_jk/dx_ik)_ij only depend on the input values of the k-th experiment. The function
-    default_jacobian_building_block implements the computation of these matrices/"building blocks".
-
-    Then, we notice that the model term values of the j-th experiment only depend on the input values of
-    the j-th experiment. Thus, to compute the partial derivative df/dx_ik we only have to compute the euclidean
-    scalar product of (K_kij)_j and (df/dy_jk)_j. The way how we built the two parts of the jacobian allows us
-    to compute this scalar product in a vectorized way for all x_ik at once, see also JacobianForLogDet.jacobian.
+    """A class implementing the evaluation of logdet(X.T@X + delta), i.e. the
+    D-optimality criterion
     """
 
-    def __init__(
-        self,
-        domain: Domain,
-        model: Formula,
-        n_experiments: int,
-        delta: float = 1e-7,
-        transform_range: Optional[Bounds] = None,
-    ) -> None:
-        super().__init__(
-            domain=domain,
-            model=model,
-            n_experiments=n_experiments,
-            delta=delta,
-            transform_range=transform_range,
-        )
-
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes the minus one times the sum of the log of the eigenvalues of X.T @ X + delta.
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            -log(det(X.T@X+delta))
-
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        return float(
-            -1
-            * torch.logdet(
-                X.detach().T @ X.detach() + self.delta * torch.eye(self.n_model_terms),
-            ),
-        )
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of minus one times the log of the determinant of X.T @ X + delta.
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            The jacobian of -log(det(X.T@X+delta)) as numpy array
-
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.logdet(X.T @ X + self.delta * torch.eye(self.n_model_terms)).backward()
-        J1 = -1 * X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments,
-            self.n_vars,
-            self.n_model_terms,
-        )
-
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
+    def _criterion(self, X: Tensor) -> Tensor:
+        return -1 * torch.logdet(X.T @ X + self.delta * torch.eye(self.n_model_terms))
 
 
 class AOptimality(ModelBasedObjective):
@@ -472,61 +271,11 @@ class AOptimality(ModelBasedObjective):
     the jacobian of tr((X.T@X + delta)^-1) instead of logdet(X.T@X + delta).
     """
 
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes the trace of the inverse of X.T @ X + delta.
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            tr((X.T@X+delta)^-1)
-
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        return float(
-            torch.trace(
-                torch.linalg.inv(
-                    X.detach().T @ X.detach()
-                    + self.delta * torch.eye(self.n_model_terms),
-                ),
-            ),
-        )
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of the trace of the inverse of X.T @ X + delta.
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            The jacobian of tr((X.T@X+delta)^-1) as numpy array
-
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.trace(
-            torch.linalg.inv(X.T @ X + self.delta * torch.eye(self.n_model_terms)),
-        ).backward()
-        J1 = X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments,
-            self.n_vars,
-            self.n_model_terms,
+    def _criterion(self, X: Tensor) -> Tensor:
+        return torch.trace(
+            torch.linalg.inv(X.T @ X + self.delta * torch.eye(self.n_model_terms))
         )
 
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
-
 
 class GOptimality(ModelBasedObjective):
     """A class implementing the evaluation of max(diag(H)) and its jacobian w.r.t. the inputs where
@@ -535,65 +284,15 @@ class GOptimality(ModelBasedObjective):
     logdet(X.T@X + delta).
     """
 
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes the maximum diagonal entry of H = X @ (X.T@X + delta)^-1 @ X.T .
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            max(diag(H))
-
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        H = (
-            X.detach()
-            @ torch.linalg.inv(
-                X.detach().T @ X.detach() + self.delta * torch.eye(self.n_model_terms),
-            )
-            @ X.detach().T
-        )
-        return float(torch.max(torch.diag(H)))
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of the maximum diagonal element of H = X @ (X.T @ X + delta)^-1 @ X.T.
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            The jacobian of max(diag(H)) as numpy array
-
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.max(
+    def _criterion(self, X: Tensor) -> Tensor:
+        return torch.max(
             torch.diag(
                 X
                 @ torch.linalg.inv(X.T @ X + self.delta * torch.eye(self.n_model_terms))
                 @ X.T,
             ),
-        ).backward()
-        J1 = X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments,
-            self.n_vars,
-            self.n_model_terms,
         )
 
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
-
 
 class EOptimality(ModelBasedObjective):
     """A class implementing the evaluation of minus one times the minimum eigenvalue of (X.T @ X + delta)
@@ -602,61 +301,11 @@ class EOptimality(ModelBasedObjective):
     logdet(X.T@X + delta).
     """
 
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes minus one times the minimum eigenvalue of (X.T@X + delta).
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            min(eigvals(X.T @ X + delta))
-
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        return -1 * float(
-            torch.min(
-                torch.linalg.eigvalsh(
-                    X.detach().T @ X.detach()
-                    + self.delta * torch.eye(self.n_model_terms),
-                ),
-            ),
-        )
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of minus one times the minimum eigenvalue of (X.T @ X + delta).
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            The jacobian of -1 * min(eigvals(X.T @ X + delta)) as numpy array
-
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.min(
+    def _criterion(self, X: Tensor) -> Tensor:
+        return -1 * torch.min(
             torch.linalg.eigvalsh(X.T @ X + self.delta * torch.eye(self.n_model_terms)),
-        ).backward()
-        J1 = -1 * X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments,
-            self.n_vars,
-            self.n_model_terms,
         )
 
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
-
 
 class KOptimality(ModelBasedObjective):
     """A class implementing the evaluation of the condition number of (X.T @ X + delta)
@@ -665,58 +314,11 @@ class KOptimality(ModelBasedObjective):
     of (X.T @ X + delta) instead of logdet(X.T@X + delta).
     """
 
-    def _evaluate(self, x: np.ndarray) -> float:
-        """Computes condition number of (X.T@X + delta).
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            cond(X.T @ X + delta)
-
-        """
-        X = self._convert_input_to_model_tensor(x, requires_grad=False)
-        return float(
-            torch.linalg.cond(
-                X.detach().T @ X.detach() + self.delta * torch.eye(self.n_model_terms),
-            ),
-        )
-
-    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Computes the jacobian of the condition number of (X.T @ X + delta).
-        Where X is the model matrix corresponding to x.
-
-        Args:
-            x (np.ndarray): values of design variables a 1d array.
-
-        Returns:
-            The jacobian of cond(X.T @ X + delta) as numpy array
-
-        """
-        # get model matrix X
-        X = self._convert_input_to_model_tensor(x, requires_grad=True)
-
-        # first part of jacobian
-        torch.linalg.cond(
+    def _criterion(self, X: Tensor) -> Tensor:
+        return torch.linalg.cond(
             X.T @ X + self.delta * torch.eye(self.n_model_terms),
-        ).backward()
-        J1 = X.grad.detach().numpy()  # type: ignore
-        J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
-            self.n_experiments,
-            self.n_vars,
-            self.n_model_terms,
         )
 
-        # second part of jacobian
-        J2 = self._model_jacobian_t(x)
-
-        # combine both parts
-        J = J1 * J2
-        J = np.sum(J, axis=-1)
-
-        return J.flatten()
-
 
 class SpaceFilling(Objective):
     def _evaluate(self, x: np.ndarray) -> float:
@@ -749,14 +351,14 @@ def get_objective_function(
     if criterion is None:
         return DOptimality(
             domain,
-            model=get_formula_from_string(domain=domain),
+            formula=get_formula_from_string(domain=domain),
             n_experiments=n_experiments,
         )
     if isinstance(criterion, DoEOptimalityCriterion):
         if isinstance(criterion, DOptimalityCriterion):
             return DOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
@@ -764,7 +366,7 @@ def get_objective_function(
         if isinstance(criterion, AOptimalityCriterion):
             return AOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
@@ -772,7 +374,7 @@ def get_objective_function(
         if isinstance(criterion, GOptimalityCriterion):
             return GOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
@@ -780,7 +382,7 @@ def get_objective_function(
         if isinstance(criterion, EOptimalityCriterion):
             return EOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
@@ -788,7 +390,7 @@ def get_objective_function(
         if isinstance(criterion, KOptimalityCriterion):
             return KOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
@@ -796,7 +398,7 @@ def get_objective_function(
         if isinstance(criterion, IOptimalityCriterion):
             return IOptimality(
                 domain,
-                model=get_formula_from_string(criterion.formula, domain),
+                formula=get_formula_from_string(criterion.formula, domain),
                 n_experiments=n_experiments,
                 delta=criterion.delta,
                 transform_range=criterion.transform_range,
diff --git a/bofire/strategies/doe/objective_base.py b/bofire/strategies/doe/objective_base.py
new file mode 100644
index 000000000..0d55eeb01
--- /dev/null
+++ b/bofire/strategies/doe/objective_base.py
@@ -0,0 +1,74 @@
+from abc import abstractmethod
+from copy import deepcopy
+from typing import Optional
+
+import numpy as np
+
+from bofire.data_models.domain.api import Domain
+from bofire.data_models.types import Bounds
+from bofire.strategies.doe.transform import IndentityTransform, MinMaxTransform
+
+
+class Objective:
+    """Base class for objectives in the context of Design of Experiments (DoE)."""
+
+    def __init__(
+        self,
+        domain: Domain,
+        n_experiments: int,
+        delta: float = 1e-7,
+        transform_range: Optional[Bounds] = None,
+    ) -> None:
+        """Args:
+        domain (Domain): A domain defining the DoE domain together with model_type.
+        model_type (str or Formula): A formula containing all model terms.
+        n_experiments (int): Number of experiments
+        delta (float): A regularization parameter for the information matrix. Default value is 1e-7.
+        transform_range (Bounds, optional): range to which the input variables are transformed before applying the objective function. Default is None.
+
+        """
+        self.domain = deepcopy(domain)
+
+        if transform_range is None:
+            self.transform = IndentityTransform()
+        else:
+            self.transform = MinMaxTransform(
+                inputs=self.domain.inputs,
+                feature_range=tuple(transform_range),  # type: ignore
+            )
+
+        self.n_experiments = n_experiments
+        self.delta = delta
+
+        self.vars = self.domain.inputs.get_keys()
+        self.n_vars = len(self.domain.inputs)
+
+    def __call__(self, x: np.ndarray) -> float:
+        return self.evaluate(x)
+
+    def evaluate(self, x: np.ndarray) -> float:
+        return self._evaluate(self.transform(x=x))
+
+    def evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
+        return self._evaluate_jacobian(self.transform(x=x)) * self.transform.jacobian(
+            x=x
+        )
+
+    def evaluate_hessian(self, x: np.ndarray) -> np.ndarray:
+        return (
+            self.transform.jacobian(x=x)[None, :]
+            * self._evaluate_hessian(self.transform(x=x))
+            * self.transform.jacobian(x=x)[:, None]
+        )
+
+    @abstractmethod
+    def _evaluate(self, x: np.ndarray) -> float:
+        pass
+
+    @abstractmethod
+    def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
+        pass
+
+    @abstractmethod
+    def _evaluate_hessian(self, x: np.ndarray) -> np.ndarray:
+        pass
diff --git a/bofire/strategies/doe/utils.py b/bofire/strategies/doe/utils.py
index aa4a6a651..95d54d915 100644
--- a/bofire/strategies/doe/utils.py
+++ b/bofire/strategies/doe/utils.py
@@ -1,11 +1,14 @@
+import importlib.util
 import sys
 from itertools import combinations
-from typing import List, Optional, Union
+from typing import List, Optional, Tuple, Union
 
 import numpy as np
 import pandas as pd
+import scipy.optimize as opt
 from formulaic import Formula
 from scipy.optimize import LinearConstraint, NonlinearConstraint
+from scipy.optimize._minimize import standardize_constraints
 
 from bofire.data_models.constraints.api import (
     Constraint,
@@ -19,9 +22,17 @@
 from bofire.data_models.domain.domain import Domain
 from bofire.data_models.features.continuous import ContinuousInput
 from bofire.data_models.strategies.api import RandomStrategy as RandomStrategyDataModel
+from bofire.strategies.doe.doe_problem import (
+    FirstOrderDoEProblem,
+    SecondOrderDoEProblem,
+)
+from bofire.strategies.doe.objective_base import Objective
 from bofire.strategies.random import RandomStrategy
 
 
+CYIPOPT_AVAILABLE = importlib.util.find_spec("cyipopt") is not None
+
+
 def get_formula_from_string(
     model_type: Union[str, Formula] = "linear",
     domain: Optional[Domain] = None,
@@ -77,6 +88,38 @@ def get_formula_from_string(
     return formula
 
 
+def convert_formula_to_string(
+    domain: Domain,
+    formula: Formula,
+) -> str:
+    """Converts a formula to a string.
+
+    Args:
+        domain (Domain): The domain that contain information about the input names.
+        formula (Formula): A formula object that should be converted to a string. If
+        the formula has both a left and right hand side, only the right hand side is
+        considered.
+
+    Returns:
+        A string representation of the formula that can be evaluated using pytorch.
+
+    """
+    if hasattr(formula, "rhs"):
+        formula = formula.rhs
+
+    term_list = [str(term) for term in list(formula)]
+
+    term_list_string = "torch.vstack(["
+    for term in term_list:
+        if term == "1":
+            term_list_string += f"torch.ones_like({domain.inputs.get_keys()[0]}), "
+        else:
+            term_list_string += term.replace(":", "*") + ", "
+    term_list_string += "]).T"
+
+    return term_list_string
+
+
 def linear_formula(
     domain: Optional[Domain],
 ) -> str:
@@ -224,10 +267,9 @@ def constraints_as_scipy_constraints(
             NonlinearInequalityConstraint,
         ):
             fun, lb, ub = get_constraint_function_and_bounds(c, domain, n_experiments)
-            if c.jacobian_expression is not None:
-                constraints.append(NonlinearConstraint(fun, lb, ub, jac=fun.jacobian))
-            else:
-                constraints.append(NonlinearConstraint(fun, lb, ub))
+            constraints.append(
+                NonlinearConstraint(fun, lb, ub, jac=fun.jacobian, hess=fun.hessian)
+            )
 
         elif isinstance(c, NChooseKConstraint):
             if ignore_nchoosek:
@@ -238,7 +280,9 @@ def constraints_as_scipy_constraints(
                     domain,
                     n_experiments,
                 )
-                constraints.append(NonlinearConstraint(fun, lb, ub, jac=fun.jacobian))
+                constraints.append(
+                    NonlinearConstraint(fun, lb, ub, jac=fun.jacobian, hess=fun.hessian)
+                )
 
         elif isinstance(c, InterpointEqualityConstraint):
             A, lb, ub = get_constraint_function_and_bounds(c, domain, n_experiments)
@@ -400,26 +444,49 @@ def __call__(self, x: np.ndarray) -> np.ndarray:
         violation[np.abs(violation) < 0] = 0
         return violation
 
-    def jacobian(self, x: np.ndarray) -> np.ndarray:
-        """Call constraint gradient with flattened numpy array."""
+    def jacobian(self, x: np.ndarray, sparse: bool = False) -> np.ndarray:
+        """Call constraint gradient with flattened numpy array.  If sparse is set to True, the output is a vector containing the entries of the sparse matrix representation of the jacobian."""
         x = pd.DataFrame(x.reshape(len(x) // self.D, self.D), columns=self.names)  # type: ignore
         gradient_compressed = self.constraint.jacobian(x).to_numpy()  # type: ignore
 
-        jacobian = np.zeros(shape=(self.n_experiments, self.D * self.n_experiments))
-        rows = np.repeat(
-            np.arange(self.n_experiments),
-            len(self.constraint_feature_indices),
-        )
         cols = np.repeat(
             self.D * np.arange(self.n_experiments),
             len(self.constraint_feature_indices),
         ).reshape((self.n_experiments, len(self.constraint_feature_indices)))
         cols = (cols + self.constraint_feature_indices).flatten()
 
+        if sparse:
+            jacobian = np.zeros(shape=(self.D * self.n_experiments))
+            jacobian[cols] = gradient_compressed.flatten()
+            return jacobian
+
+        rows = np.repeat(
+            np.arange(self.n_experiments),
+            len(self.constraint_feature_indices),
+        )
+        jacobian = np.zeros(shape=(self.n_experiments, self.D * self.n_experiments))
         jacobian[rows, cols] = gradient_compressed.flatten()
 
         return jacobian
 
+    def hessian(self, x: np.ndarray, *args):
+        """Call constraint hessian with flattened numpy array."""
+        x = pd.DataFrame(x.reshape(len(x) // self.D, self.D), columns=self.names)  # type: ignore
+        hessian_dict = self.constraint.hessian(x)  # type: ignore
+
+        hessian = np.zeros(
+            shape=(self.D * self.n_experiments, self.D * self.n_experiments)
+        )
+
+        cols, rows = np.meshgrid(
+            self.constraint_feature_indices,
+            self.constraint_feature_indices,
+        )
+        for i, idx in enumerate(hessian_dict.keys()):
+            hessian[i * self.D + cols, i * self.D + rows] = hessian_dict[idx]
+
+        return hessian
+
 
 def check_nchoosek_constraints_as_bounds(domain: Domain) -> None:
     """Checks if NChooseK constraints of domain can be formulated as bounds.
@@ -513,3 +580,66 @@ def nchoosek_constraints_as_bounds(
     bounds = [(b[0], b[1]) for b in bounds]
 
     return bounds
+
+
+def _minimize(
+    objective_function: Objective,
+    x0: np.ndarray,
+    bounds: List[Tuple[float, float]],
+    constraints: Optional[List[Union[NonlinearConstraint, LinearConstraint]]],
+    ipopt_options: dict,
+    use_hessian: bool,
+    use_cyipopt: Optional[bool] = CYIPOPT_AVAILABLE,
+) -> np.ndarray:
+    """Minimize the objective function using the given constraints and bounds.
+    Uses Ipopt if available, otherwise uses SLSQP.
+
+    Args:
+        objective_function (Objective): Objective function to minimize.
+        x0 (np.ndarray): Initial guess for the minimization.
+        bounds (List[Tuple[float, float]]): Bounds for the decision variables.
+        constraints (Optional[List[Union[NonlinearConstraint, LinearConstraint]]]): Constraints for the optimization problem.
+        ipopt_options (dict): Options for Ipopt solver. If Ipopt is not available, only the fields "max_iter" and "print_level" of this argument are used.
+        use_hessian (bool): Use hessian if set to True.
+        use_cyipopt (bool): Use cyipopt if set to True. Defaults to true if cyipopt is available.
+
+    Returns:
+        np.ndarray: The optimized design as flattened numpy array.
+    """
+    if use_cyipopt is None:
+        use_cyipopt = CYIPOPT_AVAILABLE
+
+    if use_cyipopt:
+        if use_hessian:
+            problem = SecondOrderDoEProblem(
+                doe_objective=objective_function,
+                bounds=bounds,
+                constraints=constraints,
+            )
+        else:
+            problem = FirstOrderDoEProblem(
+                doe_objective=objective_function,
+                bounds=bounds,
+                constraints=constraints,
+            )
+        for key in ipopt_options.keys():
+            problem.add_option(key, ipopt_options[key])
+
+        x, info = problem.solve(x0)
+        return x
+    else:
+        options = {}
+        if "max_iter" in ipopt_options.keys():
+            options["maxiter"] = ipopt_options["max_iter"]
+        if "print_level" in ipopt_options.keys():
+            options["disp"] = ipopt_options["print_level"]
+        result = opt.minimize(
+            fun=objective_function.evaluate,
+            x0=x0,
+            bounds=bounds,
+            options=options,
+            constraints=standardize_constraints(constraints, x0, "SLSQP"),
+            jac=objective_function.evaluate_jacobian,
+            hess=objective_function.evaluate_hessian if use_hessian else None,
+        )
+        return result.x
diff --git a/bofire/strategies/doe_strategy.py b/bofire/strategies/doe_strategy.py
index dfa3842fa..faf564dec 100644
--- a/bofire/strategies/doe_strategy.py
+++ b/bofire/strategies/doe_strategy.py
@@ -122,6 +122,8 @@ def _ask(self, candidate_count: PositiveInt) -> pd.DataFrame:  # type: ignore
                 partially_fixed_experiments=adapted_partially_fixed_candidates,
                 ipopt_options=self.data_model.ipopt_options,
                 criterion=self.data_model.criterion,
+                use_hessian=self.data_model.use_hessian,
+                use_cyipopt=self.data_model.use_cyipopt,
                 sampling=self._sampling,
             )
         # TODO adapt to when exhaustive search accepts discrete variables
@@ -139,6 +141,8 @@ def _ask(self, candidate_count: PositiveInt) -> pd.DataFrame:  # type: ignore
                 discrete_variables=new_discretes,
                 ipopt_options=self.data_model.ipopt_options,
                 criterion=self.data_model.criterion,
+                use_hessian=self.data_model.use_hessian,
+                use_cyipopt=self.data_model.use_cyipopt,
                 sampling=self._sampling,
             )
         elif self.data_model.optimization_strategy in [
@@ -156,6 +160,8 @@ def _ask(self, candidate_count: PositiveInt) -> pd.DataFrame:  # type: ignore
                 discrete_variables=new_discretes,
                 ipopt_options=self.data_model.ipopt_options,
                 criterion=self.data_model.criterion,
+                use_hessian=self.data_model.use_hessian,
+                use_cyipopt=self.data_model.use_cyipopt,
                 sampling=self._sampling,
             )
         elif self.data_model.optimization_strategy == "iterative":
@@ -181,6 +187,8 @@ def _ask(self, candidate_count: PositiveInt) -> pd.DataFrame:  # type: ignore
                     discrete_variables=new_discretes,
                     ipopt_options=self.data_model.ipopt_options,
                     criterion=self.data_model.criterion,
+                    use_hessian=self.data_model.use_hessian,
+                    use_cyipopt=self.data_model.use_cyipopt,
                     sampling=self._sampling,
                 )
                 adapted_partially_fixed_candidates = pd.concat(
diff --git a/tests/bofire/data_models/constraints/test_linear.py b/tests/bofire/data_models/constraints/test_linear.py
index be67e0ae2..fb977caa2 100644
--- a/tests/bofire/data_models/constraints/test_linear.py
+++ b/tests/bofire/data_models/constraints/test_linear.py
@@ -1,3 +1,6 @@
+import numpy as np
+import pandas as pd
+
 import tests.bofire.data_models.specs.api as specs
 from bofire.data_models.constraints.api import LinearInequalityConstraint
 
@@ -50,3 +53,44 @@ def test_as_smaller_equal():
     assert c.rhs == rhs
     assert coefficients == c.coefficients
     assert c.features == features
+
+
+def test_hessian():
+    spec = specs.constraints.valid(LinearInequalityConstraint).spec()
+    c = LinearInequalityConstraint.from_smaller_equal(
+        features=spec["features"],
+        coefficients=spec["coefficients"],
+        rhs=spec["rhs"],
+    )
+    experiments = pd.DataFrame(
+        np.random.rand(10, len(spec["features"])), columns=spec["features"]
+    )
+
+    for i, key in enumerate(c.hessian(experiments).keys()):
+        assert key == i
+    for v in c.hessian(experiments).values():
+        assert v == 0.0
+
+
+def test_jacobian():
+    spec = specs.constraints.valid(LinearInequalityConstraint).spec()
+    c = LinearInequalityConstraint.from_smaller_equal(
+        features=spec["features"],
+        coefficients=spec["coefficients"],
+        rhs=spec["rhs"],
+    )
+    experiments = pd.DataFrame(
+        np.random.rand(10, len(spec["features"])), columns=spec["features"]
+    )
+
+    res = c.jacobian(experiments)
+    assert list(res.columns) == [f"dg/d{name}" for name in spec["features"]]
+    for i, idx in enumerate(res.index):
+        assert idx == i
+    assert res.shape == (experiments.shape[0], len(spec["features"]))
+    for i in range(experiments.shape[0]):
+        np.allclose(
+            res.loc[i],
+            np.array(spec["coefficients"])
+            / np.linalg.norm(np.array(spec["coefficients"])),
+        )
diff --git a/tests/bofire/data_models/constraints/test_nonlinear.py b/tests/bofire/data_models/constraints/test_nonlinear.py
index fad3655f3..93d425b5c 100644
--- a/tests/bofire/data_models/constraints/test_nonlinear.py
+++ b/tests/bofire/data_models/constraints/test_nonlinear.py
@@ -45,9 +45,64 @@ def test_nonlinear_constraints_jacobian_expression():
     assert np.allclose(constraint2.jacobian(data), constraint0.jacobian(data))
     assert np.allclose(constraint3.jacobian(data), constraint0.jacobian(data))
 
-    with pytest.raises(ValueError):
+    with pytest.raises(
+        ValueError,
+        match="Features must be None if expression is a callable. They will be inferred from the callable.",
+    ):
         NonlinearInequalityConstraint(
             expression=lambda x1, x2, x3: x1**2 + x2**2 - x3,
             features=["x1", "x2", "x3"],
             jacobian_expression=None,
         )
+
+
+@pytest.mark.skipif(not SYMPY_AVAILABLE, reason="requires rdkit")
+@pytest.mark.skipif(not TORCH_AVAILABLE, reason="requires torch")
+def test_nonlinear_constraints_hessian_expression():
+    import torch
+
+    constraint0 = NonlinearInequalityConstraint(
+        expression="x1**3 + x2**2 - x3",
+        features=["x1", "x2", "x3"],
+    )
+    constraint1 = NonlinearInequalityConstraint(
+        expression="x1**3 + x2**2 - x3",
+        features=["x1", "x2", "x3"],
+        jacobian_expression="[3*x1**2, 2*x2, -1]",
+        hessian_expression="[[6*x1, 0, 0], [0, 2, 0], [0, 0, 0]]",
+    )
+
+    data = pd.DataFrame(np.random.rand(10, 3), columns=["x1", "x2", "x3"])
+    for i in range(10):
+        assert np.allclose(constraint0.hessian(data)[i], constraint1.hessian(data)[i])
+
+    constraint2 = NonlinearInequalityConstraint(
+        expression=lambda x1, x2, x3: x1**3 + x2**2 - x3,
+    )
+    constraint3 = NonlinearInequalityConstraint(
+        expression=lambda x1, x2, x3: x1**3 + x2**2 - x3,
+        jacobian_expression=lambda x1, x2, x3: [
+            3 * x1**2,
+            2 * x2,
+            -1.0 * torch.ones_like(x3),
+        ],
+        hessian_expression=lambda x1, x2, x3: [
+            [6 * x1, 0, 0],
+            [0, 2, 0],
+            [0, 0, 0],
+        ],
+    )
+
+    for i in range(10):
+        assert np.allclose(constraint2.hessian(data)[i], constraint0.hessian(data)[i])
+        assert np.allclose(constraint3.hessian(data)[i], constraint0.hessian(data)[i])
+
+    with pytest.raises(
+        ValueError,
+        match="Features must be None if expression is a callable. They will be inferred from the callable.",
+    ):
+        NonlinearInequalityConstraint(
+            expression=lambda x1, x2, x3: x1**2 + x2**2 - x3,
+            features=["x1", "x2", "x3"],
+            hessian_expression=None,
+        )
diff --git a/tests/bofire/data_models/specs/constraints.py b/tests/bofire/data_models/specs/constraints.py
index 0547fb57f..a3eee6a4d 100644
--- a/tests/bofire/data_models/specs/constraints.py
+++ b/tests/bofire/data_models/specs/constraints.py
@@ -47,6 +47,7 @@
     lambda: {
         "expression": "f1*f2",
         "jacobian_expression": "[f2,f1,0]",
+        "hessian_expression": "[[0,1,0],[1,0,0],[0,0,0]]",
         "features": ["f1", "f2", "f3"],
     },
 )
@@ -55,6 +56,7 @@
     lambda: {
         "expression": "f1*f2",
         "jacobian_expression": "[f2,f1,0]",
+        "hessian_expression": "[[0,1,0],[1,0,0],[0,0,0]]",
         "features": ["f1", "f2", "f3"],
     },
 )
diff --git a/tests/bofire/data_models/specs/strategies.py b/tests/bofire/data_models/specs/strategies.py
index 9f21a391c..fa8404718 100644
--- a/tests/bofire/data_models/specs/strategies.py
+++ b/tests/bofire/data_models/specs/strategies.py
@@ -235,33 +235,41 @@
             "relaxed",
             "iterative",
         ]:
-            specs.add_valid(
-                strategies.DoEStrategy,
-                lambda criterion=criterion,
-                formula=formula,
-                optimization_strategy=optimization_strategy: {
-                    "domain": domain.valid().obj().model_dump(),
-                    "optimization_strategy": optimization_strategy,
-                    "verbose": False,
-                    "seed": 42,
-                    "criterion": criterion(
-                        formula=formula, transform_range=None
-                    ).model_dump(),
-                },
-            )
-specs.add_valid(
-    strategies.DoEStrategy,
-    lambda: {
-        "domain": domain.valid().obj().dict(),
-        "optimization_strategy": "default",
-        "verbose": False,
-        "ipopt_options": {"maxiter": 200, "disp": 0},
-        "criterion": strategies.SpaceFillingCriterion(
-            sampling_fraction=0.3, transform_range=[-1, 1]
-        ).model_dump(),
-        "seed": 42,
-    },
-)
+            for use_cyipopt in [True, False, None]:
+                specs.add_valid(
+                    strategies.DoEStrategy,
+                    lambda criterion=criterion,
+                    formula=formula,
+                    optimization_strategy=optimization_strategy,
+                    use_cyipopt=use_cyipopt: {
+                        "domain": domain.valid().obj().model_dump(),
+                        "optimization_strategy": optimization_strategy,
+                        "verbose": False,
+                        "seed": 42,
+                        "criterion": criterion(
+                            formula=formula, transform_range=None
+                        ).model_dump(),
+                        "use_hessian": False,
+                        "use_cyipopt": use_cyipopt,
+                    },
+                )
+
+for use_cyipopt in [True, False, None]:
+    specs.add_valid(
+        strategies.DoEStrategy,
+        lambda use_cyipopt=use_cyipopt: {
+            "domain": domain.valid().obj().dict(),
+            "optimization_strategy": "default",
+            "verbose": False,
+            "ipopt_options": {"max_iter": 200, "print_level": 0},
+            "criterion": strategies.SpaceFillingCriterion(
+                sampling_fraction=0.3, transform_range=[-1, 1]
+            ).model_dump(),
+            "seed": 42,
+            "use_hessian": False,
+            "use_cyipopt": use_cyipopt,
+        },
+    )
 
 
 tempdomain = domain.valid().obj()
diff --git a/tests/bofire/strategies/doe/test_design.py b/tests/bofire/strategies/doe/test_design.py
index 8a8546b0c..9e4ce4311 100644
--- a/tests/bofire/strategies/doe/test_design.py
+++ b/tests/bofire/strategies/doe/test_design.py
@@ -170,7 +170,7 @@ def test_find_local_max_ipopt_mixed_results():
         domain,
         n_experiments=N,
         criterion=DOptimalityCriterion(formula="fully-quadratic"),
-        ipopt_options={"maxiter": 100},
+        ipopt_options={"max_iter": 100},
     )
     opt = np.eye(3)
     for row in A.to_numpy():
@@ -254,7 +254,7 @@ def test_find_local_max_ipopt_batch_constraint():
     result = find_local_max_ipopt(
         domain,
         criterion=DOptimalityCriterion(formula="linear"),
-        ipopt_options={"maxiter": 100},
+        ipopt_options={"max_iter": 100},
         n_experiments=30,
     )
 
@@ -380,7 +380,7 @@ def test_find_local_max_ipopt_fixed_experiments():
         domain,
         n_experiments=num_exp,
         criterion=DOptimalityCriterion(formula="fully-quadratic"),
-        ipopt_options={"maxiter": 100},
+        ipopt_options={"max_iter": 100},
         fixed_experiments=pd.DataFrame(
             [[1, 0, 0], [0, 1, 0]],
             columns=["x1", "x2", "x3"],
@@ -541,7 +541,7 @@ def test_find_local_max_ipopt_nonlinear_constraint():
         domain,
         num_exp,
         DOptimalityCriterion(formula="linear"),
-        ipopt_options={"maxiter": 100},
+        ipopt_options={"max_iter": 100},
     )
 
     assert np.allclose(domain.constraints(result), 0, atol=1e-6)
diff --git a/tests/bofire/strategies/doe/test_doe_problem.py b/tests/bofire/strategies/doe/test_doe_problem.py
new file mode 100644
index 000000000..eccc5e3df
--- /dev/null
+++ b/tests/bofire/strategies/doe/test_doe_problem.py
@@ -0,0 +1,146 @@
+import importlib.util
+
+import numpy as np
+import pytest
+
+from bofire.data_models.constraints.api import (
+    LinearInequalityConstraint,
+    NonlinearInequalityConstraint,
+)
+from bofire.data_models.domain.api import Domain
+from bofire.data_models.enum import SamplingMethodEnum
+from bofire.data_models.features.api import ContinuousInput, ContinuousOutput
+from bofire.strategies.doe.doe_problem import (
+    FirstOrderDoEProblem,
+    SecondOrderDoEProblem,
+)
+from bofire.strategies.doe.objective import DOptimalityCriterion, get_objective_function
+from bofire.strategies.doe.utils import (
+    constraints_as_scipy_constraints,
+    nchoosek_constraints_as_bounds,
+)
+
+
+CYIPOPT_AVAILABLE = importlib.util.find_spec("cyipopt") is not None
+
+
+@pytest.mark.skipif(not CYIPOPT_AVAILABLE, reason="requires cyipopt")
+def test_FirstOrderDoEProblem():
+    n_experiments = 4
+    criterion = DOptimalityCriterion(formula="linear")
+    domain = Domain.from_lists(
+        inputs=[
+            ContinuousInput(key="x1", bounds=(0, 1)),
+            ContinuousInput(key="x2", bounds=(0, 1)),
+            ContinuousInput(key="x3", bounds=(0, 1)),
+        ],
+        outputs=[ContinuousOutput(key="y")],
+        constraints=[
+            NonlinearInequalityConstraint(
+                expression="x1**2 + x2**2 - x3",
+                features=["x1", "x2", "x3"],
+                jacobian_expression="[2*x1,2*x2,-1]",
+            ),
+            LinearInequalityConstraint(
+                features=["x1", "x2", "x3"],
+                coefficients=[1, 1, 1],
+                rhs=1,
+            ),
+        ],
+    )
+
+    objective_function = get_objective_function(
+        criterion, domain=domain, n_experiments=n_experiments
+    )
+    assert objective_function is not None, "Criterion type is not supported!"
+
+    x0 = (
+        domain.inputs.sample(n=n_experiments, method=SamplingMethodEnum.UNIFORM)
+        .to_numpy()
+        .flatten()
+    )
+
+    constraints = constraints_as_scipy_constraints(
+        domain,
+        n_experiments,
+        ignore_nchoosek=True,
+    )
+
+    bounds = nchoosek_constraints_as_bounds(domain, n_experiments)
+
+    problem = FirstOrderDoEProblem(
+        doe_objective=objective_function,
+        bounds=bounds,
+        constraints=constraints,
+    )
+
+    assert len(problem.linear_constraints) == 1
+    assert len(problem.nonlinear_constraints) == 1
+    assert problem.A.shape == (4, 12)
+    assert np.allclose(problem.row, np.repeat([0, 1, 2, 3, 4, 5, 6, 7], 3))
+    assert np.allclose(problem.col, np.tile([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], 2))
+
+    with pytest.raises(AttributeError):
+        problem.hessian(x0, [0.0, 0.0], 1.0)
+    with pytest.raises(AttributeError):
+        problem.hessianstructure()
+
+
+@pytest.mark.skipif(not CYIPOPT_AVAILABLE, reason="requires cyipopt")
+def test_SecondOrderDoEProblem():
+    n_experiments = 4
+    criterion = DOptimalityCriterion(formula="linear")
+    domain = Domain.from_lists(
+        inputs=[
+            ContinuousInput(key="x1", bounds=(0, 1)),
+            ContinuousInput(key="x2", bounds=(0, 1)),
+            ContinuousInput(key="x3", bounds=(0, 1)),
+        ],
+        outputs=[ContinuousOutput(key="y")],
+        constraints=[
+            NonlinearInequalityConstraint(
+                expression="x1**2 + x2**2 - x3",
+                features=["x1", "x2", "x3"],
+                jacobian_expression="[2*x1,2*x2,-1]",
+            ),
+            LinearInequalityConstraint(
+                features=["x1", "x2", "x3"],
+                coefficients=[1, 1, 1],
+                rhs=1,
+            ),
+        ],
+    )
+
+    objective_function = get_objective_function(
+        criterion, domain=domain, n_experiments=n_experiments
+    )
+    assert objective_function is not None, "Criterion type is not supported!"
+
+    x0 = (
+        domain.inputs.sample(n=n_experiments, method=SamplingMethodEnum.UNIFORM)
+        .to_numpy()
+        .flatten()
+    )
+
+    constraints = constraints_as_scipy_constraints(
+        domain,
+        n_experiments,
+        ignore_nchoosek=True,
+    )
+
+    bounds = nchoosek_constraints_as_bounds(domain, n_experiments)
+
+    problem = SecondOrderDoEProblem(
+        doe_objective=objective_function,
+        bounds=bounds,
+        constraints=constraints,
+    )
+
+    assert len(problem.linear_constraints) == 1
+    assert len(problem.nonlinear_constraints) == 1
+    assert problem.A.shape == (4, 12)
+    assert np.allclose(problem.row, np.repeat([0, 1, 2, 3, 4, 5, 6, 7], 3))
+    assert np.allclose(problem.col, np.tile([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], 2))
+
+    problem.hessian(x0, np.array([0.0, 0.0]), 1.0)
+    problem.hessianstructure()
diff --git a/tests/bofire/strategies/doe/test_objective.py b/tests/bofire/strategies/doe/test_objective.py
index 598fde874..c6b09a30d 100644
--- a/tests/bofire/strategies/doe/test_objective.py
+++ b/tests/bofire/strategies/doe/test_objective.py
@@ -1,7 +1,6 @@
 import importlib
 
 import numpy as np
-import pandas as pd
 import pytest
 from formulaic import Formula
 
@@ -25,7 +24,6 @@
     EOptimality,
     GOptimality,
     IOptimality,
-    ModelBasedObjective,
     SpaceFilling,
     get_objective_function,
 )
@@ -35,375 +33,6 @@
 CYIPOPT_AVAILABLE = importlib.util.find_spec("cyipopt") is not None
 
 
-def test_Objective_model_jacobian_t():
-    # "small" model
-    domain = Domain.from_lists(
-        inputs=[
-            ContinuousInput(
-                key=f"x{i + 1}",
-                bounds=(0, 1),
-            )
-            for i in range(3)
-        ],
-        outputs=[ContinuousOutput(key="y")],
-    )
-
-    vars = domain.inputs.get_keys()
-    f = Formula("x1 + x2 + x3 + x1:x2 + {x3**2}")
-    x = np.array([[1, 2, 3]])
-
-    objective = ModelBasedObjective(
-        domain=domain,
-        model=f,
-        n_experiments=1,
-    )
-    model_jacobian_t = objective._model_jacobian_t
-
-    B = np.zeros(shape=(3, 6))
-    B[:, 1:4] = np.eye(3)
-    B[:, 4] = np.array([0, 0, 6])
-    B[:, 5] = np.array([2, 1, 0])
-    assert np.allclose(B, model_jacobian_t(x))
-
-    # fully quadratic model
-    f = Formula("x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + {x1**2} + {x2**2} + {x3**2}")
-    model_terms = np.array(f, dtype=str)
-    x = np.array([[1, 2, 3]])
-
-    objective = ModelBasedObjective(
-        domain=domain,
-        model=f,
-        n_experiments=1,
-    )
-    model_jacobian_t = objective._model_jacobian_t
-    B = np.zeros(shape=(3, 10))
-    B[:, 1:4] = np.eye(3)
-    B[:, 4:7] = 2 * np.diag(x[0])
-    B[:, 7:] = np.array([[2, 1, 0], [3, 0, 1], [0, 3, 2]]).T
-    B = pd.DataFrame(
-        B,
-        columns=[
-            "1",
-            "x1",
-            "x2",
-            "x3",
-            "x1 ** 2",
-            "x2 ** 2",
-            "x3 ** 2",
-            "x1:x2",
-            "x1:x3",
-            "x2:x3",
-        ],
-    )
-    B = B[model_terms].to_numpy()
-
-    assert np.allclose(B, model_jacobian_t(x)[0])
-
-    # fully cubic model
-    domain = Domain.from_lists(
-        inputs=[
-            ContinuousInput(
-                key=f"x{i + 1}",
-                bounds=(0, 1),
-            )
-            for i in range(5)
-        ],
-        outputs=[ContinuousOutput(key="y")],
-    )
-    vars = ["x1", "x2", "x3", "x4", "x5"]
-    n_vars = len(vars)
-
-    formula = ""
-    for name in vars:
-        formula += name + " + "
-
-    for name in vars:
-        formula += "{" + name + "**2} + "
-    for i in range(n_vars):
-        for j in range(i + 1, n_vars):
-            term = str(Formula(vars[j] + ":" + vars[i] + "-1")) + " + "
-            formula += term
-
-    for name in vars:
-        formula += "{" + name + "**3} + "
-    for i in range(n_vars):
-        for j in range(i + 1, n_vars):
-            for k in range(j + 1, n_vars):
-                term = (
-                    str(Formula(vars[k] + ":" + vars[j] + ":" + vars[i] + "-1")) + " + "
-                )
-                formula += term
-    f = Formula(formula[:-3])
-    x = np.array([[1, 2, 3, 4, 5]])
-    objective = ModelBasedObjective(
-        domain=domain,
-        model=f,
-        n_experiments=1,
-    )
-    model_jacobian_t = objective._model_jacobian_t
-
-    B = np.array(
-        [
-            [
-                0.0,
-                1.0,
-                2.0,
-                3.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                2.0,
-                3.0,
-                4.0,
-                5.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                6.0,
-                8.0,
-                10.0,
-                12.0,
-                15.0,
-                20.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-            ],
-            [
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                4.0,
-                12.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                0.0,
-                0.0,
-                0.0,
-                3.0,
-                4.0,
-                5.0,
-                0.0,
-                0.0,
-                0.0,
-                3.0,
-                4.0,
-                5.0,
-                0.0,
-                0.0,
-                0.0,
-                12.0,
-                15.0,
-                20.0,
-                0.0,
-            ],
-            [
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                6.0,
-                27.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                0.0,
-                0.0,
-                2.0,
-                0.0,
-                0.0,
-                4.0,
-                5.0,
-                0.0,
-                2.0,
-                0.0,
-                0.0,
-                4.0,
-                5.0,
-                0.0,
-                8.0,
-                10.0,
-                0.0,
-                20.0,
-            ],
-            [
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                8.0,
-                48.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                0.0,
-                0.0,
-                2.0,
-                0.0,
-                3.0,
-                0.0,
-                5.0,
-                0.0,
-                2.0,
-                0.0,
-                3.0,
-                0.0,
-                5.0,
-                6.0,
-                0.0,
-                10.0,
-                15.0,
-            ],
-            [
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                10.0,
-                75.0,
-                0.0,
-                0.0,
-                0.0,
-                1.0,
-                0.0,
-                0.0,
-                2.0,
-                0.0,
-                3.0,
-                4.0,
-                0.0,
-                0.0,
-                2.0,
-                0.0,
-                3.0,
-                4.0,
-                0.0,
-                6.0,
-                8.0,
-                12.0,
-            ],
-        ],
-    )
-
-    B = pd.DataFrame(
-        B,
-        columns=[
-            "1",
-            "x1",
-            "x1 ** 2",
-            "x1 ** 3",
-            "x2",
-            "x2 ** 2",
-            "x2 ** 3",
-            "x3",
-            "x3 ** 2",
-            "x3 ** 3",
-            "x4",
-            "x4 ** 2",
-            "x4 ** 3",
-            "x5",
-            "x5 ** 2",
-            "x5 ** 3",
-            "x2:x1",
-            "x3:x1",
-            "x4:x1",
-            "x5:x1",
-            "x3:x2",
-            "x4:x2",
-            "x5:x2",
-            "x4:x3",
-            "x5:x3",
-            "x5:x4",
-            "x3:x2:x1",
-            "x4:x2:x1",
-            "x5:x2:x1",
-            "x4:x3:x1",
-            "x5:x3:x1",
-            "x5:x4:x1",
-            "x4:x3:x2",
-            "x5:x3:x2",
-            "x5:x4:x2",
-            "x5:x4:x3",
-        ],
-    )
-
-    B = B[objective.model_terms].to_numpy()
-
-    assert np.allclose(B, model_jacobian_t(x)[0])
-
-
-def test_Objective_convert_input_to_model_tensor():
-    domain = Domain.from_lists(
-        inputs=[
-            ContinuousInput(
-                key=f"x{i + 1}",
-                bounds=(0, 1),
-            )
-            for i in range(3)
-        ],
-        outputs=[ContinuousOutput(key="y")],
-    )
-    model = get_formula_from_string("linear", domain=domain)
-
-    d_optimality = DOptimality(domain=domain, model=model, n_experiments=3)
-    x = np.array([1, 0, 0, 0, 2, 0, 0, 0, 3])
-    print(domain.inputs)
-    X = d_optimality._convert_input_to_model_tensor(x).detach().numpy()
-    assert np.allclose(X, np.array([[1, 1, 0, 0], [1, 0, 2, 0], [1, 0, 0, 3]]))
-
-
 def test_DOptimality_instantiation():
     # default jacobian building block
     domain = Domain.from_lists(
@@ -417,11 +46,11 @@ def test_DOptimality_instantiation():
         outputs=[ContinuousOutput(key="y")],
     )
 
-    model = Formula("x1 + x2 + x3 + x1:x2 + {x3**2}")
+    formula = Formula("x1 + x2 + x3 + x1:x2 + {x3**2}")
 
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=2,
     )
 
@@ -435,23 +64,20 @@ def test_DOptimality_instantiation():
         assert i.lower_bound == 0
     assert all(np.array(d_optimality.domain.outputs.get_keys()) == np.array(["y"]))
 
-    assert isinstance(d_optimality.model, Formula)
+    assert isinstance(d_optimality.formula, Formula)
     assert all(
-        np.array(d_optimality.model, dtype=str)
+        np.array(d_optimality.formula, dtype=str)
         == np.array(["1", "x1", "x2", "x3", "x3 ** 2", "x1:x2"]),
     )
 
-    x = np.array([[1, 2, 3], [1, 2, 3]])
-    B = np.zeros(shape=(3, 6))
-    B[:, 1:4] = np.eye(3)
-    B[:, 4] = np.array([0, 0, 6])
-    B[:, 5] = np.array([2, 1, 0])
-
-    assert np.allclose(B, d_optimality._model_jacobian_t(x))
     assert np.shape(
         d_optimality.evaluate_jacobian(np.array([[1, 1, 1], [2, 2, 2]]).flatten()),
     ) == (6,)
 
+    assert np.shape(
+        d_optimality.evaluate_hessian(np.array([[1, 1, 1], [2, 2, 2]]).flatten()),
+    ) == (6, 6)
+
     # 5th order model
     domain = Domain.from_lists(
         inputs=[
@@ -464,11 +90,11 @@ def test_DOptimality_instantiation():
         outputs=[ContinuousOutput(key="y")],
     )
 
-    model = Formula("{x1**5} + {x2**5} + {x3**5}")
+    formula = Formula("{x1**5} + {x2**5} + {x3**5}")
 
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=3,
     )
 
@@ -476,13 +102,18 @@ def test_DOptimality_instantiation():
     B = np.zeros(shape=(3, 4))
     B[:, 1:] = 5 * np.diag(x[0] ** 4)
 
-    assert np.allclose(B, d_optimality._model_jacobian_t(x))
     assert np.shape(
         d_optimality.evaluate_jacobian(
             np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]]).flatten(),
         ),
     ) == (9,)
 
+    assert np.shape(
+        d_optimality.evaluate_hessian(
+            np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]]).flatten()
+        ),
+    ) == (9, 9)
+
 
 def test_DOptimality_evaluate_jacobian():
     # n_experiment = 1, n_inputs = 2, model: x1 + x2
@@ -500,11 +131,11 @@ def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:  # type: ignore
         outputs=[ContinuousOutput(key="y")],
     )
 
-    model = Formula("x1 + x2 - 1")
+    formula = Formula("x1 + x2 - 1")
     n_experiments = 1
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=n_experiments,
         delta=1e-3,
     )
@@ -516,14 +147,14 @@ def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:  # type: ignore
             d_optimality.evaluate_jacobian(x), get_jacobian(x), rtol=1e-3
         )
 
-    # n_experiment = 1, n_inputs = 2, model: x1**2 + x2**2
+    # n_experiment = 1, n_inputs = 2, formula: x1**2 + x2**2
     def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:  # type: ignore
         return -4 * x**3 / (x[0] ** 4 + x[1] ** 4 + delta)
 
-    model = Formula("{x1**2} + {x2**2} - 1")
+    formula = Formula("{x1**2} + {x2**2} - 1")
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=n_experiments,
         delta=1e-3,
     )
@@ -534,7 +165,7 @@ def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:  # type: ignore
             d_optimality.evaluate_jacobian(x), get_jacobian(x), rtol=1e-3
         )
 
-    # n_experiment = 2, n_inputs = 2, model = x1 + x2
+    # n_experiment = 2, n_inputs = 2, formula = x1 + x2
     def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:
         X = x.reshape(2, 2)
 
@@ -574,11 +205,11 @@ def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:
 
         return y
 
-    model = Formula("x1 + x2 - 1")
+    formula = Formula("x1 + x2 - 1")
     n_experiments = 2
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=n_experiments,
         delta=1e-3,
     )
@@ -589,7 +220,7 @@ def get_jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:
             d_optimality.evaluate_jacobian(x), get_jacobian(x), rtol=1e-3
         )
 
-    # n_experiment = 2, n_inputs = 2, model = x1**2 + x2**2
+    # n_experiment = 2, n_inputs = 2, formula = x1**2 + x2**2
     def jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:
         X = x.reshape(2, 2) ** 2
 
@@ -629,10 +260,10 @@ def jacobian(x: np.ndarray, delta=1e-3) -> np.ndarray:
 
         return y
 
-    model = Formula("{x1**2} + {x2**2} - 1")
+    formula = Formula("{x1**2} + {x2**2} - 1")
     d_optimality = DOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=n_experiments,
         delta=1e-3,
     )
@@ -654,10 +285,12 @@ def test_DOptimality_evaluate():
         ],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    d_optimality = DOptimality(domain=domain, model=model, n_experiments=3)
-    x = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])
+    d_optimality = DOptimality(
+        domain=domain, formula=formula, n_experiments=3, delta=1e-7
+    )
+    x = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1], dtype=np.float64)
     assert np.allclose(d_optimality.evaluate(x), -np.log(4) - np.log(1e-7))
 
 
@@ -672,9 +305,9 @@ def test_AOptimality_evaluate():
         ],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    a_optimality = AOptimality(domain=domain, model=model, n_experiments=4)
+    a_optimality = AOptimality(domain=domain, formula=formula, n_experiments=4)
 
     x = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0])
     assert np.allclose(a_optimality.evaluate(x), 3 * 1.9999991 + 0.9999996)
@@ -685,9 +318,9 @@ def test_AOptimality_evaluate_jacobian():
         inputs=[ContinuousInput(key="x1", bounds=(0, 1))],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    a_optimality = AOptimality(domain=domain, model=model, n_experiments=2, delta=0)
+    a_optimality = AOptimality(domain=domain, formula=formula, n_experiments=2, delta=0)
 
     x = np.array([1, 0.5])
 
@@ -706,9 +339,9 @@ def test_EOptimality_evaluate():
         inputs=[ContinuousInput(key="x1", bounds=(0, 1))],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    e_optimality = EOptimality(domain=domain, model=model, n_experiments=2, delta=0)
+    e_optimality = EOptimality(domain=domain, formula=formula, n_experiments=2, delta=0)
 
     x = np.array([1, 0.5])
 
@@ -729,9 +362,9 @@ def test_EOptimality_evaluate_jacobian():
         inputs=[ContinuousInput(key="x1", bounds=(0, 1))],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    e_optimality = EOptimality(domain=domain, model=model, n_experiments=2, delta=0)
+    e_optimality = EOptimality(domain=domain, formula=formula, n_experiments=2, delta=0)
 
     x = np.array([1, 0.5])
 
@@ -754,9 +387,9 @@ def test_GOptimality_evaluate():
         inputs=[ContinuousInput(key="x1", bounds=(0, 1))],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    g_optimality = GOptimality(domain=domain, model=model, n_experiments=2, delta=0)
+    g_optimality = GOptimality(domain=domain, formula=formula, n_experiments=2, delta=0)
 
     x = np.array([1, 0.5])
 
@@ -770,9 +403,9 @@ def test_GOptimality_evaluate_jacobian():
         inputs=[ContinuousInput(key="x1", bounds=(0, 1))],
         outputs=[ContinuousOutput(key="y")],
     )
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
-    g_optimality = GOptimality(domain=domain, model=model, n_experiments=2, delta=0)
+    g_optimality = GOptimality(domain=domain, formula=formula, n_experiments=2, delta=0)
 
     x = np.array([1, 0.5])
 
@@ -875,7 +508,7 @@ def test_MinMaxTransform():
                 delta=0,
                 transform_range=None,
                 n_space_filling_points=4,
-                ipopt_options={"maxiter": 200},
+                ipopt_options={"max_iter": 200},
             ),
             domain=domain,
             n_experiments=4,
@@ -887,7 +520,7 @@ def test_MinMaxTransform():
                     delta=0,
                     transform_range=(-1.0, 1.0),
                     n_space_filling_points=4,
-                    ipopt_options={"maxiter": 200},
+                    ipopt_options={"max_iter": 200},
                 ),
                 domain=domain,
                 n_experiments=4,
@@ -902,11 +535,11 @@ def test_IOptimality_instantiation():
         outputs=[ContinuousOutput(key="y")],
     )
 
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
     i_optimality = IOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=2,
     )
     assert np.allclose(np.linspace(0, 1, 100), i_optimality.Y.to_numpy().flatten())
@@ -922,11 +555,11 @@ def test_IOptimality_instantiation():
         ],
     )
 
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
     i_optimality = IOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=2,
     )
 
@@ -944,11 +577,11 @@ def test_IOptimality_instantiation():
         ],
     )
 
-    model = get_formula_from_string("linear", domain=domain)
+    formula = get_formula_from_string("linear", domain=domain)
 
     i_optimality = IOptimality(
         domain=domain,
-        model=model,
+        formula=formula,
         n_experiments=2,
     )
 
diff --git a/tests/bofire/strategies/doe/test_utils.py b/tests/bofire/strategies/doe/test_utils.py
index dc45b8430..825ebac1c 100644
--- a/tests/bofire/strategies/doe/test_utils.py
+++ b/tests/bofire/strategies/doe/test_utils.py
@@ -1,3 +1,4 @@
+import importlib.util
 import sys
 
 import numpy as np
@@ -13,21 +14,28 @@
     NonlinearInequalityConstraint,
 )
 from bofire.data_models.domain.api import Domain
+from bofire.data_models.enum import SamplingMethodEnum
 from bofire.data_models.features.api import (
     ContinuousInput,
     ContinuousOutput,
     DiscreteInput,
 )
+from bofire.strategies.doe.objective import DOptimalityCriterion, get_objective_function
 from bofire.strategies.doe.utils import (
     ConstraintWrapper,
+    _minimize,
     check_nchoosek_constraints_as_bounds,
     constraints_as_scipy_constraints,
+    convert_formula_to_string,
     get_formula_from_string,
     n_zero_eigvals,
     nchoosek_constraints_as_bounds,
 )
 
 
+CYIPOPT_AVAILABLE = importlib.util.find_spec("cyipopt") is not None
+
+
 def get_formula_from_string_recursion_limit():
     # save recursion limit
     recursion_limit = sys.getrecursionlimit()
@@ -475,6 +483,99 @@ def test_ConstraintWrapper():
             ],
         ),
     )
+    assert np.allclose(
+        c.jacobian(x, sparse=True), np.array([2, 0, 0, 2, 1, 0, 0, 1, 6, 0, 0, 0])
+    )
+
+    assert np.allclose(
+        c.hessian(x),
+        [
+            [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2],
+        ],
+    )
+
+
+@pytest.mark.skipif(not CYIPOPT_AVAILABLE, reason="cyipopt required")
+def test_minimize():
+    # Run _minimize() with and without ipopt
+    n_experiments = 4
+    criterion = DOptimalityCriterion(formula="linear")
+    domain = Domain.from_lists(
+        inputs=[
+            ContinuousInput(key="x1", bounds=(0, 1)),
+            ContinuousInput(key="x2", bounds=(0, 1)),
+            ContinuousInput(key="x3", bounds=(0, 1)),
+        ],
+        outputs=[ContinuousOutput(key="y")],
+        constraints=[
+            LinearEqualityConstraint(
+                features=["x1", "x2", "x3"],
+                coefficients=[1, 1, 1],
+                rhs=1,
+            ),
+        ],
+    )
+
+    objective_function = get_objective_function(
+        criterion, domain=domain, n_experiments=n_experiments
+    )
+
+    x0 = (
+        domain.inputs.sample(n=n_experiments, method=SamplingMethodEnum.UNIFORM)
+        .to_numpy()
+        .flatten()
+    )
+    constraints = constraints_as_scipy_constraints(
+        domain,
+        n_experiments,
+        ignore_nchoosek=True,
+    )
+    bounds = nchoosek_constraints_as_bounds(domain, n_experiments)
+
+    result_ipopt = _minimize(
+        objective_function=objective_function,
+        x0=x0,
+        bounds=bounds,
+        constraints=constraints,
+        ipopt_options={"max_iter": 500, "print_level": 0},
+        use_hessian=False,
+        use_cyipopt=True,
+    )
+
+    result_scipy = _minimize(
+        objective_function=objective_function,
+        x0=x0,
+        bounds=bounds,
+        constraints=constraints,
+        ipopt_options={"max_iter": 500},
+        use_hessian=False,
+        use_cyipopt=False,
+    )
+
+    for i in range(n_experiments):
+        assert np.any(
+            [
+                np.allclose(result_ipopt[j], result_scipy[i])
+                for j in range(n_experiments)
+            ]
+        )
+        assert np.any(
+            [
+                np.allclose(result_scipy[j], result_ipopt[i])
+                for j in range(n_experiments)
+            ]
+        )
 
 
 def test_check_nchoosek_constraints_as_bounds():
@@ -652,3 +753,19 @@ def test_nchoosek_constraints_as_bounds():
     # assert len(bounds) == 20
     # for i in range(20):
     #     assert _bounds[i] == bounds[i]
+
+
+def test_convert_formula_to_string():
+    domain = Domain.from_lists(
+        inputs=[ContinuousInput(key=f"x{i}", bounds=(0, 1)) for i in range(3)],
+        outputs=[ContinuousOutput(key="y")],
+    )
+
+    formula = get_formula_from_string(domain=domain, model_type="fully-quadratic")
+
+    formula_str = convert_formula_to_string(domain=domain, formula=formula)
+    assert (
+        formula_str
+        == "torch.vstack([torch.ones_like(x0), x0, x1, x2, x0 ** 2, x1 ** 2, x2 ** 2,"
+        + " x0*x1, x0*x2, x1*x2, ]).T"
+    )
diff --git a/tests/bofire/strategies/test_doe.py b/tests/bofire/strategies/test_doe.py
index d342efb77..ceebd8e86 100644
--- a/tests/bofire/strategies/test_doe.py
+++ b/tests/bofire/strategies/test_doe.py
@@ -451,7 +451,12 @@ def test_functional_constraint():
     # Calculate the solid content of the formulation
     def calc_solid_content(A, B, T, W, W_T):
         # Ensure same order as in the dictionary containing the material properties
-        return np.array([A, B, T, W, W_T]).T @ (df_raw_materials["sc"].values)
+        if isinstance(A, torch.Tensor):
+            return torch.stack([A, B, T, W, W_T], 0).T @ torch.tensor(
+                df_raw_materials["sc"].values
+            )
+        else:
+            return np.array([A, B, T, W, W_T]).T @ df_raw_materials["sc"].values
 
     # Calculate the volume content of the formulation
     def calc_volume_content(A, B, T, W, W_T):
@@ -459,6 +464,11 @@ def calc_volume_content(A, B, T, W, W_T):
             A * raw_materials_data["A"][0] / raw_materials_data["A"][1]
             + B * raw_materials_data["B"][0] / raw_materials_data["B"][1]
         )
+        A = A
+        B = B
+        T = T
+        W = W
+        W_T = W_T
         volume_total = volume_solid + (1 - calc_solid_content(A, B, T, W, W_T) / 1)
         return volume_solid / volume_total
 
@@ -489,8 +499,13 @@ def calc_volume_content(A, B, T, W, W_T):
     data_model = data_models.DoEStrategy(
         domain=domain,
         criterion=DOptimalityCriterion(formula="linear"),
-        ipopt_options={"maxiter": 500},
+        ipopt_options={
+            "max_iter": 500,
+            "derivative_test": "first-order",
+            "print_level": 5,
+        },
         sampling=sampling,
+        use_cyipopt=True,
     )
     strategy = DoEStrategy(data_model=data_model)
 
diff --git a/tests/bofire/strategies/test_space_filling.py b/tests/bofire/strategies/test_space_filling.py
index 8d856b5d4..1acf45475 100644
--- a/tests/bofire/strategies/test_space_filling.py
+++ b/tests/bofire/strategies/test_space_filling.py
@@ -62,19 +62,22 @@ def test_ask(domain, num_samples):
         domain=domain,
         optimization_strategy="partially-random",
         criterion=SpaceFillingCriterion(),
-        ipopt_options={"maxiter": 300, "disp": 0},
+        ipopt_options={"max_iter": 300, "print_level": 0},
     )
     sampler = strategies.DoEStrategy(data_model=data_model)
     samples = sampler.ask(num_samples)
     assert len(samples) == num_samples
 
 
+test_ask(domain=domains[0], num_samples=1)
+
+
 def test_ask_pending_candidates():
     data_model = data_models.DoEStrategy(
         domain=domains[0],
         optimization_strategy="partially-random",
         criterion=SpaceFillingCriterion(),
-        ipopt_options={"maxiter": 300, "disp": 0},
+        ipopt_options={"max_iter": 300, "print_level": 0},
     )
     sampler = strategies.DoEStrategy(data_model=data_model)
     pending_candidates = sampler.ask(2, add_pending=True)
diff --git a/tutorials/doe/basic_examples.ipynb b/tutorials/doe/basic_examples.ipynb
index d13e9bb1a..75c447ef0 100644
--- a/tutorials/doe/basic_examples.ipynb
+++ b/tutorials/doe/basic_examples.ipynb
@@ -99,13 +99,22 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "\n",
-      "******************************************************************************\n",
-      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
-      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
-      "         For more information visit https://github.com/coin-or/Ipopt\n",
-      "******************************************************************************\n",
-      "\n"
+      "Tried to set Option: disp. It is not a valid option. Please check the list of available options.\n"
+     ]
+    },
+    {
+     "ename": "TypeError",
+     "evalue": "Error while assigning an option",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[2], line 29\u001b[0m\n\u001b[1;32m     23\u001b[0m data_model \u001b[38;5;241m=\u001b[39m DoEStrategy(\n\u001b[1;32m     24\u001b[0m     domain\u001b[38;5;241m=\u001b[39mdomain,\n\u001b[1;32m     25\u001b[0m     criterion\u001b[38;5;241m=\u001b[39mDOptimalityCriterion(formula\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlinear\u001b[39m\u001b[38;5;124m\"\u001b[39m),\n\u001b[1;32m     26\u001b[0m     ipopt_options\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdisp\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m0\u001b[39m},\n\u001b[1;32m     27\u001b[0m )\n\u001b[1;32m     28\u001b[0m strategy \u001b[38;5;241m=\u001b[39m strategies\u001b[38;5;241m.\u001b[39mmap(data_model\u001b[38;5;241m=\u001b[39mdata_model)\n\u001b[0;32m---> 29\u001b[0m candidates \u001b[38;5;241m=\u001b[39m \u001b[43mstrategy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mask\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcandidate_count\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m12\u001b[39;49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/strategy.py:128\u001b[0m, in \u001b[0;36mStrategy.ask\u001b[0;34m(self, candidate_count, add_pending, raise_validation_error)\u001b[0m\n\u001b[1;32m    123\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhas_sufficient_experiments():\n\u001b[1;32m    124\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    125\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNot enough experiments available to execute the strategy.\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    126\u001b[0m     )\n\u001b[0;32m--> 128\u001b[0m candidates \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ask\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcandidate_count\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcandidate_count\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    130\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdomain\u001b[38;5;241m.\u001b[39mvalidate_candidates(\n\u001b[1;32m    131\u001b[0m     candidates\u001b[38;5;241m=\u001b[39mcandidates,\n\u001b[1;32m    132\u001b[0m     only_inputs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m,\n\u001b[1;32m    133\u001b[0m     raise_validation_error\u001b[38;5;241m=\u001b[39mraise_validation_error,\n\u001b[1;32m    134\u001b[0m )\n\u001b[1;32m    136\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m candidate_count \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/doe_strategy.py:113\u001b[0m, in \u001b[0;36mDoEStrategy._ask\u001b[0;34m(self, candidate_count)\u001b[0m\n\u001b[1;32m    103\u001b[0m num_discrete_vars \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlen\u001b[39m(new_discretes)\n\u001b[1;32m    104\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m (\n\u001b[1;32m    105\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata_model\u001b[38;5;241m.\u001b[39moptimization_strategy \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrelaxed\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    106\u001b[0m     \u001b[38;5;129;01mor\u001b[39;00m (num_binary_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m num_discrete_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m)\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    111\u001b[0m     )\n\u001b[1;32m    112\u001b[0m ):\n\u001b[0;32m--> 113\u001b[0m     design \u001b[38;5;241m=\u001b[39m \u001b[43mfind_local_max_ipopt\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    114\u001b[0m \u001b[43m        \u001b[49m\u001b[43mnew_domain\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    115\u001b[0m \u001b[43m        \u001b[49m\u001b[43mn_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m_candidate_count\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    116\u001b[0m \u001b[43m        \u001b[49m\u001b[43mfixed_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    117\u001b[0m \u001b[43m        \u001b[49m\u001b[43mpartially_fixed_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43madapted_partially_fixed_candidates\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    118\u001b[0m \u001b[43m        \u001b[49m\u001b[43mipopt_options\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mipopt_options\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    119\u001b[0m \u001b[43m        \u001b[49m\u001b[43mcriterion\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcriterion\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    120\u001b[0m \u001b[43m        \u001b[49m\u001b[43muse_hessian\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43muse_hessian\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    121\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    122\u001b[0m \u001b[38;5;66;03m# TODO adapt to when exhaustive search accepts discrete variables\u001b[39;00m\n\u001b[1;32m    123\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m (\n\u001b[1;32m    124\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata_model\u001b[38;5;241m.\u001b[39moptimization_strategy \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexhaustive\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    125\u001b[0m     \u001b[38;5;129;01mand\u001b[39;00m num_discrete_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[1;32m    126\u001b[0m ):\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/doe/design.py:188\u001b[0m, in \u001b[0;36mfind_local_max_ipopt\u001b[0;34m(domain, n_experiments, criterion, ipopt_options, sampling, fixed_experiments, partially_fixed_experiments, use_hessian)\u001b[0m\n\u001b[1;32m    182\u001b[0m     problem \u001b[38;5;241m=\u001b[39m FirstOrderDoEProblem(\n\u001b[1;32m    183\u001b[0m         doe_objective\u001b[38;5;241m=\u001b[39mobjective_function,\n\u001b[1;32m    184\u001b[0m         bounds\u001b[38;5;241m=\u001b[39mbounds,\n\u001b[1;32m    185\u001b[0m         constraints\u001b[38;5;241m=\u001b[39mconstraints,\n\u001b[1;32m    186\u001b[0m     )\n\u001b[1;32m    187\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m key \u001b[38;5;129;01min\u001b[39;00m _ipopt_options\u001b[38;5;241m.\u001b[39mkeys():\n\u001b[0;32m--> 188\u001b[0m     \u001b[43mproblem\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_option\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m_ipopt_options\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    190\u001b[0m x, info \u001b[38;5;241m=\u001b[39m problem\u001b[38;5;241m.\u001b[39msolve(x0)\n\u001b[1;32m    192\u001b[0m design \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame(\n\u001b[1;32m    193\u001b[0m     x\u001b[38;5;241m.\u001b[39mreshape(n_experiments, \u001b[38;5;28mlen\u001b[39m(domain\u001b[38;5;241m.\u001b[39minputs)),\n\u001b[1;32m    194\u001b[0m     columns\u001b[38;5;241m=\u001b[39mdomain\u001b[38;5;241m.\u001b[39minputs\u001b[38;5;241m.\u001b[39mget_keys(),\n\u001b[1;32m    195\u001b[0m     index\u001b[38;5;241m=\u001b[39m[\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexp\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mi\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n_experiments)],\n\u001b[1;32m    196\u001b[0m )\n",
+      "File \u001b[0;32m~/anaconda3/envs/bofire/lib/python3.11/site-packages/cyipopt/cython/ipopt_wrapper.pyx:495\u001b[0m, in \u001b[0;36mipopt_wrapper.Problem.add_option\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;31mTypeError\u001b[0m: Error while assigning an option"
      ]
     }
    ],
@@ -135,7 +144,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"linear\"),\n",
-    "    ipopt_options={\"disp\": 0},\n",
+    "    ipopt_options={\"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "candidates = strategy.ask(candidate_count=12)"
@@ -260,7 +269,7 @@
     "    criterion=DOptimalityCriterion(\n",
     "        formula=\"x1 + x2 + x3 + {x1**2} + {x2**2} + {x3**2} + {x1**3} + {x2**3} + {x3**3} + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3\"\n",
     "    ),\n",
-    "    ipopt_options={\"disp\": 0},\n",
+    "    ipopt_options={\"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "candidates = strategy.ask(12)"
@@ -419,9 +428,9 @@
     "    criterion=IOptimalityCriterion(\n",
     "        formula=\"x1 + x2 + x3 + {x1**2} + {x2**2} + {x3**2} + {x1**3} + {x2**3} + {x3**3} + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3\",\n",
     "        n_space_filling_points=60,\n",
-    "        ipopt_options={\"maxiter\": 500},\n",
+    "        ipopt_options={\"max_iter\": 500},\n",
     "    ),\n",
-    "    ipopt_options={\"disp\": 0},\n",
+    "    ipopt_options={\"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "candidates = strategy.ask(12).to_numpy().T\n",
@@ -628,7 +637,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"linear\"),\n",
-    "    ipopt_options={\"maxiter\": 100, \"disp\": 0},\n",
+    "    ipopt_options={\"max_iter\": 100, \"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "result = strategy.ask(\n",
@@ -710,7 +719,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"linear\"),\n",
-    "    ipopt_options={\"maxiter\": 100, \"disp\": 0},\n",
+    "    ipopt_options={\"max_iter\": 100, \"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "result = strategy.ask(\n",
@@ -796,7 +805,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"linear\"),\n",
-    "    ipopt_options={\"maxiter\": 100, \"disp\": 0},\n",
+    "    ipopt_options={\"max_iter\": 100, \"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "result = strategy.ask(12, raise_validation_error=False)\n",
@@ -976,7 +985,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"linear\"),\n",
-    "    ipopt_options={\"maxiter\": 500, \"disp\": 0},\n",
+    "    ipopt_options={\"max_iter\": 500, \"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "result = strategy.ask(12)\n",
diff --git a/tutorials/doe/design_with_explicit_formula.ipynb b/tutorials/doe/design_with_explicit_formula.ipynb
index d42ca8b15..3657326e7 100644
--- a/tutorials/doe/design_with_explicit_formula.ipynb
+++ b/tutorials/doe/design_with_explicit_formula.ipynb
@@ -192,187 +192,33 @@
    },
    "outputs": [
     {
-     "data": {
-      "text/html": [
-       "<div>\n",
-       "<style scoped>\n",
-       "    .dataframe tbody tr th:only-of-type {\n",
-       "        vertical-align: middle;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe tbody tr th {\n",
-       "        vertical-align: top;\n",
-       "    }\n",
-       "\n",
-       "    .dataframe thead th {\n",
-       "        text-align: right;\n",
-       "    }\n",
-       "</style>\n",
-       "<table border=\"1\" class=\"dataframe\">\n",
-       "  <thead>\n",
-       "    <tr style=\"text-align: right;\">\n",
-       "      <th></th>\n",
-       "      <th>a</th>\n",
-       "      <th>b</th>\n",
-       "      <th>c</th>\n",
-       "      <th>d</th>\n",
-       "    </tr>\n",
-       "  </thead>\n",
-       "  <tbody>\n",
-       "    <tr>\n",
-       "      <th>0</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>1</th>\n",
-       "      <td>2.047832e+00</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>2</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>3</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>4</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>5</th>\n",
-       "      <td>-9.973772e-09</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>6</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>7</th>\n",
-       "      <td>-9.973772e-09</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>8</th>\n",
-       "      <td>-8.091201e-09</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>9</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>10</th>\n",
-       "      <td>-9.981833e-09</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>11</th>\n",
-       "      <td>2.907832e+00</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>12</th>\n",
-       "      <td>-9.976118e-09</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>13</th>\n",
-       "      <td>-9.906825e-09</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>199.999998</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>14</th>\n",
-       "      <td>-9.981833e-09</td>\n",
-       "      <td>40.000000</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>15</th>\n",
-       "      <td>-8.091050e-09</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>79.999999</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "    <tr>\n",
-       "      <th>16</th>\n",
-       "      <td>5.000000e+00</td>\n",
-       "      <td>800.000008</td>\n",
-       "      <td>180.000002</td>\n",
-       "      <td>800.000008</td>\n",
-       "    </tr>\n",
-       "  </tbody>\n",
-       "</table>\n",
-       "</div>"
-      ],
-      "text/plain": [
-       "               a           b           c           d\n",
-       "0   5.000000e+00   40.000000  180.000002  199.999998\n",
-       "1   2.047832e+00   40.000000  180.000002  800.000008\n",
-       "2   5.000000e+00   40.000000  180.000002  800.000008\n",
-       "3   5.000000e+00  800.000008   79.999999  199.999998\n",
-       "4   5.000000e+00  800.000008   79.999999  800.000008\n",
-       "5  -9.973772e-09   40.000000  180.000002  199.999998\n",
-       "6   5.000000e+00  800.000008  180.000002  199.999998\n",
-       "7  -9.973772e-09  800.000008  180.000002  800.000008\n",
-       "8  -8.091201e-09   40.000000   79.999999  199.999998\n",
-       "9   5.000000e+00   40.000000   79.999999  199.999998\n",
-       "10 -9.981833e-09   40.000000   79.999999  800.000008\n",
-       "11  2.907832e+00   40.000000   79.999999  800.000008\n",
-       "12 -9.976118e-09  800.000008  180.000002  199.999998\n",
-       "13 -9.906825e-09  800.000008   79.999999  199.999998\n",
-       "14 -9.981833e-09   40.000000   79.999999  800.000008\n",
-       "15 -8.091050e-09  800.000008   79.999999  800.000008\n",
-       "16  5.000000e+00  800.000008  180.000002  800.000008"
-      ]
-     },
-     "execution_count": null,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Tried to set Option: disp. It is not a valid option. Please check the list of available options.\n"
+     ]
+    },
+    {
+     "ename": "TypeError",
+     "evalue": "Error while assigning an option",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[4], line 7\u001b[0m\n\u001b[1;32m      1\u001b[0m data_model \u001b[38;5;241m=\u001b[39m DoEStrategy(\n\u001b[1;32m      2\u001b[0m     domain\u001b[38;5;241m=\u001b[39mdomain,\n\u001b[1;32m      3\u001b[0m     criterion\u001b[38;5;241m=\u001b[39mDOptimalityCriterion(formula\u001b[38;5;241m=\u001b[39mmodel_type),\n\u001b[1;32m      4\u001b[0m     ipopt_options\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmax_iter\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m100\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdisp\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m0\u001b[39m},\n\u001b[1;32m      5\u001b[0m )\n\u001b[1;32m      6\u001b[0m strategy \u001b[38;5;241m=\u001b[39m strategies\u001b[38;5;241m.\u001b[39mmap(data_model\u001b[38;5;241m=\u001b[39mdata_model)\n\u001b[0;32m----> 7\u001b[0m design \u001b[38;5;241m=\u001b[39m \u001b[43mstrategy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mask\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m17\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m      8\u001b[0m design\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/strategy.py:128\u001b[0m, in \u001b[0;36mStrategy.ask\u001b[0;34m(self, candidate_count, add_pending, raise_validation_error)\u001b[0m\n\u001b[1;32m    123\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhas_sufficient_experiments():\n\u001b[1;32m    124\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    125\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNot enough experiments available to execute the strategy.\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    126\u001b[0m     )\n\u001b[0;32m--> 128\u001b[0m candidates \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ask\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcandidate_count\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcandidate_count\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    130\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdomain\u001b[38;5;241m.\u001b[39mvalidate_candidates(\n\u001b[1;32m    131\u001b[0m     candidates\u001b[38;5;241m=\u001b[39mcandidates,\n\u001b[1;32m    132\u001b[0m     only_inputs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m,\n\u001b[1;32m    133\u001b[0m     raise_validation_error\u001b[38;5;241m=\u001b[39mraise_validation_error,\n\u001b[1;32m    134\u001b[0m )\n\u001b[1;32m    136\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m candidate_count \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/doe_strategy.py:113\u001b[0m, in \u001b[0;36mDoEStrategy._ask\u001b[0;34m(self, candidate_count)\u001b[0m\n\u001b[1;32m    103\u001b[0m num_discrete_vars \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlen\u001b[39m(new_discretes)\n\u001b[1;32m    104\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m (\n\u001b[1;32m    105\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata_model\u001b[38;5;241m.\u001b[39moptimization_strategy \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrelaxed\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    106\u001b[0m     \u001b[38;5;129;01mor\u001b[39;00m (num_binary_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m num_discrete_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m)\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    111\u001b[0m     )\n\u001b[1;32m    112\u001b[0m ):\n\u001b[0;32m--> 113\u001b[0m     design \u001b[38;5;241m=\u001b[39m \u001b[43mfind_local_max_ipopt\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    114\u001b[0m \u001b[43m        \u001b[49m\u001b[43mnew_domain\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    115\u001b[0m \u001b[43m        \u001b[49m\u001b[43mn_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m_candidate_count\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    116\u001b[0m \u001b[43m        \u001b[49m\u001b[43mfixed_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    117\u001b[0m \u001b[43m        \u001b[49m\u001b[43mpartially_fixed_experiments\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43madapted_partially_fixed_candidates\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    118\u001b[0m \u001b[43m        \u001b[49m\u001b[43mipopt_options\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mipopt_options\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    119\u001b[0m \u001b[43m        \u001b[49m\u001b[43mcriterion\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcriterion\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    120\u001b[0m \u001b[43m        \u001b[49m\u001b[43muse_hessian\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43muse_hessian\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    121\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    122\u001b[0m \u001b[38;5;66;03m# TODO adapt to when exhaustive search accepts discrete variables\u001b[39;00m\n\u001b[1;32m    123\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m (\n\u001b[1;32m    124\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata_model\u001b[38;5;241m.\u001b[39moptimization_strategy \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexhaustive\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    125\u001b[0m     \u001b[38;5;129;01mand\u001b[39;00m num_discrete_vars \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[1;32m    126\u001b[0m ):\n",
+      "File \u001b[0;32m~/Documents/temporary/bofire/bofire/strategies/doe/design.py:188\u001b[0m, in \u001b[0;36mfind_local_max_ipopt\u001b[0;34m(domain, n_experiments, criterion, ipopt_options, sampling, fixed_experiments, partially_fixed_experiments, use_hessian)\u001b[0m\n\u001b[1;32m    182\u001b[0m     problem \u001b[38;5;241m=\u001b[39m FirstOrderDoEProblem(\n\u001b[1;32m    183\u001b[0m         doe_objective\u001b[38;5;241m=\u001b[39mobjective_function,\n\u001b[1;32m    184\u001b[0m         bounds\u001b[38;5;241m=\u001b[39mbounds,\n\u001b[1;32m    185\u001b[0m         constraints\u001b[38;5;241m=\u001b[39mconstraints,\n\u001b[1;32m    186\u001b[0m     )\n\u001b[1;32m    187\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m key \u001b[38;5;129;01min\u001b[39;00m _ipopt_options\u001b[38;5;241m.\u001b[39mkeys():\n\u001b[0;32m--> 188\u001b[0m     \u001b[43mproblem\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_option\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m_ipopt_options\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    190\u001b[0m x, info \u001b[38;5;241m=\u001b[39m problem\u001b[38;5;241m.\u001b[39msolve(x0)\n\u001b[1;32m    192\u001b[0m design \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame(\n\u001b[1;32m    193\u001b[0m     x\u001b[38;5;241m.\u001b[39mreshape(n_experiments, \u001b[38;5;28mlen\u001b[39m(domain\u001b[38;5;241m.\u001b[39minputs)),\n\u001b[1;32m    194\u001b[0m     columns\u001b[38;5;241m=\u001b[39mdomain\u001b[38;5;241m.\u001b[39minputs\u001b[38;5;241m.\u001b[39mget_keys(),\n\u001b[1;32m    195\u001b[0m     index\u001b[38;5;241m=\u001b[39m[\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexp\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mi\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n_experiments)],\n\u001b[1;32m    196\u001b[0m )\n",
+      "File \u001b[0;32m~/anaconda3/envs/bofire/lib/python3.11/site-packages/cyipopt/cython/ipopt_wrapper.pyx:495\u001b[0m, in \u001b[0;36mipopt_wrapper.Problem.add_option\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;31mTypeError\u001b[0m: Error while assigning an option"
+     ]
     }
    ],
    "source": [
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=model_type),\n",
-    "    ipopt_options={\"maxiter\": 100, \"disp\": 0},\n",
+    "    ipopt_options={\"max_iter\": 100, \"print_level\": 0},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "design = strategy.ask(17)\n",
@@ -459,7 +305,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.11"
+   "version": "3.11.5"
   },
   "papermill": {
    "default_parameters": {},
diff --git a/tutorials/doe/nchoosek_constraint.ipynb b/tutorials/doe/nchoosek_constraint.ipynb
index a29ff72d3..cb296cbc0 100644
--- a/tutorials/doe/nchoosek_constraint.ipynb
+++ b/tutorials/doe/nchoosek_constraint.ipynb
@@ -315,7 +315,7 @@
     "data_model = DoEStrategy(\n",
     "    domain=domain,\n",
     "    criterion=DOptimalityCriterion(formula=\"fully-quadratic\"),\n",
-    "    ipopt_options={\"maxiter\": 500},\n",
+    "    ipopt_options={\"max_iter\": 500},\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "candidates = strategy.ask(candidate_count=12)\n",
@@ -339,7 +339,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.11"
+   "version": "3.11.5"
   },
   "papermill": {
    "default_parameters": {},
diff --git a/tutorials/doe/optimality_criteria.ipynb b/tutorials/doe/optimality_criteria.ipynb
index 2e9a4c7c2..8443d3b59 100644
--- a/tutorials/doe/optimality_criteria.ipynb
+++ b/tutorials/doe/optimality_criteria.ipynb
@@ -121,7 +121,7 @@
     "    data_model = DoEStrategy(\n",
     "        domain=domain,\n",
     "        criterion=criterion,\n",
-    "        ipopt_options={\"maxiter\": 300},\n",
+    "        ipopt_options={\"max_iter\": 300},\n",
     "    )\n",
     "    strategy = strategies.map(data_model=data_model)\n",
     "    design = strategy.ask(candidate_count=n_experiments)\n",
@@ -209,7 +209,7 @@
     "    ],\n",
     ")\n",
     "data_model = DoEStrategy(\n",
-    "    domain=domain, criterion=SpaceFillingCriterion(), ipopt_options={\"maxiter\": 500}\n",
+    "    domain=domain, criterion=SpaceFillingCriterion(), ipopt_options={\"max_iter\": 500}\n",
     ")\n",
     "strategy = strategies.map(data_model=data_model)\n",
     "X = strategy.ask(candidate_count=40).to_numpy()\n",
@@ -246,7 +246,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.19"
+   "version": "3.11.5"
   },
   "papermill": {
    "default_parameters": {},