Hypervolume is most likely the most often used set quality indicator in (evolutionary) multi-objective optimization. It may be used to compare the quality of solution sets whose images in the objective space are approximations of the true Pareto front. Although in this way we may compare two or more approximations, our knowledge is limited without the reference to the hypervolume of the true Pareto front. Calculating the true Pareto front is, however, unrealistic for most problems of practical interest. Thus in this paper, we propose an algorithm based on solving a series of scalarization problems with an exact solver, which provides consistent and unbiased estimation of hypervolume of the unknown true Pareto front, preserving linearity and differentiability of the original multi-objective problem. The estimation is associated with a confidence interval. The proposed algorithm is implemented in Julia programming language. It is numerically evaluated through several experiments using a multi-objective combinatorial optimization problem. The computational experiment confirms correctness of the proposed approach, shows that it may be applied to much larger instances than generation of the true Pareto front, and showcases illustrative applications for estimating hypervolume gaps between approximations of different nature and the true Pareto front.