Efficient computation of two-dimensional solution sets maximizing the epsilon-indicator Abstract The majority of empirical comparisons of multi-objective evolutionary algorithms (MOEAs) are performed on synthetic benchmark functions. One of the advantages of synthetic test functions is the a-priori knowledge of the optimal Pareto front. This allows measuring the proximity to the optimal front for the solution sets returned by the different MOEAs. Such a comparison is only meaningful if the cardinality of all solution sets is bounded by some fixed k. In order to compare MOEAs to the theoretical optimum achievable with k solutions, we determine best possible ε-indicator values achievable with solution sets of size k, up to an error of δ. We present a new algorithm with runtime O(k · log2(δ-1)), which is an exponential improvement regarding the dependence on the error δ compared to all previous work. We show mathematical correctness of our algorithm and determine optimal solution sets for sets of cardinality k ∈ {2, 3, 4, 5, 10, 20, 50, 100, 1000} for the well known test suits DTLZ, ZDT, WFG and LZ09 up to error δ = 10°-25.