Vehicle Platooning Optimization

Vehicle Platooning Optimization

You are an expert in coordinated truck platooning — the planning problem of routing and scheduling a fleet so that trucks share road segments at the same time, letting followers ride in the leader's slipstream and burn less fuel. This skill covers the fuel-saving objective, platoon formation on shared arcs, routing with bounded detours, scheduling with time windows and waiting, an exact MIP, a vectorized genetic algorithm for larger fleets, instance generation, and independent solution validation. Core references: Larsson, Sennton & Larson (2015), "The vehicle platooning problem: computational complexity and heuristics"; Bhoopalam, Agatz & Zuidwijk (2018), "Planning of truck platoons: a literature review and directions for future research"; Boysen, Briskorn & Schwerdfeger (2018), "The identical-path truck platooning problem". Use the framework below to take a user from a fleet description to a validated coordination plan with a defensible savings number.

Initial Assessment

Establish these facts before modeling anything:

• Fleet and network scale. How many trucks per planning run, and how large is the road network (nodes, arcs)? The pairwise platooning variables grow as O(|K|² × shared arcs); this single number decides exact vs heuristic.
• Routes: fixed or free? If every truck's path is already fixed (contracted lanes, single highway), the problem collapses to departure-time scheduling — far easier. If detours are allowed, bound them: a detour cap of 5–25% over the shortest path is typical, because detour fuel quickly eats the ~10% follower saving.
• Fuel-saving coefficients. What follower saving fraction η applies, and does the leader save anything? Field tests report follower savings of roughly 10–20% at close gaps and leader savings of a few percent (Bonnet & Fritz 2000, two electronically coupled trucks; Tsugawa 2014, Energy ITS three-truck platoon). Planning models most often use a flat η ≈ 0.10 for followers and 0 for the leader.
• Platoon formation rule. Must paired trucks enter a shared arc at exactly the same time (they waited for each other), within a tolerance δ (small en-route speed adaptation closes the gap, as in van de Hoef, Johansson & Dimarogonas 2018), or only at designated hubs? This choice shapes the synchronization constraints.
• Where is waiting allowed? At the origin only, at any intermediate node, or nowhere? Waiting at intermediate nodes is what lets trucks merge mid-route; forbid it and platooning opportunities shrink sharply.
• Time windows: hard or soft? Hard arrival deadlines kill platooning chances for tight trucks. Ask whether late arrival is forbidden or penalized, and what the per-hour cost of trip duration (driver wages, schedule slack) is — it competes directly with fuel savings.
• Platoon size cap. Most studies cap platoons at 3–5 vehicles (legal limits, braking safety, merge complexity). The cap matters: savings per arc scale with (size − 1).
• Leader compensation. If carriers differ, the leader saves nothing while followers profit — does the user need a savings-sharing scheme, or is this a single-fleet problem where only the total matters? This skill optimizes the total; flag profit allocation as a separate (cooperative game) question.
• Planning mode. One-shot offline plan, rolling horizon, or online (trucks already driving)? The MIP below is offline; rolling use re-solves with fixed past decisions.
• Solver availability. Gurobi license for the exact model and the LP schedule polish? Without it, the same model runs on HiGHS/CBC with smaller size limits.
• Data format. Real road network (OpenStreetMap via osmnx, contracted to a highway graph) or synthetic grid? Travel times constant per arc, or time-dependent? The models here assume constant arc times.
• Baseline for reporting. Savings must be reported against each truck's cheapest feasible solo path, not against the chosen (possibly detoured) routes. Fix this definition before any experiments.
• Quality requirement. Is a provably optimal coordination plan required (small instances, benchmarking) or is a good heuristic plan within minutes acceptable (operational use)?

Problem Definition and Model Landscape

Formal definition

Given a directed road network $G = (N, A)$ where arc $a \in A$ has length $d_a$ [km] and travel time $\tau_a$ [h], and a set of trucks $K$, truck $k$ has origin $o_k$, destination $s_k$, earliest departure $e_k$, and arrival deadline $\ell_k$. Each truck drives one path from a candidate set $P_k$ (the $|P_k|$ shortest paths within a detour cap $\Delta$ times the shortest distance); $d(p)$ and $\tau(p)$ are path length and drive time. Fuel cost is $c$ per km solo. When trucks travel an arc together as a platoon, one member leads and the rest follow; each follower's fuel on that arc is reduced by the fraction $\eta$. Platoons hold at most $L$ vehicles, and paired trucks must enter the arc within a tolerance $\delta \ge 0$ of each other. Trip duration (driver time) costs $\gamma$ per hour. $T$ is the planning horizon.

Decisions: which path each truck drives, when it enters every arc on that path (waiting at nodes is allowed — the slack between consecutive entries beyond drive time), and which leader-follower pairs form on which arcs.

Variables: $z_{kp} \in {0,1}$ (truck $k$ drives path $p$), $t_{ka} \ge e_k$ (entry time of $k$ into arc $a$), $w_{lka} \in {0,1}$ ($l$ leads $k$ on arc $a$), and $dep_k, arr_k$ (departure and arrival times). The objective:

