Excess Gibbs free energy models contain parameters which for a given mixture are estimated from measurements of phase-splits. Local composition models are very flexible and can accurately predict complex phase behavior. However, in many cases it has been reported that use of local composition models leads to prediction of more phase splits or more phases than measured, modeling homogeneous azeotropes as heterogeneous, etc. Here, a formulation is proposed that addresses these limitations of current parameter estimation methods. The formulation is based on a bilevel program, i.e., an optimization problem embedded in another one. Minimizing the error between model predictions and measurements gives the upper-level program. The lower-level programs are given by the minimization of the Gibbs free energy, or equivalently the satisfaction of the Gibbs tangent plane criterion. Each of the experiments is cast as a separate lower-level program. Additional requirements on the phase behavior of the system, such as enforcing the correct number of phase splits and the correct number of phases in each phase split, are similarly formulated as additional lower-level programs. Global optimization techniques are necessary even to obtain a feasible point since the lower-level programs are nonconvex. The proposed formulation is applied to problems from the literature, in which inappropriate fitting of the parameters of the non-random two-liquid (NRTL) model to experimental data has been reported to result in significant model errors, such as the prediction of an additional spurious phase split. The discussion is restricted to binary mixtures, however, the formulation can in principle be applied also to multicomponent mixtures.