1"""Stop ordering for one van: UCS, Greedy best-first and A*, plus brute force.2 3State (current stop, set of visited stops, clock time)4Action drive to an unvisited stop j, wait if early, serve it5Cost g = clock - shift start = predicted travel + waiting + service6 (service is a constant, so this is "travel + waiting" up to a constant);7 with return_to_depot the final leg back to the depot is included.8Hard mode arriving after a window closes is not allowed (successor pruned).9Soft mode late service allowed, g += LATE_PENALTY * minutes late (used only to find the10 least-late route when no on-time route exists, to locate the "first late pair").11 12Duplicate detection / dominance13 Travel times are FIFO (see travel_time.py) and waiting is allowed, so for the same14 (current stop, visited set) a state with an earlier clock (and, in soft mode, no more15 lateness) can do everything the later one can, at least as cheaply. We keep only the16 Pareto-optimal (clock, lateness) labels per (stop, visited) -- in hard mode a single label.17 18A* heuristic (admissible, proof in report section 4.2)19 LB[u,k] = min over every departure time >= now of the predicted time u -> k.20 h_sum = sum over unvisited k of ( min_{u in {cur} U unvisited, u != k} LB[u,k] + service_k )21 + min_{u in unvisited} LB[u, depot] (if returning)22 Every unvisited stop is entered by exactly one distinct leg, every leg costs at least its23 LB, waiting is >= 0, so h_sum <= true remaining cost.24 h_wait = max over unvisited k of ( max(e_k, now + min_u LB[u,k]) + service_k ) - now (+ return LB)25 Each stop k is still to be served, no earlier than its window opens: also a lower bound.26 h = max(h_sum, h_wait) is admissible (max of admissible heuristics).27 The "naive" version (cheapest incoming edge at the CURRENT hour) is NOT admissible when28 traffic gets lighter later in the day or the van must wait for a window - we keep it only29 to demonstrate that in the tests.30 A* re-opens states when a cheaper path is found, so admissibility alone gives optimality.31"""32from __future__ import annotations33 34import heapq35import itertools36import time37from dataclasses import dataclass, field38 39import numpy as np40 41LATE_PENALTY = 1000.042DEPOT = -143 44 45@dataclass46class Route:47 van: int48 sequence: list # order ids in visiting order49 arrivals: list50 begins: list # service start (after waiting)51 departures: list52 waits: list53 lateness: list54 finish: float # back at depot (or end of last service)55 cost: float # finish - start (+ penalty in soft mode)56 feasible: bool57 stats: dict = field(default_factory=dict)58 59 @property60 def n_late(self):61 return int(sum(1 for x in self.lateness if x > 1e-9))62 63 def first_late(self):64 for k, x in enumerate(self.lateness):65 if x > 1e-9:66 return k67 return None68 69 70@dataclass71class RouteProblem:72 stops: list # Order objects73 depot: int74 start: float75 tm: object # TimeModel76 return_to_depot: bool = True77 shift_end: float = float("inf")78 van: int = 079 80 81# --------------------------------------------------------------------------- simulation82def simulate(problem: RouteProblem, seq_local, travel=None) -> Route:83 """Drive a fixed sequence (local indices). `travel(i_zone, j_zone, t)` defaults to the model."""84 travel = travel or problem.tm.tt85 t, cur = problem.start, problem.depot86 arr, beg, dep, wai, late = [], [], [], [], []87 for k in seq_local:88 o = problem.stops[k]89 a = t + travel(cur, o.zone, t)90 b = max(a, o.earliest)91 arr.append(a); beg.append(b); wai.append(b - a); late.append(max(0.0, b - o.latest))92 t = b + o.service93 dep.append(t)94 cur = o.zone95 if problem.return_to_depot and seq_local:96 t = t + travel(cur, problem.depot, t)97 over = max(0.0, t - problem.shift_end)98 feasible = all(x <= 1e-9 for x in late) and over <= 1e-999 cost = t - problem.start + LATE_PENALTY * (sum(late) + over)100 return Route(problem.van, [problem.stops[k].id for k in seq_local], arr, beg, dep, wai, late,101 t, cost, feasible)102 103 104def brute_force(problem: RouteProblem, hard: bool = True) -> Route | None:105 n = len(problem.stops)106 best, t0, count = None, time.perf_counter(), 0107 for perm in itertools.permutations(range(n)):108 count += 1109 r = simulate(problem, perm)110 if hard and not r.feasible:111 continue112 if best is None or r.cost < best.cost - 1e-9:113 best = r114 if best is not None:115 best.stats = {"algo": "brute-force", "permutations": count,116 "time_s": time.perf_counter() - t0}117 return best118 119 120# --------------------------------------------------------------------------- heuristics121def make_heuristic(problem: RouteProblem, kind: str):122 stops, tm = problem.stops, problem.tm123 n = len(stops)124 zones = np.array([o.zone for o in stops], dtype=int)125 serv = np.array([o.service for o in stops], dtype=float)126 early = np.array([o.earliest for o in stops], dtype=float)127 ret = problem.return_to_depot128 129 def h(cur: int, mask: int, clock: float) -> float:130 if kind == "zero":131 return 0.0132 un = [k for k in range(n) if not mask >> k & 1]133 if not un:134 return 0.0135 cur_zone = problem.depot if cur == DEPOT else stops[cur].zone136 LB = tm.M[tm._knot(clock)[0]] if kind == "naive" else tm.lower_bound_matrix(clock)137 uz = zones[un]138 # incoming into each unvisited k from cur or another unvisited stop139 inc = LB[np.append(uz, cur_zone)][:, uz] # rows: sources (un..., cur)140 if len(un) > 1:141 sub = inc[:-1].copy()142 np.fill_diagonal(sub, np.inf) # u != k143 min_in = np.minimum(sub.min(axis=0), inc[-1])144 else:145 min_in = inc[-1]146 ret_lb = float(LB[uz, problem.depot].min()) if ret else 0.0147 if kind == "max_in":148 return float((min_in + serv[un]).max()) + ret_lb149 h_sum = float(min_in.sum() + serv[un].sum()) + ret_lb150 if kind == "naive":151 return h_sum152 h_wait = float((np.maximum(early[un], clock + min_in) + serv[un]).max()) - clock + ret_lb153 return max(h_sum, h_wait)154 155 return h156 157 158# --------------------------------------------------------------------------- best-first search159def solve_route(problem: RouteProblem, algo: str = "astar", heuristic: str = "admissible",160 hard: bool = True, trace: list | None = None, trace_limit: int = 40,161 max_expansions: int = 2_000_000) -> Route | None:162 """algo in {'ucs', 'greedy', 'astar'}. Returns the best Route or None (no feasible order)."""163 t0 = time.perf_counter()164 stops, tm, n = problem.stops, problem.tm, len(problem.stops)165 full = (1 << n) - 1166 hfun = make_heuristic(problem, "zero" if algo == "ucs" else heuristic)167 stats = {"algo": algo, "heuristic": heuristic if algo != "ucs" else "zero", "expanded": 0,168 "generated": 0, "pruned_infeasible": 0, "pruned_dominated": 0, "max_frontier": 0}169 170 # node: (cur, mask, clock, late, g, parent_id)171 nodes = []172 labels: dict = {} # (cur, mask) -> list of (clock, late, node_id)173 alive = set()174 frontier = []175 tie = itertools.count()176 177 def push(cur, mask, clock, late, g, parent, terminal=False):178 key = (cur, mask, terminal)179 lab = labels.setdefault(key, [])180 for (c2, l2, _id) in lab:181 if c2 <= clock + 1e-9 and l2 <= late + 1e-9: # an existing label dominates182 stats["pruned_dominated"] += 1183 return184 keep = []185 for (c2, l2, nid) in lab:186 if clock <= c2 + 1e-9 and late <= l2 + 1e-9:187 alive.discard(nid)188 else:189 keep.append((c2, l2, nid))190 nid = len(nodes)191 nodes.append((cur, mask, clock, late, g, parent, terminal))192 keep.append((clock, late, nid))193 labels[key] = keep194 alive.add(nid)195 hv = 0.0 if terminal else hfun(cur, mask, clock)196 pri = {"ucs": g, "astar": g + hv, "greedy": hv}[algo]197 heapq.heappush(frontier, (pri, g, next(tie), nid, hv))198 stats["generated"] += 1199 stats["max_frontier"] = max(stats["max_frontier"], len(frontier))200 201 def goal_transition(cur, mask, clock, late, g, nid):202 """All stops served: close the route (return leg, shift end)."""203 cz = problem.depot if cur == DEPOT else stops[cur].zone204 end = clock + (tm.tt(cz, problem.depot, clock) if problem.return_to_depot and n else 0.0)205 over = max(0.0, end - problem.shift_end)206 if hard and over > 1e-9:207 stats["pruned_infeasible"] += 1208 return209 push(cur, mask, end, late + over, end - problem.start + LATE_PENALTY * (late + over), nid,210 terminal=True)211 212 push(DEPOT, 0, problem.start, 0.0, 0.0, None)213 goal = None214 while frontier:215 pri, g, _, nid, hv = heapq.heappop(frontier)216 if nid not in alive:217 continue218 cur, mask, clock, late, g, parent, terminal = nodes[nid]219 if terminal:220 goal = nid221 break222 alive.discard(nid)223 stats["expanded"] += 1224 if trace is not None and len(trace) < trace_limit:225 trace.append({"step": stats["expanded"], "at": "depot" if cur == DEPOT else stops[cur].id,226 "visited": [stops[k].id for k in range(n) if mask >> k & 1],227 "clock": clock, "g": g, "h": hv, "f": g + hv})228 if stats["expanded"] > max_expansions:229 break230 if mask == full:231 goal_transition(cur, mask, clock, late, g, nid)232 continue233 cz = problem.depot if cur == DEPOT else stops[cur].zone234 for j in range(n):235 if mask >> j & 1:236 continue237 o = stops[j]238 a = clock + tm.tt(cz, o.zone, clock)239 b = max(a, o.earliest)240 lj = max(0.0, b - o.latest)241 if hard and lj > 1e-9:242 stats["pruned_infeasible"] += 1243 continue244 c2 = b + o.service245 push(j, mask | 1 << j, c2, late + lj, c2 - problem.start + LATE_PENALTY * (late + lj), nid)246 247 stats["time_s"] = time.perf_counter() - t0248 if goal is None:249 return None250 seq, k = [], nodes[goal][5]251 while k is not None and nodes[k][0] != DEPOT:252 seq.append(nodes[k][0])253 k = nodes[k][5]254 seq.reverse()255 route = simulate(problem, seq)256 route.stats = stats257 return route258 259 260def plan_van(problem: RouteProblem, algo: str = "astar", trace=None) -> Route:261 """Hard-window search first; if no on-time order exists, return the least-late route."""262 if not problem.stops:263 r = simulate(problem, [])264 r.stats = {"algo": algo, "expanded": 0, "generated": 0, "time_s": 0.0}265 return r266 r = solve_route(problem, algo=algo, hard=True, trace=trace)267 if r is not None:268 return r269 r = solve_route(problem, algo=algo, hard=False)270 r.stats["soft_fallback"] = True271 return r