A rigorous method is proposed for parameter estimation in thermodynamic models used for predicting vapor-liquid and vapor-liquid-liquid equilibrium of binary mixtures. Typically, vapor-liquid(-liquid) equilibria are modeled using asymmetric models, i.e., an excess Gibbs free energy model for the liquid phase and an equation of state for the vapor phase. Parameter estimation for asymmetric models is a nontrivial task. The existence of spurious additional liquid splits predicted by the models, but not observed experimentally, must be excluded and the capability of the models to predict homogeneous or heterogeneous azeotropes accurately is crucial. The proposed formulation for parameter estimation in this class of nonideal phase equilibria is based on bilevel optimization with nonconvex lower-level programs. The bilevel formulation proposed can reliably estimate model parameter values that can be employed to correlate the correct (experimentally measured) number of phases, phase splits and azeotropes and can distinguish between heterogeneous and homogeneous azeotropes. Using the bilevel formulation multiple minimization problems are solved globally: the upper level program (minimization of the discrepancy between measured and predicted values for the species compositions in each phase) that yields the best possible thermodynamic description of the systems studied, and several lower-level programs per experiment that enforce thermodynamic stability and correlation of the correct number and type of phases and phase splits. Examples from the literature are given, in which conventional parameter estimation methods yield significant qualitative and quantitative errors in the prediction of the phase behavior using the NRTL, UNIQUAC and Wilson models.