The activity coefficient models, when based on parameters obtained through binary phase equilibrium information result in multiple sets that may either not be able to predict ternary and higher order azeotropes, or computed values are significantly different from experimental data, particularly at thermodynamic landmarks. Sensitivity of azeotropic states to model parameters and operating conditions is imbedded into a homotopy continuation formulation, providing a functionally sound basis for prediction of such states and their corresponding sensitivity. Furthermore, the validity of parameter set depends on the range of operating variables. The methodology allows systematic variation of model parameters and selection of a parameter set that allows computation of thermodynamic landmarks, such as azeotropy. The systematic analysis of the effects of parametric perturbation on the prediction of azeotropes and heteroazeotropes also provide useful insights into capabilities and limitations of phase equilibrium models. The technique is also useful in developing the effect of system variables on azeotropy conditions. We demonstrated the utility of our approach through modeling benzene-cyclohexane-n-propanol, ethyl acetate-ethanol-water and water-formic acid-acetic acid system. (c) 2005 Elsevier B.V. All rights reserved.