diff --git a/bench/solve.py b/bench/solve.py
index 78a11f7..0d9e42c 100644
--- a/bench/solve.py
+++ b/bench/solve.py
@@ -60,6 +60,24 @@ def write_solution(where: Path, data, result):
         fh.write(f"Cost: {round(result.cost(), 2)}\n")
 
 
+def pi(stats, bks_value: int) -> float:
+    """
+    Computes the primal integral over the given statistics, using the provided
+    best-known solution value.
+    """
+    if len(stats.data) == 0:
+        return 100
+
+    bks_value = min(bks_value, stats.data[-1].best_cost)
+    best_values = np.array([datum.best_cost for datum in stats], dtype=float)
+    gaps = (best_values - bks_value) / best_values
+
+    is_feas = np.array([datum.best_feas for datum in stats], dtype=bool)
+    gaps[~is_feas] = 1
+
+    return 100 * np.sum(gaps * stats.runtimes) / sum(stats.runtimes)
+
+
 class SolveResult(NamedTuple):
     """
     Named tuple to store the results of a single solver run.
@@ -79,6 +97,16 @@ class SolveResult(NamedTuple):
     gap
         The gap to the best-known solution if there is one, otherwise
         ``float('nan')``.
+    primal_integral
+        The primal integral of the solver run if a best-known solution is
+        provided and statistics are collected. Otherwise, ``float('nan')``.
+        See [1]_ for details.
+
+    References
+    ----------
+    .. [1] Berthold, T. (2013). Measuring the impact of primal heuristics.
+            *Operations Research Letters*, 41(6): 611-614.
+            https://doi.org/10.1016/j.orl.2013.08.007.
     """
 
     instance: str
@@ -87,6 +115,7 @@ class SolveResult(NamedTuple):
     num_iterations: int
     runtime: float
     gap: float
+    primal_integral: float
 
 
 def _solve(
@@ -189,11 +218,14 @@ def _solve(
         write_solution(sol_dir / (instance_name + ".sol"), data, result)
 
     gap = float("nan")
+    primal_integral = float("nan")
+
     if bks_loc:
         sol = read_solution(bks_loc, data)
         cost_eval = CostEvaluator([0] * data.num_load_dimensions, 0, 0)
         bks = cost_eval.cost(sol)
         gap = 100 * (result.cost() - bks) / bks
+        primal_integral = pi(result.stats, bks)
 
     return SolveResult(
         instance_name,
@@ -202,6 +234,7 @@ def _solve(
         result.num_iterations,
         round(result.runtime, 3),
         round(gap, 2),
+        round(primal_integral, 2),
     )
 
 
@@ -246,15 +279,24 @@ def benchmark(
         ("iters", int),
         ("time", float),
         ("gap", float),
+        ("pi", float),
     ]
 
     data = np.asarray(res, dtype=dtypes)
-    headers = ["Instance", "OK", "Obj.", "Iters. (#)", "Time (s)", "Gap (%)"]
+    headers = [
+        "Instance",
+        "OK",
+        "Obj.",
+        "Iters. (#)",
+        "Time (s)",
+        "Gap (%)",
+        "PI (%)",
+    ]
 
-    exclude_gap = solutions is None
-    if exclude_gap:
+    exclude_headers = solutions is None
+    if exclude_headers:
         data = data[["inst", "ok", "obj", "iters", "time"]]
-        headers = headers[:-1]
+        headers = headers[:-2]
 
     print("\n", tabulate(headers, data), "\n", sep="")
     print(f"     Avg. objective: {data['obj'].mean():.0f}")
@@ -262,8 +304,9 @@ def benchmark(
     print(f"      Avg. run-time: {data['time'].mean():.2f}s")
     print(f"       Total not OK: {np.count_nonzero(data['ok'] == 'N')}")
 
-    if not exclude_gap:
+    if not exclude_headers:
         print(f"           Avg. gap: {data['gap'].mean():.2f}%")
+        print(f"            Avg. PI: {data['pi'].mean():.2f}%")
 
 
 def setup_parser(subparser):
@@ -280,8 +323,8 @@ def setup_parser(subparser):
 
     msg = """
     Optional paths to best-known solutions in VRPLIB format, used to calculate
-    gaps. If provided, it must match the number of instances. Instances and
-    solutions are paired in the given order.
+    gaps and primal integrals. If provided, it must match the number of
+    instances. Instances and solutions are paired in the given order.
     """
     parser.add_argument("--solutions", nargs="+", type=Path, help=msg)
 
diff --git a/pyproject.toml b/pyproject.toml
index 778d2cf..05716af 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ readme = "README.md"
 
 
 [tool.poetry.dependencies]
-python = ">=3.10,<4.0"
+python = ">=3.11"
 numpy = [
     # Numpy 1.26 is the first version of numpy that supports Python 3.12+.
     { version = ">=1.15.2", python = "<3.12" },