The fitting of physical models is often done only using a single target observable. However, when multiple targets are considered, the fitting procedure becomes cumbersome, there being no easy way to quantify the robustness of the model for all different observables. Here, we illustrate that one can jointly search for the best model for each desired observable through multi-objective optimization. To do so we construct the Pareto front to study if there exists a set of parameters of the model that can jointly describe multiple, or all, observables. To alleviate the computational cost, the predicted error for each targeted objective is approximated with a Gaussian process model, as it is commonly done in the Bayesian optimization framework. We applied this methodology to improve three different models used in the simulation of stationary state $cis-trans$ photoisomerization of retinal in rhodopsin. Optimization was done with respect to different experimental measurements, including emission spectra, peak absorption frequencies for the $cis$ and $trans$ conformers, and the energy storage.