Several approaches have been proposed for the computation of solutions to the phase and chemical equilibrium problem when the problem is posed as the minimization of the Gibbs free energy function. None of them can guarantee convergence to the true optimal solution, and are highly dependent on the supplied initial point. Convergence to local solutions often occurs, yielding incorrect phase and component distributions. This work examines the problems when the liquid phase is adequately modeled by the Non-Random Two Liquid (NRTL) activity coefficient expression and the vapor phase is assumed to be ideal. The contribution of the proposed approach is two-fold. Firstly, a novel and important property of the Gibbs free energy expression involving the NRTL equation is provided. It is subsequently shown that by introducing new variables, the problem can then be transformed into one where a biconvex objective function is minimized over a set of bilinear constraints. Secondly, the Global OPtimization (GOP) algorithm is used to exploit these induced properties of the formulation to guarantee convergence to an epsilon-global solution, regardless of the starting point. A geometrical interpretation is provided for a selected smaller, but challenging, example. Numerous examples are presented which demonstrate the broad applicability of the proposed approach.