Many researchers in industry and academia are showing an increasing interest in the definition of fuel surrogates for Computational Fluid Dynamics simulation applications. This need is mainly driven by the necessity of the engine research community to anticipate the effects of new gasoline formulations and combustion modes (e.g., Homogeneous Charge Compression Ignition, Spark Assisted Compression Ignition) to meet future emission regulations. Since those solutions strongly rely on the tailored mixture distribution, the simulation and accurate prediction of the mixture formation will be mandatory. Focusing purely on the definition of surrogates to emulate liquid phase and liquid-vapor equilibrium of gasolines, the following target properties are considered in this work: density, Reid vapor pressure, chemical macro-composition and volatility. A set of robust algorithms has been developed for the prediction of volatility and Reid vapor pressure. A Bayesian optimization algorithm based on a customized merit function has been developed to allow for the efficient definition of surrogate formulations from a palette of 15 pure compounds. The developed methodology has been applied on different real gasolines from literature in order to identify their optima surrogates. Furthermore, the ‘unicity’ of the surrogate composition is discussed by comparing the optimum solution with the most different one available in the pool of equivalent-valuable solutions. The proposed methodology has proven the potential to formulate surrogates characterized by an overall good agreement with the target properties of the experimental gasolines (max relative error below 10%, average relative error around 3%). In particular, the shape and the end-tails of the distillation curve are well captured. Furthermore, an accurate prediction of key chemical macro-components such as ethanol and aromatics and their influence on evaporative behavior is achieved. The study of the ‘unicity’ of the surrogate composition has revealed that (i) the unicity is strongly correlated with the accuracy and that (ii) both ‘unicity’ and accuracy of the prediction are very sensitive to the high presence of aromatics.