Repository: PyVRP/PyVRP (python)

Issue:
Infeasible result on small instance with optional client groups for multiple time windows
#### Describe the bug

Hello, i'm interfacing PyVRP and create some dummy test cases to verify its behavior and the correctness of the interface.

I've defined a problem with 2 "real clients", 3 timewindows each. The timewindows and the duration matrix are defined such that only one client is feasible.

Both clients and groups are marked "required: false".

The solver returns the following:

>  A penalty parameter has reached its maximum value. This means PyVRP
> struggles to find a feasible solution for the instance that's being
> solved, either because the instance has no feasible solution, or it
> is very hard to find one. Check the instance carefully to determine
> if a feasible solution exists.

I would expect a feasible solution with 1 client active in the solution. It works well with 2 timewindows, but struggles with 3.

#### Instance details

```json
{
	"depots": [{
		"x": 0,
		"y": 0,
		"tw_early": 0,
		"tw_late": 9223372036854775807,
		"name": "point_0"
	}, {
		"x": 0,
		"y": 0,
		"tw_early": 0,
		"tw_late": 9223372036854775807,
		"name": "_null_store"
	}], 
        "clients": [{
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 3,
		"tw_late": 4,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_1_tw0",
		"group": 0
	}, {
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 7,
		"tw_late": 8,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_1_tw1",
		"group": 0
	}, {
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 11,
		"tw_late": 12,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_1_tw2",
		"group": 0
	}, {
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 5,
		"tw_late": 6,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_2_tw0",
		"group": 1
	}, {
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 10,
		"tw_late": 11,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_2_tw1",
		"group": 1
	}, {
		"x": 0,
		"y": 0,
		"delivery": [],
		"pickup": [],
		"service_duration": 0,
		"tw_early": 15,
		"tw_late": 16,
		"release_time": 0,
		"prize": 429496729,
		"required": false,
		"name": "service_2_tw2",
		"group": 1
	}], 
        "vehicle_types": [{
		"num_available": 1,
		"capacity": [],
		"start_depot": 0,
		"end_depot": 1,
		"fixed_cost": 0,
		"tw_early": 0,
		"tw_late": 9223372036854775807,
		"max_duration": 9223372036854775807,
		"max_distance": 9223372036854775807,
		"unit_distance_cost": 0,
		"unit_duration_cost": 1,
		"profile": 0,
		"start_late": null,
		"reload_depots": [],
		"max_reloads": 9223372036854775807,
		"name": "vehicle_0"
	}], 
        "distance_matrices": [
		[
			[0, 0, 9, 9, 9, 9, 9, 9],
			[0, 0, 0, 0, 0, 0, 0, 0],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 9, 9, 9, 0, 0, 0],
			[9, 0, 9, 9, 9, 0, 0, 0],
			[9, 0, 9, 9, 9, 0, 0, 0]
		]
	], 
        "duration_matrices": [
		[
			[0, 0, 9, 9, 9, 9, 9, 9],
			[0, 0, 0, 0, 0, 0, 0, 0],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 0, 0, 0, 9, 9, 9],
			[9, 0, 9, 9, 9, 0, 0, 0],
			[9, 0, 9, 9, 9, 0, 0, 0],
			[9, 0, 9, 9, 9, 0, 0, 0]
		]
	], 
        "groups": [{
		"clients": [2, 3, 4],
		"required": false
	}, {
		"clients": [5, 6, 7],
		"required": false
	}]
}
```
INSTALLED VERSIONS
------------------
pyvrp: 0.11.1
numpy: 1.26.3
matplotlib: 3.10.3
vrplib: 1.5.1
tqdm: 4.67.1
tomli: 2.2.1
Python: 3.10.12

#### Expected Output
Solution results
INFO: ================
INFO: # routes: 1
INFO: # trips: 1
INFO: # clients: 1
objective: xxxxxxxxx
INFO: distance: 9
INFO: duration: 9


Pull request description:
Make penalty term bounds configurable via PenaltyParams
This PR:

- [x] Closes #836.
- [x] Adds and passes relevant tests.
- [x] Has well-formatted code and documentation.

<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>


Patch:
diff --git a/pyvrp/PenaltyManager.py b/pyvrp/PenaltyManager.py
index b65fdbe..76f5b42 100644
--- a/pyvrp/PenaltyManager.py
+++ b/pyvrp/PenaltyManager.py
@@ -20,8 +20,7 @@ class PenaltyParams:
     repair_booster
         A repair booster value :math:`r \\ge 1`. This value is used to
         temporarily multiply the current penalty terms, to force feasibility.
