Penalty per load dimension (PyVRP#713)

* Store a vector for load penalties in CostEvaluator

* Update constructor arguments, fix tests

* Update PenaltyManager and tests

* Update notebook

* Fix Result calculation

* Complete implementation

* Add another test for penalty computations with multiple load dimensions

Author: N-Wouda

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)

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()"

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 @@ def init_from(
         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 @@ def _compute(self, penalty: float, feas_percentage: float) -> float:
         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 @@ def _compute(self, penalty: float, feas_percentage: float) -> float:
             """
             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 @@ def register(self, sol: Solution):
         """
         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 @@ def cost_evaluator(self) -> CostEvaluator:
         """
         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)

pyvrp/Result.py

@@ -48,7 +48,8 @@ def cost(self) -> float:
         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:
         """

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: ...

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)
 {

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 @@ class CostEvaluator
                    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)