The ReaxFF incorporates complex functions with associated parameters in order to describe the inter and intra-atomic interactions in materials systems. A typical ReaxFF force field consists of approximately a hundred parameters per element type. During the development of a force field for a molecular system of interest, these parameters are optimized to reproduce reference values with reasonable accuracy, which is a multi-objective optimization problem. A standard approach is to reduce multi-objective optimization problem to a one-dimensional optimization problem by defining the objective function in the certain form. However, the optimization problem is non-trivial and challenging due to the high dimensionality of the input space and the difficulty of fitting the objective function well with a surrogate model. In this talk, we propose an approach that is based on a combination of modern machine learning approaches and careful ``statistical thinking'' to achieve the target.