One of the most important applications of thermodynamics is the accurate prediction of fluid phase equilibria problems related to real chemical engineering processes. Various equations of state as well as activity coefficient models have been developed for such calculations with many interaction, size, and randomness parameters, which should be optimized based on powerful and effective computational methods. Leading to globally optimal values, genetic algorithm (GA) as a powerful and effective tool can be used for prediction of the interaction parameters of thermodynamic models in complex liquid-liquid equilibrium (LLE) systems. It requires only lower and upper bounds for the interaction parameters and the necessary initial guesses are produced automatically. In the present work, based on the GA method, a global optimization procedure is introduced for calculation of the binary interaction and nonrandomness parameters of NRTL, NRTL1, and NRTL2 activity coefficient models for 20 ternary aromatic extraction systems containing 16 different ionic liquids at various temperatures. The values of the parameters along with the root-mean-square deviations (rmsd) are reported. The results, in terms of rmsd for NRTL, NRTL1, and NRTL2 models, are very satisfactory, with global values of 0.0031, 0.0020, and 0.0053 for 187 tie-lines respectively. The obtained rmsd values for the NRTL model using the GA method are better than those reported in the literature. The rmsd results for the three studied models show that NRTL1 can handle the LLE calculations with more accuracy than the original NRTL and NRTL2 activity coefficient models.