The determination of the correct number of equilibrium phases and their corresponding compositions at fixed temperature and pressure (TP flash) is studied. The novel aspects of this work center around unique initialization strategies and successive quadratic programming (SQP) enhancements that include the use of (I) only binary tangent plane analyses; (2) the determination of all partially miscible binary pairs and a dominant immiscible pair; (3) novel relative solubility calculations based on component activities and double tangency separation; (4) least squares solutions to compute phase fraction estimates; (5) a variety of algorithmic features that dynamically trap difficulties such as compositions well below machine accuracy, and trivial and collapsed solutions and (6) a posteriori testing of phase and solution stability. The overall algorithmic framework is one based on using a combination of binary tangent plane analyses, bubble point calculations and dimensionless Gibbs free energy minimizations. Binary tangent plane analyses are used to identify all immiscible or partially miscible binary pairs and to avoid dimensionality difficulties associated with locating all stationary points in the tangent plane distance function in the full composition space. The proposed approach consists of solving a sequence of subproblems (i.e. LE, LLE, VLLE,...) until the global minimum dimensionless Gibbs free energy (G/RT) is found. Maximum information from binary tangent plane analyses and previously solved subproblems are used to generate initial values for the next subproblem. The concept of relative solubilities is introduced and used to initialize phase compositions in all LLE calculations (i.e. phase split or flash). All completely miscible component relative solubilities are calculated using component activities while those for immiscible or partially miscible components are initialized using double tangency separation. Phase fractions are initialized using a least-square solution to the set of component mass balances. All subproblems are formulated in terms of component hows and solved using a full space SQP method using a modified Broyden-Fletcher-Goldfarb-Shanno (BFGS) update of the Lagrangian Hessian matrix. The proposed algorithm was tested within the Aspen Plus process simulator using a Variety of physical properties options. Twenty six multicomponent mixtures including some four-phase (VLLLE) emulsion polymerization problems were used to test the proposed algorithm. All problems were easily solved and clearly demonstrate the capabilities of the present multiphase TP flash model. (C) 2000 Elsevier Science Ltd. All rights reserved.