Large neighbourhood search (ruin and recreate)

Deep dive

How it works

Parts:

• destroy operators — random, worst (largest detour saving, randomised), related (Shaw: items similar in place and demand)
• repair operators — greedy cheapest insertion; regret-k (insert first the item that loses most if not placed now)
• selector — roulette over weights, re-weighted at segment ends
• acceptance — hill climbing, simulated annealing or record-to-record

1. Pick a destroy/repair pair with probability proportional to its weight.
2. Draw q (here 3 to 10 of 60), copy the current solution, remove q items.
3. Re-insert them; a new route opens if nothing fits.
4. Score the pair: new best 33, better than current 9, accepted worse 3 (my values).
5. Accept if better, else with probability exp(-delta/T); cool T geometrically.
6. Every 100 iterations: w = (1-r)*w + r*score/uses, r = 0.2; reset scores.

 current --destroy(q)--> partial --repair--> candidate --accept?--> current
    ^                                           |                      |
    +---- weights <-- score(pair) <-------------+       best <---------+

Acceptance. Hill climbing stalls early; annealing accepts worse candidates with falling probability; record-to-record accepts anything within a fraction of the best. Destruction degree is drawn from a range at each iteration to mix small and large steps.

Inside exact solvers. RINS (Danna, Rothberg, Le Pape 2005) fixes variables where the incumbent and the LP relaxation agree and solves the sub-MIP. CP-SAT runs use_rins_lns plus neighbourhood generators (random variables or constraints, variable/constraint graph, scheduling windows, routing paths) in parallel workers; difficulty adapts (lns_initial_difficulty, default 0.5).

Code

Capacitated routing, 60 customers, capacity 40. Destroy: random or worst; repair: cheapest insertion, as-removed or heaviest-first order; four pairs.

