Repository: PyVRP/PyVRP
Pull request #713 (2024-12-06 10:33:28)

Penalty per load dimension
This PR generalises the load penalty system from a single penalty for all loads, to one per load dimension. Part of #458.

<details>

**Notes**:

Please read our [contributing guidelines](https://pyvrp.org/dev/contributing.html) first.
In particular:

- You must add tests when making code changes.
  This keeps the code coverage level up, and helps ensure the changes work as intended.
- When fixing a bug, you must add a test that would produce the bug in the master branch, and then show that it is fixed with the new code. 
- New code additions must be well formatted. Changes should pass the pre-commit workflow, which you can set up locally using [pre-commit](https://pre-commit.com/#intro). 
- Docstring additions must render correctly, including escapes and LaTeX.
- Finally, it is essential that all contributions in this PR are license-compatible with PyVRP's MIT license.
  Please check that this PR can be included into PyVRP under the MIT license.

</details>


codspeed-hq[bot]: ## [CodSpeed Performance Report](https://codspeed.io/PyVRP/PyVRP/branches/penalty-per-load-dim)

### Merging #713 will **not alter performance**

<sub>Comparing <code>penalty-per-load-dim</code> (5ad9a56) with <code>main</code> (0534f00)</sub>



### Summary

`✅ 42` untouched benchmarks  




codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PyVRP/PyVRP/pull/713?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 96.39%. Comparing base [(`0534f00`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/0534f008cd263c995131f98a484f01cc8e19a52e?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) to head [(`5ad9a56`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/5ad9a562002bace61e90bf1002dc96d520bb098d?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #713      +/-   ##
==========================================
- Coverage   96.73%   96.39%   -0.34%     
==========================================
  Files          76       76              
  Lines        3738     3722      -16     
  Branches      991      963      -28     
==========================================
- Hits         3616     3588      -28     
- Misses         14       15       +1     
- Partials      108      119      +11     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PyVRP/PyVRP/pull/713?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).

N-Wouda: @wouterkool @leonlan could you review this?

Patch:
diff --git a/examples/using_pyvrp_components.ipynb b/examples/using_pyvrp_components.ipynb
index d3c1b78..04850df 100644
--- a/examples/using_pyvrp_components.ipynb
+++ b/examples/using_pyvrp_components.ipynb
@@ -137,7 +137,11 @@
    "source": [
     "from pyvrp import Solution, CostEvaluator\n",
     "\n",
-    "cost_evaluator = CostEvaluator(load_penalty=20, tw_penalty=20, dist_penalty=0)\n",
+    "cost_evaluator = CostEvaluator(\n",
+    "    load_penalties=[20],\n",
+    "    tw_penalty=20,\n",
+    "    dist_penalty=0,\n",
+    ")\n",
     "sol = Solution.make_random(INSTANCE, rng)"
    ]
   },
@@ -197,7 +201,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cost_evaluator = CostEvaluator(200, 200, 0)\n",
+    "cost_evaluator = CostEvaluator([200], 200, 0)\n",
     "new_sol = ls.search(sol, cost_evaluator)\n",
     "\n",
     "assert new_sol.is_feasible()"
diff --git a/pyvrp/PenaltyManager.py b/pyvrp/PenaltyManager.py
index a6123b0..18f6c9f 100644
--- a/pyvrp/PenaltyManager.py
+++ b/pyvrp/PenaltyManager.py
@@ -107,13 +107,12 @@ class PenaltyManager:
 
     Parameters
     ----------
+    initial_penalties
+        Initial penalty values for units of load (idx 0), duration (1), and
+        distance (2) violations. These values are clipped to the range
+        ``[MIN_PENALTY, MAX_PENALTY]``.
     params
         PenaltyManager parameters. If not provided, a default will be used.
-    initial_penalties
-        Initial penalty values for unit load (idx 0), duration (1), and
-        distance (2) violations. Defaults to ``(20, 6, 6)`` for backwards
-        compatibility. These values are clipped to the range ``[MIN_PENALTY,
-        MAX_PENALTY]``.
     """
 
     MIN_PENALTY: float = 0.1
@@ -122,22 +121,31 @@ class PenaltyManager:
 
     def __init__(
         self,
+        initial_penalties: tuple[list[float], float, float],
         params: PenaltyParams = PenaltyParams(),
-        initial_penalties: tuple[float, float, float] = (20, 6, 6),
     ):
         self._params = params
         self._penalties = np.clip(
-            initial_penalties,
+            initial_penalties[0] + list(initial_penalties[1:]),
             self.MIN_PENALTY,
             self.MAX_PENALTY,
         )
 
+        # Tracks recent feasibilities for each penalty dimension.
         self._feas_lists: list[list[bool]] = [
-            [],  # tracks recent load feasibility
-            [],  # track recent time feasibility
-            [],  # track recent distance feasibility
+            [] for _ in range(len(self._penalties))
         ]
 
+    def penalties(self) -> tuple[list[float], float, float]:
+        """
+        Returns the current penalty values.
+        """
+        return (
+            self._penalties[:-2].tolist(),  # loads
+            self._penalties[-2],  # duration
+            self._penalties[-1],  # distance
+        )
+
     @classmethod
     def init_from(
         cls,
@@ -170,18 +178,18 @@ class PenaltyManager:
         avg_distance = np.minimum.reduce(distances).mean()
         avg_duration = np.minimum.reduce(durations).mean()
 
-        avg_load = 0
+        avg_load = np.zeros((data.num_load_dimensions,))
         if data.num_clients != 0 and data.num_load_dimensions != 0:
             pickups = np.array([c.pickup for c in data.clients()])
             deliveries = np.array([c.delivery for c in data.clients()])
-            avg_load = np.maximum(pickups, deliveries).mean()
+            avg_load = np.maximum(pickups, deliveries).mean(axis=0)
 
         # Initial penalty parameters are meant to weigh an average increase
         # in the relevant value by the same amount as the average edge cost.
-        init_load = avg_cost / max(avg_load, 1)
+        init_load = avg_cost / np.maximum(avg_load, 1)
         init_tw = avg_cost / max(avg_duration, 1)
         init_dist = avg_cost / max(avg_distance, 1)
-        return cls(params, (init_load, init_tw, init_dist))
+        return cls((init_load.tolist(), init_tw, init_dist), params)
 
     def _compute(self, penalty: float, feas_percentage: float) -> float:
         # Computes and returns the new penalty value, given the current value
@@ -196,9 +204,7 @@ class PenaltyManager:
         else:
             new_penalty = self._params.penalty_decrease * penalty
 
-        clipped = np.clip(new_penalty, self.MIN_PENALTY, self.MAX_PENALTY)
-
-        if clipped == self.MAX_PENALTY:
+        if new_penalty >= self.MAX_PENALTY:
             msg = """
             A penalty parameter has reached its maximum value. This means PyVRP
             struggles to find a feasible solution for the instance that's being
@@ -208,7 +214,7 @@ class PenaltyManager:
             """
             warn(msg, PenaltyBoundWarning)
 
-        return clipped
+        return np.clip(new_penalty, self.MIN_PENALTY, self.MAX_PENALTY)
 
     def _register(self, feas_list: list[bool], penalty: float, is_feas: bool):
         feas_list.append(is_feas)
@@ -224,13 +230,13 @@ class PenaltyManager:
         """
         Registers the feasibility dimensions of the given solution.
         """
-        args = [
-            not sol.has_excess_load(),
+        is_feasible = [
+            *[excess == 0 for excess in sol.excess_load()],
             not sol.has_time_warp(),
             not sol.has_excess_distance(),
         ]
 
-        for idx, is_feas in enumerate(args):
+        for idx, is_feas in enumerate(is_feasible):
             feas_list = self._feas_lists[idx]
             penalty = self._penalties[idx]
             self._penalties[idx] = self._register(feas_list, penalty, is_feas)
@@ -239,10 +245,12 @@ class PenaltyManager:
         """
         Get a cost evaluator using the current penalty values.
         """
-        return CostEvaluator(*self._penalties)
+        *loads, tw, dist = self._penalties
+        return CostEvaluator(loads, tw, dist)
 
     def booster_cost_evaluator(self) -> CostEvaluator:
         """
         Get a cost evaluator using the boosted current penalty values.
         """
-        return CostEvaluator(*(self._penalties * self._params.repair_booster))
+        *loads, tw, dist = self._penalties * self._params.repair_booster
+        return CostEvaluator(loads, tw, dist)
diff --git a/pyvrp/Result.py b/pyvrp/Result.py
index 2284315..c60ef0c 100644
--- a/pyvrp/Result.py
+++ b/pyvrp/Result.py
@@ -48,7 +48,8 @@ class Result:
         if not self.best.is_feasible():
             return math.inf
 
-        return CostEvaluator(0, 0, 0).cost(self.best)
+        num_load_dims = len(self.best.excess_load())
+        return CostEvaluator([0] * num_load_dims, 0, 0).cost(self.best)
 
     def is_feasible(self) -> bool:
         """
diff --git a/pyvrp/_pyvrp.pyi b/pyvrp/_pyvrp.pyi
index b14b599..530aad8 100644
--- a/pyvrp/_pyvrp.pyi
+++ b/pyvrp/_pyvrp.pyi
@@ -5,11 +5,13 @@ import numpy as np
 class CostEvaluator:
     def __init__(
         self,
-        load_penalty: float,
+        load_penalties: list[float],
         tw_penalty: float,
         dist_penalty: float,
     ) -> None: ...
-    def load_penalty(self, load: int, capacity: int) -> int: ...
+    def load_penalty(
+        self, load: int, capacity: int, dimension: int
+    ) -> int: ...
     def tw_penalty(self, time_warp: int) -> int: ...
     def dist_penalty(self, distance: int, max_distance: int) -> int: ...
     def penalised_cost(self, solution: Solution) -> int: ...
diff --git a/pyvrp/cpp/CostEvaluator.cpp b/pyvrp/cpp/CostEvaluator.cpp
index 8e8a1d8..c13cd78 100644
--- a/pyvrp/cpp/CostEvaluator.cpp
+++ b/pyvrp/cpp/CostEvaluator.cpp
@@ -2,10 +2,10 @@
 
 using pyvrp::CostEvaluator;
 
-CostEvaluator::CostEvaluator(double loadPenalty,
+CostEvaluator::CostEvaluator(std::vector<double> loadPenalties,
                              double twPenalty,
                              double distPenalty)
-    : loadPenalty_(loadPenalty),
+    : loadPenalties_(std::move(loadPenalties)),
       twPenalty_(twPenalty),
       distPenalty_(distPenalty)
 {
diff --git a/pyvrp/cpp/CostEvaluator.h b/pyvrp/cpp/CostEvaluator.h
index 3631187..e0db887 100644
--- a/pyvrp/cpp/CostEvaluator.h
+++ b/pyvrp/cpp/CostEvaluator.h
@@ -4,8 +4,10 @@
 #include "Measure.h"
 #include "Solution.h"
 
+#include <cassert>
 #include <concepts>
 #include <limits>
+#include <vector>
 
 namespace pyvrp
 {
@@ -42,7 +44,11 @@ concept DeltaCostEvaluatable = requires(T arg, size_t dimension) {
 };
 
 /**
- * CostEvaluator(load_penalty: float, tw_penalty: float, dist_penalty: float)
+ * CostEvaluator(
+ *     load_penalties: list[float],
+ *     tw_penalty: float,
+ *     dist_penalty: float,
+ * )
  *
  * Creates a CostEvaluator instance.
  *
@@ -51,8 +57,9 @@ concept DeltaCostEvaluatable = requires(T arg, size_t dimension) {
  *
  * Parameters
  * ----------
- * load_penalty
- *    The penalty for each unit of excess load over the vehicle capacity.
+ * load_penalties
+ *    The penalty terms (one for each load dimension) for each unit of load in
+ *    excess of the vehicle capacity.
  * tw_penalty
  *    The penalty for each unit of time warp.
  * dist_penalty
@@ -61,18 +68,28 @@ concept DeltaCostEvaluatable = requires(T arg, size_t dimension) {
  */
 class CostEvaluator
 {
-    double loadPenalty_;
+    std::vector<double> loadPenalties_;  // per load dimension
     double twPenalty_;
     double distPenalty_;
 
+    /**
+     * Computes the cost penalty incurred from the given excess loads. This is
+     * a convenient shorthand for calling ``loadPenalty`` for each dimension.
+     */
+    [[nodiscard]] inline Cost
+    excessLoadPenalties(std::vector<Load> const &excessLoads) const;
+
 public:
-    CostEvaluator(double loadPenalty, double twPenalty, double distPenalty);
+    CostEvaluator(std::vector<double> loadPenalties,
+                  double twPenalty,
+                  double distPenalty);
 
     /**
      * Computes the total excess load penalty for the given load and vehicle
-     * capacity.
+     * capacity, and dimension.
      */
-    [[nodiscard]] inline Cost loadPenalty(Load load, Load capacity) const;
+    [[nodiscard]] inline Cost
+    loadPenalty(Load load, Load capacity, size_t dimension) const;
 
     /**
      * Computes the time warp penalty for the given time warp.
@@ -167,10 +184,25 @@ public:
                    T<vArgs...> const &vProposal) const;
 };
 
-Cost CostEvaluator::loadPenalty(Load load, Load capacity) const
+Cost CostEvaluator::excessLoadPenalties(
+    std::vector<Load> const &excessLoads) const
 {
+    assert(excessLoads.size() == loadPenalties_.size());
+
+    Cost cost = 0;
+    for (size_t dim = 0; dim != loadPenalties_.size(); ++dim)
+        cost += loadPenalties_[dim] * excessLoads[dim].get();
+
+    return cost;
+}
+
+Cost CostEvaluator::loadPenalty(Load load,
+                                Load capacity,
+                                size_t dimension) const
+{
+    assert(dimension < loadPenalties_.size());
     auto const excessLoad = std::max<Load>(load - capacity, 0);
-    return static_cast<Cost>(excessLoad.get() * loadPenalty_);
+    return static_cast<Cost>(excessLoad.get() * loadPenalties_[dimension]);
 }
 
 Cost CostEvaluator::twPenalty([[maybe_unused]] Duration timeWarp) const
@@ -192,13 +224,11 @@ template <CostEvaluatable T>
 Cost CostEvaluator::penalisedCost(T const &arg) const
 {
     // Standard objective plus infeasibility-related penalty terms.
-    auto cost = arg.distanceCost() + arg.durationCost()
-                + (!arg.empty() ? arg.fixedVehicleCost() : 0)
-                + twPenalty(arg.timeWarp())
-                + distPenalty(arg.excessDistance(), 0);
-
-    for (auto const excess : arg.excessLoad())
-        cost += loadPenalty(excess, 0);
+    auto const cost = arg.distanceCost() + arg.durationCost()
+                      + (!arg.empty() ? arg.fixedVehicleCost() : 0)
+                      + excessLoadPenalties(arg.excessLoad())
+                      + twPenalty(arg.timeWarp())
+                      + distPenalty(arg.excessDistance(), 0);
 
     if constexpr (PrizeCostEvaluatable<T>)
         return cost + arg.uncollectedPrizes();
@@ -228,8 +258,7 @@ bool CostEvaluator::deltaCost(Cost &out, T<Args...> const &proposal) const
     out -= distPenalty(route->distance(), route->maxDistance());
 
     if constexpr (!skipLoad)
-        for (auto const excess : route->excessLoad())
-            out -= loadPenalty(excess, 0);
+        out -= excessLoadPenalties(route->excessLoad());
 
     out -= route->durationCost();
     out -= twPenalty(route->timeWarp());
@@ -246,7 +275,8 @@ bool CostEvaluator::deltaCost(Cost &out, T<Args...> const &proposal) const
     {
         auto const &capacity = route->capacity();
         for (size_t dim = 0; dim != capacity.size(); ++dim)
-            out += loadPenalty(proposal.loadSegment(dim).load(), capacity[dim]);
+            out += loadPenalty(
+                proposal.loadSegment(dim).load(), capacity[dim], dim);
     }
 
     auto const duration = proposal.durationSegment();
@@ -279,11 +309,8 @@ bool CostEvaluator::deltaCost(Cost &out,
 
     if constexpr (!skipLoad)
     {
-        for (auto const excess : uRoute->excessLoad())
-            out -= loadPenalty(excess, 0);
-
-        for (auto const excess : vRoute->excessLoad())
-            out -= loadPenalty(excess, 0);
+        out -= excessLoadPenalties(uRoute->excessLoad());
+        out -= excessLoadPenalties(vRoute->excessLoad());
     }
 
     out -= uRoute->durationCost();
@@ -308,13 +335,13 @@ bool CostEvaluator::deltaCost(Cost &out,
     {
         auto const &uCapacity = uRoute->capacity();
         for (size_t dim = 0; dim != uCapacity.size(); ++dim)
-            out += loadPenalty(uProposal.loadSegment(dim).load(),
-                               uCapacity[dim]);
+            out += loadPenalty(
+                uProposal.loadSegment(dim).load(), uCapacity[dim], dim);
 
         auto const &vCapacity = vRoute->capacity();
         for (size_t dim = 0; dim != vCapacity.size(); ++dim)
-            out += loadPenalty(vProposal.loadSegment(dim).load(),
-                               vCapacity[dim]);
+            out += loadPenalty(
+                vProposal.loadSegment(dim).load(), vCapacity[dim], dim);
     }
 
     if constexpr (!exact)
diff --git a/pyvrp/cpp/bindings.cpp b/pyvrp/cpp/bindings.cpp
index ca1cf79..701b55b 100644
--- a/pyvrp/cpp/bindings.cpp
+++ b/pyvrp/cpp/bindings.cpp
@@ -728,14 +728,15 @@ PYBIND11_MODULE(_pyvrp, m)
              });
 
     py::class_<CostEvaluator>(m, "CostEvaluator", DOC(pyvrp, CostEvaluator))
-        .def(py::init<double, double, double>(),
-             py::arg("load_penalty"),
+        .def(py::init<std::vector<double>, double, double>(),
+             py::arg("load_penalties"),
              py::arg("tw_penalty"),
              py::arg("dist_penalty"))
         .def("load_penalty",
              &CostEvaluator::loadPenalty,
              py::arg("load"),
              py::arg("capacity"),
+             py::arg("dimension"),
              DOC(pyvrp, CostEvaluator, loadPenalty))
         .def("tw_penalty",
              &CostEvaluator::twPenalty,
diff --git a/pyvrp/cpp/search/SwapStar.cpp b/pyvrp/cpp/search/SwapStar.cpp
index 8b88f0a..b98ade2 100644
--- a/pyvrp/cpp/search/SwapStar.cpp
+++ b/pyvrp/cpp/search/SwapStar.cpp
@@ -82,21 +82,21 @@ Cost SwapStar::deltaLoadCost(Route::Node *U,
     // added as well. The following addresses this issue with an approximation,
     // which is inexact when there are both pickups and deliveries in the data.
     // So it's pretty rough but fast and seems to mostly work well enough.
-    Cost delta = 0;
+    Cost cost = 0;
     for (size_t dim = 0; dim != data.numLoadDimensions(); ++dim)
     {
-        auto const loadDiff
+        auto const delta
             = std::max(uClient.delivery[dim], uClient.pickup[dim])
               - std::max(vClient.delivery[dim], vClient.pickup[dim]);
 
-        delta += costEvaluator.loadPenalty(uLoad[dim] - loadDiff, uCap[dim]);
-        delta -= costEvaluator.loadPenalty(uLoad[dim], uCap[dim]);
+        cost += costEvaluator.loadPenalty(uLoad[dim] - delta, uCap[dim], dim);
+        cost -= costEvaluator.loadPenalty(uLoad[dim], uCap[dim], dim);
 
-        delta += costEvaluator.loadPenalty(vLoad[dim] + loadDiff, vCap[dim]);
-        delta -= costEvaluator.loadPenalty(vLoad[dim], vCap[dim]);
+        cost += costEvaluator.loadPenalty(vLoad[dim] + delta, vCap[dim], dim);
+        cost -= costEvaluator.loadPenalty(vLoad[dim], vCap[dim], dim);
     }
 
-    return delta;
+    return cost;
 }
 
 SwapStar::InsertPoint SwapStar::bestInsertPoint(


Test patch:
diff --git a/benchmarks/search/test_LocalSearch.py b/benchmarks/search/test_LocalSearch.py
index cd51f86..441a6b1 100644
--- a/benchmarks/search/test_LocalSearch.py
+++ b/benchmarks/search/test_LocalSearch.py
@@ -26,7 +26,7 @@ def test_all_operators(instance, benchmark, request):
         ls.add_route_operator(route_op(data))
 
     sol = Solution.make_random(data, rng)
-    cost_evaluator = CostEvaluator(20, 6, 6)
+    cost_evaluator = CostEvaluator([20], 6, 6)
     benchmark(ls, sol, cost_evaluator)
 
 
@@ -42,7 +42,7 @@ def test_each_node_operator(node_op, instance, benchmark, request):
     ls.add_node_operator(node_op(data))
 
     sol = Solution.make_random(data, rng)
-    cost_evaluator = CostEvaluator(20, 6, 6)
+    cost_evaluator = CostEvaluator([20], 6, 6)
     benchmark(ls, sol, cost_evaluator)
 
 
@@ -58,5 +58,5 @@ def test_each_route_operator(route_op, instance, benchmark, request):
     ls.add_route_operator(route_op(data))
 
     sol = Solution.make_random(data, rng)
-    cost_evaluator = CostEvaluator(20, 6, 6)
+    cost_evaluator = CostEvaluator([20], 6, 6)
     benchmark(ls, sol, cost_evaluator)
diff --git a/tests/crossover/test_ordered_crossover.py b/tests/crossover/test_ordered_crossover.py
index 913a1dc..ec90b4a 100644
--- a/tests/crossover/test_ordered_crossover.py
+++ b/tests/crossover/test_ordered_crossover.py
@@ -10,7 +10,7 @@ def test_raises_when_not_tsp(ok_small):
     Tests that the ordered crossover (OX) operator raises when used on a data
     instance that is not a TSP.
     """
-    cost_eval = CostEvaluator(20, 6, 0)
+    cost_eval = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     sol = Solution(ok_small, [[1, 2], [3, 4]])
 
@@ -77,7 +77,7 @@ def test_empty_solution(prize_collecting):
     """
     data = prize_collecting.replace(vehicle_types=[VehicleType(capacity=[0])])
 
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     empty = Solution(data, [])
diff --git a/tests/crossover/test_selective_route_exchange.py b/tests/crossover/test_selective_route_exchange.py
index f3641d6..3f47584 100644
--- a/tests/crossover/test_selective_route_exchange.py
+++ b/tests/crossover/test_selective_route_exchange.py
@@ -26,7 +26,7 @@ def test_same_parents_same_offspring(ok_small):
     Tests that SREX produces identical offspring when both parents are the
     same.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     solution = Solution(ok_small, [[1, 2], [3, 4]])
@@ -41,7 +41,7 @@ def test_srex_empty_solution(prize_collecting):
     empty. This can occur during prize collecting, and in that case there are
     no routes to exchange.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     empty = Solution(prize_collecting, [])
@@ -71,7 +71,7 @@ def test_raise_invalid_arguments(ok_small, idx1, idx2, num_moved_routes):
     """
     Tests that SREX raises when some of the arguments are invalid.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     sol1 = Solution(ok_small, [[1], [2], [3, 4]])
     sol2 = Solution(ok_small, [[1, 2, 3, 4]])
@@ -91,7 +91,7 @@ def test_srex_move_all_routes(ok_small):
     Tests if SREX produces an offspring that is identical to the second parent
     when all routes are replaced during crossover.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     sol1 = Solution(ok_small, [[1], [2], [3, 4]])
     sol2 = Solution(ok_small, [[1, 2], [3], [4]])
@@ -110,7 +110,7 @@ def test_srex_sorts_routes(ok_small):
     routes, SREX should be invariant to route permutations and always produce
     the exact same offspring.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     routes1 = [[1], [2], [3, 4]]
     routes2 = [[1, 2], [3], [4]]
@@ -134,7 +134,7 @@ def test_srex_does_not_reinsert_unplanned_clients(ok_small):
     case, some clients will no longer in the solution after crossover: SREX
     does not reinsert those, but instead leaves that to the search method.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     # Note that when routes are exchanged, some customers will necessarily be
     # removed from the solution because the routes are overlapping. The
@@ -157,7 +157,7 @@ def test_srex_changed_start_indices(ok_small):
     """
     Tests the case where the initial start indices are changed in SREX.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     sol1 = Solution(ok_small, [[4], [1, 2, 3]])
     sol2 = Solution(ok_small, [[3], [1, 2, 4]])
@@ -187,7 +187,7 @@ def test_srex_heterogeneous_changed_start_indices(ok_small):
             VehicleType(1, capacity=[20]),
         ]
     )
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     # We create the routes sorted by angle such that SREX sorting doesn't
     # affect them
@@ -228,7 +228,7 @@ def test_srex_a_right_move(ok_small):
     Tests the case where the initial start indices are changed by moving the
     A index to the right.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     sol1 = Solution(ok_small, [[4], [2], [1, 3]])
     sol2 = Solution(ok_small, [[3], [2], [4, 1]])
@@ -272,7 +272,7 @@ def test_srex_a_left_move(ok_small):
     Tests the case where the initial start indices are changed by moving to
     A index to the left.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     # We create the routes sorted by angle such that SREX sorting doesn't
     # affect them.
@@ -289,7 +289,7 @@ def test_srex_b_left_move(ok_small):
     Tests the case where the initial start indices are changed by moving the
     B index to the left.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     # We create the routes sorted by angle such that SREX sorting doesn't
     # affect them.
@@ -305,7 +305,7 @@ def test_srex_b_right_move(ok_small):
     Tests the case where the initial start indices are changed by moving the
     B index to the right.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     # We create the routes sorted by angle such that SREX sorting doesn't
     # affect them.
@@ -323,7 +323,7 @@ def test_srex_warns_for_tsp_instances(pr107):
     since in those cases SREX cannot do anything but return one of the two
     parent solutions.
     """
-    cost_eval = CostEvaluator(20, 6, 0)
+    cost_eval = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     sol1 = Solution.make_random(pr107, rng)
@@ -340,7 +340,7 @@ def test_srex_tsp_instance_returns_a_parent_solution(pr107):
     in one the return of one of the two parent solutions. No new solution can
     be created in this case, and SREX is not the right tool to use!
     """
-    cost_eval = CostEvaluator(20, 6, 0)
+    cost_eval = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     sol1 = Solution.make_random(pr107, rng)
diff --git a/tests/plotting/test_plotting.py b/tests/plotting/test_plotting.py
index d615a40..b704dc1 100644
--- a/tests/plotting/test_plotting.py
+++ b/tests/plotting/test_plotting.py
@@ -78,7 +78,7 @@ def test_plot_result():
 
     data = read("data/RC208.vrp", round_func="trunc")
     bks = read_solution("data/RC208.sol")
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     params = PopulationParams()
diff --git a/tests/repair/test_greedy_repair.py b/tests/repair/test_greedy_repair.py
index bf2ad51..167651f 100644
--- a/tests/repair/test_greedy_repair.py
+++ b/tests/repair/test_greedy_repair.py
@@ -11,7 +11,7 @@ def test_raises_given_no_routes_and_unplanned_clients(ok_small):
     to insert unplanned clients into. Greedy repair does not create new routes,
     so this is an impossible situation.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # This call should not raise since unplanned is empty: there are no routes
     # to insert into, which is OK since we have nothing to insert.
@@ -27,7 +27,7 @@ def test_insert_into_empty_route(ok_small):
     Although greedy repair does not create *new* routes, existing empty routes
     will be used if they're available.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # We want to insert client one into an empty route. That should result in
     # a single route with just client 1.
@@ -41,7 +41,7 @@ def test_empty_routes_or_unplanned_is_a_no_op(ok_small):
     If there are no routes, or no unplanned clients, then the returned routes
     should be the same as those given as an argument.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # When unplanned is empty, there is nothing for greedy repair to do, so it
     # should return the exact same routes it received.
@@ -66,7 +66,7 @@ def test_after_depot(ok_small):
     """
     Tests moves where it is optimal to insert directly after the depot.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # We want to insert client 4 into the following route. It is optimal to do
     # so directly after the depot, just before client 3.
@@ -86,7 +86,7 @@ def test_OkSmall(ok_small):
     """
     Tests greedy repair on a small instance.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # We want to insert 1 and 4 into these routes. Both 1 and 4 are close to
     # 3, so it would be cheapest to insert these into the second route.
@@ -116,7 +116,7 @@ def test_RC208(rc208, seed: int):
     routes = [[idx + 1] for idx in range(rc208.num_vehicles)]
     to_repair = Solution(rc208, routes).routes()
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     unplanned = list(range(rc208.num_vehicles + 1, rc208.num_locations))
 
     # Greedily repair by inserting all clients that are not already in the
diff --git a/tests/repair/test_nearest_route_insert.py b/tests/repair/test_nearest_route_insert.py
index ca614ad..acc9fb9 100644
--- a/tests/repair/test_nearest_route_insert.py
+++ b/tests/repair/test_nearest_route_insert.py
@@ -11,7 +11,7 @@ def test_raises_given_no_routes_and_unplanned_clients(ok_small):
     unplanned clients into. The operator does not create new routes, so this is
     an impossible situation.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # This call should not raise since unplanned is empty: there are no routes
     # to insert into, which is OK since we have nothing to insert.
@@ -27,7 +27,7 @@ def test_insert_into_empty_route(ok_small):
     Although nearest route insert does not create *new* routes, existing empty
     routes will be used if they're the only ones available.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # We want to insert client one into an empty route. That should result in
     # a single route with just client 1.
@@ -41,7 +41,7 @@ def test_empty_routes_or_unplanned_is_a_no_op(ok_small):
     If there are no routes, or no unplanned clients, then the returned routes
     should be the same as those given as an argument.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # When unplanned is empty, there is nothing for the operator to do, so it
     # should return the exact same routes it received.
@@ -69,7 +69,7 @@ def test_OkSmall(ok_small):
     """
     Tests nearest route insert on a small instance.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     routes = Solution(ok_small, [[2], [3]]).routes()
     unplanned = [1, 4]
@@ -99,7 +99,7 @@ def test_RC208(rc208, seed: int):
     # same number of vehicles, we initialise dummy routes.
     routes = [Route(rc208, [idx + 1], 0) for idx in range(rc208.num_vehicles)]
     unplanned = list(range(rc208.num_vehicles + 1, rc208.num_clients + 1))
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # Repair inserting all clients that are not already in the dummy routes.
     nearest = nearest_route_insert(routes, unplanned, rc208, cost_eval)
diff --git a/tests/search/test_Exchange.py b/tests/search/test_Exchange.py
index dcf62cf..b20dec6 100644
--- a/tests/search/test_Exchange.py
+++ b/tests/search/test_Exchange.py
@@ -46,7 +46,7 @@ def test_swap_single_route_stays_single_route(rc208, operator):
     only move within the same route, so they can never find a solution that has
     more than one route.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=rc208.num_clients)
@@ -70,7 +70,7 @@ def test_relocate_uses_empty_routes(rc208, operator):
     Unlike the swapping exchange operators, relocate should be able to relocate
     clients to empty routes if that is an improvement.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=rc208.num_clients)
@@ -106,7 +106,7 @@ def test_cannot_exchange_when_parts_overlap_with_depot(ok_small, operator):
     nodes that does not contain the depot (because the routes are very short),
     then no exchange is possible.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=ok_small.num_clients)
@@ -125,7 +125,7 @@ def test_cannot_exchange_when_segments_overlap(ok_small, operator):
     (3, 2)- and (3, 3)-exchange cannot exchange anything on a length-four
     single route solution: there's always overlap between the segments.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=ok_small.num_clients)
@@ -143,7 +143,7 @@ def test_cannot_swap_adjacent_segments(ok_small):
     (2, 2)-exchange on a single route cannot swap adjacent segments, since
     that's already covered by (2, 0)-exchange.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=ok_small.num_clients)
@@ -164,7 +164,7 @@ def test_swap_between_routes_OkSmall(ok_small):
     On the OkSmall example, (2, 1)-exchange should be able to swap parts of a
     two route solution, resulting in something better.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     nb_params = NeighbourhoodParams(nb_granular=ok_small.num_clients)
@@ -201,7 +201,7 @@ def test_relocate_after_depot_should_work(ok_small):
     # This solution can be improved by moving 3 into its own route, that is,
     # inserting it after the depot of an empty route. Before the bug was fixed,
     # (1, 0)-exchange never performed this move.
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     assert_(route1[3] is nodes[-1])
     assert_(op.evaluate(nodes[-1], route2[0], cost_evaluator) < 0)
 
@@ -262,7 +262,7 @@ def test_relocate_only_happens_when_distance_and_duration_allow_it():
     assert_(distance_optimal.distance() < duration_optimal.distance())
     assert_(duration_optimal.time_warp() < distance_optimal.time_warp())
 
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([], 1, 0)
     rng = RandomNumberGenerator(seed=42)
     ls = LocalSearch(data, rng, compute_neighbours(data))
     ls.add_node_operator(Exchange10(data))
@@ -280,7 +280,7 @@ def test_relocate_to_heterogeneous_empty_route(ok_small):
     data = ok_small.replace(vehicle_types=vehicle_types)
 
     # Use a huge cost for load penalties to make other aspects irrelevant
-    cost_evaluator = CostEvaluator(100_000, 6, 0)
+    cost_evaluator = CostEvaluator([100_000], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     # This is a non-empty neighbourhood (so LS does not complain), but the only
@@ -345,7 +345,7 @@ def test_relocate_fixed_vehicle_cost(ok_small, op, base_cost, fixed_cost):
     # First route is not empty, second route is. The operator evaluates moving
     # some nodes to the second route, which would use both of them. That should
     # add to the fixed vehicle cost.
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(
         op.evaluate(route1[1], route2[0], cost_eval), base_cost + fixed_cost
     )
@@ -390,7 +390,7 @@ def test_exchange_with_max_duration_constraint(ok_small, op, max_dur, cost):
     assert_equal(route1.duration(), 5_229)
     assert_equal(route2.duration(), 5_814)
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route1[1], route2[1], cost_eval), cost)
 
 
@@ -429,7 +429,7 @@ def test_within_route_simultaneous_pickup_and_delivery(operator):
     # excess load. For (1, 1)-exchange, we evaluate swapping 1 and 3, which
     # would also resolve the excess load: the important bit is that we visit 3
     # before 1.
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route[1], route[3], cost_eval), -5)
 
 
@@ -461,7 +461,7 @@ def test_relocate_max_distance(ok_small, max_distance: int, expected: int):
     assert_equal(route1.distance(), 5_501)
     assert_equal(route1.excess_distance(), max(5_501 - max_distance, 0))
 
-    cost_eval = CostEvaluator(0, 0, 10)
+    cost_eval = CostEvaluator([0], 0, 10)
     op = Exchange10(data)
 
     # Moving client #2 from route1 to route2 does not improve the overall
@@ -520,7 +520,7 @@ def test_swap_max_distance(ok_small, max_distance: int, expected: int):
     assert_equal(route2.distance(), 3_994)
     assert_equal(route2.excess_distance(), max(3_994 - max_distance, 0))
 
-    cost_eval = CostEvaluator(0, 0, 10)
+    cost_eval = CostEvaluator([0], 0, 10)
     op = Exchange11(data)
 
     # Swapping client #2 in route1 and client #3 in route2 improves the overall
@@ -564,7 +564,7 @@ def test_swap_with_different_profiles(ok_small_two_profiles):
     route2.update()
 
     op = Exchange11(data)
-    cost_eval = CostEvaluator(0, 0, 0)  # all zero so no costs from penalties
+    cost_eval = CostEvaluator([0], 0, 0)  # all zero so no costs from penalties
 
     # This move evaluates swapping loc 3 and 4 between routes. The cost delta
     # is as follows, taking into account the different profiles' distances.
diff --git a/tests/search/test_LocalSearch.py b/tests/search/test_LocalSearch.py
index dadd572..8926ba6 100644
--- a/tests/search/test_LocalSearch.py
+++ b/tests/search/test_LocalSearch.py
@@ -30,7 +30,7 @@ def test_local_search_returns_same_solution_with_empty_neighbourhood(ok_small):
     neighbourhood definition, they cannot do anything with an empty
     neighbourhood.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     neighbours = [[] for _ in range(ok_small.num_locations)]
@@ -161,14 +161,14 @@ def test_reoptimize_changed_objective_timewarp_OkSmall(ok_small):
 
     # With 0 timewarp penalty, the solution should not change since
     # the solution [2, 1, 3, 4] has larger distance.
-    improved_sol = ls.search(sol, CostEvaluator(0, 0, 0))
+    improved_sol = ls.search(sol, CostEvaluator([0], 0, 0))
     assert_equal(sol, improved_sol)
 
     # Now doing it again with a large TW penalty, we must find the alternative
     # solution
     # (previously this was not the case since due to caching the current TW was
     # computed as being zero, causing the move to be evaluated as worse)
-    cost_evaluator_tw = CostEvaluator(0, 1000, 0)
+    cost_evaluator_tw = CostEvaluator([0], 1000, 0)
     improved_sol = ls.search(sol, cost_evaluator_tw)
     improved_cost = cost_evaluator_tw.penalised_cost(improved_sol)
     assert_(improved_cost < cost_evaluator_tw.penalised_cost(sol))
@@ -179,7 +179,7 @@ def test_prize_collecting(prize_collecting):
     Tests that local search works on a small prize-collecting instance.
     """
     rng = RandomNumberGenerator(seed=42)
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([1], 1, 0)
 
     sol = Solution.make_random(prize_collecting, rng)
     sol_cost = cost_evaluator.penalised_cost(sol)
@@ -207,7 +207,7 @@ def test_cpp_shuffle_results_in_different_solution(rc208):
     ls.add_node_operator(Exchange10(rc208))
     ls.add_node_operator(Exchange11(rc208))
 
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([1], 1, 0)
     sol = Solution.make_random(rc208, rng)
 
     # LocalSearch::search is deterministic, so two calls with the same base
@@ -237,7 +237,7 @@ def test_vehicle_types_are_preserved_for_locally_optimal_solutions(rc208):
     ls.add_node_operator(Exchange10(rc208))
     ls.add_node_operator(Exchange11(rc208))
 
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([1], 1, 0)
     sol = Solution.make_random(rc208, rng)
 
     improved = ls.search(sol, cost_evaluator)
@@ -281,7 +281,7 @@ def test_bugfix_vehicle_type_offsets(ok_small):
     ls = cpp_LocalSearch(data, compute_neighbours(data))
     ls.add_node_operator(Exchange10(data))
 
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([1], 1, 0)
 
     current = Solution(data, [Route(data, [1, 3], 1), Route(data, [2, 4], 1)])
     current_cost = cost_evaluator.penalised_cost(current)
@@ -306,7 +306,7 @@ def test_intensify_overlap_tolerance(method, rc208):
 
     func = getattr(ls, method)
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     sol = Solution.make_random(rc208, rng)
 
     # Overlap tolerance is zero, so no routes should have overlap and thus
@@ -333,7 +333,7 @@ def test_intensify_overlap_tolerance_raises_outside_unit_interval(rc208, tol):
     ls = LocalSearch(rc208, rng, neighbours)
     ls.add_route_operator(SwapStar(rc208))
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     sol = Solution.make_random(rc208, rng)
 
     with assert_raises(RuntimeError):  # each tolerance value is outside [0, 1]
@@ -348,7 +348,7 @@ def test_no_op_results_in_same_solution(ok_small):
     """
     rng = RandomNumberGenerator(seed=42)
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     sol = Solution.make_random(ok_small, rng)
 
     # Empty local search does not actually search anything, so it should return
@@ -370,7 +370,7 @@ def test_intensify_can_improve_solution_further(rc208):
     ls.add_node_operator(Exchange11(rc208))
     ls.add_route_operator(SwapStar(rc208))
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # The following solution is locally optimal w.r.t. the node operators. This
     # solution cannot be improved further by repeated calls to ``search()``.
@@ -404,7 +404,7 @@ def test_local_search_completes_incomplete_solutions(ok_small_prizes):
     ls = LocalSearch(ok_small_prizes, rng, compute_neighbours(ok_small_prizes))
     ls.add_node_operator(Exchange10(ok_small_prizes))
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     sol = Solution(ok_small_prizes, [[2], [3, 4]])
     assert_(not sol.is_complete())  # 1 is required but not visited
 
@@ -442,7 +442,7 @@ def test_local_search_does_not_remove_required_clients():
     # Test that the improved solution contains the first client, but removes
     # the second. The first client is required, so could not be removed, but
     # the second could and that is an improving move.
-    cost_eval = CostEvaluator(100, 100, 0)
+    cost_eval = CostEvaluator([100], 100, 0)
     new_sol = ls.search(sol, cost_eval)
     assert_equal(new_sol.num_clients(), 1)
     assert_(new_sol.is_complete())
@@ -466,7 +466,7 @@ def test_mutually_exclusive_group(gtsp):
     ls.add_node_operator(Exchange10(gtsp))
 
     sol = Solution.make_random(gtsp, rng)
-    cost_eval = CostEvaluator(20, 6, 0)
+    cost_eval = CostEvaluator([20], 6, 0)
     improved = ls(sol, cost_eval)
 
     assert_(not sol.is_group_feasible())
@@ -493,7 +493,7 @@ def test_mutually_exclusive_group_not_in_solution(
     sol = Solution(ok_small_mutually_exclusive_groups, [[4]])
     assert_(not sol.is_group_feasible())
 
-    improved = ls(sol, CostEvaluator(20, 6, 0))
+    improved = ls(sol, CostEvaluator([20], 6, 0))
     assert_(improved.is_group_feasible())
 
 
@@ -510,7 +510,7 @@ def test_swap_if_improving_mutually_exclusive_group(
     ls = LocalSearch(ok_small_mutually_exclusive_groups, rng, neighbours)
     ls.add_node_operator(Exchange10(ok_small_mutually_exclusive_groups))
 
-    cost_eval = CostEvaluator(20, 6, 0)
+    cost_eval = CostEvaluator([20], 6, 0)
     sol = Solution(ok_small_mutually_exclusive_groups, [[1, 4]])
     improved = ls(sol, cost_eval)
     assert_(cost_eval.penalised_cost(improved) < cost_eval.penalised_cost(sol))
diff --git a/tests/search/test_SwapRoutes.py b/tests/search/test_SwapRoutes.py
index 3dd5682..59df761 100644
--- a/tests/search/test_SwapRoutes.py
+++ b/tests/search/test_SwapRoutes.py
@@ -63,7 +63,7 @@ def test_evaluate_same_vehicle_type(ok_small):
     route2.update()
 
     op = SwapRoutes(ok_small)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route1, route2, cost_eval), 0)
 
 
@@ -77,7 +77,7 @@ def test_same_route(ok_small):
     route.update()
 
     op = SwapRoutes(ok_small)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route, route, cost_eval), 0)
 
 
@@ -102,7 +102,7 @@ def test_evaluate_empty_routes(ok_small):
     route2.update()
 
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # Vehicle types are no longer the same, but one of the routes is empty.
     # That situation is not currently handled.
@@ -143,7 +143,7 @@ def test_evaluate_capacity_differences(ok_small):
     assert_equal(route2.load(), [3])
 
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(40, 1, 0)
+    cost_eval = CostEvaluator([40], 1, 0)
 
     # Swapping the route plans should alleviate the excess load, since the load
     # of 15 on route1 is below route2's capacity, and similarly for route2's
@@ -194,7 +194,7 @@ def test_evaluate_shift_time_window_differences(ok_small):
     # route. Thus, we should have that swapping the vehicle types results in
     # a lower cost, due to decreased time warp on the routes.
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_(op.evaluate(route1, route2, cost_eval) < 0)
 
 
@@ -240,7 +240,7 @@ def test_evaluate_shift_latest_start_differences(ok_small):
     # Swapping the routes results in a reduction of 1000 units of duration,
     # since the wait before the first client is reduced.
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route1, route2, cost_eval), -1_000)
 
 
@@ -280,7 +280,7 @@ def test_evaluate_max_duration_constraints(ok_small):
     # Swapping the routes results in a reduction of 5'332 - 5'323 = 9 units of
     # time warp.
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(op.evaluate(route1, route2, cost_eval), -9)
 
 
@@ -318,7 +318,7 @@ def test_evaluate_with_different_depots():
     route2.update()
 
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([], 1, 0)
 
     # The routes each cost 16 distance which is not as efficient as swapping
     # them, as that would reduce each route's cost to 4, for an improvement
@@ -346,7 +346,7 @@ def test_different_objectives(ok_small_multi_depot):
 
     data = ok_small_multi_depot.replace(vehicle_types=vehicle_types)
     op = SwapRoutes(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # First route is first depot -> second client -> first depot.
     route1 = Route(data, idx=0, vehicle_type=0)
diff --git a/tests/search/test_SwapStar.py b/tests/search/test_SwapStar.py
index c7712cd..60f6299 100644
--- a/tests/search/test_SwapStar.py
+++ b/tests/search/test_SwapStar.py
@@ -28,7 +28,7 @@ def test_swap_star_identifies_additional_moves_over_regular_swap(rc208):
     Thus, SWAP* should still be able to identify additional improving moves
     after (1, 1)-exchange gets stuck.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     # For a fair comparison we should not hamper the node operator with
@@ -62,7 +62,7 @@ def test_swap_star_on_RC208_instance(rc208, seed: int):
     """
     Evaluate SWAP* on the RC208 instance, over a few seeds.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=seed)
 
     ls = LocalSearch(rc208, rng, compute_neighbours(rc208))
@@ -123,7 +123,7 @@ def test_swap_star_can_swap_in_place():
     route2.append(nodes[3])
     route2.update()
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([], 1, 0)
     swap_star = SwapStar(data)
 
     # Best is to exchange clients 1 and 3. The cost delta is all distance: it
@@ -180,7 +180,7 @@ def test_wrong_load_calculation_bug():
     route2.append(nodes[4])
     route2.update()
 
-    cost_eval = CostEvaluator(1_000, 1, 0)
+    cost_eval = CostEvaluator([1_000], 1, 0)
     swap_star = SwapStar(data)
 
     # Optimal is 0 -> 3 -> 1 -> 0 and 0 -> 4 -> 2 -> 0. This exchanges four
@@ -215,7 +215,7 @@ def test_max_distance(ok_small):
 
         return route1, route2
 
-    cost_eval = CostEvaluator(0, 0, 10)  # only max distance is penalised
+    cost_eval = CostEvaluator([0], 0, 10)  # only max distance is penalised
 
     route1, route2 = make_routes(ok_small)
     assert_equal(route1.distance(), 6_220)
diff --git a/tests/search/test_SwapTails.py b/tests/search/test_SwapTails.py
index 7b38669..9c7d480 100644
--- a/tests/search/test_SwapTails.py
+++ b/tests/search/test_SwapTails.py
@@ -37,7 +37,7 @@ def test_OkSmall_multiple_vehicle_types(
     """
     data = ok_small.replace(vehicle_types=vehicle_types)
 
-    cost_evaluator = CostEvaluator(10_000, 6, 0)  # large load penalty
+    cost_evaluator = CostEvaluator([10_000], 6, 0)  # large load penalty
     rng = RandomNumberGenerator(seed=42)
 
     neighbours: list[list[int]] = [[], [2], [], [], []]  # only 1 -> 2
@@ -89,7 +89,7 @@ def test_move_involving_empty_routes():
     route2.update()  # depot -> depot
 
     op = SwapTails(data)
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([], 0, 0)
 
     # This move does not change the route structure, so the delta cost is 0.
     assert_equal(op.evaluate(route1[2], route2[0], cost_eval), 0)
@@ -164,7 +164,7 @@ def test_move_involving_multiple_depots():
     assert_equal(route2.distance(), 16)
 
     op = SwapTails(data)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([], 1, 0)
 
     assert_equal(op.evaluate(route1[1], route2[1], cost_eval), 0)  # no-op
 
@@ -194,7 +194,7 @@ def test_move_with_different_profiles(ok_small_two_profiles):
     route2.update()
 
     op = SwapTails(data)
-    cost_eval = CostEvaluator(0, 0, 0)  # all zero so no costs from penalties
+    cost_eval = CostEvaluator([0], 0, 0)  # all zero so no costs from penalties
 
     # First route has profile 0, and its distance is thus computed using the
     # first distance matrix.
diff --git a/tests/search/test_primitives.py b/tests/search/test_primitives.py
index f400638..d062754 100644
--- a/tests/search/test_primitives.py
+++ b/tests/search/test_primitives.py
@@ -16,7 +16,7 @@ def test_insert_cost_zero_when_not_allowed(ok_small):
     Tests that insert_cost() returns zero when a move is not possible. This is
     the only case where it shortcuts to return a non-negative delta cost.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     route = Route(ok_small, idx=0, vehicle_type=0)
     for loc in [1, 2]:
@@ -35,7 +35,7 @@ def test_insert_cost(ok_small):
     """
     Tests that the insert_cost() method works correctly on a few basic cases.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     route = Route(ok_small, idx=0, vehicle_type=0)
     for loc in [1, 2]:
@@ -84,7 +84,7 @@ def test_insert_cost_between_different_depots(ok_small_multi_depot):
     route = Route(data, idx=0, vehicle_type=0)
     route.update()
 
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([0], 0, 0)
     dist_mat = data.distance_matrix(0)
 
     delta = dist_mat[0, 2] + dist_mat[2, 1] - dist_mat[0, 1]
@@ -96,7 +96,7 @@ def test_remove_cost_zero_when_not_allowed(ok_small):
     Tests that remove_cost() returns zero when a move is not possible. This is
     the only case where it shortcuts to return a non-negative delta cost.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     route = Route(ok_small, idx=0, vehicle_type=0)
     for loc in [1, 2]:
@@ -115,7 +115,7 @@ def test_remove(ok_small):
     """
     Tests that the remove_cost() method works correctly on a few basic cases.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     route = Route(ok_small, idx=0, vehicle_type=0)
     for loc in [1, 2]:
@@ -135,7 +135,7 @@ def test_insert_fixed_vehicle_cost():
     """
     Tests that insert_cost() adds the fixed vehicle cost if the route is empty.
     """
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([], 0, 0)
     data = ProblemData(
         clients=[Client(x=1, y=1), Client(x=1, y=0)],
         depots=[Depot(x=0, y=0)],
@@ -162,7 +162,7 @@ def test_remove_fixed_vehicle_cost():
     Tests that remove_cost() subtracts the fixed vehicle cost if the route will
     be left empty.
     """
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([], 0, 0)
     data = ProblemData(
         clients=[Client(x=1, y=1), Client(x=1, y=0)],
         depots=[Depot(x=0, y=0)],
@@ -194,7 +194,7 @@ def test_inplace_cost_zero_when_shortcutting_on_guard_clauses(ok_small):
     clauses: either when the first node is in a route, or the second node is
     not.
     """
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     route = Route(ok_small, idx=0, vehicle_type=0)
     node1 = Node(loc=1)
     node2 = Node(loc=2)
@@ -238,5 +238,5 @@ def test_inplace_cost_delta_distance_computation(ok_small):
     #   dist(0, 1) + dist(1, 2) = 3536
     #   dist(0, 3) + dist(3, 2) = 2578
     #                     delta = -958
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([0], 0, 0)
     assert_equal(inplace_cost(node3, node1, ok_small, cost_eval), -958)
diff --git a/tests/stop/test_FirstFeasible.py b/tests/stop/test_FirstFeasible.py
index 1782cc0..edb5769 100644
--- a/tests/stop/test_FirstFeasible.py
+++ b/tests/stop/test_FirstFeasible.py
@@ -15,6 +15,6 @@ def test_stops_on_first_feasible_solution(ok_small):
     assert_(feas.is_feasible())
 
     stop = FirstFeasible()
-    cost_eval = CostEvaluator(0, 0, 0)
+    cost_eval = CostEvaluator([0], 0, 0)
     assert_(not stop(cost_eval.cost(infeas)))
     assert_(stop(cost_eval.cost(feas)))
diff --git a/tests/test_CostEvaluator.py b/tests/test_CostEvaluator.py
index 69daa09..99a08a4 100644
--- a/tests/test_CostEvaluator.py
+++ b/tests/test_CostEvaluator.py
@@ -1,6 +1,6 @@
 import numpy as np
+import pytest
 from numpy.testing import assert_, assert_equal
-from pytest import mark
 
 from pyvrp import (
     Client,
@@ -13,58 +13,63 @@ from pyvrp import (
 )
 
 
-def test_load_penalty():
+@pytest.mark.parametrize("load_penalty", (2, 4))
+def test_load_penalty(load_penalty: float):
     """
     This test asserts that load penalty computations are correct.
     """
-    cost_evaluator = CostEvaluator(2, 1, 0)
+    pen = load_penalty
+    cost_eval = CostEvaluator([pen], 1, 0)
 
-    assert_equal(cost_evaluator.load_penalty(0, 1), 0)  # below capacity
-    assert_equal(cost_evaluator.load_penalty(1, 1), 0)  # at capacity
+    assert_equal(cost_eval.load_penalty(0, 1, 0), 0)  # below capacity
+    assert_equal(cost_eval.load_penalty(1, 1, 0), 0)  # at capacity
+    assert_equal(cost_eval.load_penalty(2, 1, 0), 1 * pen)  # 1 above cap
+    assert_equal(cost_eval.load_penalty(3, 1, 0), 2 * pen)  # 2 above cap
 
-    # Penalty per unit excess capacity is 2
-    # 1 unit above capacity
-    assert_equal(cost_evaluator.load_penalty(2, 1), 2)
-    # 2 units above capacity
-    assert_equal(cost_evaluator.load_penalty(3, 1), 4)
 
-    # Penalty per unit excess capacity is 4
-    cost_evaluator = CostEvaluator(4, 1, 0)
+def test_load_penalty_multiple_dimensions():
+    """
+    Tests that the load penalty computation uses the correct penalty value for
+    each load dimension.
+    """
+    load_penalties = [1, 2]
+    cost_eval = CostEvaluator(load_penalties, 0, 0)
 
-    # 1 unit above capacity
-    assert_equal(cost_evaluator.load_penalty(2, 1), 4)
-    # 2 units above capacity
-    assert_equal(cost_evaluator.load_penalty(3, 1), 8)
+    for dim, pen in enumerate(load_penalties):
+        assert_equal(cost_eval.load_penalty(0, 1, dim), 0)  # below capacity
+        assert_equal(cost_eval.load_penalty(1, 1, dim), 0)  # at capacity
+        assert_equal(cost_eval.load_penalty(2, 1, dim), 1 * pen)  # 1 above cap
+        assert_equal(cost_eval.load_penalty(3, 1, dim), 2 * pen)  # 2 above cap
 
 
-@mark.parametrize("cap", [5, 15, 29, 51, 103])
+@pytest.mark.parametrize("cap", [5, 15, 29, 51, 103])
 def test_load_penalty_always_zero_when_below_capacity(cap: int):
     """
     This test asserts that load penalties are only applied to excess load, that
     is, load in excess of the vehicle's capacity.
     """
     penalty = 2
-    cost_eval = CostEvaluator(penalty, 1, 0)
+    cost_eval = CostEvaluator([penalty], 1, 0)
 
-    assert_equal(cost_eval.load_penalty(0, cap), 0)  # below cap
-    assert_equal(cost_eval.load_penalty(cap - 1, cap), 0)
-    assert_equal(cost_eval.load_penalty(cap, cap), 0)  # at cap
-    assert_equal(cost_eval.load_penalty(cap + 1, cap), penalty)  # above cap
-    assert_equal(cost_eval.load_penalty(cap + 2, cap), 2 * penalty)
+    assert_equal(cost_eval.load_penalty(0, cap, 0), 0)  # below cap
+    assert_equal(cost_eval.load_penalty(cap - 1, cap, 0), 0)
+    assert_equal(cost_eval.load_penalty(cap, cap, 0), 0)  # at cap
+    assert_equal(cost_eval.load_penalty(cap + 1, cap, 0), penalty)  # above cap
+    assert_equal(cost_eval.load_penalty(cap + 2, cap, 0), 2 * penalty)
 
 
 def test_tw_penalty():
     """
     This test asserts that time window penalty computations are correct.
     """
-    cost_evaluator = CostEvaluator(1, 2, 0)
+    cost_evaluator = CostEvaluator([1], 2, 0)
 
     # Penalty per unit time warp is 2.
     assert_equal(cost_evaluator.tw_penalty(0), 0)
     assert_equal(cost_evaluator.tw_penalty(1), 2)
     assert_equal(cost_evaluator.tw_penalty(2), 4)
 
-    cost_evaluator = CostEvaluator(1, 4, 0)
+    cost_evaluator = CostEvaluator([1], 4, 0)
 
     # Penalty per unit excess capacity is now 4.
     assert_equal(cost_evaluator.tw_penalty(0), 0)
@@ -76,7 +81,7 @@ def test_dist_penalty():
     """
     This test asserts that excess distance penalty computations are correct.
     """
-    cost_eval = CostEvaluator(1, 1, 2)
+    cost_eval = CostEvaluator([1], 1, 2)
 
     # Penalty per unit excess distance is 2.
     assert_equal(cost_eval.dist_penalty(-1, 0), 0)
@@ -84,7 +89,7 @@ def test_dist_penalty():
     assert_equal(cost_eval.dist_penalty(1, 0), 2)
     assert_equal(cost_eval.dist_penalty(2, 0), 4)
 
-    cost_eval = CostEvaluator(1, 1, 4)
+    cost_eval = CostEvaluator([1], 1, 4)
 
     # Penalty per unit excess capacity is now 4.
     assert_equal(cost_eval.dist_penalty(-1, 0), 0)
@@ -99,8 +104,8 @@ def test_cost(ok_small):
     solutions, and is a large value (representing infinity) for infeasible
     solutions.
     """
-    default_cost_evaluator = CostEvaluator(0, 0, 0)
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    default_cost_evaluator = CostEvaluator([0], 0, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
 
     feas_sol = Solution(ok_small, [[1, 2], [3], [4]])  # feasible solution
     distance = feas_sol.distance()
@@ -122,7 +127,7 @@ def test_cost_with_prizes(prize_collecting):
     plus a prize term.
     """
     data = prize_collecting
-    cost_evaluator = CostEvaluator(1, 1, 0)
+    cost_evaluator = CostEvaluator([1], 1, 0)
 
     sol = Solution(data, [[1, 2], [3, 4, 5]])
     cost = cost_evaluator.cost(sol)
@@ -146,8 +151,8 @@ def test_penalised_cost(ok_small):
     """
     penalty_capacity = 20
     penalty_tw = 6
-    default_evaluator = CostEvaluator(0, 0, 0)
-    cost_evaluator = CostEvaluator(penalty_capacity, penalty_tw, 0)
+    default_evaluator = CostEvaluator([0], 0, 0)
+    cost_evaluator = CostEvaluator([penalty_capacity], penalty_tw, 0)
 
     feas = Solution(ok_small, [[1, 2], [3], [4]])
     assert_(feas.is_feasible())
@@ -192,7 +197,7 @@ def test_excess_distance_penalised_cost(ok_small):
     assert_equal(routes[0].excess_distance(), 501)
     assert_equal(routes[1].excess_distance(), 0)
 
-    cost_eval = CostEvaluator(0, 0, 10)
+    cost_eval = CostEvaluator([0], 0, 10)
     assert_equal(cost_eval.penalised_cost(sol), 5501 + 4224 + 10 * 501)
 
 
@@ -227,12 +232,11 @@ def test_excess_load_penalised_cost():
     assert_equal(routes[1].excess_load(), [2, 1])
     assert_equal(sol.excess_load(), [3, 1])
 
-    # Both dimensions are penalised.
-    cost_eval = CostEvaluator(10, 0, 0)
+    cost_eval = CostEvaluator([10, 10], 0, 0)
     assert_equal(cost_eval.penalised_cost(sol), 10 * (1 + 2) + 10 * (0 + 1))
 
 
-@mark.parametrize(
+@pytest.mark.parametrize(
     ("assignment", "expected"), [((0, 0), 0), ((0, 1), 10), ((1, 1), 20)]
 )
 def test_cost_with_fixed_vehicle_cost(
@@ -257,7 +261,7 @@ def test_cost_with_fixed_vehicle_cost(
     ]
 
     sol = Solution(data, routes)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     # Solution is feasible, so penalised cost and regular cost are equal. Both
     # should contain the fixed vehicle cost.
@@ -282,7 +286,7 @@ def test_unit_distance_duration_cost(ok_small):
     assert_equal(sol.distance(), 5_501 + 4_224)
     assert_equal(sol.duration(), 6_221 + 5_004)
 
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
     assert_equal(sol.distance_cost(), 31_729)
     assert_equal(sol.duration_cost(), 31_241)
     assert_equal(cost_eval.penalised_cost(sol), 31_729 + 31_241)
diff --git a/tests/test_GeneticAlgorithm.py b/tests/test_GeneticAlgorithm.py
index 97932a7..fa09347 100644
--- a/tests/test_GeneticAlgorithm.py
+++ b/tests/test_GeneticAlgorithm.py
@@ -65,7 +65,7 @@ def test_raises_when_no_initial_solutions(rc208):
     Tests that GeneticAlgorithm raises when no initial solutions are provided,
     since that is insufficient to do crossover.
     """
-    pen_manager = PenaltyManager()
+    pen_manager = PenaltyManager(initial_penalties=([20], 6, 6))
     rng = RandomNumberGenerator(seed=42)
     ls = LocalSearch(rc208, rng, compute_neighbours(rc208))
 
@@ -86,7 +86,7 @@ def test_initial_solutions_added_when_running(rc208):
     Tests that GeneticAlgorithm adds initial solutions to the population
     when running the algorithm.
     """
-    pm = PenaltyManager()
+    pm = PenaltyManager(initial_penalties=([20], 6, 6))
     rng = RandomNumberGenerator(seed=42)
     pop = Population(bpd)
     ls = LocalSearch(rc208, rng, compute_neighbours(rc208))
@@ -107,7 +107,7 @@ def test_initial_solutions_added_when_restarting(rc208):
     Tests that GeneticAlgorithm clears the population and adds the initial
     solutions when restarting.
     """
-    pm = PenaltyManager()
+    pm = PenaltyManager(initial_penalties=([20], 6, 6))
     rng = RandomNumberGenerator(seed=42)
     pop = Population(bpd)
 
@@ -144,7 +144,7 @@ def test_best_solution_improves_with_more_iterations(rc208):
     execution.
     """
     rng = RandomNumberGenerator(seed=42)
-    pm = PenaltyManager()
+    pm = PenaltyManager(initial_penalties=([20], 6, 6))
     pop_params = PopulationParams()
     pop = Population(bpd, params=pop_params)
     init = [
@@ -173,7 +173,7 @@ def test_best_initial_solution(rc208):
     the best found solution.
     """
     rng = RandomNumberGenerator(seed=42)
-    pm = PenaltyManager()
+    pm = PenaltyManager(initial_penalties=([20], 6, 6))
     pop = Population(bpd)
 
     init = [Solution.make_random(rc208, rng) for _ in range(24)]
@@ -197,7 +197,7 @@ def test_infeasible_offspring_is_repaired(rc208):
     """
     bks = Solution(rc208, read_solution("data/RC208.sol"))
 
-    pm = PenaltyManager()
+    pm = PenaltyManager(initial_penalties=([20], 6, 6))
     rng = RandomNumberGenerator(seed=42)
     pop = Population(bpd)
 
@@ -229,7 +229,7 @@ def test_never_repairs_when_zero_repair_probability(rc208):
     probability parameter is set to zero.
     """
     rng = RandomNumberGenerator(seed=42)
-    pm = PenaltyManager(PenaltyParams(repair_booster=10))
+    pm = PenaltyManager(([20], 6, 6), PenaltyParams(repair_booster=10))
     pop = Population(bpd)
 
     ls = LocalSearch(rc208, rng, compute_neighbours(rc208))
@@ -246,7 +246,7 @@ def test_never_repairs_when_zero_repair_probability(rc208):
     # Now we patch the penalty manager: when asked for a booster cost evaluator
     # (as used during repair), this will now raise a runtime error. Since the
     # repair probability is still 100%, this should certainly raise.
-    def raise_when_called(self):
+    def raise_when_called(_):
         raise RuntimeError
 
     pm.booster_cost_evaluator = MethodType(raise_when_called, pm)
diff --git a/tests/test_PenaltyManager.py b/tests/test_PenaltyManager.py
index d139a53..8fc6eb2 100644
--- a/tests/test_PenaltyManager.py
+++ b/tests/test_PenaltyManager.py
@@ -58,12 +58,12 @@ def test_repair_booster():
     penalises constraint violations much more severely.
     """
     params = PenaltyParams(5, 1, 1, 1, 1)
-    pm = PenaltyManager(params, initial_penalties=(1, 1, 1))
+    pm = PenaltyManager(([1], 1, 1), params)
 
     cost_evaluator = pm.cost_evaluator()
 
     assert_equal(cost_evaluator.tw_penalty(1), 1)
-    assert_equal(cost_evaluator.load_penalty(2, 1), 1)  # 1 unit above capacity
+    assert_equal(cost_evaluator.load_penalty(2, 1, 0), 1)  # 1 unit above cap
 
     # With the booster, the penalty values are multiplied by the
     # repairBooster term (x5 in this case).
@@ -71,12 +71,12 @@ def test_repair_booster():
     assert_equal(booster.tw_penalty(1), 5)
     assert_equal(booster.tw_penalty(2), 10)
 
-    assert_equal(booster.load_penalty(2, 1), 5)  # 1 unit above capacity
-    assert_equal(booster.load_penalty(3, 1), 10)  # 2 units above capacity
+    assert_equal(booster.load_penalty(2, 1, 0), 5)  # 1 unit above capacity
+    assert_equal(booster.load_penalty(3, 1, 0), 10)  # 2 units above capacity
 
     # Test that using booster did not affect normal cost_evaluator.
     assert_equal(cost_evaluator.tw_penalty(1), 1)
-    assert_equal(cost_evaluator.load_penalty(2, 1), 1)  # 1 unit above capacity
+    assert_equal(cost_evaluator.load_penalty(2, 1, 0), 1)  # 1 unit above cap
 
 
 def test_load_penalty_update_increase(ok_small):
@@ -86,10 +86,10 @@ def test_load_penalty_update_increase(ok_small):
     """
     num_registrations = 4
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 1, 1))
+    pm = PenaltyManager(([1], 1, 1), params)
 
     # Within bandwidth, so penalty should not change.
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 1)
 
     feas = Solution(ok_small, [[1, 2]])
     infeas = Solution(ok_small, [[1, 2, 3]])
@@ -99,23 +99,23 @@ def test_load_penalty_update_increase(ok_small):
 
     for sol in [feas, infeas, feas, infeas]:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 1)
 
     # Below targetFeasible, so should increase the loadPenalty to 1.1 due to
     # penaltyIncrease, and then int((2 - 1) * 1.1) = 1.
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 1)
 
     # Now we start from a much bigger initial loadPenalty. Here we want the
     # penalty to increase by 10% due to penaltyIncrease = 1.1.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(100, 1, 1))
+    pm = PenaltyManager(([100], 1, 1), params)
 
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 100)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 100)
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 110)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 110)
 
 
 def test_load_penalty_update_decrease(ok_small):
@@ -125,7 +125,7 @@ def test_load_penalty_update_decrease(ok_small):
     """
     num_registrations = 4
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(4, 1, 1))
+    pm = PenaltyManager(([4], 1, 1), params)
 
     feas = Solution(ok_small, [[1, 2]])
     infeas = Solution(ok_small, [[1, 2, 3]])
@@ -134,36 +134,36 @@ def test_load_penalty_update_decrease(ok_small):
     assert_(infeas.has_excess_load())
 
     # Within bandwidth, so penalty should not change.
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 4)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 4)
     for sol in [feas, infeas, feas, infeas]:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 4)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 4)
 
     # Above targetFeasible, so should decrease the loadPenalty to 90%, and -1
     # from the bounds check. So 0.9 * 4 = 3.6, and int(3.6) = 3.
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 3)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 3)
 
     # Now we start from a much bigger initial loadPenalty. Here we want the
     # penalty to decrease by 10% due to penaltyDecrease = 0.9, and -1 due to
     # double -> int.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(100, 1, 1))
+    pm = PenaltyManager(([100], 1, 1), params)
 
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 100)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 100)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 90)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 90)
 
     # Test that the penalty cannot decrease beyond MIN_PENALTY.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(0.1, 1, 1))
+    pm = PenaltyManager(([0.1], 1, 1), params)
 
-    assert_equal(pm.cost_evaluator().load_penalty(11, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(11, 1, 0), 1)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(11, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(11, 1, 0), 1)
 
 
 def test_time_warp_penalty_update_increase(ok_small):
@@ -173,7 +173,7 @@ def test_time_warp_penalty_update_increase(ok_small):
     """
     num_registrations = 4
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 1, 1))
+    pm = PenaltyManager(([1], 1, 1), params)
 
     feas = Solution(ok_small, [[1, 2]])
     infeas = Solution(ok_small, [[1, 2, 3]])
@@ -197,7 +197,7 @@ def test_time_warp_penalty_update_increase(ok_small):
     # penalty to increase by 10% due to penaltyIncrease = 1.1, and +1 due
     # to double -> int.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 100, 1))
+    pm = PenaltyManager(([1], 100, 1), params)
 
     assert_equal(pm.cost_evaluator().tw_penalty(1), 100)
     for sol in [infeas] * num_registrations:
@@ -213,7 +213,7 @@ def test_time_warp_penalty_update_decrease(ok_small):
     """
     num_registrations = 4
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 4, 1))
+    pm = PenaltyManager(([1], 4, 1), params)
 
     feas = Solution(ok_small, [[1, 2]])
     infeas = Solution(ok_small, [[1, 2, 3]])
@@ -237,7 +237,7 @@ def test_time_warp_penalty_update_decrease(ok_small):
     # penalty to decrease by 10% due to penaltyDecrease = 0.9, and -1 due
     # to double -> int.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 100, 1))
+    pm = PenaltyManager(([1], 100, 1), params)
 
     assert_equal(pm.cost_evaluator().tw_penalty(1), 100)
     for sol in [feas] * num_registrations:
@@ -246,7 +246,7 @@ def test_time_warp_penalty_update_decrease(ok_small):
 
     # Test that the penalty cannot decrease beyond MIN_PENALTY.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 0.1, 1))
+    pm = PenaltyManager(([1], 0.1, 1), params)
 
     assert_equal(pm.cost_evaluator().tw_penalty(10), 1)
     for sol in [feas] * num_registrations:
@@ -264,7 +264,7 @@ def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
 
     num_registrations = 4
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(4, 4, 4))
+    pm = PenaltyManager(([4], 4, 4), params)
 
     feas = Solution(data, [[1, 2], [3, 4]])
     infeas = Solution(data, [[1, 2, 3, 4]])
@@ -274,7 +274,7 @@ def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
 
     # Both have four initial penalty, and vehicle capacity is one.
     assert_equal(pm.cost_evaluator().tw_penalty(1), 4)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 4)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 4)
     assert_equal(pm.cost_evaluator().dist_penalty(2, 1), 4)
 
     # Register three times. We need at least four registrations before the
@@ -282,13 +282,13 @@ def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
     for sol in [feas, infeas, feas]:
         pm.register(sol)
         assert_equal(pm.cost_evaluator().tw_penalty(1), 4)
-        assert_equal(pm.cost_evaluator().load_penalty(2, 1), 4)
+        assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 4)
         assert_equal(pm.cost_evaluator().dist_penalty(1, 0), 4)
 
     # Register a fourth time. Now the penalties should change. Since there are
     # more feasible registrations than desired, the penalties should decrease.
     pm.register(feas)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 3)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 3)
     assert_equal(pm.cost_evaluator().tw_penalty(1), 3)
     assert_equal(pm.cost_evaluator().dist_penalty(1, 0), 3)
 
@@ -304,7 +304,7 @@ def test_max_min_penalty(ok_small):
         penalty_decrease=0,
         penalty_increase=2,
     )
-    pm = PenaltyManager(params, (20, PenaltyManager.MAX_PENALTY, 6))
+    pm = PenaltyManager(([20], PenaltyManager.MAX_PENALTY, 6), params)
 
     # Initial penalty is MAX_PENALTY, so one unit of time warp should be
     # penalised by that value.
@@ -337,8 +337,9 @@ def test_warns_max_penalty_value(ok_small):
     This typically indicates a data issue that PyVRP is struggling with.
     """
     params = PenaltyParams(solutions_between_updates=1)
-    initial = (1, PenaltyManager.MAX_PENALTY, 1)
-    pm = PenaltyManager(params, initial_penalties=initial)
+    initial = ([1], PenaltyManager.MAX_PENALTY, 1)
+    pm = PenaltyManager(initial, params)
+    assert_equal(pm.penalties(), initial)
 
     infeas = Solution(ok_small, [[1, 2, 3, 4]])
     assert_(infeas.has_time_warp())
@@ -357,7 +358,7 @@ def test_init_from_load_penalty_value(ok_small):
 
     avg_cost = ok_small.distance_matrix(0).mean()
     avg_load = np.mean([c.delivery for c in ok_small.clients()])
-    assert_equal(cost_eval.load_penalty(1, 0), int(avg_cost / avg_load))
+    assert_equal(cost_eval.load_penalty(1, 0, 0), int(avg_cost / avg_load))
 
 
 def test_init_from_tw_penalty_value(ok_small):
@@ -402,13 +403,31 @@ def test_init_clips_penalties():
     Tests that the initial penalty values are clipped to the [MIN_PENALTY,
     MAX_PENALTY] range.
     """
-    penalties = (0, PenaltyManager.MAX_PENALTY + 1, 2)
+    penalties = ([0], PenaltyManager.MAX_PENALTY + 1, 2)
     pm = PenaltyManager(initial_penalties=penalties)
 
     cost_eval = pm.cost_evaluator()
     assert_allclose(  # MIN
-        cost_eval.load_penalty(10, 0),
+        cost_eval.load_penalty(10, 0, 0),
         10 * PenaltyManager.MIN_PENALTY,
     )
     assert_equal(cost_eval.tw_penalty(1), PenaltyManager.MAX_PENALTY)  # MAX
     assert_equal(cost_eval.dist_penalty(1, 0), 2)  # already OK, so unchanged
+
+
+def test_init_from_multiple_load_penalties(ok_small_multiple_load):
+    """
+    Tests that init_from correctly sets up multiple, different load penalties,
+    one for each load dimension.
+    """
+    pm = PenaltyManager.init_from(ok_small_multiple_load)
+    load_penalties, *_ = pm.penalties()
+    assert_equal(
+        len(load_penalties),
+        ok_small_multiple_load.num_load_dimensions,
+    )
+
+    # The first load dimension has 18 total demand. The second 5. The ratio of
+    # the load penalties should reflect this difference.
+    load_penalty1, load_penalty2 = load_penalties
+    assert_allclose(load_penalty1 / load_penalty2, 5 / 18)
diff --git a/tests/test_Population.py b/tests/test_Population.py
index ee1632d..ad905bd 100644
--- a/tests/test_Population.py
+++ b/tests/test_Population.py
@@ -110,7 +110,7 @@ def test_add_triggers_purge(ok_small):
     survivor selection, that is, a purge that reduces the relevant population
     back to minimum size.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     params = PopulationParams()
@@ -156,7 +156,7 @@ def test_select_returns_same_parents_if_no_other_option(ok_small):
     Tests that the ``select()`` method tries to return two different parents,
     but will return the same solution twice if there is no other option.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=2_147_483_647)
 
     params = PopulationParams(min_pop_size=0)
@@ -197,7 +197,7 @@ def test_pop_is_empty_with_zero_min_pop_size_and_generation_size(ok_small):
     generation size is set to 0: adding a new solution then immediately
     triggers survivor selection.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=12)
 
     params = PopulationParams(min_pop_size=0, generation_size=0)
@@ -220,7 +220,7 @@ def test_elite_solutions_are_not_purged(rc208, nb_elite: int):
     they should be given special treatment - are not purged during survivor
     selection.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     params = PopulationParams(nb_elite=nb_elite)
     rng = RandomNumberGenerator(seed=42)
 
@@ -259,7 +259,7 @@ def test_tournament_ranks_by_fitness(rc208, k: int):
     Tests that the tournament-based parent selection on average returns
     solutions about as often as their relative fitness value would suggest.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(bpd)
 
@@ -307,7 +307,7 @@ def test_tournament_raises_for_invalid_k(rc208, k: int):
     """
     Tests that k >= 0 is required for tournament-based selection.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     pop = Population(bpd)
@@ -323,7 +323,7 @@ def test_purge_removes_duplicates(rc208):
     Tests that purging/survivor selection first removes duplicate solutions,
     before purging by (biased) fitness.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     params = PopulationParams(min_pop_size=5, generation_size=20)
     rng = RandomNumberGenerator(seed=42)
 
@@ -358,7 +358,7 @@ def test_clear(rc208):
     """
     Tests that clearing the population reduces its size to zero.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
 
     pop = Population(bpd)
diff --git a/tests/test_ProgressPrinter.py b/tests/test_ProgressPrinter.py
index 4257b22..8d2dd26 100644
--- a/tests/test_ProgressPrinter.py
+++ b/tests/test_ProgressPrinter.py
@@ -75,7 +75,7 @@ def test_iteration(ok_small, capsys):
     """
     pop = Population(bpd)
     rng = RandomNumberGenerator(seed=42)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     for _ in range(10):
         pop.add(Solution.make_random(ok_small, rng), cost_eval)
@@ -123,7 +123,7 @@ def test_should_print_false_no_output(ok_small, capsys):
     """
     pop = Population(bpd)
     rng = RandomNumberGenerator(seed=42)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     for _ in range(10):
         pop.add(Solution.make_random(ok_small, rng), cost_eval)
@@ -163,7 +163,7 @@ def test_print_dash_when_subpopulation_is_empty(ok_small, routes, capsys):
     Tests that a "-" is printed as cost if one of the subpopulations is empty.
     """
     pop = Population(bpd)
-    cost_eval = CostEvaluator(1, 1, 0)
+    cost_eval = CostEvaluator([1], 1, 0)
 
     solution = Solution(ok_small, routes)
     pop.add(solution, cost_eval)
diff --git a/tests/test_Result.py b/tests/test_Result.py
index ccb7a07..793a11b 100644
--- a/tests/test_Result.py
+++ b/tests/test_Result.py
@@ -27,7 +27,7 @@ def test_fields_are_correctly_set(ok_small, routes, num_iterations, runtime):
     assert_allclose(res.runtime, runtime)
 
     if sol.is_feasible():
-        assert_equal(res.cost(), CostEvaluator(0, 0, 0).cost(sol))
+        assert_equal(res.cost(), CostEvaluator([0], 0, 0).cost(sol))
     else:
         assert_equal(res.cost(), math.inf)
 
@@ -54,7 +54,7 @@ def test_num_iterations(ok_small, num_iterations: int):
     """
     Tests access to the ``num_iterations`` property.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(broken_pairs_distance)
     stats = Statistics()
@@ -81,7 +81,7 @@ def test_summary_contains_essential_information(ok_small, routes):
 
     # Test that feasibility status and solution cost are presented.
     if sol.is_feasible():
-        cost = CostEvaluator(0, 0, 0).cost(sol)
+        cost = CostEvaluator([0], 0, 0).cost(sol)
         assert_(str(cost) in summary)
     else:
         assert_("INFEASIBLE" in summary)
@@ -120,7 +120,7 @@ def test_result_can_be_pickled(ok_small, num_iterations: int):
     Tests that a ``Result`` object can be pickled: it can be serialised and
     unserialised. This is useful for e.g. storing results to disk.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(broken_pairs_distance)
     stats = Statistics()
diff --git a/tests/test_Statistics.py b/tests/test_Statistics.py
index 0ceb50e..cccbdd2 100644
--- a/tests/test_Statistics.py
+++ b/tests/test_Statistics.py
@@ -16,7 +16,7 @@ def test_csv_serialises_correctly(ok_small, tmp_path):
     Tests that writing a CSV of a ``Statistics`` object and then reading that
     CSV again returns the same object.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(broken_pairs_distance)
 
@@ -45,7 +45,7 @@ def test_collect_a_data_point_per_iteration(ok_small, num_iterations: int):
     Tests that the statistics object collects on feasible and infeasible data
     point every time ``collect_from`` is called.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(broken_pairs_distance)
 
@@ -65,7 +65,7 @@ def test_eq(ok_small, num_iterations: int):
     """
     Tests the equality operator.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=42)
     pop = Population(broken_pairs_distance)
 
@@ -93,7 +93,7 @@ def test_more_eq():
     """
     Tests the equality operator for the same population trajectory.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     pop = Population(broken_pairs_distance)
     assert_equal(len(pop), 0)
 
@@ -126,7 +126,7 @@ def test_not_collecting():
     stats = Statistics(collect_stats=False)
     assert_(not stats.is_collecting())
 
-    cost_eval = CostEvaluator(1, 1, 1)
+    cost_eval = CostEvaluator([1], 1, 1)
     pop = Population(broken_pairs_distance)
     stats.collect_from(pop, cost_eval)
     stats.collect_from(pop, cost_eval)
diff --git a/tests/test_SubPopulation.py b/tests/test_SubPopulation.py
index 90f4e4d..ae54f99 100644
--- a/tests/test_SubPopulation.py
+++ b/tests/test_SubPopulation.py
@@ -18,7 +18,7 @@ def test_avg_distance_closest_is_same_up_to_nb_close(rc208, nb_close: int):
     Tests that the average distance of a solution to other solutions only looks
     at the nearest ``nb_close`` solutions.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=5)
 
     params = PopulationParams(
@@ -55,7 +55,7 @@ def test_avg_distance_closest_for_single_route_solutions(rc208):
     Tests that the closest computations are exactly right for a simple, single
     route solution where it's easy to reason about what's going on.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     params = PopulationParams(min_pop_size=0, nb_close=10)
 
     subpop = SubPopulation(bpd, params)
@@ -93,7 +93,7 @@ def test_fitness_is_purely_based_on_cost_when_only_elites(rc208):
     Tests than when all solutions are considered elite, the fitness values
     are completely determines by the solutions' cost.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=51)
     params = PopulationParams(nb_elite=25, min_pop_size=25)
     subpop = SubPopulation(bpd, params)
@@ -128,7 +128,7 @@ def test_fitness_is_average_of_cost_and_diversity_when_no_elites(rc208):
     When there are no elite solutions, the fitness ranking averages the cost
     and diversity rank.
     """
-    cost_evaluator = CostEvaluator(20, 6, 0)
+    cost_evaluator = CostEvaluator([20], 6, 0)
     rng = RandomNumberGenerator(seed=52)
     params = PopulationParams(nb_elite=0, min_pop_size=25)
     subpop = SubPopulation(bpd, params)