diff --git a/NAMESPACE b/NAMESPACE
index 0f2cc50..8b35d48 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -74,6 +74,7 @@ export(setupParentSelector)
 export(setupRecombinator)
 export(setupStoppingConditions)
 export(setupSurvivalSelector)
+export(smsemoa)
 export(which.dominated)
 export(which.nondominated)
 import(BBmisc)
diff --git a/R/emoa.sms-emoa.R b/R/emoa.sms-emoa.R
new file mode 100644
index 0000000..233c8f4
--- /dev/null
+++ b/R/emoa.sms-emoa.R
@@ -0,0 +1,100 @@
+#' @title
+#'   Implementation of the SMS-EMOA by Emmerich et al.
+#'
+#' @note
+#'   This helper function hides the regular \pkg{ecr} interface and offers a more
+#'   R like interface of this state of the art EMOA while still being quite adaptable.
+#'
+#' @param task [\code{ecr_optimization_task} | \code{smoof_function}]\cr
+#'   Optimization task or objective function of type \code{smoof_function}.
+#' @param n.population [\code{integer(1)}]\cr
+#'   Population size. Default is \code{100}.
+#' @param n.offspring [\code{integer(1)}]\cr
+#'   Offspring size, i.e., number of individuals generated by variation operators
+#'   in each iteration. Default is \code{n.population}.
+#' @param ref.point [\code{numeric}]\cr
+#'   Reference point for the hypervolume computation. Default is (11, ..., 11)'
+#'   with the corresponding dimension.
+#' @template arg_parent_selector
+#' @template arg_mutator
+#' @template arg_recombinator
+#' @param max.iter [\code{integer(1)}]\cr
+#'   Maximal number of iterations. Default ist \code{100L}.
+#' @param max.evals [\code{integer(1)}]\cr
+#'   Maximal number of iterations/generations. Default is \code{Inf}.
+#' @param max.time [\code{integer(1)}]\cr
+#'   Time budget in seconds. Default ist \code{Inf}.
+#' @return [\code{ecr_ecr_multi_objective_result}]
+#' @export
+smsemoa = function(
+  task,
+  n.population = 100L, n.offspring = n.population,
+  ref.point = NULL,
+  parent.selector = makeSimpleSelector(),
+  mutator = makeGaussMutator(),
+  recombinator = makeCrossoverRecombinator(),
+  max.iter = 100L,
+  max.evals = NULL,
+  max.time = NULL) {
+
+  hypervolumeSelector = makeSelector(
+    selector = function(population, storage, n.select, control) {
+      fitness = population$fitness
+      population = population$individuals
+
+      # do non-dominated sorting
+      nds.res = doNondominatedSorting(fitness)
+      ranks = nds.res$ranks
+      idx.max = which(ranks == max(ranks))
+
+      # there is exactly one individual that is "maximally" dominated
+      if (length(idx.max) == 1L) {
+        return(makePopulation(population[-idx.max], fitness[, -idx.max, drop = FALSE]))
+      }
+
+      # compute exclusive hypervolume contributions and remove the one with the smallest
+      hvctrbs = computeHypervolumeContribution(fitness[, idx.max, drop = FALSE], ref.point = control$ref.point)
+      die.idx = idx.max[getMinIndex(hvctrbs)]
+
+      return(makePopulation(population[-die.idx], fitness[, -die.idx, drop = FALSE]))
+    },
+    supported.objectives = "multi-objective",
+    name = "Hypervolume contribution selector",
+    description = "description"
+  )
+
+  if (isSmoofFunction(task)) {
+    task = makeOptimizationTask(task)
+  }
+  assertClass(task, "ecr_optimization_task")
+
+  # SMS-EMOA control object
+  ctrl = setupECRControl(
+    n.population = n.population,
+    n.offspring = n.offspring,
+    representation = "float",
+    monitor = makeConsoleMonitor(),
+    stopping.conditions = list(
+      makeMaximumEvaluationsStoppingCondition(max.evals),
+      makeMaximumTimeStoppingCondition(max.time),
+      makeMaximumIterationsStoppingCondition(max.iter)
+    )
+  )
+  ctrl = setupEvolutionaryOperators(
+    ctrl,
+    parent.selector = parent.selector,
+    recombinator = recombinator,
+    mutator = mutator,
+    survival.selector = hypervolumeSelector
+  )
+
+  ctrl$ref.point = ref.point
+  if (is.null(ref.point)) {
+    ctrl$ref.point = rep(11, task$n.objectives)
+  }
+  if (length(ctrl$ref.point) != task$n.objectives) {
+    stopf("Reference point ref.point needs to have as many components as objectives.")
+  }
+
+  return(doTheEvolution(task, ctrl))
+}
diff --git a/tests/testthat/test_emoa.R b/tests/testthat/test_emoa.R
index 2a0297a..bc01e45 100644
--- a/tests/testthat/test_emoa.R
+++ b/tests/testthat/test_emoa.R
@@ -1,16 +1,44 @@
 context("Evolutionary Multi-Objective Algorithms")
 
 test_that("preimplemented EMOAs work well", {
-  fn = smoof::makeZDT1Function(dimensions = 2L)
-  res = nsga2(
-    makeOptimizationTask(fn),
-    n.population = 10L,
-    n.offspring = 3L,
-    max.evals = 100L
+
+  #FIXME: is something like this in BBmisc? Needed that multiple times now.
+  printList = function(l) {
+    ns = names(l)
+    pairs = sapply(ns, function(n) {
+      paste0(n, l[[n]], sep = "=")
+    })
+    collapse(pairs, ", ")
+  }
+
+  expect_is_pareto_approximation = function(pf, n.obj, algo, prob, algo.pars) {
+    info.suffix = sprintf("Algo '%s' failed on problem '%s' with params '%s'",
+      algo, prob, printList(algo.pars))
+    expect_equal(nrow(pf), length(which.nondominated(t(pf))), info = paste0(info.suffix, "Not all returned points are nondominated."))
+    expect_true(all(is.numeric(pf)), info = paste0(info.suffix, "Not all returned points are numeric."))
+    expect_equal(ncol(pf), n.obj, info = paste0(info.suffix, "Number of columns is not equal to the number of objectives."))
+  }
+
+  fns = list(
+    zdt1 = smoof::makeZDT1Function(dimensions = 2L),
+    zdt2 = smoof::makeZDT2Function(dimensions = 2L),
+    zdt3 = smoof::makeZDT3Function(dimensions = 2L)
   )
-  pf = res$pareto.front
-  print(pf)
-  expect_equal(nrow(pf), length(which.nondominated(t(pf))))
-  expect_true(all(is.numeric(pf)))
-  expect_equal(ncol(pf), getNumberOfObjectives(fn))
+  max.evals = 100L
+
+  for (emoa in c("nsga2", "smsemoa")) {
+    for (n.pop in c(5, 10, 15)) {
+      for (fn in names(fns)) {
+        res = do.call(emoa, list(
+          task = makeOptimizationTask(fns[[fn]]),
+          n.population = n.pop,
+          n.offspring = 5L,
+          max.evals = max.evals)
+        )
+        expect_is_pareto_approximation(res$pareto.front, 2L, emoa, fn,
+          list(n.pop = n.pop, n.offspring = 5L, max.evals = max.evals)
+        )
+      }
+    }
+  }
 })