diff --git a/desdeo/problem/gurobipy_evaluator.py b/desdeo/problem/gurobipy_evaluator.py
index 744ccb0b0..968e71bc3 100644
--- a/desdeo/problem/gurobipy_evaluator.py
+++ b/desdeo/problem/gurobipy_evaluator.py
@@ -33,7 +33,8 @@ class GurobipyEvaluatorWarning(UserWarning):
 class GurobipyEvaluator:
     """Defines as evaluator that transforms an instance of Problem into a GurobipyModel."""
 
-    # gp.Model does not support these, so the evaluator will handle them
+    variables: dict[str, gp.Var | gp.MVar]
+    constraints: dict[str, gp.Constr]
     objective_functions: dict[str, gp.Var | gp.MVar | gp.LinExpr | gp.QuadExpr | gp.MLinExpr | gp.MQuadExpr]
     scalarizations: dict[str, gp.Var | gp.MVar | gp.LinExpr | gp.QuadExpr | gp.MLinExpr | gp.MQuadExpr]
     extra_functions: dict[
@@ -54,13 +55,14 @@ def __init__(self, problem: Problem):
         self.scalarizations = {}
         self.extra_functions = {}
         self.constants = {}
-        self.mvars = {}
+        self.variables = {}
+        self.constraints = {}
 
         # set the parser
         self.parse = MathParser(to_format=FormatEnum.gurobipy).parse
 
         # Add variables
-        self.model = self.init_variables(problem)
+        self.variables = self.init_variables(problem)
 
         # Add constants, if any
         if problem.constants is not None:
@@ -75,7 +77,7 @@ def __init__(self, problem: Problem):
 
         # Add constraints, if any
         if problem.constraints is not None:
-            self.model = self.init_constraints(problem)
+            self.constraints = self.init_constraints(problem)
 
         # Add scalarization functions, if any
         if problem.scalarization_funcs is not None:
@@ -83,7 +85,7 @@ def __init__(self, problem: Problem):
 
         self.problem = problem
 
-    def init_variables(self, problem: Problem) -> gp.Model:
+    def init_variables(self, problem: Problem) -> dict[str, gp.Var | gp.MVar]:
         """Add variables to the GurobipyModel.
 
         Args:
@@ -94,8 +96,9 @@ def init_variables(self, problem: Problem) -> gp.Model:
                 I.e., the variables are of a non supported type.
 
         Returns:
-            GurobipyModel: the GurobipyModel with the variables added as attributes.
+            dict[str, gp.Var | gp.MVar]: the variables added to the model.
         """
+        variables = {}
         for var in problem.variables:
             if isinstance(var, Variable):
                 # handle regular variables
@@ -121,6 +124,7 @@ def init_variables(self, problem: Problem) -> gp.Model:
                 # set the initial value, if one has been defined
                 if var.initial_value is not None:
                     gvar.setAttr("Start", var.initial_value)
+                variables[var.symbol] = gvar
 
             elif isinstance(var, TensorVariable):
                 # handle tensor variables, i.e., vectors etc..
@@ -160,12 +164,12 @@ def init_variables(self, problem: Problem) -> gp.Model:
                 # set the initial value, if one has been defined
                 if var.initial_values is not None:
                     gvar.setAttr("Start", np.array(var.get_initial_values()))
-                self.mvars[var.symbol] = gvar
+                variables[var.symbol] = gvar
 
         # update the model before returning, so that other expressions can reference the variables
         self.model.update()
 
-        return self.model
+        return variables
 
     def init_constants(self, problem: Problem) -> dict[str, int | float | list[int] | list[float]]:
         """Add constants to a GurobipyEvaluator.
@@ -276,6 +280,7 @@ def init_constraints(self, problem: Problem) -> gp.Model:
         Returns:
             GurobipyModel: the GurobipyModel with the constraint expressions added.
         """
+        constraints = {}
         for cons in problem.constraints:
             gp_expr = self.parse(cons.func, callback=self.get_expression_by_name)
 
@@ -289,10 +294,10 @@ def init_constraints(self, problem: Problem) -> gp.Model:
                     msg = f"Constraint type of {con_type} not supported. Must be one of {ConstraintTypeEnum}."
                     raise GurobipyEvaluatorError(msg)
 
-            self.model.addConstr(gp_expr, name=cons.symbol)
+            constraints[cons.symbol] = self.model.addConstr(gp_expr, name=cons.symbol)
 
         self.model.update()
-        return self.model
+        return constraints
 
     def init_scalarizations(
         self, problem: Problem
@@ -345,6 +350,7 @@ def add_constraint(self, constraint: Constraint) -> gp.Constr:
 
         return_cons = self.model.addConstr(gp_expr, name=constraint.symbol)
         self.model.update()
+        self.constraints[constraint.symbol] = return_cons
         return return_cons
 
     def add_objective(self, obj: Objective):
@@ -461,14 +467,16 @@ def add_variable(self, var: Variable | TensorVariable) -> gp.Var | gp.MVar:
             # set the initial value, if one has been defined
             if var.initial_values is not None:
                 gvar.setAttr("Start", np.array(var.get_initial_values()))
-            self.mvars[var.symbol] = gvar
 
         self.model.update()
+        self.variables[var.symbol] = gvar
         return gvar
 
     def get_expression_by_name(
         self, name: str
-    ) -> gp.Var | gp.MVar | gp.LinExpr | gp.QuadExpr | gp.MLinExpr | gp.MQuadExpr | gp.GenExpr | int | float:
+    ) -> (
+        gp.Var | gp.MVar | gp.LinExpr | gp.QuadExpr | gp.MLinExpr | gp.MQuadExpr | gp.GenExpr | int | float | np.ndarray
+    ):
         """Returns a gurobipy expression corresponding to the name.
 
         Only looks for variables, objective functions, scalarizations, extra functions, and constants.
@@ -480,20 +488,22 @@ def get_expression_by_name(
         Returns:
             gurobipy expression: A mathematical expression that gp.Model can use either as a constraint or an objective
         """
-        expression = self.model.getVarByName(name)
-        if expression is None:
-            # check if an MVar by checking gurobi.MVars stored in the evaluator directly,
-            # which results in terms multiplied by zero being removed from the equations:
-            if name in self.mvars:
-                expression = self.mvars[name]
-            elif name in self.objective_functions:
-                expression = self.objective_functions[name]
-            elif name in self.scalarizations:
-                expression = self.scalarizations[name]
-            elif name in self.extra_functions:
-                expression = self.extra_functions[name]
-            elif name in self.constants:
+        if name in self.variables:
+            expression = self.variables[name]
+        elif name in self.objective_functions:
+            expression = self.objective_functions[name]
+        elif name in self.scalarizations:
+            expression = self.scalarizations[name]
+        elif name in self.extra_functions:
+            expression = self.extra_functions[name]
+        elif name in self.constants:
+            if isinstance(self.constants[name], list):
+                expression = np.array(self.constants[name])
+            else:
                 expression = self.constants[name]
+        else:
+            msg = f"No expression with name {name} found in the gurobipy model."
+            raise GurobipyEvaluatorError(msg)
         return expression
 
     def get_values(self) -> dict[str, float | int | bool | list[float] | list[int]]:
@@ -507,12 +517,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int]]:
         result_dict = {}
 
         for var in self.problem.variables:
-            # if var is type MVar, get the values of MVar
-            if var.symbol in self.mvars:
-                result_dict[var.symbol] = self.mvars[var.symbol].getAttr(gp.GRB.Attr.X)
-            else:
-                result_dict[var.symbol] = self.model.getVarByName(var.symbol).getAttr(gp.GRB.Attr.X)
-
+            result_dict[var.symbol] = self.variables[var.symbol].getAttr(gp.GRB.Attr.X)
         for obj in self.problem.objectives:
             result_dict[obj.symbol] = self.objective_functions[obj.symbol].getValue()
 
@@ -526,7 +531,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int]]:
 
         if self.problem.constraints is not None:
             for const in self.problem.constraints:
-                result_dict[const.symbol] = -self.model.getConstrByName(const.symbol).getAttr("Slack")
+                result_dict[const.symbol] = -self.constraints[const.symbol].getAttr("Slack")
 
         if self.problem.scalarization_funcs is not None:
             for scal in self.problem.scalarization_funcs:
@@ -544,7 +549,8 @@ def remove_constraint(self, symbol: str):
         Args:
             symbol (str): a str representing the symbol of the constraint to be removed.
         """
-        self.model.remove(self.model.getConstrByName(symbol))
+        self.model.remove(self.constraints[symbol])
+        self.constraints.pop(symbol)
         self.model.update()
 
     def remove_variable(self, symbol: str):
@@ -557,18 +563,15 @@ def remove_variable(self, symbol: str):
         Args:
             symbol (str): a str representing the symbol of the variable to be removed.
         """
