diff --git a/examples/btsp_kroab.rs b/examples/btsp_kroab.rs
new file mode 100644
index 0000000..3cf8906
--- /dev/null
+++ b/examples/btsp_kroab.rs
@@ -0,0 +1,187 @@
+//! Bi-objective TSP using NSGA-II on the **Kroak/Krobk** instance family
+//! (Lust & Teghem, 2010).
+//!
+//! Two TSP instances over the **same** set of cities define two distance
+//! matrices A and B; the search trades off tour length under A versus tour
+//! length under B. This is the canonical multi-objective combinatorial
+//! benchmark, and it gives a rich Pareto front because the geographies
+//! disagree.
+//!
+//! The instance embedded here is **KroAB-25**: the first 25 cities of
+//! TSPLIB KroA100 and KroB100 (both EUC_2D). Same city *indices*, two
+//! coordinate listings.
+//!
+//! - **Algorithm**: [`Nsga2`].
+//! - **Variation**: [`EdgeRecombinationCrossover`] (the gold-standard TSP
+//!   crossover) piped into [`InversionMutation`] via [`CompositeVariation`].
+//! - **Initializer**: [`ShuffledPermutation`].
+//! - **Encoding**: strict permutation of `[0..25)`.
+//!
+//! Sources:
+//! - TSPLIB95 KroA100 / KroB100 (Reinelt, 1991).
+//! - Lust & Teghem (2010), "The Multiobjective Traveling Salesman Problem:
+//!   A Survey and a New Approach."
+//!
+//! Run with:
+//!
+//! ```bash
+//! cargo run --release --example btsp_kroab
+//! ```
+
+use heuropt::metrics::hypervolume_2d;
+use heuropt::prelude::*;
+
+/// First 25 cities of TSPLIB KroA100 (EUC_2D).
+const KROA_25: [(f64, f64); 25] = [
+    (1380.0, 939.0), (2848.0, 96.0),  (3510.0, 1671.0), (457.0, 334.0),
+    (3888.0, 666.0), (984.0, 965.0),  (2721.0, 1482.0), (1286.0, 525.0),
+    (2716.0, 1432.0),(738.0, 1325.0), (1251.0, 1832.0), (2728.0, 1698.0),
+    (3815.0, 169.0), (3683.0, 1533.0),(1247.0, 1945.0), (123.0, 862.0),
+    (1234.0, 1946.0),(252.0, 1240.0), (611.0, 673.0),   (2576.0, 1676.0),
+    (928.0, 1700.0), (53.0, 857.0),   (1807.0, 1711.0), (274.0, 1420.0),
+    (2574.0, 946.0),
+];
+
+/// First 25 cities of TSPLIB KroB100 (EUC_2D).
+const KROB_25: [(f64, f64); 25] = [
+    (3140.0, 1401.0),(556.0, 1056.0), (3675.0, 1522.0), (1182.0, 1853.0),
+    (3595.0, 1340.0),(1936.0, 953.0), (2722.0, 1311.0), (2839.0, 2055.0),
+    (2253.0, 1242.0),(3142.0, 1591.0),(627.0, 1336.0),  (936.0, 211.0),
+    (4014.0, 471.0), (1376.0, 1452.0),(3289.0, 593.0),  (1453.0, 67.0),
+    (1014.0, 1944.0),(2811.0, 1080.0),(3010.0, 1290.0), (1817.0, 1517.0),
+    (510.0, 458.0),  (1717.0, 1693.0),(1252.0, 1633.0), (1693.0, 1374.0),
+    (539.0, 1378.0),
+];
+
+const N_CITIES: usize = 25;
+
+/// TSPLIB EUC_2D distance: rounded Euclidean.
+fn euc2d_matrix(coords: &[(f64, f64)]) -> Vec<Vec<f64>> {
+    let n = coords.len();
+    let mut d = vec![vec![0.0_f64; n]; n];
+    for i in 0..n {
+        for j in (i + 1)..n {
+            let dx = coords[i].0 - coords[j].0;
+            let dy = coords[i].1 - coords[j].1;
+            let dij = (dx * dx + dy * dy).sqrt().round();
+            d[i][j] = dij;
+            d[j][i] = dij;
+        }
+    }
+    d
+}
+
+struct BTspKroAB {
+    dist_a: Vec<Vec<f64>>,
+    dist_b: Vec<Vec<f64>>,
+}
+
+impl BTspKroAB {
+    fn new() -> Self {
+        Self {
+            dist_a: euc2d_matrix(&KROA_25),
+            dist_b: euc2d_matrix(&KROB_25),
+        }
+    }
+
+    fn tour_length(d: &[Vec<f64>], tour: &[usize]) -> f64 {
+        let n = tour.len();
+        let mut total = 0.0;
+        for i in 0..n {
+            total += d[tour[i]][tour[(i + 1) % n]];
+        }
+        total
+    }
+}
+
+impl Problem for BTspKroAB {
+    type Decision = Vec<usize>;
+
+    fn objectives(&self) -> ObjectiveSpace {
+        ObjectiveSpace::new(vec![
+            Objective::minimize("length_A"),
+            Objective::minimize("length_B"),
+        ])
+    }
+
+    fn evaluate(&self, tour: &Vec<usize>) -> Evaluation {
+        Evaluation::new(vec![
+            Self::tour_length(&self.dist_a, tour),
+            Self::tour_length(&self.dist_b, tour),
+        ])
+    }
+
+    fn decision_schema(&self) -> Vec<DecisionVariable> {
+        (0..N_CITIES)
+            .map(|k| DecisionVariable::new(format!("tour_position_{k}")))
+            .collect()
+    }
+}
+
+fn main() {
+    let problem = BTspKroAB::new();
+
+    let mut optimizer = Nsga2::new(
+        Nsga2Config {
+            population_size: 200,
+            generations: 600,
+            seed: 11,
+        },
+        ShuffledPermutation { n: N_CITIES },
+        CompositeVariation {
+            crossover: EdgeRecombinationCrossover,
+            mutation: InversionMutation,
+        },
+    );
+    let result = optimizer.run(&problem);
+
+    println!("bTSP KroAB-25 — bi-objective TSP via NSGA-II");
+    println!("Source: TSPLIB95 KroA100/KroB100 (first 25 cities), Lust & Teghem bTSP family");
+    println!();
+    println!("Total evaluations: {}", result.evaluations);
+    println!("Pareto-front size: {}", result.pareto_front.len());
+    println!();
+
+    let mut front: Vec<&Candidate<Vec<usize>>> = result.pareto_front.iter().collect();
+    front.sort_by(|a, b| {
+        a.evaluation.objectives[0]
+            .partial_cmp(&b.evaluation.objectives[0])
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+
+    // Print a spread sample of the front (no more than 12 rows).
+    let stride = (front.len() / 12).max(1);
+    println!("    length_A     length_B");
+    let mut printed = 0_usize;
+    for (i, c) in front.iter().enumerate() {
+        if i % stride == 0 || i + 1 == front.len() {
+            let o = &c.evaluation.objectives;
+            println!("    {:>8.0}     {:>8.0}", o[0], o[1]);
+            printed += 1;
+            if printed >= 12 {
+                break;
+            }
+        }
+    }
+    println!();
+
+    if let (Some(corner_a), Some(corner_b)) = (front.first(), front.last()) {
+        println!(
+            "A-corner: A={:.0}, B={:.0}",
+            corner_a.evaluation.objectives[0], corner_a.evaluation.objectives[1]
+        );
+        println!(
+            "B-corner: A={:.0}, B={:.0}",
+            corner_b.evaluation.objectives[0], corner_b.evaluation.objectives[1]
+        );
+    }
+
+    // Hypervolume vs. a generous reference point. Pick a reference well past
+    // the worst values likely to appear so different runs can be compared.
+    let ref_point = [40_000.0, 40_000.0];
+    let owned: Vec<Candidate<Vec<usize>>> = result.pareto_front.to_vec();
+    let hv = hypervolume_2d(&owned, &problem.objectives(), ref_point);
+    println!();
+    println!("Hypervolume vs. reference ({}, {}): {:.0}",
+        ref_point[0], ref_point[1], hv);
+}
diff --git a/examples/mo_jss_la01.rs b/examples/mo_jss_la01.rs
new file mode 100644
index 0000000..73932a1
--- /dev/null
+++ b/examples/mo_jss_la01.rs
@@ -0,0 +1,269 @@
+//! 3-objective Job-Shop Scheduling on Lawrence's LA01 instance, solved with
+//! NSGA-III (the many-objective successor to NSGA-II).
+//!
+//! - **Benchmark**: Lawrence LA01 (1984), 10 jobs × 5 machines, 50 operations
+//!   total. Each operation has a fixed machine and processing time;
+//!   operations within a job run in order. Data taken from the OR-Library /
+//!   JSPLIB la01 instance file.
+//! - **Three objectives** (this example):
+//!   - f₁ = makespan
+//!   - f₂ = total flow time Σⱼ Cⱼ
+//!   - f₃ = total tardiness Σⱼ max(0, Cⱼ − dⱼ), with synthetic due dates
+//!     dⱼ = 1.3 × (sum of processing times of job j)
+//! - **Algorithm**: [`Nsga3`] — designed for ≥ 3 objectives (NSGA-II's
+//!   crowding distance degrades in higher dim).
+//! - **Encoding**: operation-based string of length 50.
+//! - **Variation**: a local POX (multiset-preserving) crossover piped through
+//!   a small randomly-chosen mutation that alternates between
+//!   [`InsertionMutation`] and [`ScrambleMutation`]. Strict-permutation
+//!   crossovers cannot be used on multiset encodings.
+//! - **Initializer**: [`ShuffledMultisetPermutation`].
+//!
+//! Sources:
+//! - Lawrence (1984), thesis benchmark instances.
+//! - OR-Library / JSPLIB LA01 instance file.
+//! - Deb & Jain (2014), "An evolutionary many-objective optimization
+//!   algorithm using reference-point based non-dominated sorting approach,
+//!   Part I" — NSGA-III.
+//!
+//! Run with:
+//!
+//! ```bash
+//! cargo run --release --example mo_jss_la01
+//! ```
+
+use heuropt::prelude::*;
+use rand::Rng as _;
+
+const N_JOBS: usize = 10;
+const N_MACHINES: usize = 5;
+
+/// LA01 routing — machine id for the k-th operation of job j.
+const LA01_MACHINE: [[usize; N_MACHINES]; N_JOBS] = [
+    [1, 0, 4, 3, 2],
+    [0, 3, 4, 2, 1],
+    [3, 4, 1, 2, 0],
+    [1, 0, 4, 2, 3],
+    [0, 3, 2, 1, 4],
+    [1, 2, 4, 0, 3],
+    [3, 4, 1, 2, 0],
+    [2, 0, 1, 3, 4],
+    [3, 1, 4, 0, 2],
+    [4, 3, 1, 2, 0],
+];
+
+/// LA01 processing times — duration of the k-th operation of job j.
+const LA01_TIME: [[f64; N_MACHINES]; N_JOBS] = [
+    [21.0, 53.0, 95.0, 55.0, 34.0],
+    [21.0, 52.0, 16.0, 26.0, 71.0],
+    [39.0, 98.0, 42.0, 31.0, 12.0],
+    [77.0, 55.0, 79.0, 66.0, 77.0],
+    [83.0, 34.0, 64.0, 19.0, 37.0],
+    [54.0, 43.0, 79.0, 92.0, 62.0],
+    [69.0, 77.0, 87.0, 87.0, 93.0],
+    [38.0, 60.0, 41.0, 24.0, 66.0],
+    [17.0, 49.0, 25.0, 44.0, 98.0],
+    [77.0, 79.0, 43.0, 75.0, 96.0],
+];
+
+/// Synthetic due dates: 1.3 × total processing time of each job.
+fn due_dates() -> [f64; N_JOBS] {
+    let mut d = [0.0_f64; N_JOBS];
+    for (j, row) in LA01_TIME.iter().enumerate() {
+        d[j] = 1.3 * row.iter().sum::<f64>();
+    }
+    d
+}
+
+struct La01ThreeObjective {
+    due: [f64; N_JOBS],
+}
+
+impl Problem for La01ThreeObjective {
+    type Decision = Vec<usize>;
+
+    fn objectives(&self) -> ObjectiveSpace {
+        ObjectiveSpace::new(vec![
+            Objective::minimize("makespan"),
+            Objective::minimize("total_flow_time"),
+            Objective::minimize("total_tardiness"),
+        ])
+    }
+
+    fn evaluate(&self, schedule: &Vec<usize>) -> Evaluation {
+        let mut job_next = [0_usize; N_JOBS];
+        let mut job_clock = [0.0_f64; N_JOBS];
+        let mut machine_clock = [0.0_f64; N_MACHINES];
+        for &job in schedule {
+            let k = job_next[job];
+            let m = LA01_MACHINE[job][k];
+            let t = LA01_TIME[job][k];
+            let start = job_clock[job].max(machine_clock[m]);
+            let end = start + t;
+            job_clock[job] = end;
+            machine_clock[m] = end;
+            job_next[job] = k + 1;
+        }
+        let makespan = machine_clock.iter().cloned().fold(0.0_f64, f64::max);
+        let flow_time: f64 = job_clock.iter().sum();
+        let tardiness: f64 = job_clock
+            .iter()
+            .zip(self.due.iter())
+            .map(|(&c, &d)| (c - d).max(0.0))
+            .sum();
+        Evaluation::new(vec![makespan, flow_time, tardiness])
+    }
+
+    fn decision_schema(&self) -> Vec<DecisionVariable> {
+        (0..N_JOBS * N_MACHINES)
+            .map(|k| DecisionVariable::new(format!("op_slot_{k}")))
+            .collect()
+    }
+}
+
+/// POX — multiset-preserving crossover for operation-string encodings.
+/// (Identical in spirit to the one in `jss_ft06_bi.rs`; copied locally so
+/// each example stays self-contained.)
+#[derive(Debug, Clone, Copy, Default)]
+struct PrecedenceOrderCrossover;
+
+impl Variation<Vec<usize>> for PrecedenceOrderCrossover {
+    fn vary(&mut self, parents: &[Vec<usize>], rng: &mut Rng) -> Vec<Vec<usize>> {
+        assert!(parents.len() >= 2, "POX requires 2 parents");
+        let p1 = &parents[0];
+        let p2 = &parents[1];
+        let mut in_j1 = [false; N_JOBS];
+        loop {
+            for slot in &mut in_j1 {
+                *slot = rng.random_bool(0.5);
+            }
+            let n_in_j1 = in_j1.iter().filter(|&&b| b).count();
+            if n_in_j1 > 0 && n_in_j1 < N_JOBS {
+                break;
+            }
+        }
+        vec![pox_child(p1, p2, &in_j1), pox_child(p2, p1, &in_j1)]
+    }
+}
+
+fn pox_child(donor: &[usize], filler: &[usize], in_donor_set: &[bool]) -> Vec<usize> {
+    let n = donor.len();
+    let mut child = vec![usize::MAX; n];
+    for k in 0..n {
+        if in_donor_set[donor[k]] {
+            child[k] = donor[k];
+        }
+    }
+    let mut fill_idx = 0;
+    for &v in filler {
+        if !in_donor_set[v] {
+            while fill_idx < n && child[fill_idx] != usize::MAX {
+                fill_idx += 1;
+            }
+            child[fill_idx] = v;
+            fill_idx += 1;
+        }
+    }
+    child
+}
+
+/// Per-call random choice between Insertion and Scramble. Both preserve the
+/// multiset; flipping a coin gives the schedule access to two complementary
+/// neighborhood moves.
+#[derive(Debug, Clone, Copy, Default)]
+struct InsertionOrScramble;
+
+impl Variation<Vec<usize>> for InsertionOrScramble {
+    fn vary(&mut self, parents: &[Vec<usize>], rng: &mut Rng) -> Vec<Vec<usize>> {
+        if rng.random_bool(0.5) {
+            InsertionMutation.vary(parents, rng)
+        } else {
+            ScrambleMutation.vary(parents, rng)
+        }
+    }
+}
+
+fn main() {
+    let problem = La01ThreeObjective { due: due_dates() };
+
+    let mut optimizer = Nsga3::new(
+        Nsga3Config {
+            population_size: 120,
+            generations: 600,
+            reference_divisions: 12,
+            seed: 9,
+        },
+        ShuffledMultisetPermutation::new(vec![N_MACHINES; N_JOBS]),
+        CompositeVariation {
+            crossover: PrecedenceOrderCrossover,
+            mutation: InsertionOrScramble,
+        },
+    );
+    let result = optimizer.run(&problem);
+
+    println!("LA01 — 3-objective JSS via NSGA-III");
+    println!("Source: Lawrence (1984), OR-Library la01 instance");
+    println!();
+    println!("Objectives: f1 = makespan,  f2 = total flow time,  f3 = total tardiness");
+    println!("Due dates:  dⱼ = 1.3 × Σ(processing times of job j)");
+    println!();
+    println!("Total evaluations: {}", result.evaluations);
+    println!("Pareto-front size: {}", result.pareto_front.len());
+    println!();
+
+    // Sort by makespan and print up to 12 well-spaced rows.
+    let mut front: Vec<&Candidate<Vec<usize>>> = result.pareto_front.iter().collect();
+    front.sort_by(|a, b| {
+        a.evaluation.objectives[0]
+            .partial_cmp(&b.evaluation.objectives[0])
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+
+    let stride = (front.len() / 12).max(1);
+    println!("    f1 makespan   f2 flow time   f3 tardiness");
+    let mut printed = 0_usize;
+    for (i, c) in front.iter().enumerate() {
+        if i % stride == 0 || i + 1 == front.len() {
+            let o = &c.evaluation.objectives;
+            println!("    {:>11.0}   {:>12.0}   {:>11.0}", o[0], o[1], o[2]);
+            printed += 1;
+            if printed >= 12 {
+                break;
+            }
+        }
+    }
+    println!();
+
+    if let (Some(corner_ms), Some(corner_ft), Some(corner_td)) = (
+        front.first(),
+        front.iter().min_by(|a, b| {
+            a.evaluation.objectives[1]
+                .partial_cmp(&b.evaluation.objectives[1])
+                .unwrap_or(std::cmp::Ordering::Equal)
+        }),
+        front.iter().min_by(|a, b| {
+            a.evaluation.objectives[2]
+                .partial_cmp(&b.evaluation.objectives[2])
+                .unwrap_or(std::cmp::Ordering::Equal)
+        }),
+    ) {
+        println!(
+            "Makespan corner:  f1={:.0}, f2={:.0}, f3={:.0}",
+            corner_ms.evaluation.objectives[0],
+            corner_ms.evaluation.objectives[1],
+            corner_ms.evaluation.objectives[2],
+        );
+        println!(
+            "Flow-time corner: f1={:.0}, f2={:.0}, f3={:.0}",
+            corner_ft.evaluation.objectives[0],
+            corner_ft.evaluation.objectives[1],
+            corner_ft.evaluation.objectives[2],
+        );
+        println!(
+            "Tardiness corner: f1={:.0}, f2={:.0}, f3={:.0}",
+            corner_td.evaluation.objectives[0],
+            corner_td.evaluation.objectives[1],
+            corner_td.evaluation.objectives[2],
+        );
+    }
+}
diff --git a/examples/mo_knapsack.rs b/examples/mo_knapsack.rs
new file mode 100644
index 0000000..567cf28
--- /dev/null
+++ b/examples/mo_knapsack.rs
@@ -0,0 +1,206 @@
+//! Bi-objective 0/1 knapsack — Zitzler & Thiele's textbook multi-objective
+//! combinatorial benchmark, solved with NSGA-II.
+//!
+//! - **Benchmark family**: Zitzler & Thiele (1999) bi-objective knapsack.
+//!   Each item has two profit values and a single weight; a single capacity
+//!   constraint. We use a 30-item instance with values drawn from the same
+//!   U(10, 100) distribution scheme as the published instances, embedded as
+//!   `const` tables so the example stays self-contained.
+//! - **Algorithm**: [`Nsga2`].
+//! - **Decision**: `Vec<bool>` of length 30 (take / leave each item).
+//! - **Variation**: a local one-point crossover (binary GAs' workhorse) piped
+//!   into [`BitFlipMutation`] via [`CompositeVariation`]. **A future PR could
+//!   lift `OnePointCrossover` / `UniformCrossover` into the library proper**
+//!   so users don't need to roll their own.
+//! - **Initializer**: a tiny local `RandomBinary` (one-liner; would be a
+//!   reasonable library addition too).
+//! - **Constraint handling**: weight overruns are penalized in both
+//!   objectives by `-large * overrun`. With the penalty dominating profit
+//!   range, the Pareto front is composed entirely of feasible solutions
+//!   (standard heuristic-MO practice).
+//!
+//! Sources:
+//! - Zitzler & Thiele (1999), "Multiobjective evolutionary algorithms: A
+//!   comparative case study and the Strength Pareto approach."
+//! - Deb (2001), "Multi-Objective Optimization Using Evolutionary Algorithms"
+//!   for the standard penalty-based MO constraint handling.
+//!
+//! Run with:
+//!
+//! ```bash
+//! cargo run --release --example mo_knapsack
+//! ```
+
+use heuropt::metrics::hypervolume_2d;
+use heuropt::prelude::*;
+use rand::Rng as _;
+
+const N_ITEMS: usize = 30;
+
+/// Profit vector A (one of two objectives), U(10, 100) style.
+const PROFITS_A: [f64; N_ITEMS] = [
+    61.0, 17.0, 92.0, 49.0, 73.0, 28.0, 84.0, 36.0, 55.0, 78.0,
+    23.0, 91.0, 12.0, 67.0, 45.0, 58.0, 33.0, 71.0, 14.0, 26.0,
+    87.0, 42.0, 19.0, 65.0, 30.0, 51.0, 79.0, 22.0, 47.0, 88.0,
+];
+
+/// Profit vector B (the other objective). Intentionally anti-correlated with
+/// A on many items so the Pareto front spans a wide trade-off.
+const PROFITS_B: [f64; N_ITEMS] = [
+    24.0, 81.0, 16.0, 67.0, 29.0, 73.0, 41.0, 60.0, 52.0, 19.0,
+    77.0, 34.0, 95.0, 22.0, 71.0, 88.0, 56.0, 27.0, 64.0, 90.0,
+    18.0, 43.0, 79.0, 31.0, 85.0, 25.0, 38.0, 92.0, 70.0, 13.0,
+];
+
+/// Item weights.
+const WEIGHTS: [f64; N_ITEMS] = [
+    35.0, 58.0, 22.0, 71.0, 14.0, 86.0, 31.0, 53.0, 78.0, 19.0,
+    44.0, 16.0, 67.0, 88.0, 25.0, 51.0, 33.0, 74.0, 12.0, 47.0,
+    63.0, 28.0, 91.0, 36.0, 55.0, 17.0, 82.0, 41.0, 24.0, 68.0,
+];
+
+/// Capacity = roughly half the total weight (standard Zitzler-Thiele convention).
+fn capacity() -> f64 {
+    0.5 * WEIGHTS.iter().sum::<f64>()
+}
+
+struct BiKnapsack {
+    cap: f64,
+}
+
+impl Problem for BiKnapsack {
+    type Decision = Vec<bool>;
+
+    fn objectives(&self) -> ObjectiveSpace {
+        ObjectiveSpace::new(vec![
+            Objective::maximize("profit_A"),
+            Objective::maximize("profit_B"),
+        ])
+    }
+
+    fn evaluate(&self, take: &Vec<bool>) -> Evaluation {
+        let (pa, pb, w) = take.iter().enumerate().fold(
+            (0.0_f64, 0.0_f64, 0.0_f64),
+            |(pa, pb, w), (i, &t)| {
+                if t {
+                    (pa + PROFITS_A[i], pb + PROFITS_B[i], w + WEIGHTS[i])
+                } else {
+                    (pa, pb, w)
+                }
+            },
+        );
+        // Penalty: large coefficient on weight overrun, applied to both objectives.
+        let overrun = (w - self.cap).max(0.0);
+        let penalty = 1000.0 * overrun;
+        Evaluation::new(vec![pa - penalty, pb - penalty])
+    }
+
+    fn decision_schema(&self) -> Vec<DecisionVariable> {
+        (0..N_ITEMS)
+            .map(|i| DecisionVariable::new(format!("item_take_{i}")))
+            .collect()
+    }
+}
+
+/// Random binary initializer — each bit is 50/50 independently.
+#[derive(Debug, Clone, Copy)]
+struct RandomBinary {
+    n: usize,
+}
+
+impl Initializer<Vec<bool>> for RandomBinary {
+    fn initialize(&mut self, size: usize, rng: &mut Rng) -> Vec<Vec<bool>> {
+        (0..size)
+            .map(|_| (0..self.n).map(|_| rng.random_bool(0.5)).collect())
+            .collect()
+    }
+}
+
+/// One-point crossover for binary chromosomes.
+#[derive(Debug, Clone, Copy, Default)]
+struct OnePointCrossoverBool;
+
+impl Variation<Vec<bool>> for OnePointCrossoverBool {
+    fn vary(&mut self, parents: &[Vec<bool>], rng: &mut Rng) -> Vec<Vec<bool>> {
+        assert!(parents.len() >= 2, "OnePointCrossoverBool requires 2 parents");
+        let p1 = &parents[0];
+        let p2 = &parents[1];
+        assert_eq!(p1.len(), p2.len(), "parent lengths differ");
+        let n = p1.len();
+        if n < 2 {
+            return vec![p1.clone(), p2.clone()];
+        }
+        let cut = rng.random_range(1..n);
+        let mut c1 = Vec::with_capacity(n);
+        let mut c2 = Vec::with_capacity(n);
+        c1.extend_from_slice(&p1[..cut]);
+        c1.extend_from_slice(&p2[cut..]);
+        c2.extend_from_slice(&p2[..cut]);
+        c2.extend_from_slice(&p1[cut..]);
+        vec![c1, c2]
+    }
+}
+
+fn main() {
+    let cap = capacity();
+    let problem = BiKnapsack { cap };
+
+    let mut optimizer = Nsga2::new(
+        Nsga2Config {
+            population_size: 120,
+            generations: 400,
+            seed: 19,
+        },
+        RandomBinary { n: N_ITEMS },
+        CompositeVariation {
+            crossover: OnePointCrossoverBool,
+            mutation: BitFlipMutation { probability: 1.0 / N_ITEMS as f64 },
+        },
+    );
+    let result = optimizer.run(&problem);
+
+    println!("Bi-objective 0/1 knapsack — Zitzler–Thiele style, 30 items");
+    println!("Capacity = {:.0} (≈ half of total weight {:.0})",
+        cap, WEIGHTS.iter().sum::<f64>());
+    println!();
+    println!("Total evaluations: {}", result.evaluations);
+    println!("Pareto-front size: {}", result.pareto_front.len());
+    println!();
+
+    // Sort by profit_A descending for display, dedupe by integer-rounded objective values.
+    let mut front: Vec<&Candidate<Vec<bool>>> = result.pareto_front.iter().collect();
+    front.sort_by(|a, b| {
+        b.evaluation.objectives[0]
+            .partial_cmp(&a.evaluation.objectives[0])
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+    let mut seen: Vec<(i64, i64)> = Vec::new();
+    println!("    profit_A     profit_B   weight");
+    for c in &front {
+        let o = &c.evaluation.objectives;
+        let key = (o[0] as i64, o[1] as i64);
+        if seen.contains(&key) {
+            continue;
+        }
+        seen.push(key);
+        let w: f64 = c
+            .decision
+            .iter()
+            .enumerate()
+            .filter(|&(_, &t)| t)
+            .map(|(i, _)| WEIGHTS[i])
+            .sum();
+        println!("    {:>8.0}     {:>8.0}   {:>6.0}", o[0], o[1], w);
+    }
+    println!("    ({} unique objective-space points)", seen.len());
+
+    // Hypervolume against a reference point of (0, 0): since these are
+    // maximization objectives, we transform to minimization by negation in
+    // the metric — hypervolume_2d uses ObjectiveSpace::as_minimization() so
+    // it Just Works.
+    let ref_point = [0.0, 0.0];
+    let owned: Vec<Candidate<Vec<bool>>> = result.pareto_front.to_vec();
+    let hv = hypervolume_2d(&owned, &problem.objectives(), ref_point);
+    println!();
+    println!("Hypervolume vs. reference (profit_A=0, profit_B=0): {:.0}", hv);
+}