routik_solver/

anneal.rs

//! Simulated-annealing metaheuristic over the local-search neighbourhoods.
//!
//! Starting from the Clarke-Wright construction, the search repeatedly proposes
//! a random move, accepts every improving move and — with a temperature-decaying
//! probability — some worsening ones, so it can climb out of local optima. Only
//! capacity- and time-window-feasible moves are ever applied, so the returned
//! solution is feasible by construction. Cooling is geometric; the run stops on
//! an iteration and/or wall-clock budget.
//!
//! Determinism: all randomness comes from a seeded [`ChaCha8Rng`], so the same
//! `(problem, matrix, config)` always yields the same solution.

use crate::clarke_wright::{SolveError, clarke_wright};
use crate::matrix::CostMatrix;
use crate::model::{Problem, Solution};
use crate::objective::Objective;
use crate::search::State;
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use std::time::{Duration, Instant};

/// How long the search may run: capped by iterations, wall-clock, or both
/// (whichever bound is hit first). At least one bound should be set.
#[derive(Debug, Clone, Copy)]
pub struct Budget {
    pub max_iterations: Option<u64>,
    pub max_time: Option<Duration>,
}

impl Budget {
    /// Stop after `n` iterations.
    #[must_use]
    pub fn iterations(n: u64) -> Self {
        Self {
            max_iterations: Some(n),
            max_time: None,
        }
    }

    /// Stop after `d` of wall-clock time.
    #[must_use]
    pub fn time(d: Duration) -> Self {
        Self {
            max_iterations: None,
            max_time: Some(d),
        }
    }
}

impl Budget {
    /// Divide the budget into `n` equal slices (one per restart). At least one
    /// iteration / a non-zero duration per slice.
    fn split(&self, n: usize) -> Self {
        let n = n.max(1) as u64;
        Self {
            max_iterations: self.max_iterations.map(|m| (m / n).max(1)),
            max_time: self.max_time.map(|d| d / n as u32),
        }
    }
}

impl Default for Budget {
    fn default() -> Self {
        Self::iterations(500_000)
    }
}

/// Knobs for [`solve`]. [`Default`] minimises distance with a 100k-iteration
/// budget and an auto-calibrated temperature schedule.
#[derive(Debug, Clone, Copy)]
pub struct SolverConfig {
    /// Seed for the internal RNG — fixes the whole run.
    pub seed: u64,
    /// What to minimise.
    pub objective: Objective,
    /// Total run length, split evenly across [`Self::restarts`].
    pub budget: Budget,
    /// Number of independent multi-start cool-downs; the best is returned.
    /// Diversifies the search so it doesn't get trapped in one local optimum.
    pub restarts: usize,
    /// Starting temperature; `None` auto-calibrates from sampled move deltas.
    pub initial_temperature: Option<f64>,
    /// Ending temperature; `None` uses `initial / 1000`.
    pub final_temperature: Option<f64>,
}

impl Default for SolverConfig {
    fn default() -> Self {
        Self {
            seed: 0x5EED,
            objective: Objective::default(),
            budget: Budget::default(),
            restarts: 10,
            initial_temperature: None,
            final_temperature: None,
        }
    }
}

/// How often to poll the wall clock (polling every iteration is wasteful).
const TIME_CHECK_STRIDE: u64 = 2048;

/// Rolled-up run statistics from [`solve_with_stats`], for surfacing solver
/// metadata (e.g. in an API response). Carries no solution data.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct SolveStats {
    /// Total simulated-annealing iterations summed across every restart.
    pub iterations: u64,
    /// Number of multi-start restarts actually performed.
    pub restarts: usize,
}

/// Optimise a problem with simulated annealing, returning a feasible solution.
///
/// The construction baseline is Clarke-Wright; the returned solution is never
/// worse than it under the configured [`Objective`].
///
/// # Errors
/// Propagates [`SolveError`] from construction (no vehicles, an unroutable
/// stop, or more routes than vehicles).
pub fn solve<M: CostMatrix>(
    problem: &Problem,
    matrix: &M,
    config: &SolverConfig,
) -> Result<Solution, SolveError> {
    solve_with_stats(problem, matrix, config).map(|(solution, _)| solution)
}