-        if symbol in self.mvars:
-            self.model.remove(self.mvars[symbol])
-            self.mvars.pop(symbol)
-        else:
-            self.model.remove(self.model.getVarByName(symbol))
+        self.model.remove(self.variables[symbol])
+        self.variables.pop(symbol)
         self.model.update()
 
-    def set_optimization_target(self, target: str, maximize: bool = False):  # noqa: FBT001, FBT002
+    def set_optimization_target(self, target: str, maximize: bool = False):
         """Sets a minimization objective to match the target objective or scalarization of the gurobipy model.
 
         Args:
-            target (str): an str representing a symbol. Needs to match an objective function or scaralization
+            target (str): an str representing a symbol. Needs to match an objective function or scalarization
             function already found in the model.
             maximize (bool): If true, the target function is maximized instead of minimized
 
diff --git a/desdeo/problem/testproblems/__init__.py b/desdeo/problem/testproblems/__init__.py
index 87f5d8390..0f0fcd3f0 100644
--- a/desdeo/problem/testproblems/__init__.py
+++ b/desdeo/problem/testproblems/__init__.py
@@ -29,6 +29,7 @@
     "river_pollution_problem_discrete",
     "river_pollution_scenario",
     "rocket_injector_design",
+    "simple_constrained_quadratic_tensor_test_problem",
     "simple_data_problem",
     "simple_integer_test_problem",
     "simple_knapsack",
@@ -74,6 +75,7 @@
 )
 from .rocket_injector_design_problem import rocket_injector_design
 from .simple_problem import (
+    simple_constrained_quadratic_tensor_test_problem,
     simple_data_problem,
     simple_integer_test_problem,
     simple_linear_test_problem,
diff --git a/desdeo/problem/testproblems/simple_problem.py b/desdeo/problem/testproblems/simple_problem.py
index cf23e3138..bb5bf5a6d 100644
--- a/desdeo/problem/testproblems/simple_problem.py
+++ b/desdeo/problem/testproblems/simple_problem.py
@@ -1,3 +1,5 @@
+"""Defines simple test problems for testing purposes."""
+
 from desdeo.problem.schema import (
     Constant,
     Constraint,
@@ -7,10 +9,13 @@
     Objective,
     ObjectiveTypeEnum,
     Problem,
+    TensorConstant,
+    TensorVariable,
     Variable,
     VariableTypeEnum,
 )
 
+
 def simple_test_problem() -> Problem:
     """Defines a simple problem suitable for testing purposes."""
     variables = [
@@ -131,8 +136,8 @@ def simple_data_problem() -> Problem:
             func=None,
             objective_type=ObjectiveTypeEnum.data_based,
             maximize=i == 1,
-            ideal=3000 if i == 1 else -60.0 if i == 3 else 0,
-            nadir=0 if i == 1 else 15 - 2.0 if i == 3 else 15,
+            ideal=3000 if i == 1 else -60.0 if i == 3 else 0,  # noqa: PLR2004
+            nadir=0 if i == 1 else 15 - 2.0 if i == 3 else 15,  # noqa: PLR2004
         )
         for i in range(1, n_objectives + 1)
     ]
@@ -188,6 +193,66 @@ def simple_linear_test_problem() -> Problem:
     )
 
 
+def simple_constrained_quadratic_tensor_test_problem() -> Problem:
+    """Defines a simple constrained quadratic problem with tensor variables, suitable for testing purposes."""
+    xvar = TensorVariable(
+        name="X",
+        symbol="X",
+        variable_type=VariableTypeEnum.real,
+        shape=[
+            2,
+        ],
+        initial_values=[1, 1],
+        lowerbounds=[-10, -10],
+        upperbounds=[10, 10],
+    )
+
+    mmult = TensorConstant(
+        name="Mmult",
+        symbol="A",
+        shape=[2, 2],
+        values=[[1, 0.5], [0.5, 1]],
+    )
+
+    bvector = TensorConstant(
+        name="mcon",
+        symbol="b",
+        shape=[
+            2,
+        ],
+        values=[1, 1],
+    )
+
+    cons = Constraint(
+        name="cons",
+        symbol="cons",
+        cons_type=ConstraintTypeEnum.LTE,
+        func="b-A@X",
+        is_linear=True,
+        is_convex=True,
+        is_twice_differentiable=True,
+    )
+
+    obj = Objective(
+        name="f_1",
+        symbol="f_1",
+        func="-0.5*X@X",  # this is equivalent to 0.5*X.T@X, since X is a 1D tensor variable
+        maximize=True,
+        is_linear=False,
+        is_convex=True,
+        is_twice_differentiable=True,
+    )
+
+    return Problem(
+        name="Simple constrained quadratic tensor test problem.",
+        description="A simple problem for testing purposes.",
+        variables=[xvar],
+        constants=[mmult, bvector],
+        constraints=[cons],
+        objectives=[obj],
+    )
+
+
 def simple_scenario_test_problem():
     """Returns a simple, scenario-based multiobjective optimization test problem."""
     constants = [Constant(name="c_1", symbol="c_1", value=3)]
diff --git a/desdeo/tools/gurobipy_solver_interfaces.py b/desdeo/tools/gurobipy_solver_interfaces.py
index 615c9ae00..56a5702e2 100644
--- a/desdeo/tools/gurobipy_solver_interfaces.py
+++ b/desdeo/tools/gurobipy_solver_interfaces.py
@@ -103,6 +103,10 @@ def __init__(self, problem: Problem, options: dict[str, any] | None = None):
         if options is not None:
             for key, value in options.items():
                 self.evaluator.model.setParam(key, value)
+        else:
+            # Set some default parameters that are good for most problems.
+            self.evaluator.model.setParam("OutputFlag", 0)  # Suppress Gurobi output
+            self.evaluator.model.setParam("LogToConsole", 0)  # Suppress Gurobi logging to console
 
     def solve(self, target: str) -> SolverResults:
         """Solve the problem for the given target.
