Liquid-liquid extraction is an important separation technology in chemical processing industries. Due to its low energy consumption, extraction remains a competitive and sustainable option when distillation is not feasible or too expensive.Generally, there are three types of methods to design multistage liquid extractors: (1) graphical [1,2], (2) simulation-based , and (3) optimization-based . Graphical methods are based on the graphic representation of the liquid-liquid equilibrium (LLE) relationship and material and energy balances. Simulation-based methods use iterative procedures to simulate the process and find the desired configuration. Finally, optimization-based methods identify a configuration by solving a mathematical programming model which comprises sets of constraints, representing: the mass and energy balances; equilibrium relationships; and an objective function representing the design objective. Both graphical and simulation methods are relatively simple to implement; however, they are not flexible in adopting different design objectives. Furthermore, they cannot be readily integrated in optimization models for the synthesis of larger systems. However, graphical methods have also two important advantages: (1) they are intuitive and offer significant insights, and (2) they can readily represent “complex thermodynamics” graphically. On the other hand, while optimization methods are powerful and can be used to design large systems, they are computationally expensive due to the use of thermodynamic activity coefficient models (e.g. NRTL, UNIQUAC, and UNIFAC). Accordingly, the aim of this work is to develop an optimization model while exploiting the advantages of graphical methods. Specifically, we propose a mixed integer nonlinear programming (MINLP) model for a multistage liquid extractor design using some concepts from the modified McCabe-Thiele graphical method . In this graphical method, we use operating and equilibrium curves to determine the required number of equilibriums stages.Mathematically, we represent the operating curve as material balances and solubility equations and the equilibrium curve as tie-lines equations. Solubility and tie-lines equations are commonly obtained using complex thermodynamics equations or empiric correlations. However, we can approximate these equations using piece-wise linear functions thus avoiding high nonlinearities. Furthermore, we use a superstructure representation for the extractor where weintroduce binary variables to indicate the stages that are used in the optimal configuration. In the modified McCabe-Thiele method, we have a fixed solvent flow rate and need to determine the number of stages. With the proposed model, we can set this flow rate as a decision variable and have different objective functions, such as minimizing cost which depends on both the number of stages and the solvent flow rate. Moreover, the proposed model is extended to account for the optimal design of configurations with reflux, where the final extract stream is stripped of its solvent component and reenters the cascade as a reflux therefore allowing the extract product to be richer. For this configuration, additional binary variables are employed to choose the optimal feed location since it is no longer at the first stage. In addition, we consider some special cases for our model, such as dilute systems, systems with insoluble solvents, and non-ideal stages. Finally, we demonstrate the applicability of the proposed model through several illustrative cases.