diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index ff3803a..b2840c1 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -15,6 +15,7 @@ pub mod paes;
 pub(crate) mod parallel_eval;
 pub mod particle_swarm;
 pub mod random_search;
+pub mod rvea;
 pub mod simulated_annealing;
 pub mod sms_emoa;
 pub mod spea2;
@@ -35,6 +36,7 @@ pub use nsga3::*;
 pub use paes::*;
 pub use particle_swarm::*;
 pub use random_search::*;
+pub use rvea::*;
 pub use simulated_annealing::*;
 pub use sms_emoa::*;
 pub use spea2::*;
diff --git a/src/algorithms/rvea.rs b/src/algorithms/rvea.rs
new file mode 100644
index 0000000..7a1b821
--- /dev/null
+++ b/src/algorithms/rvea.rs
@@ -0,0 +1,328 @@
+//! `Rvea` — Cheng, Jin, Olhofer & Sendhoff 2016 Reference Vector-guided EA.
+
+use rand::Rng as _;
+
+use crate::algorithms::parallel_eval::evaluate_batch;
+use crate::core::candidate::Candidate;
+use crate::core::objective::ObjectiveSpace;
+use crate::core::population::Population;
+use crate::core::problem::Problem;
+use crate::core::result::OptimizationResult;
+use crate::core::rng::rng_from_seed;
+use crate::pareto::front::{best_candidate, pareto_front};
+use crate::pareto::reference_points::das_dennis;
+use crate::traits::{Initializer, Optimizer, Variation};
+
+/// Configuration for [`Rvea`].
+#[derive(Debug, Clone)]
+pub struct RveaConfig {
+    /// Constant population size.
+    pub population_size: usize,
+    /// Number of generations.
+    pub generations: usize,
+    /// Number of divisions `H` for Das–Dennis reference vectors. Pop size
+    /// should be roughly `binomial(H + M − 1, M − 1)`.
+    pub reference_divisions: usize,
+    /// Penalty exponent `α`. The paper recommends 2.0.
+    pub alpha: f64,
+    /// Seed for the deterministic RNG.
+    pub seed: u64,
+}
+
+impl Default for RveaConfig {
+    fn default() -> Self {
+        Self {
+            population_size: 100,
+            generations: 250,
+            reference_divisions: 12,
+            alpha: 2.0,
+            seed: 42,
+        }
+    }
+}
+
+/// Reference Vector-guided Evolutionary Algorithm.
+#[derive(Debug, Clone)]
+pub struct Rvea<I, V> {
+    /// Algorithm configuration.
+    pub config: RveaConfig,
+    /// Initial-decision sampler.
+    pub initializer: I,
+    /// Offspring-producing variation operator.
+    pub variation: V,
+}
+
+impl<I, V> Rvea<I, V> {
+    /// Construct an `Rvea`.
+    pub fn new(config: RveaConfig, initializer: I, variation: V) -> Self {
+        Self { config, initializer, variation }
+    }
+}
+
+impl<P, I, V> Optimizer<P> for Rvea<I, V>
+where
+    P: Problem + Sync,
+    P::Decision: Send,
+    I: Initializer<P::Decision>,
+    V: Variation<P::Decision>,
+{
+    fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
+        assert!(self.config.population_size > 0, "Rvea population_size must be > 0");
+        let n = self.config.population_size;
+        let objectives = problem.objectives();
+        let m = objectives.len();
+        // Reference vectors normalized to unit norm.
+        let raw_refs = das_dennis(m, self.config.reference_divisions);
+        let references: Vec<Vec<f64>> = raw_refs.into_iter().map(unit_normalize).collect();
+        assert!(!references.is_empty(), "Rvea: no reference vectors generated");
+
+        // Smallest angle between any two reference vectors — used to scale
+        // the APD penalty term.
+        let theta_max = smallest_neighbor_angle(&references);
+        let mut rng = rng_from_seed(self.config.seed);
+
+        let initial_decisions = self.initializer.initialize(n, &mut rng);
+        let mut population: Vec<Candidate<P::Decision>> =
+            evaluate_batch(problem, initial_decisions);
+        let mut evaluations = population.len();
+
+        for gen_idx in 0..self.config.generations {
+            // Phase 1: random parent selection + variation.
+            let mut offspring_decisions: Vec<P::Decision> = Vec::with_capacity(n);
+            while offspring_decisions.len() < n {
+                let p1 = rng.random_range(0..population.len());
+                let p2 = rng.random_range(0..population.len());
+                let parents =
+                    vec![population[p1].decision.clone(), population[p2].decision.clone()];
+                let children = self.variation.vary(&parents, &mut rng);
+                assert!(!children.is_empty(), "Rvea variation returned no children");
+                for child in children {
+                    if offspring_decisions.len() >= n {
+                        break;
+                    }
+                    offspring_decisions.push(child);
+                }
+            }
+            let offspring = evaluate_batch(problem, offspring_decisions);
+            evaluations += offspring.len();
+
+            // Combine + APD-based survival.
+            let mut combined: Vec<Candidate<P::Decision>> = Vec::with_capacity(2 * n);
+            combined.extend(population);
+            combined.extend(offspring);
+
+            // Ideal point z*.
+            let m_dim = m;
+            let mut ideal = vec![f64::INFINITY; m_dim];
+            for c in &combined {
+                let oriented = objectives.as_minimization(&c.evaluation.objectives);
+                for (k, v) in oriented.iter().enumerate() {
+                    if *v < ideal[k] {
+                        ideal[k] = *v;
+                    }
+                }
+            }
+            // Translate.
+            let translated: Vec<Vec<f64>> = combined
+                .iter()
+                .map(|c| {
+                    let oriented = objectives.as_minimization(&c.evaluation.objectives);
+                    oriented.iter().enumerate().map(|(k, v)| v - ideal[k]).collect()
+                })
+                .collect();
+
+            // Associate each member with its closest-angle reference vector.
+            let mut assoc: Vec<usize> = vec![0; combined.len()];
+            let mut angles: Vec<f64> = vec![0.0; combined.len()];
+            for (i, t) in translated.iter().enumerate() {
+                let (best_ref, best_angle) = closest_reference(t, &references);
+                assoc[i] = best_ref;
+                angles[i] = best_angle;
+            }
+
+            // For each occupied reference vector, keep the member with the
+            // smallest APD score.
+            let alpha_t = (gen_idx as f64 / (self.config.generations as f64).max(1.0))
+                .powf(self.config.alpha);
+            let mut keep: Vec<Option<(usize, f64)>> = vec![None; references.len()];
+            for i in 0..combined.len() {
+                let r = assoc[i];
+                let length: f64 = translated[i].iter().map(|v| v * v).sum::<f64>().sqrt();
+                let theta_max_safe = theta_max.max(1e-12);
+                let penalty = 1.0 + (m_dim as f64) * alpha_t * (angles[i] / theta_max_safe);
+                let apd = penalty * length;
+                match keep[r] {
+                    None => keep[r] = Some((i, apd)),
+                    Some((_, current)) if apd < current => keep[r] = Some((i, apd)),
+                    _ => {}
+                }
+            }
+
+            let mut next: Vec<Candidate<P::Decision>> =
+                keep.into_iter().flatten().map(|(i, _)| combined[i].clone()).collect();
+            // If we ended up with fewer than n (some references unfilled),
+            // backfill with the lowest-APD remaining candidates.
+            if next.len() < n {
+                let mut all_apds: Vec<(usize, f64)> = (0..combined.len())
+                    .map(|i| {
+                        let length: f64 =
+                            translated[i].iter().map(|v| v * v).sum::<f64>().sqrt();
+                        let theta_max_safe = theta_max.max(1e-12);
+                        let penalty = 1.0 + (m_dim as f64) * alpha_t * (angles[i] / theta_max_safe);
+                        (i, penalty * length)
+                    })
+                    .collect();
+                all_apds
+                    .sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
+                for (i, _) in all_apds {
+                    if next.len() >= n {
+                        break;
+                    }
+                    if !next
+                        .iter()
+                        .any(|c| std::ptr::eq(c as *const _, &combined[i] as *const _))
+                    {
+                        next.push(combined[i].clone());
+                    }
+                }
+            }
+            // If too many (only possible if the reference set has > n
+            // vectors), truncate by APD.
+            if next.len() > n {
+                next.truncate(n);
+            }
+            population = next;
+        }
+
+        let front = pareto_front(&population, &objectives);
+        let best = best_candidate(&population, &objectives);
+        OptimizationResult::new(
+            Population::new(population),
+            front,
+            best,
+            evaluations,
+            self.config.generations,
+        )
+    }
+}
+
+fn unit_normalize(mut v: Vec<f64>) -> Vec<f64> {
+    let n: f64 = v.iter().map(|x| x * x).sum::<f64>().sqrt();
+    if n > 1e-12 {
+        for x in v.iter_mut() {
+            *x /= n;
+        }
+    }
+    v
+}
+
+fn closest_reference(point: &[f64], references: &[Vec<f64>]) -> (usize, f64) {
+    let length: f64 = point.iter().map(|v| v * v).sum::<f64>().sqrt().max(1e-12);
+    let mut best = 0;
+    let mut best_angle = f64::INFINITY;
+    for (i, r) in references.iter().enumerate() {
+        let dot: f64 = point.iter().zip(r.iter()).map(|(a, b)| a * b).sum();
+        let cosine = (dot / length).clamp(-1.0, 1.0);
+        let angle = cosine.acos();
+        if angle < best_angle {
+            best_angle = angle;
+            best = i;
+        }
+    }
+    (best, best_angle)
+}
+
+fn smallest_neighbor_angle(references: &[Vec<f64>]) -> f64 {
+    let mut min_angle = f64::INFINITY;
+    for i in 0..references.len() {
+        for j in (i + 1)..references.len() {
+            let dot: f64 = references[i]
+                .iter()
+                .zip(references[j].iter())
+                .map(|(a, b)| a * b)
+                .sum();
+            let angle = dot.clamp(-1.0, 1.0).acos();
+            if angle < min_angle {
+                min_angle = angle;
+            }
+        }
+    }
+    if !min_angle.is_finite() { std::f64::consts::FRAC_PI_4 } else { min_angle }
+}
+
+#[allow(unused_imports)]
+use crate::core::objective::Objective;
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::operators::{
+        CompositeVariation, PolynomialMutation, RealBounds, SimulatedBinaryCrossover,
+    };
+    use crate::tests_support::SchafferN1;
+
+    fn make_optimizer(
+        seed: u64,
+    ) -> Rvea<RealBounds, CompositeVariation<SimulatedBinaryCrossover, PolynomialMutation>> {
+        let bounds = vec![(-5.0, 5.0)];
+        let initializer = RealBounds::new(bounds.clone());
+        let variation = CompositeVariation {
+            crossover: SimulatedBinaryCrossover::new(bounds.clone(), 15.0, 0.5),
+            mutation: PolynomialMutation::new(bounds, 20.0, 1.0),
+        };
+        Rvea::new(
+            RveaConfig {
+                population_size: 20,
+                generations: 15,
+                reference_divisions: 19,
+                alpha: 2.0,
+                seed,
+            },
+            initializer,
+            variation,
+        )
+    }
+
+    #[test]
+    fn produces_pareto_front() {
+        let mut opt = make_optimizer(1);
+        let r = opt.run(&SchafferN1);
+        assert!(!r.pareto_front.is_empty());
+    }
+
+    #[test]
+    fn deterministic_with_same_seed() {
+        let mut a = make_optimizer(99);
+        let mut b = make_optimizer(99);
+        let ra = a.run(&SchafferN1);
+        let rb = b.run(&SchafferN1);
+        let oa: Vec<Vec<f64>> =
+            ra.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect();
+        let ob: Vec<Vec<f64>> =
+            rb.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect();
+        assert_eq!(oa, ob);
+    }
+
+    #[test]
+    #[should_panic(expected = "population_size must be > 0")]
+    fn zero_population_size_panics() {
+        let bounds = vec![(0.0, 1.0)];
+        let initializer = RealBounds::new(bounds.clone());
+        let variation = CompositeVariation {
+            crossover: SimulatedBinaryCrossover::new(bounds.clone(), 15.0, 0.5),
+            mutation: PolynomialMutation::new(bounds, 20.0, 1.0),
+        };
+        let mut opt = Rvea::new(
+            RveaConfig {
+                population_size: 0,
+                generations: 1,
+                reference_divisions: 5,
+                alpha: 2.0,
+                seed: 0,
+            },
+            initializer,
+            variation,
+        );
+        let _ = opt.run(&SchafferN1);
+    }
+}
diff --git a/src/prelude.rs b/src/prelude.rs
index 656cefe..4720687 100644
--- a/src/prelude.rs
+++ b/src/prelude.rs
@@ -26,7 +26,8 @@ pub use crate::algorithms::{
     DifferentialEvolutionConfig,
     GeneticAlgorithm, GeneticAlgorithmConfig, HillClimber, HillClimberConfig, Hype,
     HypeConfig, Ibea, IbeaConfig, Moead, MoeadConfig, Mopso, MopsoConfig, Nsga2, Nsga2Config, Nsga3, Nsga3Config, Paes, PaesConfig, ParticleSwarm,
-    ParticleSwarmConfig, RandomSearch, RandomSearchConfig, SimulatedAnnealing,
+    ParticleSwarmConfig, RandomSearch, RandomSearchConfig, Rvea, RveaConfig,
+    SimulatedAnnealing,
     SimulatedAnnealingConfig, SmsEmoa, SmsEmoaConfig, Spea2, Spea2Config, TabuSearch,
     TabuSearchConfig, Umda,
     UmdaConfig,