A systematic, formal approach to optimization of force field parameters for molecular simulations is presented. The procedure is based on response surface mapping methodology that allows simultaneous parameter optimization against multiple property targets while constraining the number of required computationally expensive numerical experiments. The approach was implemented for prediction of vapor-liquid equilibrium properties of alkanes, alkenes, and their fluorinated derivatives via Monte Carlo molecular simulations. To further reduce computational costs, a bootstrap procedure that involves a sequence of parameter optimization for four pairs of compounds (ethane and propane, ethene and propene, perfluoroethane and perfluoropropane, and 2,3,3,3-tetrafluoropropene and (E)-1,3,3,3-tetrafluoropropene) was used. The results of simulations utilizing the optimized force field parameters agree well with associated reference equations of state.