/// Like [`solve`], but also returns [`SolveStats`] about the run.
///
/// # Errors
/// Propagates [`SolveError`] from construction (no vehicles, an unroutable
/// stop, or more routes than vehicles).
pub fn solve_with_stats<M: CostMatrix>(
    problem: &Problem,
    matrix: &M,
    config: &SolverConfig,
) -> Result<(Solution, SolveStats), SolveError> {
    let initial = clarke_wright(problem, matrix)?;

    // Local search relocates/swaps stops *between* routes but never drops or adds
    // one, so the set of unassigned stops is invariant across the whole search —
    // capture it now and stamp it back onto the optimised solution at the end.
    let unassigned = initial.unassigned.clone();

    // Nothing to move around with 0 or 1 stop.
    if problem.stops.len() <= 1 {
        return Ok((
            initial,
            SolveStats {
                iterations: 0,
                restarts: 0,
            },
        ));
    }

    let capacity = problem
        .vehicles
        .iter()
        .map(|v| v.capacity)
        .min()
        .unwrap_or(0);
    let obj = &config.objective;
    let restarts = config.restarts.max(1);
    let per_start = config.budget.split(restarts);

    // Multi-start: each restart cools fully from a fresh construction under its
    // own RNG stream, so the runs both intensify (a full cool each) and
    // diversify (different trajectories). Keeping the best is far more robust
    // than a single cool, which can get trapped on clustered instances.
    let mut best_routes: Option<Vec<Vec<usize>>> = None;
    let mut best_cost = f64::INFINITY;
    let mut iterations: u64 = 0;

    for s in 0..restarts {
        let seed = config
            .seed
            .wrapping_add((s as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15));
        let mut rng = ChaCha8Rng::seed_from_u64(seed);
        let mut state = State::from_solution(problem, matrix, capacity, &initial);
        iterations += anneal_once(&mut state, &mut rng, config, &per_start, obj);
        let cost = state.cost(obj);
        if cost < best_cost {
            best_cost = cost;
            best_routes = Some(state.snapshot());
        }
    }

    let mut state = State::from_solution(problem, matrix, capacity, &initial);
    if let Some(routes) = best_routes {
        state.restore(routes);
    }
    let stats = SolveStats {
        iterations,
        restarts,
    };
    let mut solution = state.to_solution();
    solution.feasible = solution.feasible && unassigned.is_empty();
    solution.unassigned = unassigned;
    Ok((solution, stats))
}

/// One simulated-annealing cool-down. Leaves `state` at the best solution it
/// found (always feasible, since only feasible moves are applied) and returns
/// the number of iterations it ran.
fn anneal_once<M: CostMatrix>(
    state: &mut State<M>,
    rng: &mut ChaCha8Rng,
    config: &SolverConfig,
    budget: &Budget,
    obj: &Objective,
) -> u64 {
    let t0 = config
        .initial_temperature
        .unwrap_or_else(|| auto_initial_temperature(state, obj, rng));
    let t_end = config.final_temperature.unwrap_or(t0 / 1000.0).max(1e-9);
    let alpha = cooling_rate(t0, t_end, budget.max_iterations);

    let mut temperature = t0;
    let mut best_routes = state.snapshot();
    let mut best_cost = state.cost(obj);

    let start = Instant::now();
    let mut iter: u64 = 0;
    loop {
        if budget.max_iterations.is_some_and(|m| iter >= m) {
            break;
        }
        if iter.is_multiple_of(TIME_CHECK_STRIDE)
            && budget
                .max_time
                .is_some_and(|limit| start.elapsed() >= limit)
        {
            break;
        }

        if let Some(mv) = state.random_move(rng) {
            let delta = state.cost_delta(&mv, obj);
            let accept = delta <= 0.0 || rng.random::<f64>() < (-delta / temperature).exp();
            if accept && state.try_apply(&mv) {
                let cost = state.cost(obj);
                if cost < best_cost {
                    best_cost = cost;
                    best_routes = state.snapshot();
                }
            }
        }

        temperature = (temperature * alpha).max(t_end);
        iter += 1;
    }

    state.restore(best_routes);
    iter
}

/// Calibrate the starting temperature so a move that worsens the objective by
/// the average sampled magnitude is accepted with probability ~0.5.
fn auto_initial_temperature<M: CostMatrix>(
    state: &State<M>,
    obj: &Objective,
    rng: &mut ChaCha8Rng,
) -> f64 {
    const SAMPLES: usize = 200;
    let mut sum = 0.0;
    let mut count = 0u32;
    for _ in 0..SAMPLES {
        if let Some(mv) = state.random_move(rng) {
            sum += state.cost_delta(&mv, obj).abs();
            count += 1;
        }
    }
    if count == 0 {
        return 1.0;
    }
    let mean = sum / f64::from(count);
    (mean / std::f64::consts::LN_2).max(1e-6)
}