$$ \min ;; c \sum_{k \in K} \sum_{p \in P_k} d(p), z_{kp} ;-; \eta, c \sum_{(l,k,a)} d_a, w_{lka} ;+; \gamma \sum_{k \in K} \left( arr_k - dep_k \right) $$

subject to seven constraint families:

1. Path selection. $\sum_{p \in P_k} z_{kp} = 1$ for every truck $k$.
2. Time propagation. For consecutive arcs $a \to b$ on path $p$: $t_{kb} \ge t_{ka} + \tau_a - M(1 - z_{kp})$. The inequality's slack is waiting time at the intermediate node.
3. Time windows. $t_{ka} \ge e_k$ (variable lower bound) and $arr_k \le \ell_k$ (variable upper bound), with $arr_k \ge t_{k,\mathrm{last}(p)} + \tau_{\mathrm{last}(p)} - M(1 - z_{kp})$ linking arrival to the chosen path.
4. Route linking. $w_{lka} \le \sum_{p \in P_k : a \in p} z_{kp}$ and $w_{lka} \le \sum_{q \in P_l : a \in q} z_{lq}$ — pairing only on arcs both trucks actually drive.
5. One leader per follower. $\sum_{l \ne k} w_{lka} \le 1$ for every $(k, a)$.
6. Two-level platoons and size cap. $\sum_{k} w_{lka} \le (L-1)\bigl(1 - \sum_{m} w_{mla}\bigr)$: a truck that follows on arc $a$ leads nobody there, and a leader has at most $L-1$ followers. This is the standard pairwise leader-follower structure of Larsson, Sennton & Larson (2015).
7. Synchronization. $|t_{ka} - t_{la}| \le \delta + M(1 - w_{lka})$: paired trucks enter the arc together (within tolerance).

A platoon of size $m$ on arc $a$ saves $(m-1),\eta, c, d_a$; the leader pays full price (leader savings are an extension — see Advanced Techniques). Model size: $\sum_k |P_k|$ path variables, $\sum_k |A_k|$ time variables ($A_k$ = arcs on any candidate path of $k$), and up to $|K|(|K|-1) \times$ (average shared arcs) pairing variables. The pairing variables are the bottleneck.

Complexity

The platooning problem is NP-hard, and remains NP-hard even without deadlines (the "unlimited platooning problem"); a constant number of trucks admits polynomial algorithms (Larsson, Sennton & Larson 2015). The identical-path special case — all trucks drive the same route, choose departure times only — has polynomially solvable variants and its own structure (Boysen, Briskorn & Schwerdfeger 2018). In practice, exact MIPs of the form above are reported solvable for tens of trucks on metropolitan-scale networks (Larson, Munson & Sokolov 2016, coordinated platoon routing in a metropolitan network); beyond that, heuristics take over.

Variant taxonomy

Decision guidance

• Routes fixed and identical → identical-path machinery. Do not deploy the general MIP on it.
• ≤ ~15–25 trucks, candidate paths ≤ 5, planning offline → exact MIP. Provable optimality, and the instance sizes where the pairwise model stays tractable.
• Hundreds of trucks or online re-planning → GA + LP schedule polish (below), or a matheuristic that runs the MIP on truck clusters.
• Sparse overlap (trucks rarely share corridors) → filter first. Compute pairwise corridor overlap and time-window compatibility; most $w$ variables can be pruned before any model is built.
• Need savings-sharing across carriers → optimize total savings first, then allocate; do not bake allocation into the routing model.

Instance Generator and Data Model

Synthetic grid networks with random arc lengths reproduce the structure that matters — multiple O-D pairs whose shortest paths overlap on central corridors — while staying fully seeded and reproducible. Candidate paths come from Yen's k-shortest-paths (networkx.shortest_simple_paths) truncated by the detour cap.

import networkx as nx
import numpy as np
from dataclasses import dataclass


@dataclass(frozen=True)
class Truck:
    """One truck: origin/destination nodes and a departure/arrival time window."""
    origin: int
    dest: int
    earliest_dep: float   # earliest departure time from origin [h]
    latest_arr: float     # hard arrival deadline at destination [h]


@dataclass
class PlatooningInstance:
    """Road network, trucks, candidate paths, and cost parameters."""
    graph: nx.DiGraph                       # arcs carry 'dist' [km] and 'time' [h]
    trucks: list[Truck]
    paths: list[list[tuple[int, ...]]]      # paths[k][p] = node sequence of a candidate path
    fuel_per_km: float = 0.55               # solo fuel cost [EUR/km]
    eta_follow: float = 0.10                # follower fuel-saving fraction
    wait_cost: float = 2.0                  # cost per hour of trip duration [EUR/h]
    max_platoon: int = 4                    # vehicles per platoon (1 leader + followers)
    sync_tol: float = 0.0                   # entry-time tolerance for platooning [h]

    def path_arcs(self, k: int, p: int) -> list[tuple[int, int]]:
        """Arc sequence of candidate path p of truck k."""
        nodes = self.paths[k][p]
        return list(zip(nodes[:-1], nodes[1:]))

    def arc_dist(self, a: tuple[int, int]) -> float:
        """Arc length in km."""
        return self.graph.edges[a]["dist"]

    def arc_time(self, a: tuple[int, int]) -> float:
        """Arc travel time in hours."""
        return self.graph.edges[a]["time"]


