Approximate computing has gained high attention in various applications that can benefit from a reduction in costs by lowering the accuracy. In this paper we present an optimization approach for functional approximation of Binary Decision Diagrams (BDDs) which are known for their widespread applications in electronic design automation and formal verification. We propose a three-objective epsilon-preferred evolutionary algorithm with the first objective set to the BDD size which is given higher priority to the two other objectives set to errors caused by approximation. This is highly demanded by the application to ensure that the minimum size for the approximated BDD is accessible when the error metrics meet certain threshold values. While BDD size minimization is guaranteed by incorporating priority, the use of epsilon in the proposed approach ensures to guide the search towards desired error values in parallel. Experiments confirm the efficiency of the proposed approach by a size improvement of 64.24% at a fair cost of 3.86% inaccuracy on average.