Skip to content

Commit 7294e65

Browse files
daniprecdaniprec
authored
Merge pull request #62 from IEResearchDatalab/feat/swopp3-final-results
feat: final SWOPP3 result improvements (hourly ERA5, weather constraints, EDT land penalty)
2 parents 2d29822 + a989627 commit 7294e65

37 files changed

Lines changed: 4252 additions & 163 deletions

‎codabench/scoring_program/scoring.py‎

Lines changed: 279 additions & 53 deletions
Large diffs are not rendered by default.

‎docs/swopp3_sweep_results.md‎

Lines changed: 178 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,178 @@
1+
# SWOPP3 Parameter Sweep Results
2+
3+
**Date:** 2026-03-19
4+
**Branch:** `feat/swopp3-final-results`
5+
**Route:** Pacific noWPS (no waypoint system)
6+
**Operational constraints:** TWS ≤ 20 m/s, Hs ≤ 7.0 m
7+
8+
---
9+
10+
## How the Weather Penalty Works
11+
12+
The smooth weather penalty for a batch of routes is:
13+
14+
$$
15+
\text{penalty\_total} = \texttt{wpw} \cdot \sum_{i} \texttt{sharpness} \cdot \max(0,\; x_i - \text{limit})^2
16+
$$
17+
18+
where $x_i$ is TWS or Hs at segment midpoint $i$.
19+
20+
Two parameters control the penalty magnitude:
21+
22+
| Parameter | Symbol | Role |
23+
|-----------|--------|------|
24+
| `weather_penalty_weight` (wpw) | outer multiplier | **Global scale** — multiplies the entire penalty term before it is added to the energy cost. Controls how much the penalty matters *relative to fuel cost*. |
25+
| `sharpness` | inner multiplier | **Violation sensitivity** — multiplies each individual squared excess *before* summation. Controls how steeply the penalty ramps up *per segment* as conditions worsen beyond the threshold. |
26+
27+
They are mathematically interchangeable in a single-field scenario
28+
($\texttt{wpw} \cdot \texttt{sharpness}$ acts as one effective multiplier),
29+
but they decouple when both wind and wave penalties are active: `sharpness`
30+
scales each field's excess identically (same ramp steepness for TWS and Hs),
31+
while `wpw` scales the combined total against the fuel cost.
32+
33+
**In practice for this sweep:**
34+
- Higher **wpw** → optimizer sees weather avoidance as more important than fuel savings → routes detour more.
35+
- Higher **sharpness** → small exceedances above the limit are punished more aggressively → optimizer reacts earlier to marginal violations.
36+
37+
---
38+
39+
## Stage A — Catastrophic Detour Departures
40+
41+
**Goal:** Find penalty parameters that eliminate the massive energy
42+
over-spend (+30 to +99%) seen on stormy Pacific crossings.
43+
44+
**Setup:**
45+
- 10 departures with worst delta%: [24, 337, 338, 339, 340, 341, 348, 349, 362, 364]
46+
- Grid: `wpw ∈ {5, 10, 20, 40, 100}` × `sharpness ∈ {1, 2, 5}` × `σ₀ ∈ {0.1, 0.2}`
47+
- 30 configs × 10 departures = **300 runs** completed in 622 s (~2 s/run)
48+
49+
### Summary by Configuration
50+
51+
| wpw | sharp | σ₀ | mean Δ% | #viol TWS | #viol Hs | #any viol | mean TWS | mean Hs |
52+
|----:|------:|---:|--------:|----------:|---------:|----------:|---------:|--------:|
53+
| 5 | 1 | 0.2 | **-5.0** | 5 | 10 | 10 | 19.5 | 7.92 |
54+
| 5 | 2 | 0.2 | -3.0 | 5 | 9 | 9 | 19.4 | 7.72 |
55+
| 10 | 1 | 0.2 | -3.0 | 3 | 9 | 9 | 19.4 | 7.72 |
56+
| 5 | 1 | 0.1 | -2.6 | 5 | 9 | 9 | 19.9 | 8.11 |
57+
| 10 | 1 | 0.1 | 0.0 | 5 | 8 | 9 | 19.7 | 7.66 |
58+
| 5 | 2 | 0.1 | 0.3 | 4 | 8 | 8 | 19.7 | 7.52 |
59+
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
60+
| 100 | 5 | 0.1 | **+24.6** | 1 | 7 | 7 | 19.0 | 7.22 |
61+
| 100 | 2 | 0.1 | +25.4 | 2 | 7 | 7 | 18.7 | 7.24 |
62+
63+
### Best Config vs Baseline — Per Departure
64+
65+
| Dep | Baseline (wpw=100,s=5,σ=0.1) Δ% | Baseline viol | Best (wpw=5,s=1,σ=0.2) Δ% | Best viol |
66+
|----:|---:|---|---:|---|
67+
| 24 | +28.2 | Hs | -21.6 | Hs |
68+
| 337 | +18.2 | TWS+Hs | -5.7 | TWS+Hs |
69+
| 338 | +34.9 | Hs | +4.6 | Hs |
70+
| 339 | +35.2 | Hs | +12.0 | TWS+Hs |
71+
| 340 | +98.7 | Hs | +10.6 | TWS+Hs |
72+
| 341 | +66.2 | Hs | +16.9 | Hs |
73+
| 348 | +18.1 | **none** | +7.0 | Hs |
74+
| 349 | +8.2 | Hs | +5.5 | Hs |
75+
| 362 | -30.1 | none | -38.7 | TWS+Hs |
76+
| 364 | -31.5 | none | -40.2 | TWS+Hs |
77+
78+
### Key Observations
79+
80+
1. **Lower wpw → dramatically lower energy cost.** The best config
81+
(wpw=5, sharpness=1, σ₀=0.2) averages −5.0% vs GC, compared to
82+
+24.6% for baseline. That is a **29.6 pp improvement**.
83+
84+
2. **σ₀=0.2 consistently outperforms σ₀=0.1** — the larger initial
85+
CMA-ES step size helps exploration in these stormy departures.
86+
87+
3. **Weather violations are largely unavoidable for these 10 departures.**
88+
Even the strongest penalty (wpw=100) still has 7/10 Hs violations.
89+
The optimizer is forced through heavy weather regardless — the
90+
penalty just makes it take longer, costlier detours that still
91+
violate Hs.
92+
93+
4. **TWS violations increase with weaker penalty.** baseline has 1 TWS
94+
violation; wpw=5 has 5. If TWS compliance is critical, wpw=10
95+
(3 violations) is a safer middle ground.
96+
97+
5. **Dep 348 is the cautionary case**: baseline keeps it violation-free,
98+
but wpw=5 introduces an Hs violation. The weak penalty trades
99+
compliance for efficiency.
100+
101+
---
102+
103+
## Stage B — Zero-Delta (GC-Sticking) Departures
104+
105+
**Goal:** Find CMA-ES exploration parameters that improve routes for
106+
departures currently converging exactly to the great circle (0% delta).
107+
108+
**Setup:**
109+
- 20 sampled zero-delta departures: [44, 45, 46, 111, 113, 116, 118, 125,
110+
143, 159, 160, 169, 174, 177, 185, 186, 232, 233, 239, 274]
111+
- Grid: `σ₀ ∈ {0.1, 0.3, 0.5}` × `popsize ∈ {200, 400}` × `maxfevals ∈ {25k, 50k}` × `K ∈ {10, 15}`
112+
- 24 configs × 20 departures = **480 runs** completed in 1294 s (~2.7 s/run)
113+
114+
### Summary by Configuration
115+
116+
| σ₀ | pop | maxfevals | K | mean Δ% | median Δ% | min Δ% | #any viol |
117+
|---:|----:|----------:|--:|--------:|----------:|-------:|----------:|
118+
| 0.3 | 200 | 50000 | 15 | **-0.7** | 0.0 | -7.6 | 0 |
119+
| 0.5 | 200 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
120+
| 0.3 | 400 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
121+
| 0.5 | 400 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
122+
| 0.3 | 200 | 25000 | 15 | -0.7 | 0.0 | -7.5 | 0 |
123+
| 0.5 | 200 | 25000 | 15 | -0.7 | 0.0 | -7.3 | 0 |
124+
| 0.1 | 200 | 25000 | 15 | -0.4 | 0.0 | -7.6 | 0 |
125+
| ... | ... | ... | ... | ... | ... | ... | ... |
126+
| 0.1 | 200 | 25000 | 10 | **0.0** | 0.0 | 0.0 | 0 |
127+
| 0.1 | 400 | 25000 | 10 | 0.0 | 0.0 | 0.0 | 0 |
128+
129+
### Key Observations
130+
131+
1. **K=15 is the dominant factor.** Every K=15 config beats every K=10
132+
config. More Bézier control points give the curve freedom to deviate
133+
from GC where beneficial.
134+
135+
2. **σ₀ ≥ 0.3 gives a small extra boost** (−0.7% vs −0.4% for σ₀=0.1
136+
with K=15). Larger initial step helps explore alternatives to the
137+
GC.
138+
139+
3. **popsize and maxfevals do not matter** for these calm departures.
140+
200 pop / 25k fevals performs identically to 400 pop / 50k fevals.
141+
142+
4. **Zero weather violations across all 480 runs.** These departures
143+
sail through benign weather regardless of parameters.
144+
145+
5. **Only 2–3 departures respond to tuning** (deps 111, 169 show ~7%
146+
improvements; the rest stay at 0.0%). Most zero-delta departures are
147+
genuinely GC-optimal.
148+
149+
---
150+
151+
## Combined Recommendations
152+
153+
### Parameter Selection
154+
155+
| Parameter | Recommended | Rationale |
156+
|-----------|-------------|-----------|
157+
| **K** | **15** | Consistently improves zero-delta departures; no downside |
158+
| **σ₀** | **0.3** | Best tradeoff: helps Stage A detours and Stage B exploration |
159+
| **popsize** | **200** | 400 shows no benefit; saves compute |
160+
| **maxfevals** | **25000** | 50k shows no benefit on these departures |
161+
| **wpw** | **10** | Compromise: mean Δ% ≈ 0% (vs +24.6% baseline), TWS violations drop from 5 to 3 |
162+
| **sharpness** | **1** | Lower sharpness reduces over-reaction to marginal exceedances |
163+
164+
### Tradeoff: Energy vs Compliance
165+
166+
The Stage A results reveal a **fundamental tension**:
167+
168+
- These 10 departures cross inherently stormy regions where Hs > 7 m
169+
is unavoidable (7+ violations even at wpw=100).
170+
- Strong penalties force expensive detours that still violate Hs.
171+
- Weak penalties accept similar violation levels but at much lower
172+
energy cost.
173+
174+
**The weather penalty cannot solve the compliance problem for these
175+
departures** — it can only control how much extra fuel the ship burns
176+
while still violating constraints. This suggests that for truly stormy
177+
crossings, a **departure delay or alternative port** strategy may be
178+
more effective than in-route weather avoidance.