def generate_instance(n_rows: int, n_cols: int, n_trucks: int, n_paths: int,
                      detour_cap: float, seed: int) -> PlatooningInstance:
    """Random grid road network with seeded trucks and k-shortest candidate paths.

    Every truck's shortest path is deadline-feasible by construction (slack >= 1.2x
    the shortest drive time), so the instance always has a solo-feasible solution.
    """
    rng = np.random.default_rng(seed)
    g = nx.DiGraph()
    for r in range(n_rows):
        for c in range(n_cols):
            u = r * n_cols + c
            for dr, dc in ((0, 1), (1, 0)):
                r2, c2 = r + dr, c + dc
                if r2 < n_rows and c2 < n_cols:
                    v = r2 * n_cols + c2
                    dist = float(rng.uniform(30.0, 80.0))               # km
                    for s, t2 in ((u, v), (v, u)):
                        g.add_edge(s, t2, dist=dist, time=dist / 80.0)  # 80 km/h
    horizon = 24.0
    trucks: list[Truck] = []
    paths: list[list[tuple[int, ...]]] = []
    while len(trucks) < n_trucks:
        o, d = rng.choice(n_rows * n_cols, size=2, replace=False)
        sp_time = nx.shortest_path_length(g, int(o), int(d), weight="time")
        if sp_time < 0.5:                       # skip trivial neighbor-to-neighbor trips
            continue
        e_k = float(rng.uniform(0.0, 4.0))
        slack = float(rng.uniform(1.2, 1.6))
        trucks.append(Truck(int(o), int(d), e_k, min(horizon, e_k + slack * sp_time)))
        sp_dist = nx.shortest_path_length(g, int(o), int(d), weight="dist")
        cand: list[tuple[int, ...]] = []
        for path in nx.shortest_simple_paths(g, int(o), int(d), weight="dist"):
            length = sum(g.edges[a]["dist"] for a in zip(path[:-1], path[1:]))
            if length > detour_cap * sp_dist or len(cand) == n_paths:
                break
            cand.append(tuple(path))
        paths.append(cand)
    return PlatooningInstance(g, trucks, paths)


inst = generate_instance(n_rows=4, n_cols=4, n_trucks=6, n_paths=3,
                         detour_cap=1.25, seed=42)
print(len(inst.trucks), sum(len(p) for p in inst.paths))
# Expected: 6 trucks and roughly 12-18 candidate paths in total; the exact output is
# identical on every run with seed=42.

Scale the generator by raising n_rows/n_cols and n_trucks; hardness grows with truck density per corridor (more overlap → more pairing decisions) and with tighter detour_cap/deadline slack. For real road networks, replace the grid loop with an osmnx-extracted highway graph and keep everything else unchanged.

Exact MIP in Gurobi

The model follows the formulation above. Each constraint family lives in its own builder function so families can be unit-tested, relaxed, or swapped independently — the construction discipline of milp-modeling-gurobi. Big-M values derive from the horizon; see Advanced Techniques for per-pair tightening, which is the single most effective model improvement.

import gurobipy as gp
from gurobipy import GRB

# Uses PlatooningInstance and generate_instance from the instance-generator block.


def truck_arcs(inst: PlatooningInstance, k: int) -> set[tuple[int, int]]:
    """All arcs that appear on at least one candidate path of truck k."""
    return {a for p in range(len(inst.paths[k])) for a in inst.path_arcs(k, p)}


def shared_arc_pairs(inst: PlatooningInstance) -> dict[tuple[int, int], list[tuple[int, int]]]:
    """Ordered truck pairs (l, k), l != k, mapped to the arcs both trucks can traverse."""
    arcs = [truck_arcs(inst, k) for k in range(len(inst.trucks))]
    out: dict[tuple[int, int], list[tuple[int, int]]] = {}
    for l in range(len(inst.trucks)):
        for k in range(len(inst.trucks)):
            if l != k and (common := sorted(arcs[l] & arcs[k])):
                out[(l, k)] = common
    return out


def add_path_selection_constraints(model: gp.Model, v: dict, inst: PlatooningInstance) -> None:
    """Each truck commits to exactly one candidate path."""
    for k in range(len(inst.trucks)):
        model.addConstr(
            gp.quicksum(v["z"][k, p] for p in range(len(inst.paths[k]))) == 1,
            name=f"one_path[{k}]")


