This paper aims to enhance the accuracy of the Pareto front estimation model as an aggregative representation of the non-dominated solutions and proposes an algorithm named the Pareto front Model Optimization Algorithm (PFMOA). The typical output of multi-objective optimization is a set of non-dominated solutions approximating the Pareto front, which is the optimal trade-off between objective function values. The more non-dominated solutions there are, the more accurately the Pareto front can be approximated. However, especially in real-world problems, there is often a limitation on increasing the number of solutions due to the time required to execute objective functions. For the issue, a Pareto front estimation method interpolates between a limited number of non-dominated solutions to represent changes in objective function values even in regions where non-dominated solutions are not actually obtained. The proposed PFMOA enhances the accuracy of the Pareto front estimation model while evaluating new solutions. PFMOA generates the estimated Pareto front based on a known non-dominated solution set using Kriging. PFMOA focuses on the point with the lowest confidence levels on the estimated Pareto front and evaluates the corresponding point on the estimated Pareto set. PFMOA also generates new solutions using evolutionary variation. PFMOA stochastically switches these two model-based and evolutionary-based solution generation methods. If the new solution is non-dominated, it is included in the known non-dominated solution set, and the process is repeated to enhance the accuracy of the Pareto front estimation model. The effectiveness of the proposed PFMOA is verified using DTLZ1-3, and WFG4 problems. The results show that, in all cases, the accuracy of the Pareto front approximation by the Pareto front estimation model is higher than that by the obtained solution set itself. Additionally, the combination of model-based and evolutionary-based solution generation is bene...