The refining industry is facing increased volatility in the global market such as lower economic margins. In this situation, refinery production planning becomes a substantial tool to utilize all potentials to extend the economic margin to its maximum level. Classical refinery production planning models employ input-output relationships for different refinery units. However, these input-output relationships usually do not include decision making variables related to the units’ operation (pressure, temperature, etc.). In addition, they fail to capture the nonlinear nature of the refinery units and the inaccuracy caused by these models may reduce the overall profitability or compromise product quality. The refinery business involves tasks spanning from unloading crude oil at the refinery’s front end, generating utilities for refinery operations, blending different crudes and procuring crude blends to intermediate products, blending intermediates into final products and distributing refined products from refinery to distribution centers by means like pipelines. The traditional approach for the planning problem of these sub-systems is a sequential hierarchical approach. In reality, refinery processes are tightly interconnected and any attempt to solve this problem in a sequential push or pull manner may result in infeasible or suboptimal solutions. Eventually, the oil industry is subject to uncertainties such as unstable prices and unpredictable product demands. The existence of uncertainty in such parameters can seriously affect the optimization results and lead to inefficient operations. Nevertheless, there are challenges associated with the integrated optimization approach which arise from the difficulty of obtaining a solution with reasonable quality and time. The corresponding integrated model of the refinery operations with nonlinear unit correlations under uncertainty results in a computationally intensive nonconvex MINLP (Mixed Integer Nonlinear Programming) model due to the presence of bilinear, signomial, exponential or logarithmic terms in some of the mass balance, product yield and quality constraints. A further complexity is introduced in attempting to obtain a globally optimal solution. The standard methods for solving this complex refinery wide optimization problem may fail to converge to a solution or lead to sub-optimal solutions. Despite these complexities, this research presents a novel optimization heuristic for the integrated supply chain model of a refinery network. The production planning model of the refinery employs two different nonlinear models for the distillation unit utilizing correlations based on the Geddes fractionation index (FI) and a swing-cut (SW) data-based model and empirical nonlinear models for the remaining refinery units. The main contribution of this research is the development of an aggregation /disaggregation methodology based on lumped variable linearization (LVL) and normalized multiparametric disaggregation (NMDT) techniques through a two-level optimization algorithm and obtaining ε-global optimal solutions. To model the uncertainties, three modeling schemes such as robust programming, fuzzy possibilistic programming and two-stage stochastic programming with financial risk management have been utilized. The results indicated the effectiveness of the nonlinear models in terms of improving the overall refinery profit margin over a linear input-output model. In addition, they substantiated the huge advantage of a nonlinear distillation model utilizing correlations based on the Geddes fractionation index (FI) over the swing-cut (SW) data-based model. The results on the integration of the refinery supply chain echelons through a deterministic refinery model reinforced that the integrated model when solved by the proposed heuristic can provide significantly better solutions than its sequential counterpart in terms of both economic and operational objectives. Eventually, the results of the stochastic refinery model represented significant economical and operational differences between the outcomes of three uncertainty quantification methodologies discussed earlier while highlighting an obvious advantage of the robust optimization scheme.