def add_schedule_constraints(model: gp.Model, v: dict, inst: PlatooningInstance,
                             horizon: float) -> None:
    """Entry-time propagation along the chosen path; the inequality slack is waiting.

    dep/arr are linked one-sided only: the duration term in the objective pushes
    dep up and arr down, so both bind exactly at any solver incumbent.
    """
    z, t, dep, arr = v["z"], v["t"], v["dep"], v["arr"]
    for k in range(len(inst.trucks)):
        for p in range(len(inst.paths[k])):
            arcs = inst.path_arcs(k, p)
            for a, b in zip(arcs[:-1], arcs[1:]):
                big_m = horizon + inst.arc_time(a)
                model.addConstr(
                    t[k, b] >= t[k, a] + inst.arc_time(a) - big_m * (1 - z[k, p]),
                    name=f"prop[{k},{p},{a[0]}_{a[1]}]")
            first, last = arcs[0], arcs[-1]
            model.addConstr(dep[k] <= t[k, first] + horizon * (1 - z[k, p]),
                            name=f"dep_link[{k},{p}]")
            model.addConstr(
                arr[k] >= t[k, last] + inst.arc_time(last)
                - (horizon + inst.arc_time(last)) * (1 - z[k, p]),
                name=f"arr_link[{k},{p}]")


def add_platoon_pairing_constraints(model: gp.Model, v: dict, inst: PlatooningInstance) -> None:
    """Route linking, one leader per follower, two-level platoons, and the size cap."""
    z, w = v["z"], v["w"]
    leaders_of: dict[tuple, list] = {}   # (k, a) -> [w[l,k,a] over leaders l]
    out_of: dict[tuple, list] = {}       # (l, a) -> [w[l,k,a] over followers k]
    for (l, k, a), var in w.items():
        p] for p in range(len(inst.paths[k]))
                           if a in inst.path_arcs(k, p))
        q] for q in range(len(inst.paths[l]))
                           if a in inst.path_arcs(l, q))
        model.addConstr(var <= on_k, name=f"route_f[{l},{k},{a[0]}_{a[1]}]")
        model.addConstr(var <= on_l, name=f"route_l[{l},{k},{a[0]}_{a[1]}]")
        leaders_of.setdefault((k, a), []).append(var)
        out_of.setdefault((l, a), []).append(var)
    for (k, a), vars_in in leaders_of.items():
        model.addConstr(gp.quicksum(vars_in) <= 1,
                        name=f"one_leader[{k},{a[0]}_{a[1]}]")
    for (l, a), vars_out in out_of.items():
        follows_l = gp.quicksum(leaders_of.get((l, a), []))
        model.addConstr(
            gp.quicksum(vars_out) <= (inst.max_platoon - 1) * (1 - follows_l),
            name=f"leader_cap[{l},{a[0]}_{a[1]}]")


def add_synchronization_constraints(model: gp.Model, v: dict, inst: PlatooningInstance,
                                    horizon: float) -> None:
    """Paired trucks enter the shared arc within sync_tol of each other (big-M)."""
    t, w = v["t"], v["w"]
    for (l, k, a), var in w.items():
        model.addConstr(t[k, a] - t[l, a] <= inst.sync_tol + horizon * (1 - var),
                        name=f"sync_p[{l},{k},{a[0]}_{a[1]}]")
        model.addConstr(t[l, a] - t[k, a] <= inst.sync_tol + horizon * (1 - var),
                        name=f"sync_m[{l},{k},{a[0]}_{a[1]}]")


def build_platooning_model(inst: PlatooningInstance,
                           horizon: float) -> tuple[gp.Model, dict]:
    """Assemble variables, all four constraint families, and the objective."""
    model = gp.Model("platooning")
    K = range(len(inst.trucks))
    z = {(k, p): model.addVar(vtype=GRB.BINARY, name=f"z[{k},{p}]")
         for k in K for p in range(len(inst.paths[k]))}
    t = {(k, a): model.addVar(lb=inst.trucks[k].earliest_dep, ub=horizon,
                              name=f"t[{k},{a[0]}_{a[1]}]")
         for k in K for a in truck_arcs(inst, k)}
    w = {(l, k, a): model.addVar(vtype=GRB.BINARY, name=f"w[{l},{k},{a[0]}_{a[1]}]")
         for (l, k), common in shared_arc_pairs(inst).items() for a in common}
    dep = {k: model.addVar(lb=inst.trucks[k].earliest_dep, ub=horizon,
                           name=f"dep[{k}]") for k in K}
    arr = {k: model.addVar(ub=inst.trucks[k].latest_arr, name=f"arr[{k}]") for k in K}
    v = {"z": z, "t": t, "w": w, "dep": dep, "arr": arr}
    add_path_selection_constraints(model, v, inst)
    add_schedule_constraints(model, v, inst, horizon)
    add_platoon_pairing_constraints(model, v, inst)
    add_synchronization_constraints(model, v, inst, horizon)
    base = gp.quicksum(
        inst.fuel_per_km * sum(inst.arc_dist(a) for a in inst.path_arcs(k, p)) * z[k, p]
        for k in K for p in range(len(inst.paths[k])))
    savings = gp.quicksum(inst.eta_follow * inst.fuel_per_km * inst.arc_dist(a) * var
                          for (l, k, a), var in w.items())
    duration = gp.quicksum(arr[k] - dep[k] for k in K)
    model.setObjective(base - savings + inst.wait_cost * duration, GRB.MINIMIZE)
    return model, v