-        See also
-        :meth:`~pyvrp.PenaltyManager.PenaltyManager.booster_cost_evaluator`.
+        See also :meth:`~PenaltyManager.booster_cost_evaluator`.
     solutions_between_updates
         Number of feasibility registrations between penalty value updates. The
         penalty manager updates the penalty terms every once in a while based
@@ -47,6 +46,19 @@ class PenaltyParams:
         penalty terms are decreased. This ensures a balanced population, with a
         fraction :math:`p_f` feasible and a fraction :math:`1 - p_f` infeasible
         solutions.
+    feas_tolerance
+        Deviation tolerance (in :math:`[0, 1]`) between actual and target
+        percentage of feasible solutions between updates. If the deviation is
+        smaller than this tolerance, the penalty terms are not updated. See
+        also ``target_feasible`` and ``solutions_between_updates``.
+    min_penalty
+        Minimum penalty term value. Must not be negative.
+    max_penalty
+        Maximum penalty term value. Must not be negative.
+
+        .. warning::
+           Setting a (too) large maximum penalty value may cause integer
+           overflow in PyVRP's native extensions.
 
     Attributes
     ----------
@@ -66,6 +78,12 @@ class PenaltyParams:
     target_feasible
         Target percentage :math:`p_f \\in [0, 1]` of feasible registrations
         in the last ``solutions_between_updates`` registrations.
+    feas_tolerance
+        Deviation tolerance for ``target_feasible``.
+    min_penalty
+        Minimum penalty term value.
+    max_penalty
+        Maximum penalty term value.
     """
 
     repair_booster: int = 12
@@ -73,6 +91,9 @@ class PenaltyParams:
     penalty_increase: float = 1.34
     penalty_decrease: float = 0.32
     target_feasible: float = 0.43
+    feas_tolerance: float = 0.05
+    min_penalty: float = 0.1
+    max_penalty: float = 100_000.0
 
     def __post_init__(self):
         if not self.repair_booster >= 1:
@@ -90,6 +111,15 @@ class PenaltyParams:
         if not (0.0 <= self.target_feasible <= 1.0):
             raise ValueError("Expected target_feasible in [0, 1].")
 
+        if not (0.0 <= self.feas_tolerance <= 1.0):
+            raise ValueError("Expected feas_tolerance in [0, 1].")
+
+        if self.min_penalty < 0:
+            raise ValueError("Expected min_penalty >= 0.")
+
+        if self.max_penalty < self.min_penalty:
+            raise ValueError("Expected max_penalty >= min_penalty.")
+
 
 class PenaltyManager:
     """
@@ -110,15 +140,12 @@ class PenaltyManager:
     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]``.
+        [:attr:`~pyvrp.PenaltyManager.PenaltyParams.min_penalty`,
+        :attr:`~pyvrp.PenaltyManager.PenaltyParams.max_penalty`].
     params
         PenaltyManager parameters. If not provided, a default will be used.
     """
 
-    MIN_PENALTY: float = 0.1
-    MAX_PENALTY: float = 100_000.0
-    FEAS_TOL: float = 0.05
-
     def __init__(
         self,
         initial_penalties: tuple[list[float], float, float],
@@ -127,8 +154,8 @@ class PenaltyManager:
         self._params = params
         self._penalties = np.clip(
             initial_penalties[0] + list(initial_penalties[1:]),
-            self.MIN_PENALTY,
-            self.MAX_PENALTY,
+            params.min_penalty,
+            params.max_penalty,
         )
 
         # Tracks recent feasibilities for each penalty dimension.
@@ -209,7 +236,7 @@ class PenaltyManager:
         # and the percentage of feasible solutions since the last update.
         diff = self._params.target_feasible - feas_percentage
 
-        if abs(diff) < self.FEAS_TOL:
+        if abs(diff) < self._params.feas_tolerance:
             return penalty
 
         if diff > 0:
@@ -217,17 +244,21 @@ class PenaltyManager:
         else:
             new_penalty = self._params.penalty_decrease * penalty
 
-        if new_penalty >= self.MAX_PENALTY:
+        if new_penalty >= self._params.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
-            solved, either because the instance has no feasible solution, or it
-            is very hard to find one. Check the instance carefully to determine
-            if a feasible solution exists.
+            struggles to find a feasible solution for this instance, either
+            because the instance has no feasible solution, or it is hard to
+            find one - possibly due to large data scaling differences. Check
+            the instance carefully to determine if a feasible solution exists.
             """
             warn(msg, PenaltyBoundWarning)
 