struct Rng(u64);
impl Rng {
    fn next(&mut self) -> u64 { // xorshift64*
        self.0 ^= self.0 >> 12; self.0 ^= self.0 << 25; self.0 ^= self.0 >> 27;
        self.0.wrapping_mul(0x2545F4914F6CDD1D)
    }
    fn f(&mut self) -> f64 { (self.next() >> 11) as f64 / (1u64 << 53) as f64 }
    fn below(&mut self, n: usize) -> usize { (self.f() * n as f64) as usize }
}
struct Inst { xy: Vec<(f64, f64)>, dem: Vec<u32>, cap: u32 } // node 0 = depot
impl Inst {
    fn d(&self, a: usize, b: usize) -> f64 { (self.xy[a].0 - self.xy[b].0).hypot(self.xy[a].1 - self.xy[b].1) }
    fn route_cost(&self, r: &[usize]) -> f64 {
        let p = std::iter::once(&0).chain(r).zip(r.iter().chain(&[0]));
        p.map(|(&a, &b)| self.d(a, b)).sum()
    }
    fn cost(&self, s: &[Vec<usize>]) -> f64 { s.iter().map(|r| self.route_cost(r)).sum() }
    fn load(&self, r: &[usize]) -> u32 { r.iter().map(|&n| self.dem[n]).sum() }
}
fn nearest_neighbour(ins: &Inst) -> Vec<Vec<usize>> {
    let mut left: Vec<usize> = (1..ins.xy.len()).collect();
    let mut sol = vec![];
    while !left.is_empty() {
        let (mut route, mut load, mut at) = (vec![], 0, 0usize);
        while let Some((i, &n)) = left.iter().enumerate()
            .filter(|(_, &n)| load + ins.dem[n] <= ins.cap)
            .min_by(|a, b| ins.d(at, *a.1).total_cmp(&ins.d(at, *b.1))) {
            load += ins.dem[n]; route.push(n); at = n; left.swap_remove(i);
        }
        sol.push(route);
    }
    sol
}
// Destroy: returns the removed customers. kind 0 = random, 1 = worst (largest detour saving, randomised).
fn destroy(ins: &Inst, s: &mut Vec<Vec<usize>>, q: usize, kind: usize, rng: &mut Rng) -> Vec<usize> {
    let mut out = vec![];
    for _ in 0..q {
        let mut cand = vec![]; // (saving, route, pos)
        for (ri, r) in s.iter().enumerate() {
            for (pi, &n) in r.iter().enumerate() {
                let p = if pi == 0 { 0 } else { r[pi - 1] };
                let nx = if pi + 1 == r.len() { 0 } else { r[pi + 1] };
                cand.push((ins.d(p, n) + ins.d(n, nx) - ins.d(p, nx), ri, pi));
            }
        }
        if kind == 1 { cand.sort_by(|a, b| b.0.total_cmp(&a.0)); }
        let k = if kind == 0 { rng.below(cand.len()) } else { (rng.f().powi(3) * cand.len() as f64) as usize };
        out.push(s[cand[k].1].remove(cand[k].2));
    }
    s.retain(|r| !r.is_empty());
    out
}
// Repair: cheapest insertion; kind 1 inserts the heaviest demand first.
fn repair(ins: &Inst, s: &mut Vec<Vec<usize>>, mut rem: Vec<usize>, kind: usize) {
    if kind == 1 { rem.sort_by_key(|&n| std::cmp::Reverse(ins.dem[n])); }
    for n in rem {
        let mut best = (f64::MAX, usize::MAX, 0);
        for (ri, r) in s.iter().enumerate() {
            if ins.load(r) + ins.dem[n] > ins.cap { continue; }
            for pi in 0..=r.len() {
                let p = if pi == 0 { 0 } else { r[pi - 1] };
                let nx = if pi == r.len() { 0 } else { r[pi] };
                let delta = ins.d(p, n) + ins.d(n, nx) - ins.d(p, nx);
                if delta < best.0 { best = (delta, ri, pi); }
            }
        }
        if best.1 == usize::MAX { s.push(vec![n]); } else { s[best.1].insert(best.2, n); }
    }
}
fn run(ins: &Inst, seed: u64, iters: usize) -> (f64, Vec<Vec<usize>>, Vec<f64>) {
    let mut rng = Rng(seed);
    let mut cur = nearest_neighbour(ins);
    let c0 = ins.cost(&cur);
    let (mut cc, mut bc, mut best) = (c0, c0, cur.clone());
    let mut t = 0.05 * c0 / std::f64::consts::LN_2; // 5 % worse accepted with p = 0.5 at the start
    let cool = (0.0005f64).powf(1.0 / iters as f64); // ends at 0.05 % of T0
    let (mut w, mut score, mut used) = ([1.0f64; 4], [0.0f64; 4], [0u32; 4]); // pair = 2*destroy + repair
    for it in 0..iters {
        let tot: f64 = w.iter().sum();
        let (mut x, mut pair) = (rng.f() * tot, 3);
        for (i, wi) in w.iter().enumerate() { if x < *wi { pair = i; break; } x -= wi; }
        let q = 3 + rng.below(8); // destroy 3..10 of 60 customers
        let mut cand = cur.clone();
        let rem = destroy(ins, &mut cand, q, pair / 2, &mut rng);
        repair(ins, &mut cand, rem, pair % 2);
        let c = ins.cost(&cand);
        used[pair] += 1;
        if c < bc - 1e-9 { best = cand.clone(); bc = c; score[pair] += 33.0; }
        else if c < cc - 1e-9 { score[pair] += 9.0; }
        if c < cc || rng.f() < ((cc - c) / t).exp() {
            if c >= cc { score[pair] += 3.0; }
            cur = cand; cc = c;
        }
        t *= cool;
        if (it + 1) % 100 == 0 { // segment end: w = (1-r) w + r score/uses, r = 0.2
            for i in 0..4 {
                if used[i] > 0 { w[i] = (0.8 * w[i] + 0.2 * score[i] / used[i] as f64).max(0.05); }
                score[i] = 0.0; used[i] = 0;
            }
        }
    }
    (bc, best, w.to_vec())
}
fn main() {
    let mut rng = Rng(42);
    let (mut xy, mut dem) = (vec![(50.0, 50.0)], vec![0]);
    for _ in 0..60 { xy.push((rng.f() * 100.0, rng.f() * 100.0)); dem.push(1 + rng.below(9) as u32); }
    let ins = Inst { xy, dem, cap: 40 };
    let start = ins.cost(&nearest_neighbour(&ins));
    let t0 = std::time::Instant::now();
    let (c1, sol, w) = run(&ins, 7, 20_000);
    let el = t0.elapsed();
    let (c2, _, _) = run(&ins, 7, 20_000);
    let mut seen: Vec<usize> = sol.iter().flatten().copied().collect();
    seen.sort();
    assert!(seen == (1..61).collect::<Vec<_>>() && sol.iter().all(|r| ins.load(r) <= ins.cap));
    assert!(c1 < 0.9 * start && c1 == c2);
    assert!(w.iter().any(|x| (x - 1.0).abs() > 0.05));
    println!("start {start:.1} final {c1:.1} routes {} iters 20000 time {el:?} weights {w:.2?}", sol.len());
}