Solving and extracting a structured solution — routes, explicit arc entry times, and the leader-follower pairs — keeps everything downstream (validation, plots, KPI tables) independent of gurobipy objects:

import gurobipy as gp
from gurobipy import GRB

# Uses generate_instance and build_platooning_model from the blocks above.


def solve_platooning(inst: PlatooningInstance, horizon: float = 24.0,
                     time_limit: float = 60.0) -> dict | None:
    """Solve the platooning MIP; extract routes, entry times, and platoon pairs."""
    model, v = build_platooning_model(inst, horizon)
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    model.Params.MIPGap = 1e-4
    model.optimize()
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        return None
    K = range(len(inst.trucks))
    route = {k: next(p for p in range(len(inst.paths[k])) if v["z"][k, p].X > 0.5)
             for k in K}
    entry = {(k, a): v["t"][k, a].X for k in K for a in inst.path_arcs(k, route[k])}
    pairs = [(l, k, a) for (l, k, a), var in v["w"].items() if var.X > 0.5]
    return {"objective": model.ObjVal, "bound": model.ObjBound,
            "route": route, "entry": entry, "pairs": pairs}


inst = generate_instance(n_rows=4, n_cols=4, n_trucks=6, n_paths=3,
                         detour_cap=1.25, seed=42)
solo_fuel = sum(
    min(inst.fuel_per_km * sum(inst.arc_dist(a) for a in inst.path_arcs(k, p))
        for p in range(len(inst.paths[k])))
    for k in range(len(inst.trucks)))
sol = solve_platooning(inst)
print(round(sol["objective"], 1), round(solo_fuel, 1), len(sol["pairs"]))
# Expected: optimal in well under a minute at this size. The fuel part of the objective
# sits a few percent below solo_fuel, and pairs lists leader-follower matches on shared
# corridor arcs (empty only if no two time windows overlap on a common arc).

With sync_tol = 0 the model forces exact simultaneous entry — achievable because trucks can wait at nodes. Deadline-infeasible candidate paths need no special handling: choosing one makes arr exceed its upper bound, so the solver simply never selects it.

Genetic Algorithm for Larger Fleets

Beyond ~20–30 trucks the pairwise MIP stalls, but the problem has a natural compact genome: per truck, a path gene (index into the deadline-feasible candidate paths) and a departure key in [0, 1] mapped to the feasible departure window of that path. Trucks then drive at constant speed with no intermediate waiting; platoons form wherever entry times on a shared arc fall within sync_tol, interpreted as bounded en-route speed adaptation (van de Hoef et al. 2018). The genome encodes no infeasible schedule — the key mapping makes every individual deadline-feasible by construction. Operator theory (selection pressure, crossover/mutation choices, population sizing) is in genetic-algorithms; this section only supplies the platooning decoder and fitness.

The fitness is fully vectorized across the population: one lexsort over (individual, arc, entry time) groups co-located trucks for all individuals at once. Grouping uses chain semantics (consecutive gaps ≤ tolerance), which slightly overestimates savings versus the leader-anchored rule the validator enforces — so the final answer is always the decoded and validated objective, never the raw fitness.

import numpy as np

# Uses PlatooningInstance and generate_instance from the instance-generator block.


def precompute_tables(inst: PlatooningInstance) -> dict:
    """Pad per-(truck, path) arc ids, entry offsets, and totals into numpy tables.

    Keeps only deadline-feasible candidate paths; the shortest path of every truck
    is feasible by instance construction, so each truck retains at least one path.
    """
    arc_list = sorted(inst.graph.edges)
    arc_index = {a: i for i, a in enumerate(arc_list)}
    arc_dist = np.array([inst.graph.edges[a]["dist"] for a in arc_list])
    K = len(inst.trucks)
    feas: list[list[int]] = []
    for k in range(K):
        keep = [p for p in range(len(inst.paths[k]))
                if inst.trucks[k].earliest_dep
                + sum(inst.arc_time(a) for a in inst.path_arcs(k, p))
                <= inst.trucks[k].latest_arr]
        feas.append(keep)
    n_path = max(len(f) for f in feas)
    max_len = max(len(inst.path_arcs(k, p)) for k in range(K) for p in feas[k])
    arc_ids = -np.ones((K, n_path, max_len), dtype=np.int64)
    offsets = np.zeros((K, n_path, max_len))
    tot_time = np.zeros((K, n_path))
    tot_dist = np.zeros((K, n_path))
    dep_lo = np.zeros((K, n_path))
    dep_span = np.zeros((K, n_path))
    for k in range(K):
        for j, p in enumerate(feas[k]):
            t = 0.0
            for s, a in enumerate(inst.path_arcs(k, p)):
                arc_ids[k, j, s] = arc_index[a]
                offsets[k, j, s] = t
                t += inst.arc_time(a)
            tot_time[k, j] = t
            tot_dist[k, j] = sum(inst.arc_dist(a) for a in inst.path_arcs(k, p))
            dep_lo[k, j] = inst.trucks[k].earliest_dep
            dep_span[k, j] = max(0.0, inst.trucks[k].latest_arr - t - dep_lo[k, j])
    return {"arc_ids": arc_ids, "offsets": offsets, "tot_time": tot_time,
            "tot_dist": tot_dist, "dep_lo": dep_lo, "dep_span": dep_span,
            "n_paths": np.array([len(f) for f in feas]), "arc_dist": arc_dist,
            "feas": feas}


