Predictions of the solubility of flavonoids in a large variety of ionic liquids (ILs) with over 1800 available structures were examined based on COSMO-RS computation. The results show that the solubilities of flavonoids are strongly anion-dependent. Experimental measurement of the solubilities of esculin and rutin in 12 ILs with varying anions and cations show that predicted and experimental results generally have a good agreement. Based on the sound physical basis of COSMO-RS, the solubility changes of flavonoids were quantitatively associated with solvation interactions and structural characteristics of ILs. COSMO-RS derived parameters, i.e. misfit, H-bonding and van der Waals interaction energy, are shown to be capable of characterizing the complicated multiple interactions in the IL system effectively. H-bonding interaction is the most dominant interaction for ILs (followed by misfit and van der Waals interactions) to determine the solubility of flavonoids, and the anionic part has greater effect on the overall H-bonding capability of the IL. Based on basicity of anions, ILs were categorized into 3 groups, corresponding to the classification of the solubility of flavonoid. COSMO sigma-moment descriptors, which roughly denote the characteristic properties of the ILs, might be of general value to have a fast estimation for the solubilities of flavonoids as well as those compounds with massive moieties as H-bonding donors. The results obtained in this work may be important for achieving an improved understanding of IL solvations and the tailoring of the desired structures of ILs used as the media for efficient enzymatic esterification of flavonoids.