diff --git a/tests/test_gurobipy_solver.py b/tests/test_gurobipy_solver.py
index 378489721..9ce1de416 100644
--- a/tests/test_gurobipy_solver.py
+++ b/tests/test_gurobipy_solver.py
@@ -13,7 +13,12 @@
     Variable,
     VariableTypeEnum,
 )
-from desdeo.problem.testproblems import simple_knapsack_vectors, simple_linear_test_problem
+from desdeo.problem.gurobipy_evaluator import GurobipyEvaluatorError
+from desdeo.problem.testproblems import (
+    simple_constrained_quadratic_tensor_test_problem,
+    simple_knapsack_vectors,
+    simple_linear_test_problem,
+)
 from desdeo.tools import GurobipySolver, PersistentGurobipySolver
 
 
@@ -36,7 +41,7 @@ def test_gurobipy_solver():
 @pytest.mark.slow
 @pytest.mark.gurobipy
 def test_gurobipy_persistent_solver():
-    """Tests the bonmin solver."""
+    """Tests the gurobipy solver."""
     problem = simple_linear_test_problem()
     solver = PersistentGurobipySolver(problem)
 
@@ -73,7 +78,8 @@ def test_gurobipy_persistent_solver():
     assert solver.evaluator.model.getConstrByName("c_test") is None
 
     solver.remove_variable("y")
-    assert solver.evaluator.get_expression_by_name("y") is None
+    with pytest.raises(GurobipyEvaluatorError):
+        solver.evaluator.get_expression_by_name("y")
 
     # Check that the solver can still solve the original problem
     # after removing the added variables and constraints
@@ -97,7 +103,8 @@ def test_gurobipy_persistent_solver():
     assert isinstance(solver.evaluator.get_expression_by_name("y"), gp.MVar)
 
     solver.remove_variable("y")
-    assert solver.evaluator.get_expression_by_name("y") is None
+    with pytest.raises(GurobipyEvaluatorError):
+        solver.evaluator.get_expression_by_name("y")
 
 
 @pytest.mark.slow
@@ -124,3 +131,19 @@ def test_gurobipy_solver_with_tensors():
     assert np.allclose(xs["X"], [0.0, 0.0, 1.0, 0.0])
     assert np.isclose(ys["f_1"], 6.0)
     assert np.isclose(ys["f_2"], 7.0)
+
+
+@pytest.mark.slow
+@pytest.mark.gurobipy
+def test_gurobipy_solver_qp_with_tensors():
+    """Test gurobipy solver with a quadratic problem with TensorVariables."""
+    problem = simple_constrained_quadratic_tensor_test_problem()
+    solver = GurobipySolver(problem)
+
+    results = solver.solve("f_1_min")
+
+    assert results.success
+    xs, ys = results.optimal_variables, results.optimal_objectives
+
+    assert np.allclose(xs["X"], [2 / 3, 2 / 3])
+    assert np.isclose(ys["f_1"], -((2 / 3) ** 2))