def fitness(pop_path: np.ndarray, pop_key: np.ndarray, tab: dict,
            inst: PlatooningInstance) -> np.ndarray:
    """Vectorized population cost: base fuel + duration cost - platoon savings."""
    P, K = pop_path.shape
    kk = np.arange(K)
    dep = tab["dep_lo"][kk, pop_path] + pop_key * tab["dep_span"][kk, pop_path]
    base = (inst.fuel_per_km * tab["tot_dist"][kk, pop_path]).sum(axis=1)
    dur = tab["tot_time"][kk, pop_path].sum(axis=1)
    arcs = tab["arc_ids"][kk, pop_path]                       # (P, K, L)
    times = dep[:, :, None] + tab["offsets"][kk, pop_path]    # (P, K, L)
    valid = arcs >= 0
    ind = np.broadcast_to(np.arange(P)[:, None, None], arcs.shape)[valid]
    a_f, t_f = arcs[valid], times[valid]
    order = np.lexsort((t_f, a_f, ind))                       # by individual, arc, time
    ind, a_f, t_f = ind[order], a_f[order], t_f[order]
    brk = np.ones(a_f.size, dtype=bool)
    brk[1:] = ((ind[1:] != ind[:-1]) | (a_f[1:] != a_f[:-1])
               | (t_f[1:] - t_f[:-1] > inst.sync_tol))
    gid = np.cumsum(brk) - 1
    sizes = np.bincount(gid)
    followers = sizes - np.ceil(sizes / inst.max_platoon).astype(np.int64)
    starts = np.flatnonzero(brk)
    sav = followers * inst.eta_follow * inst.fuel_per_km * tab["arc_dist"][a_f[starts]]
    savings = np.bincount(ind[starts], weights=sav, minlength=P)
    return base + inst.wait_cost * dur - savings


def decode_individual(path_idx: np.ndarray, key: np.ndarray, tab: dict,
                      inst: PlatooningInstance) -> dict:
    """Explicit schedule + leader-anchored platoon pairs for one individual."""
    route: dict[int, int] = {}
    entry: dict[tuple[int, tuple[int, int]], float] = {}
    by_arc: dict[tuple[int, int], list[tuple[float, int]]] = {}
    for k in range(len(inst.trucks)):
        j = int(path_idx[k])
        route[k] = tab["feas"][k][j]
        t_cur = float(tab["dep_lo"][k, j] + key[k] * tab["dep_span"][k, j])
        for a in inst.path_arcs(k, route[k]):
            entry[(k, a)] = t_cur
            by_arc.setdefault(a, []).append((t_cur, k))
            t_cur += inst.arc_time(a)
    pairs: list[tuple[int, int, tuple[int, int]]] = []
    for a, members in by_arc.items():
        members.sort()
        i = 0
        while i < len(members):
            t0, leader = members[i]
            j2 = i + 1
            while (j2 < len(members) and members[j2][0] - t0 <= inst.sync_tol
                   and j2 - i < inst.max_platoon):
                pairs.append((leader, members[j2][1], a))
                j2 += 1
            i = j2
    return {"route": route, "entry": entry, "pairs": pairs}


inst = generate_instance(n_rows=4, n_cols=4, n_trucks=6, n_paths=3,
                         detour_cap=1.25, seed=42)
inst.sync_tol = 0.1
tab = precompute_tables(inst)
rng = np.random.default_rng(0)
pop_path = rng.integers(0, tab["n_paths"], size=(32, len(inst.trucks)))
pop_key = rng.random((32, len(inst.trucks)))
print(round(float(fitness(pop_path, pop_key, tab, inst).min()), 1))
# Expected: one float — the best cost in a random population of 32; identical on every
# run with these seeds, and above (worse than) the MIP optimum on the same instance.

The GA engine itself is a standard elitist generational loop over the two gene matrices. Path and departure genes of the same truck travel together through crossover, since a departure key is only meaningful relative to its path's time window:

import numpy as np

# Uses fitness, precompute_tables, decode_individual, generate_instance from above.