/// Per-iteration geometric cooling factor taking `t0` down to `t_end` over the
/// iteration budget; a slow default when only a time budget is set.
fn cooling_rate(t0: f64, t_end: f64, max_iterations: Option<u64>) -> f64 {
    match max_iterations {
        Some(iters) if iters > 1 => (t_end / t0).powf(1.0 / iters as f64),
        _ => 0.99997,
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::matrix::EuclideanMatrix;
    use crate::model::{Coord, Stop, StopId, TimeWindow, Vehicle, VehicleId};

    /// Four spatial clusters of stops with generous windows: there is a clearly
    /// better routing than Clarke-Wright tends to build, so SA should improve.
    fn clustered_problem() -> Problem {
        let centers = [(0.0, 0.0), (100.0, 0.0), (0.0, 100.0), (100.0, 100.0)];
        let mut stops = Vec::new();
        let mut id = 1u32;
        for (cx, cy) in centers {
            for k in 0..5 {
                let off = f64::from(k) * 2.0;
                stops.push(Stop {
                    id: StopId(id),
                    coord: Coord::new(cy + off, cx + off), // lat, lon
                    demand: 5,
                    time_window: Some(TimeWindow {
                        start: 0.0,
                        end: 100_000.0,
                    }),
                    service_time: 1.0,
                });
                id += 1;
            }
        }
        let vehicles = (1..=8u32)
            .map(|i| Vehicle {
                id: VehicleId(i),
                capacity: 30,
            })
            .collect();
        Problem {
            depot: Coord::new(50.0, 50.0),
            stops,
            vehicles,
            depot_window: Some(TimeWindow {
                start: 0.0,
                end: 1_000_000.0,
            }),
        }
    }

    /// A single-vehicle, no-window instance over scattered points — a TSP where
    /// the Clarke-Wright tour has crossings that 2-opt/relocate can iron out, so
    /// SA reliably improves on the construction.
    fn scattered_problem() -> Problem {
        let pts = [
            (20.0, 20.0),
            (80.0, 25.0),
            (15.0, 70.0),
            (70.0, 80.0),
            (45.0, 10.0),
            (90.0, 55.0),
            (30.0, 45.0),
            (60.0, 30.0),
            (10.0, 40.0),
            (75.0, 65.0),
            (40.0, 85.0),
            (55.0, 60.0),
        ];
        let stops = pts
            .iter()
            .enumerate()
            .map(|(i, &(x, y))| Stop {
                id: StopId((i + 1) as u32),
                coord: Coord::new(y, x),
                demand: 1,
                time_window: None,
                service_time: 0.0,
            })
            .collect();
        Problem {
            depot: Coord::new(50.0, 50.0),
            stops,
            vehicles: vec![Vehicle {
                id: VehicleId(1),
                capacity: 1000,
            }],
            depot_window: None,
        }
    }

    fn all_feasible(problem: &Problem, sol: &Solution) -> bool {
        let served: usize = sol.routes.iter().map(|r| r.stop_ids.len()).sum();
        sol.feasible && served == problem.stops.len()
    }

    #[test]
    fn same_seed_is_deterministic() {
        let p = clustered_problem();
        let m = EuclideanMatrix::from_problem(&p);
        let cfg = SolverConfig {
            seed: 123,
            budget: Budget::iterations(20_000),
            ..Default::default()
        };
        let a = solve(&p, &m, &cfg).expect("feasible");
        let b = solve(&p, &m, &cfg).expect("feasible");
        assert_eq!(a.total_distance, b.total_distance);
        let route_a: Vec<_> = a.routes.iter().map(|r| &r.stop_ids).collect();
        let route_b: Vec<_> = b.routes.iter().map(|r| &r.stop_ids).collect();
        assert_eq!(route_a, route_b);
    }

    #[test]
    fn different_seeds_can_diverge_but_stay_feasible() {
        let p = clustered_problem();
        let m = EuclideanMatrix::from_problem(&p);
        for seed in [1u64, 2, 99] {
            let cfg = SolverConfig {
                seed,
                budget: Budget::iterations(20_000),
                ..Default::default()
            };
            let sol = solve(&p, &m, &cfg).expect("feasible");
            assert!(all_feasible(&p, &sol), "seed {seed} produced infeasible");
        }
    }

    #[test]
    fn improves_on_clarke_wright() {
        let p = scattered_problem();
        let m = EuclideanMatrix::from_problem(&p);
        let cw = clarke_wright(&p, &m).expect("feasible");
        let cfg = SolverConfig {
            seed: 7,
            budget: Budget::iterations(40_000),
            ..Default::default()
        };
        let sol = solve(&p, &m, &cfg).expect("feasible");
        assert!(all_feasible(&p, &sol));
        assert!(
            sol.total_distance < cw.total_distance,
            "SA {} did not improve on CW {}",
            sol.total_distance,
            cw.total_distance
        );
    }

    #[test]
    fn respects_tight_time_windows() {
        // Stops on a line, each only serviceable in a narrow staggered window,
        // forcing a specific visiting order.
        let stops = (1..=6u32)
            .map(|i| Stop {
                id: StopId(i),
                coord: Coord::new(0.0, f64::from(i) * 10.0),
                demand: 1,
                time_window: Some(TimeWindow {
                    start: f64::from(i) * 10.0 - 2.0,
                    end: f64::from(i) * 10.0 + 2.0,
                }),
                service_time: 0.5,
            })
            .collect();
        let p = Problem {
            depot: Coord::new(0.0, 0.0),
            stops,
            vehicles: (1..=3u32)
                .map(|i| Vehicle {
                    id: VehicleId(i),
                    capacity: 100,
                })
                .collect(),
            depot_window: Some(TimeWindow {
                start: 0.0,
                end: 1000.0,
            }),
        };
        let m = EuclideanMatrix::from_problem(&p);
        let cfg = SolverConfig {
            seed: 5,
            budget: Budget::iterations(20_000),
            ..Default::default()
        };
        let sol = solve(&p, &m, &cfg).expect("feasible");
        assert!(all_feasible(&p, &sol));
    }
}