‎routetools/cmaes.py‎

Lines changed: 78 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
from routetools.land import Land
1818
from routetools.vectorfield import vectorfield_fourvortices
1919
from routetools.weather import weather_penalty as _weather_penalty
20+
from routetools.weather import weather_penalty_smooth as _weather_penalty_smooth
2021

2122

2223
@jit # type: ignore[misc]
@@ -288,7 +289,11 @@ def _cma_evolution_strategy(
288289
]
289290
| None = None,
290291
penalty: float = 1e10,
292+
land_distance_weight: float = 0.0,
293+
land_distance_epsilon: float = 1.0,
291294
weather_penalty_weight: float = 0.0,
295+
weather_penalty_type: str = "hard",
296+
weather_penalty_sharpness: float = 5.0,
292297
tws_limit: float = 20.0,
293298
hs_limit: float = 7.0,
294299
travel_stw: float | None = None,
@@ -310,21 +315,25 @@ def _cma_evolution_strategy(
310315
cost_fn: Callable[[jnp.ndarray], jnp.ndarray] | None = None,
311316
land_margin: int = 0,
312317
verbose: bool = True,
318+
bounds: list[list[float]] | None = None,
313319
**kwargs: dict[str, Any],
314320
) -> cma.CMAEvolutionStrategy:
315321
curve: jnp.ndarray
316322
# Initialize the optimizer
323+
inopts: dict[str, Any] = {
324+
"popsize": popsize,
325+
"tolfun": tolfun,
326+
"maxfevals": maxfevals,
327+
"seed": seed,
328+
"CSA_dampfac": damping, # v positive multiplier for step-size damping
329+
}
330+
if bounds is not None:
331+
inopts["bounds"] = bounds
332+
inopts |= kwargs
317333
es = cma.CMAEvolutionStrategy(
318334
x0,
319335
sigma0,
320-
inopts={
321-
"popsize": popsize,
322-
"tolfun": tolfun,
323-
"maxfevals": maxfevals,
324-
"seed": seed,
325-
"CSA_dampfac": damping, # v positive multiplier for step-size damping
326-
}
327-
| kwargs,
336+
inopts=inopts,
328337
)
329338
# Check if the land penalization is consistent
330339
if land is not None:
@@ -379,18 +388,33 @@ def _cma_evolution_strategy(
379388
# toward fewer land points.
380389
cost = jnp.where(has_land, penalty + land_count, cost)
381390

391+
# Smooth distance-to-land penalty via EDT
392+
if land is not None and land_distance_weight > 0:
393+
cost += land.distance_penalty(
394+
curve, weight=land_distance_weight, epsilon=land_distance_epsilon
395+
)
396+
382397
# Weather constraint penalization
383398
if weather_penalty_weight > 0 and (
384399
windfield is not None or wavefield is not None
385400
):
386-
cost += _weather_penalty(
387-
curve,
401+
_wp_fn = (
402+
_weather_penalty_smooth
403+
if weather_penalty_type == "smooth"
404+
else _weather_penalty
405+
)
406+
_wp_kwargs: dict[str, Any] = dict(
388407
windfield=windfield,
389408
wavefield=wavefield,
390409
tws_limit=tws_limit,
391410
hs_limit=hs_limit,
392411
penalty=weather_penalty_weight,
412+
travel_time=travel_time,
413+
time_offset=time_offset,
393414
)
415+
if weather_penalty_type == "smooth":
416+
_wp_kwargs["sharpness"] = weather_penalty_sharpness
417+
cost += _wp_fn(curve, **_wp_kwargs)
394418

395419
# Replace the worst solutions with the best found so far
396420
if keep_top > 0 and es.countiter > 1:
@@ -435,7 +459,11 @@ def optimize(
435459
]
436460
| None = None,
437461
penalty: float = 1e10,
462+
land_distance_weight: float = 0.0,
463+
land_distance_epsilon: float = 1.0,
438464
weather_penalty_weight: float = 0.0,
465+
weather_penalty_type: str = "hard",
466+
weather_penalty_sharpness: float = 5.0,
439467
tws_limit: float = 20.0,
440468
hs_limit: float = 7.0,
441469
travel_stw: float | None = None,
@@ -458,6 +486,7 @@ def optimize(
458486
cost_fn: Callable[[jnp.ndarray], jnp.ndarray] | None = None,
459487
land_margin: int = 0,
460488
verbose: bool = True,
489+
bounds: list[list[float]] | None = None,
461490
) -> tuple[jnp.ndarray, dict[str, Any]]:
462491
"""
463492
Solve the vessel routing problem for a given vector field.
@@ -488,9 +517,17 @@ def optimize(
488517
penalty : float, optional
489518
Large penalty applied to routes that intersect land (death-penalty
490519
scheme), by default 1e10
520+
land_distance_weight : float, optional
521+
Weight for the smooth distance-to-land penalty via EDT.
522+
Set to 0 (default) to disable.
523+
land_distance_epsilon : float, optional
524+
Regularisation constant for the EDT penalty (default 1.0).
491525
weather_penalty_weight : float, optional
492526
Penalty weight for weather constraint violations (TWS, Hs).
493527
Set to 0 (default) to disable weather penalties.
528+
weather_penalty_type : str, optional
529+
``"hard"`` (step function, default) or ``"smooth"`` (squared-ReLU
530+
ramp from :func:`weather_penalty_smooth`).
494531
tws_limit : float, optional
495532
Maximum allowed true wind speed in m/s, by default 20.0
496533
hs_limit : float, optional
@@ -529,6 +566,10 @@ def optimize(
529566
Random seed for reproducibility. By default jnp.nan
530567
verbose : bool, optional
531568
By default True
569+
bounds : list[list[float]] | None, optional
570+
Per-dimension ``[lower, upper]`` bounds for CMA-ES control-point
571+
parameters. Each list has length ``2*(K-2)``. ``None`` disables
572+
bounds (default).
532573
533574
Returns
534575
-------
@@ -594,7 +635,11 @@ def optimize(
594635
wavefield=wavefield,
595636
windfield=windfield,
596637
penalty=penalty,
638+
land_distance_weight=land_distance_weight,
639+
land_distance_epsilon=land_distance_epsilon,
597640
weather_penalty_weight=weather_penalty_weight,
641+
weather_penalty_type=weather_penalty_type,
642+
weather_penalty_sharpness=weather_penalty_sharpness,
598643
tws_limit=tws_limit,
599644
hs_limit=hs_limit,
600645
travel_stw=travel_stw,
@@ -616,6 +661,7 @@ def optimize(
616661
cost_fn=cost_fn,
617662
land_margin=land_margin,
618663
verbose=verbose,
664+
bounds=bounds,
619665
)
620666
time_end = time.time()
621667
if verbose:
@@ -660,14 +706,23 @@ def optimize(
660706
if weather_penalty_weight > 0 and (
661707
windfield is not None or wavefield is not None
662708
):
663-
cost_initial += _weather_penalty(
664-
curve0[jnp.newaxis, :, :],
709+
_wp_fn_c0 = (
710+
_weather_penalty_smooth
711+
if weather_penalty_type == "smooth"
712+
else _weather_penalty
713+
)
714+
_wp_kwargs_c0: dict[str, Any] = dict(
665715
windfield=windfield,
666716
wavefield=wavefield,
667717
tws_limit=tws_limit,
668718
hs_limit=hs_limit,
669719
penalty=weather_penalty_weight,
670-
).item()
720+
travel_time=travel_time,
721+
time_offset=time_offset,
722+
)
723+
if weather_penalty_type == "smooth":
724+
_wp_kwargs_c0["sharpness"] = weather_penalty_sharpness
725+
cost_initial += _wp_fn_c0(curve0[jnp.newaxis, :, :], **_wp_kwargs_c0).item()
671726
if cost_initial < cost_best:
672727
warnings.warn(
673728
"[WARNING] The optimized curve has a higher cost "
@@ -706,7 +761,11 @@ def optimize_with_increasing_penalization(
706761
penalty_init: float = 0,
707762
penalty_increment: float = 10,
708763
maxiter: int = 10,
764+
land_distance_weight: float = 0.0,
765+
land_distance_epsilon: float = 1.0,
709766
weather_penalty_weight: float = 0.0,
767+
weather_penalty_type: str = "hard",
768+
weather_penalty_sharpness: float = 5.0,
710769
tws_limit: float = 20.0,
711770
hs_limit: float = 7.0,
712771
travel_stw: float | None = None,
@@ -761,6 +820,8 @@ def optimize_with_increasing_penalization(
761820
weather_penalty_weight : float, optional
762821
Penalty weight for weather constraint violations (TWS, Hs).
763822
Set to 0 (default) to disable weather penalties.
823+
weather_penalty_type : str, optional
824+
``"hard"`` (step function, default) or ``"smooth"`` (squared-ReLU).
764825
tws_limit : float, optional
765826
Maximum allowed true wind speed in m/s, by default 20.0
766827
hs_limit : float, optional
@@ -823,7 +884,11 @@ def optimize_with_increasing_penalization(
823884
wavefield=wavefield,
824885
windfield=windfield,
825886
penalty=penalty,
887+
land_distance_weight=land_distance_weight,
888+
land_distance_epsilon=land_distance_epsilon,
826889
weather_penalty_weight=weather_penalty_weight,
890+
weather_penalty_type=weather_penalty_type,
891+
weather_penalty_sharpness=weather_penalty_sharpness,
827892
tws_limit=tws_limit,
828893
hs_limit=hs_limit,
829894
travel_stw=travel_stw,

0 commit comments

Comments
 (0)