Run: start 1425.8, final 993.9, 9 routes, 20,000 iterations, about 0.13 s (release build, one seeded instance). Weights end at [1.89, 2.14, 2.02, 1.97] (start 1.0), no pair dominates. Asserts: feasibility, 10 % improvement, same cost twice, changed weights.

Complexity and constants

Per iteration O(q*n): removal scans all n positions per item; insertion tries every position of every feasible route; full cost recomputation is O(n). Measured: about 6 us per iteration at n = 60. Delta evaluation removes the dependence on n.

Efficiency

• Cost: about 6 us per iteration at n = 60 (measured, one thread, compute bound).
• Levers:
  • delta cost instead of full recomputation — removes the O(n) pass — always
  • copy only touched routes or undo by journal (journaled-state-revert) — cuts the clone — large n
  • neighbour lists in insertion (k nearest positions only) — scan n to k per item (standard in routing codes) — n over a few hundred
  • parallel LNS workers sharing the best solution, as CP-SAT does — spare cores
• Naive → standard → state of the art: full exact re-solve · exponential; 2-opt/or-opt only · local optima; ALNS with annealing and delta evaluation · usual routing baseline.

Implementation notes

Cost drivers are cloning the solution and recomputing cost; use flat Vecs and reuse buffers. One seeded PRNG (PCG or xoshiro) passed everywhere gives reproducible runs. Set T0 relatively (here 5 % worse accepted with p = 0.5) so it transfers across scales.

Pitfalls

• Weights collapse onto one pair; floor them (0.05 here) and keep r small.
• Pure greedy repair repeats itself; randomise removal or insertion order.
• Compare with a tuned 2-opt/or-opt at equal time, not only with the start.

Variants and alternatives

In production

• OR-Tools routing: LNS operators (path, TSP, inactive, cheapest-insertion variants), switchable in the search parameters.
• OR-Tools CP-SAT: LNS workers beside complete search.
• MIP solvers: RINS-type neighbourhoods as primal heuristics (MIP).
• Python ALNS package: roulette-wheel selection; hill-climbing, record-to-record and annealing acceptance.

Implementations

Sources

• Shaw, Using Constraint Programming and Local Search Methods to Solve Vehicle Routing Problems, CP 1998, LNCS 1520, pp. 417-431. paper
• Ropke, Pisinger, An Adaptive Large Neighborhood Search Heuristic for the Pickup and Delivery Problem with Time Windows, Transportation Science 40(4), 2006, pp. 455-472. paper
• Danna, Rothberg, Le Pape, Exploring relaxation induced neighborhoods to improve MIP solutions, Mathematical Programming 102, 2005, pp. 71-90. paper
• OR-Tools source: ortools/sat/cp_model_lns.h, sat_parameters.proto. code