def run_ga(inst: PlatooningInstance, tab: dict, pop_size: int = 200,
           generations: int = 300, pc: float = 0.9, pm: float = 0.15,
           seed: int = 0) -> tuple[np.ndarray, np.ndarray, float]:
    """Elitist GA over (path index, departure key) genomes.

    See genetic-algorithms for operator and parameter theory; everything
    problem-specific lives in fitness() and decode_individual().
    """
    rng = np.random.default_rng(seed)
    K = len(inst.trucks)
    high = tab["n_paths"]                       # per-truck feasible path counts
    pop_path = rng.integers(0, high, size=(pop_size, K))
    pop_key = rng.random((pop_size, K))
    fit = fitness(pop_path, pop_key, tab, inst)
    for _ in range(generations):
        # tournament selection, size 2 (vectorized)
        cand = rng.integers(0, pop_size, size=(pop_size, 2))
        winners = np.where(fit[cand[:, 0]] <= fit[cand[:, 1]], cand[:, 0], cand[:, 1])
        par_path, par_key = pop_path[winners], pop_key[winners]
        # uniform crossover on truck columns; path + key genes travel together
        mate = rng.permutation(pop_size)
        mask = (rng.random((pop_size, K)) < 0.5) & (rng.random((pop_size, 1)) < pc)
        child_path = np.where(mask, par_path[mate], par_path)
        child_key = np.where(mask, par_key[mate], par_key)
        # mutation: resample path gene; Gaussian creep on the departure key
        mut = rng.random((pop_size, K)) < pm
        child_path = np.where(mut, rng.integers(0, high, size=(pop_size, K)), child_path)
        creep = rng.random((pop_size, K)) < pm
        child_key = np.clip(
            np.where(creep, child_key + rng.normal(0.0, 0.15, (pop_size, K)), child_key),
            0.0, 1.0)
        child_fit = fitness(child_path, child_key, tab, inst)
        # elitism: best parent replaces the worst child
        best, worst = int(np.argmin(fit)), int(np.argmax(child_fit))
        child_path[worst], child_key[worst] = pop_path[best], pop_key[best]
        child_fit[worst] = fit[best]
        pop_path, pop_key, fit = child_path, child_key, child_fit
    best = int(np.argmin(fit))
    return pop_path[best], pop_key[best], float(fit[best])


inst = generate_instance(n_rows=6, n_cols=6, n_trucks=20, n_paths=3,
                         detour_cap=1.25, seed=7)
inst.sync_tol = 0.1            # bounded en-route speed adaptation closes <= 6 min gaps
tab = precompute_tables(inst)
best_path, best_key, best_fit = run_ga(inst, tab, seed=1)
ga_sol = decode_individual(best_path, best_key, tab, inst)
print(round(best_fit, 1), len(ga_sol["pairs"]))
# Expected: fitness improves over the random-population best from the previous block's
# setup; pairs concentrate on central grid corridors. Deterministic for fixed seeds.

Parameter guidance for this specific GA:

Independent Solution Validation

Never trust the producing code. One validator checks any solution — MIP or decoded GA — against the raw instance: routes are candidate paths, schedules propagate with non-negative waiting, time windows hold, every platoon pair is physically consistent (shared arc, synchronized entry, one leader per follower, leaders never follow, size cap), and the objective recomputed from scratch matches.

import math

# Uses PlatooningInstance from the instance-generator block; validates solution dicts
# in the {"route", "entry", "pairs"} format produced by both the MIP and the GA decode.


def validate_solution(inst: PlatooningInstance, sol: dict) -> tuple[bool, list[str], float]:
    """Independent feasibility check + objective recomputation from raw instance data."""
    errors: list[str] = []
    eps = 1e-6
    route, entry, pairs = sol["route"], sol["entry"], sol["pairs"]
    for k in range(len(inst.trucks)):
        if route[k] not in range(len(inst.paths[k])):
            errors.append(f"truck {k}: route index {route[k]} is not a candidate path")
            continue
        arcs = inst.path_arcs(k, route[k])
        if entry[(k, arcs[0])] < inst.trucks[k].earliest_dep - eps:
            errors.append(f"truck {k}: departs before earliest departure")
        for a, b in zip(arcs[:-1], arcs[1:]):
            if entry[(k, b)] < entry[(k, a)] + inst.arc_time(a) - eps:
                errors.append(f"truck {k}: negative waiting before arc {b}")
        last = arcs[-1]
        if entry[(k, last)] + inst.arc_time(last) > inst.trucks[k].latest_arr + eps:
            errors.append(f"truck {k}: misses arrival deadline")
    leaders_of: dict[tuple, list[int]] = {}
    follower_count: dict[tuple, int] = {}
    for l, k, a in pairs:
        for m in (l, k):
            if a not in inst.path_arcs(m, route[m]):
                errors.append(f"pair ({l},{k},{a}): truck {m} does not traverse the arc")
        gap = abs(entry.get((k, a), math.inf) - entry.get((l, a), -math.inf))
        if gap > inst.sync_tol + eps:
            errors.append(f"pair ({l},{k},{a}): entry gap {gap:.4f} beyond tolerance")
        leaders_of.setdefault((k, a), []).append(l)
        follower_count[(l, a)] = follower_count.get((l, a), 0) + 1
    for (k, a), ls in leaders_of.items():
        if len(ls) > 1:
            errors.append(f"truck {k} has {len(ls)} leaders on arc {a}")
        for l in ls:
            if (l, a) in leaders_of:
                errors.append(f"leader {l} is itself a follower on arc {a}")
    for (l, a), cnt in follower_count.items():
        if cnt > inst.max_platoon - 1:
            errors.append(f"leader {l} exceeds the platoon size cap on arc {a}")
    base = sum(inst.fuel_per_km * inst.arc_dist(a)
               for k in range(len(inst.trucks)) for a in inst.path_arcs(k, route[k]))
    savings = sum(inst.eta_follow * inst.fuel_per_km * inst.arc_dist(a)
                  for (_, _, a) in pairs)
    duration = 0.0
    for k in range(len(inst.trucks)):
        arcs = inst.path_arcs(k, route[k])
        duration += (entry[(k, arcs[-1])] + inst.arc_time(arcs[-1])
                     - entry[(k, arcs[0])])
    objective = base - savings + inst.wait_cost * duration
    return (not errors), errors, objective