-        return np.clip(new_penalty, self.MIN_PENALTY, self.MAX_PENALTY)
+        return np.clip(
+            new_penalty,
+            self._params.min_penalty,
+            self._params.max_penalty,
+        )
 
     def _register(self, feas_list: list[bool], penalty: float, is_feas: bool):
         feas_list.append(is_feas)


Test patch:
diff --git a/tests/test_PenaltyManager.py b/tests/test_PenaltyManager.py
index 798959e..2103177 100644
--- a/tests/test_PenaltyManager.py
+++ b/tests/test_PenaltyManager.py
@@ -19,17 +19,24 @@ from pyvrp.exceptions import PenaltyBoundWarning
         "penalty_increase",
         "penalty_decrease",
         "target_feasible",
+        "feas_tolerance",
+        "min_penalty",
+        "max_penalty",
     ),
     [
-        (1, -1, 1, 0.5, 0.5),  # -1 solutions between updates
-        (1, 0, 1, 0.5, 0.5),  # 0 solutions between updates
-        (1, 1, -1.0, 0.5, 0.5),  # -1 penalty increase
-        (1, 1, 0.5, 0.5, 0.5),  # 0.5 penalty increase
-        (1, 1, 1.5, -1, 0.5),  # -1 penalty decrease
-        (1, 1, 1.5, 2, 0.5),  # 2 penalty decrease
-        (1, 1, 1, 1, -1),  # -1 target feasible
-        (1, 1, 1, 1, 2),  # 2 target feasible
-        (0, 1, 1, 1, 1),  # 0 repair booster
+        (0, 1, 1, 1, 1, 0, 0, 0),  # 0 repair booster
+        (1, -1, 1, 0.5, 0.5, 0, 0, 0),  # -1 solutions between updates
+        (1, 0, 1, 0.5, 0.5, 0, 0, 0),  # 0 solutions between updates
+        (1, 1, -1.0, 0.5, 0.5, 0, 0, 0),  # -1 penalty increase
+        (1, 1, 0.5, 0.5, 0.5, 0, 0, 0),  # 0.5 penalty increase
+        (1, 1, 1.5, -1, 0.5, 0, 0, 0),  # -1 penalty decrease
+        (1, 1, 1.5, 2, 0.5, 0, 0, 0),  # 2 penalty decrease
+        (1, 1, 1, 1, -1, 0, 0, 0),  # -1 target feasible
+        (1, 1, 1, 1, 2, 0, 0, 0),  # 2 target feasible
+        (1, 1, 1.5, 0.5, 0.5, -1, 0, 0),  # -1 feas tolerance
+        (1, 1, 1.5, 0.5, 0.5, 2, 0, 0),  # 2 feas tolerance
+        (1, 1, 1.5, 0.5, 0.5, 0, -1, 0),  # -1 min_penalty
+        (1, 1, 1.5, 0.5, 0.5, 0, 2, 1),  # max_penalty < min_penalty
     ],
 )
 def test_constructor_throws_when_arguments_invalid(
@@ -38,6 +45,9 @@ def test_constructor_throws_when_arguments_invalid(
     penalty_increase: float,
     penalty_decrease: float,
     target_feasible: float,
+    feas_tolerance: float,
+    min_penalty: float,
+    max_penalty: float,
 ):
     """
     Tests that invalid arguments are not accepted.
@@ -49,6 +59,9 @@ def test_constructor_throws_when_arguments_invalid(
             penalty_increase,
             penalty_decrease,
             target_feasible,
+            feas_tolerance,
+            min_penalty,
+            max_penalty,
         )
 
 
@@ -156,7 +169,7 @@ def test_load_penalty_update_decrease(ok_small):
         pm.register(sol)
     assert_equal(pm.cost_evaluator().load_penalty(2, 1, 0), 90)
 
-    # Test that the penalty cannot decrease beyond MIN_PENALTY.
+    # Test that the penalty cannot decrease beyond min_penalty.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(([0.1], 1, 1), params)
 
@@ -244,7 +257,7 @@ def test_time_warp_penalty_update_decrease(ok_small):
         pm.register(sol)
     assert_equal(pm.cost_evaluator().tw_penalty(1), 90)
 
-    # Test that the penalty cannot decrease beyond MIN_PENALTY.
+    # Test that the penalty cannot decrease beyond min_penalty.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(([1], 0.1, 1), params)
 
@@ -296,7 +309,7 @@ def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
 @pytest.mark.filterwarnings("ignore::pyvrp.exceptions.PenaltyBoundWarning")
 def test_max_min_penalty(ok_small):
     """
-    Tests that penalty parameters are clipped to [MIN_PENALTY, MAX_PENALTY]
+    Tests that penalty parameters are clipped to [min_penalty, max_penalty]
     when updating their values.
     """
     params = PenaltyParams(
@@ -304,11 +317,11 @@ def test_max_min_penalty(ok_small):
         penalty_decrease=0,
         penalty_increase=2,
     )
-    pm = PenaltyManager(([20], PenaltyManager.MAX_PENALTY, 6), params)
+    pm = PenaltyManager(([20], params.max_penalty, 6), params)
 
-    # Initial penalty is MAX_PENALTY, so one unit of time warp should be
+    # Initial penalty is max_penalty, so one unit of time warp should be
     # penalised by that value.
-    assert_equal(pm.cost_evaluator().tw_penalty(1), PenaltyManager.MAX_PENALTY)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), params.max_penalty)
 
     infeas = Solution(ok_small, [[1, 2, 3, 4]])
     assert_(infeas.has_time_warp())
@@ -317,27 +330,27 @@ def test_max_min_penalty(ok_small):
     # up by two times due to the penalty_increase parameter. But it's already
     # at the upper limit, and can thus not increase further.
     pm.register(infeas)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), PenaltyManager.MAX_PENALTY)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), params.max_penalty)
 
     feas = Solution(ok_small, [[1, 2], [3, 4]])
     assert_(not feas.has_time_warp())
 
     # But when we register a feasible solution, the time warp penalty parameter
     # should drop to zero due to the penalty_decrease parameter. But penalty
-    # parameters cannot drop below MIN_PENALTY.
+    # parameters cannot drop below min_penalty.
     pm.register(feas)
     assert_allclose(
-        pm.cost_evaluator().tw_penalty(10), 10 * PenaltyManager.MIN_PENALTY
+        pm.cost_evaluator().tw_penalty(10), 10 * params.min_penalty
     )
 
 
 def test_warns_max_penalty_value(ok_small):
     """
-    Tests that a penalty parameter clipped to MAX_PENALTY raises a warning.
+    Tests that a penalty parameter clipped to max_penalty raises a warning.
     This typically indicates a data issue that PyVRP is struggling with.
     """
     params = PenaltyParams(solutions_between_updates=1)
-    initial = ([1], PenaltyManager.MAX_PENALTY, 1)
+    initial = ([1], params.max_penalty, 1)
     pm = PenaltyManager(initial, params)
     assert_equal(pm.penalties(), initial)
 
@@ -400,18 +413,19 @@ def test_init_from_different_unit_costs(ok_small):
 
 def test_init_clips_penalties():
     """
-    Tests that the initial penalty values are clipped to the [MIN_PENALTY,
-    MAX_PENALTY] range.
+    Tests that the initial penalty values are clipped to the [min_penalty,
+    max_penalty] range.
     """
-    penalties = ([0], PenaltyManager.MAX_PENALTY + 1, 2)
-    pm = PenaltyManager(initial_penalties=penalties)
+    params = PenaltyParams()
+    penalties = ([0], params.max_penalty + 1, 2)
+    pm = PenaltyManager(initial_penalties=penalties, params=params)
 
     cost_eval = pm.cost_evaluator()
     assert_allclose(  # MIN
         cost_eval.load_penalty(10, 0, 0),
-        10 * PenaltyManager.MIN_PENALTY,
+        10 * params.min_penalty,
     )
-    assert_equal(cost_eval.tw_penalty(1), PenaltyManager.MAX_PENALTY)  # MAX
+    assert_equal(cost_eval.tw_penalty(1), params.max_penalty)  # MAX
     assert_equal(cost_eval.dist_penalty(1, 0), 2)  # already OK, so unchanged