diff --git a/src/poli/objective_repository/toy_continuous_problem/definitions.py b/src/poli/objective_repository/toy_continuous_problem/definitions.py
index 9b7a4c2b..4e530d8f 100644
--- a/src/poli/objective_repository/toy_continuous_problem/definitions.py
+++ b/src/poli/objective_repository/toy_continuous_problem/definitions.py
@@ -430,6 +430,33 @@ def levy(x: np.ndarray):
     return -(term1 + term2 + term3)
 
 
+def himmelblau(x: np.ndarray):
+    """
+    Compute the Himmelblau function.
+
+    Parameters
+    ----------
+    x (np.ndarray): A 2D numpy array of shape [b, d] where 'b' is the batch size and 'd' is the dimensionality.
+
+    Returns
+    -------
+    np.ndarray: The value of the Himmelblau function at each point in 'x'. Shape is [b,].
+
+    References
+    ----------
+    [1] Surjanovic, S. and Bingham, D. Virtual Library of Simulation Experiments:
+    Test Functions and Datasets. [https://www.sfu.ca/~ssurjano/optimization.html]
+    """
+    assert len(x.shape) == 2, "Input x must be a 2D array."
+    d = x.shape[1]
+    assert d == 2, "Dimensionality must be 2."
+
+    x1 = x[:, 0]
+    x2 = x[:, 1]
+
+    return -((x1**2 + x2 - 11) ** 2 + (x1 + x2**2 - 7) ** 2)
+
+
 if __name__ == "__main__":
     b = branin_2d
     maximal_b = b(np.array([[-np.pi, 12.275], [np.pi, 2.275], [9.42478, 2.475]]))
diff --git a/src/poli/objective_repository/toy_continuous_problem/toy_continuous_problem.py b/src/poli/objective_repository/toy_continuous_problem/toy_continuous_problem.py
index 04d132ac..83173ba6 100644
--- a/src/poli/objective_repository/toy_continuous_problem/toy_continuous_problem.py
+++ b/src/poli/objective_repository/toy_continuous_problem/toy_continuous_problem.py
@@ -26,6 +26,7 @@
     easom,
     egg_holder,
     hartmann_6d,
+    himmelblau,
     levy,
     rosenbrock,
     shifted_sphere,
@@ -55,6 +56,7 @@
     "branin_2d",
     "rosenbrock",
     "levy",
+    "himmelblau",
 ]
 TWO_DIMENSIONAL_PROBLEMS = [
     "easom",
@@ -62,6 +64,7 @@
     "egg_holder",
     "camelback_2d",
     "branin_2d",
+    "himmelblau",
 ]
 SIX_DIMENSIONAL_PROBLEMS = ["hartmann_6d"]
 
@@ -248,6 +251,12 @@ def __init__(
             self.optima_location = np.array([1.0] * n_dims)
             self.solution_length = n_dims
             self.x0 = np.array([[0.0] * self.solution_length])
+        elif name == "himmelblau":
+            self.function = himmelblau
+            self.limits = [-5.0, 5.0]
+            self.optima_location = np.array([3.0, 2.0])
+            self.solution_length = n_dims
+            self.x0 = np.array([[0.0] * self.solution_length])
         else:
             raise ValueError(f"Expected {name} to be one of {POSSIBLE_FUNCTIONS}")