inst = generate_instance(n_rows=4, n_cols=4, n_trucks=6, n_paths=3,
                         detour_cap=1.25, seed=42)
mip_sol = solve_platooning(inst)
ok, problems, obj = validate_solution(inst, mip_sol)
print(ok, len(problems), round(obj - mip_sol["objective"], 6))
# Expected: True 0 0.0 (up to ~1e-6 solver tolerance). Any nonzero difference or any
# listed problem means a modeling or extraction bug, not a bad instance.

Run the validator on every solution that leaves the pipeline. For GA solutions it also serves as the semantic bridge: the decoded leader-anchored pairs are exactly what it checks, so a validated GA objective is comparable with a validated MIP objective under the same sync_tol. Cross-validate the two methods on small instances: the MIP optimum is a regression anchor the GA must approach.

Advanced Techniques

Tightening and pruning the synchronization model

The horizon big-M is correct but weak. Compute, per truck and arc, the earliest possible entry $est_{ka} = e_k + \mathrm{SP}\tau(o_k, \mathrm{tail}(a))$ and the latest useful entry $lst{ka} = \ell_k - \tau_a - \mathrm{SP}\tau(\mathrm{head}(a), s_k)$, where $\mathrm{SP}\tau$ is the shortest drive time (a handful of Dijkstra runs — see network-flow-optimization). Then: (1) prune every pairing variable $w_{lka}$ whose entry intervals $[est, lst]$ do not overlap within $\delta$ — on instances with spread-out departures this removes the majority of pairing variables before the model exists; (2) tighten the remaining sync big-Ms to $M_{lka} = \max(lst_{ka} - est_{la},, lst_{la} - est_{ka})$; (3) tighten the time-variable bounds themselves to $[est_{ka}, lst_{ka}]$. Gurobi indicator constraints (addGenConstrIndicator) are the robust alternative when horizons are long and big-Ms stay large despite tightening — slightly slower per node, immune to big-M numerics.

Leader savings and platoon-size-dependent factors

To give the leader a saving fraction $\eta_L$ (aerodynamic measurements support a few percent), add binaries $g_{la}$ ("$l$ leads at least one truck on $a$") with $g_{la} \le \sum_k w_{lka}$ and the term $-\eta_L, c, d_a, g_{la}$ in the objective; the negative coefficient pushes $g$ up, so only the upper link is needed. Size-dependent follower savings (middle trucks save more than the tail truck) need position modeling — usually not worth it at planning level; calibrate a flat per-follower $\eta$ from the aerodynamic curve instead, and state that simplification in the report.

Maximal shared segments instead of single arcs

Platooning decisions per arc multiply variables and let the solver form physically silly one-arc platoons. Precompute, per truck pair, the maximal contiguous shared subpaths of their candidate paths and define one pairing variable per (pair, segment). Synchronization is then one constraint per segment (entry into its first arc), savings aggregate the segment length, and a minimum-segment-length filter (say ≥ 30 km) encodes that forming a platoon has a fixed coordination overhead. This typically shrinks the pairing model by 3–10x and improves solution realism at the same time.

Matheuristic: routes and pairs from the GA, exact timing from an LP

The GA's tolerance-based pairs are promises, not schedules. With routes and pairs fixed, finding exactly synchronized entry times is a pure LP — small, always fast, and it certifies the promise or exposes it as over-constrained (then drop the cheapest-saving pair and retry). This route-first, schedule-second split is the workhorse decomposition of the platooning literature (see the Bhoopalam et al. 2018 review) and an instance of the fix-and-optimize pattern in matheuristics:

import gurobipy as gp
from gurobipy import GRB

# Uses PlatooningInstance from the instance-generator block.


def polish_schedule(inst: PlatooningInstance, route: dict,
      

…(truncated)