ABSTRACT

The process of crude distillation impacts the economy of any refinery in a considerable manner. Therefore, it is necessary to improve it taking good advantage of the available infrastructure, generating products that conform to the specifications without violating the equipment operating constraints or plant restrictions at industrial units. The objective of this paper is to present the development of an optimization model for a Crude Distillation Unit (CDU) system at a ECOPETROL S.A. refinery in Barrancabermeja, involving the typical restrictions (flow according to pipeline capacity, pumps, distillation columns, etc) and a restriction that has not been included in bibliographic reports for this type of models: the heat integration of streams from Atmospheric Distillation Towers (ADTs) and Vacuum Distillation Towers (VDT) with the heat exchanger networks for crude pre-heating. On the other hand, ADTs were modeled with Metamodels in function of column temperatures and pressures, pumparounds flows and return temperatures, stripping steam flows, Jet EBP ASTM D-86 and Diesel EBP ASTM D-86. Pre-heating trains were modeled with mass and energy balances, and design equation of each heat exchanger. The optimization model is NLP, maximizing the system profit. This model was implemented in GAMSide 22,2 using the CONOPT solver and it found new operating points with better economic results than those obtained with the normal operation in the real plants. It predicted optimum operation conditions of 3 ADTs for constant composition crude and calculated the yields and properties of atmospheric products, additional to temperatures and duties of 27 Crude Oil exchangers.

Keywords: crude distillation units system, CDU, optimization, NLP, metamodels, operational constraints, energy constraints (restrictions), heat exchanger network.

RESUMEN

La destilación de crudo es un proceso que impacta considerablemente la economía de cualquier refinería y por tanto es necesario mejorarlo aprovechando la infraestructura instalada, generando productos dentro de especificaciones, sin violar las diferentes ventanas operativas o restricciones de planta que existen en las unidades industriales. El objetivo de este artículo es presentar el desarrollo de un modelo de optimización de un Sistema de Unidades de Destilación de Crudo (UDCs) perteneciente a la refinería de Barrancabermeja de ECOPETROL S.A. que involucra las restricciones de planta típicas (flujo por capacidad de tuberías, bombas, torres de destilación, etc.) y una restricción aún no contemplada en la literatura para modelos de este tipo: la integración energética de las corrientes de las Torres de Destilación Atmosféricas (TDAs) y las Torres de Destilación al Vacío (TDVs) con la red de intercambiadores de calor que precalientan el crudo. Por otro lado, las TDAs fueron modeladas a través de Metamodelos en función de las temperaturas y presiones de las torres, los flujos y temperaturas de retorno de los pumparounds, los flujos de vapor de despojo, el PFE de la ASTM D-86 para el Jet y el Diesel. Los trenes de precalentamiento fueron modelados con balances de masa y energía, y la ecuación de diseño de cada intercambiador. El modelo de optimización es de tipo NLP, maximizando la utilidad del sistema. Este modelo se implementó en GAMSide 22.2 usando el solver CONOPT y predice nuevos puntos de operación óptimos con mejores resultados económicos que aquellos obtenidos con la operación normal en las plantas reales. El modelo calculó las condiciones de operación óptimas de 3 TDAs para un crudo de composición constante y calculó los rendimientos y propiedades de los productos atmosféricos, adicionalmente a las temperaturas y duties de 27 intercambiadores de crudo.

Palabras Clave: sistema de unidades de destilación de crudo, UDC, optimización, PNL, metamodelos, restricciones operacionales, restricción energética, red de intercambiadores de calor.

RESUMEN

A destilação de cru é um processo que impacta consideravelmente a economia de qualquer refinaria e portanto é necessário melhorálo aproveitando a infraestrutura instalada, gerando produtos dentro de especificações, sem violar as diferentes janelas operativas ou restrições de planta que existem nas unidades industriais. O objetivo deste artigo é apresentar o desenvolvimento de um modelo de otimização de um Sistema de Unidades de Destilação de Cru (UDCs) pertencente à refinaria de Barrancabermeja da ECOPETROL S.A. que envolve as restrições de planta típicas (fluxo por capacidade de tubulações, bombas, torres de destilação, etc.) e uma restrição ainda não contemplada na literatura para modelos deste tipo: a integração energética das correntes das Torres de Destilação Atmosféricas (TDAs) e as Torres de Destilação a Vácuo (TDVs) com a rede de intercambiadores de calor que pré-aquecem o cru. Por outro lado, as TDAs foram modeladas através de Metamodelos em função das temperaturas e pressões das torres, os fluxos e temperaturas de retorno dos pumparounds, os fluxos de vapor de despojo, o PFE da ASTM D-86 para o Jet e o Diesel. Os trens de pré-aquecimento foram modelados com balanços de massa e energia, e a equação de desenho de cada intercambiador. O modelo de otimização é de tipo NLP, maximizando a utilidade do sistema. Este modelo foi implementado em GAMSide 22.2 usando o solver CONOPT e prediz novos pontos de operação ótimos com melhores resultados econômicos que aqueles obtidos com a operação normal nas plantas reais. O modelo calculou as condições de operação ótimas de 3 TDAs para um cru de composição constante e calculou os rendimentos e propriedades dos produtos atmosféricos, adicionalmente às temperaturas e duties de 27 intercambiadores de cru.

Palavras Chave: sistema de unidades de destilação de cru, UDCotimização, PNL, metamodelos, restrições operacionais, restrição energética, rede de intercambiadores de calor.

INTRODUCTION

Crude Distillation Units are key process plants in a petroleum refinery as they produce intermediate streams that are used in downstream process units. Changes in these units have a great impact on product yield and quality and, therefore, it is recommended to operate these units at optimal conditions from technical and economical points of view; that means operating conditions such as temperatures, pressures and flows of the units that maximize their economic performance (increasing product yield), subject to their real physical restrictions and their design capabilities. Mathematical modelling has become very common to develop these optimization studies.

There are several strategies that can be used to perform mathematical optimization processes. Linear Programming (LP) is generally used in production planning and programming due to their simplicity. However, its precision has been the subject of discussion for decades due to the simplified linear formulation of non-linear processes that can lead to non real solutions. These linear models usually do not consider interactions among CDUs occurring when these units share pre-heating sections or exchange streams among them (Wenkai, Chi-Wai, & AnXue, 2005), (Zhang, Zhu, & Towler, 2000, 2001). One of the most important strengths of linear models is that they always converge in reasonable calculation times. Another optimization strategy consists in the utilization of rigorous models, based on thermodynamic principles linked to an optimizer (Basak, 2002).  Even though this type of optimization begins with a precise CDU representation, it exhibits serious problems for its practical implementation since it does not ensure convergence and the calculation time increases significantly for very complex models. Furthermore, synchronization problems might occur when the model and the optimizer do not use the same platform. In refinery planning, using rigorous process models entails serious complications because these models require long solution times and frequently hide critical characteristics and parameters in order to increase profits (Hartmann, 2001).

There are other optimization strategies combining LP models with rigorous models, as it is reported by Zhang (2000), who independently performs non-linear optimization of each process, such as CDUs, in order to minimize operation costs by using fixed information from the LP model results. Due to the deficiencies in the above mentioned strategy, Zhang (2000) proposes an alternative where process optimization results are compared to the results from LP model. If they match with each other and do not violate any equipment operating constraints, implementation of these results is considered feasible. If this is not the case, process level results are used to feedback and update production yields of units in the LP model. This is a cyclic procedure that might take long execution times when using optimization with rigorous process models. 

Furthermore, Zhang (2001) proposes another strategy in order to take advantage of synergy among processes. The idea applies again the LP model structure, although this time proposes the inclusion of non-linear process models within the global model through linearization techniques. As a consequence, representation capabi-lity of each unit within the global model improves. This concept was implemented in some commercial software such as Aspen-PIMS, thus creating an evolution of such programs.

The last strategy mentioned here was proposed by the Brazilian researchers (Pinto, Joly & Moro, 2000) who have shown their interest in the improvement of planning models by enhancing process representation using non-linear optimization. They propose the inclusion of process complexity within the model, in order to improve the accuracy of the results. In other words, they consider the non-linear representation of units such as CDU within the macromodel in order to take advantage of the interaction among processes and provide operational guidelines to the units through a simple but robust model.

This paper describes a new strategy based on metamodelling to optimize complex systems and presents the simultaneous optimization of three crude oil distillation towers with their heat integration (crude oil preheating trains) as an example of application. Metamodels are approximated models from rigorous simulation that can be used on a routine basis due to their quick implementation, execution, and convergence, surpassing the rigorous models (Palmer & Realff, 2002a and 2002b).  Crude oil preheating trains were represented by rigorous exchanger models. Metamodels were developed from previous PRO/II rigorous models of the industrial units. Metamodels of the distillation towers and rigorous models of the heat exchange trains were integrated in the GAMSide 22.2 - based optimization algorithm that maximizes the CDU system profits, finding optimum operation conditions for each atmospheric distillation tower.

METHODOLOGY

The system to optimize comprises three crude oil distillation units (CDU1, CDU2 and CDU3) and these crude units are composed by three atmospheric distillation towers with different design characteristics (ADT1, ADT2 and ADT3), three atmospheric furnaces (AF1, AF2 and AF3), two vacuum towers (VDT1 and VDT2) and two heat exchanger networks (NET1 and NET2). Figure 1 shows a simplified diagram of the CDUs studied in this work.

CDU1 and CDU2 share the network NET1 conformed by 17 exchangers and the NET2 is composed of 10 heat exchangers. The vacuum residue stream from CDU3 interacts with CDU2 through the NET1 network preheating the Crude Oil that feeds this unit.

PRO/II Model Construction 

A rigorous model using the commercial simulator PRO/II ® version 8.0 was constructed for each CDU. As thermodynamic packages, the Grayson-Streed and Braun K10 (BK10) equations of state were used in the atmospheric and vacuum towers respectively. Feedstock stream was characterized using the TBP curve, API gravity curves, sulfur content, neutralization number (NN), Conradson Carbon content (CCR) and lightends content (C₁, C₂, C₃, iC₄, nC₄, iC₅, nC₅, hexane, H₂O, CO₂, CO, N₂ and H₂). PRO/II models were tuned with the information from industrial runs. Regarding ADTs, the set of manipulated variables is conformed by the ope-ration variables (Crude flow, input and output furnace temperatures, temperature and pressure on the top of the column, pumparounds flows and return temperatures, stripping steam flows, condenser temperature and pressure). A quality control variable of the JET and DIESEL products is also included: End Boiling Point (EBP) ASTM D-86 distillation curve.

Construction of Atmospheric Tower Metamodels

For the purpose of this research, non-linear semi-rigorous models (Metamodels) of each atmospheric tower were constructed. The methodology used to develop the Metamodels includes:1) Selection of input variables, 2) Design of Experiments (Latin Hypercube), 3) Generation of data using PRO/II models, 4) Parame-trization of Metamodels using multivariate stepwise regression, 5) Validation and evaluation of their significance. Output variables in Metamodels were flows, temperatures, and properties of all the streams leaving the towers (products and pumparounds toward NET1 and NET2). The properties calculated by Metamodels were API gravity, ASTM D-86 curve (10, 30, 50, 70 and 90%) and calorific capacities. Pumparound duty Metamodels of ADT1, ADT2 and ADT3 were created, which were linked to the NET1 and NET2.

The input variables set is composed for the PRO/II manipulated variables and changes depending on the design of the tower (see Table 1).

The data collection plan or Design of Experiments (DoE) used was the sampling process by Latin Hypercube (LHS: Latin Hypercube Sample) which provides unique values for each point and exhibits better dispersion than other sampling procedures such as the random and grid sampling (Palmer & Realff, 2002a and 2002b).

The Metamodels were a second-order polynomial function with binary interaction:

Where ŷ(X) was the output variable vector (flows, temperatures, properties, etc.), X is the normalized input variable vector,  (feed flow, operating conditions, and quality properties) with i=1,..,13 for ADT1, i=1,..,138, for ADT2 and i=1,..,17 for ADT3. Each output varia-ble has a Metamodel with characteristic parameters according to its behavior in regards to the input varia-bles. For ADT1, 105 parameters were calculated; for ADT2, 190 parameters and for ADT3, 171 parameters were found.

Furthermore, 5000 sampling points were used in Metamodels construction (N=5000). These points represented 5000 simulation runs for each distillation tower (ADT1, ADT2 and ADT3) using the configured and tuned rigorous models in PRO/II. Due to the random nature of LHS for each independent variable and to the complexity of the distillation process, 100% convergence of the model is not attained. Therefore, only converged data was used in Metamodels construction.

Parameters for each Metamodel were estimated based on the matrix of experimentation results in PRO/II, using multivariate stepwise regression. During regression, input variables were normalized between -1 and 1, and output variable results (ŷ(x) were in the real magnitude corresponding to each variable.

A typical residual analysis, an escalated residual analysis (standardized and studentized) and a PRESS analysis were used for the detection outliers within the Metamodel construction set. In order to measure the Metamodel prediction capability, the R²_Prediction (Montgomery, 2001) was calculated. Each metamodel was evaluated by a validation set with 1.000 simulated points in PRO/II to verify their prediction quality before using them in the optimization model. Figure 2a shows the correlations between calculated volumetric flow for the ADT2 products with Metamodels and the results found in PRO/II. Figure 2b shows the error distribution of product flow Metamodels for ADT1 and ADT2.

Atmospheric tower Metamodels showed a maximum average error of 2,5% with a standard deviation of 0,41% for flows, while for the rest of the properties the maximum average error was 2,0% with standard deviation of 0,75%.

Construction of Heat Exchanger Integrated Networks

NET1 and NET2 of the CDU system were modeled through energy balances by Equation 2 on the hot side, by Equation 3 on the cold side and by the design equation expressed in Equation 4.

Q_hx,u is the duty of the hx exchanger belonging to the set of exchangers corresponding to the CDU, u belon-ging to the U set (CDU1, CDU2 and CDU3).  F_h,hx,u is the mass flow of hot stream (h), entering the hx exchanger, of the unit u. Th_in,hx,u is the input temperature (in) of h that enters the hx exchanger of the unit u. Th_{_out,hx,u} is the output temperature (out) of h that enters the hx exchanger of unit u. The same naming system is used for the cold fluid (c). Uhx,u y Ahx,u are global heat transfer coefficient (BTU/day ft² °F) and transfer area in ft² (see Table 3) of exchanger hx belonging to unit u.

The values for h and c in Equation 5 were determined by correlations (API Technical Data Book, 2005) in function of mean input temperature to the exchanger of c and h () and the constants A₁, A₂ and A₃ (functions of the ASTM D-86 curve, Watson factor and Specific Gravity of each stream).

The path flow of each fluid through the whole set of exchangers was considered in network configuration, including recycles and bypass from real plants.

Determination of CDUs Operating Constraints

There are different operating constraints within an industrial plant that ensure equipment integrity. In the optimization model, flow restrictions were included based on pump capacity, input and output temperature restrictions of exchangers to ensure their useful life, and restrictions at furnace duties (Tables 2 and 3). As product quality restrictions for JET and Diesel, maximum EBP values of 570°F and 760°F were considered respectively.

Identification of the Optimization Model Objective

The objective of the optimization model was maximize the CDU system profit (Pr) that considers distillation towers, furnaces, and heat integration as operating restriction for a Crude feedstock with constant composition. Pr of the CDU system is defined as the income due to the sale of products, minus the raw material costs, minus the operating cost for each CDU:

The income due to product sale is defined in Pa for all products. This is the set of atmospheric pro-ducts (Naphtha, Jet, Diesel, AGO, and Reduced Crude -REDC) from all the atmospheric towers Ua (ADT1, ADT2 and ADT3). The raw material cost is defined for all crudes entering to each atmospheric tower Ua. The operation cost is defined considering all the resources consumed in the process, R, such as fuel, compressed air, fresh water, industrial water, steam and electricity.

The stream price data and the operation costs in dollars per barrel (US$/bl) are listed in Tables 4 and 5, respectively. The manipulated variables in maximization were the Metamodel input variables shown in Table 1.

Selection of Solver in GAMS and Initial Point

Non-linear models created with GAMSide version 22,2 (Brooke, Kendrick, & Meeraus,1992) must be solved with Non-Linear Programming algorithms (NLP). This research used the CONOPT solver since it correctly manages correctly models with many non-linear restrictions and a moderate number of degrees of freedom. These are characteristics of the problem analyzed here. The Initial Estimation (IE) was a set of operating conditions close to an operation of the industrial unit. Therefore, an operation with low yield of REDC and maximum Jet or Diesel production were selected.

Running the Optimization Model and Results

Once the objective function and the model restrictions were configured, and the solver in GAMS was specified, the next step was to solve the optimization problem. The results of the model included the profit of the system and of each CDU. Regarding ADT1, ADT2 y ADT3, they included: operating variables, Crude flow, products yields, flow, temperatures and properties of both products and pumparounds. The following is included for NET1 and NET2: Duties, temperatures and properties of Crude and hot streams at the entry and exit of each exchanger. The furnace duty was calculated.

Simulation of the Optimum Point Generated by GAMS

An algorithm is designed to find the feasible optimum point (Figure 3). This algorithm contemplates n simulations in PRO/II, including the first run to find the IE of optimizer and n-1 optimization runs in the GAMS model that generates the operational and quality conditions called the Optimum Feasible Point.

The algorithm in Figure 3 starts from available runs plants in CDUs and/or converged points of individual models in PRO/II of ADT1, ADT2 and ADT3 such as IE, which enters the CDU metamodel optimizer in GAMS. Next the run is completed and the Optimum Variables (V_opt) are obtained; and they are compared to the Initial Estimation Variables (V_IE). The cycle finishes when these values are close enough to each other (<TOL).

Optimum Point Sensitivity Analysis

Due to the fact that it is a NLP optimization, the new operational point found by GAMS is a local optimum point. Several IE were tested in order to figure out the contour of the objective function. The optimum points were compared and the best one was selected as the most profitable of the system.

MATHEMATICAL MODEL OPTIMIZATION

An optimization of the model, maximizing Pr in Equation 6, was developed. This model was subject to semi-rigorous process models (Metamodels), restrictions such as product specifications, operational variable boundaries, limits of other equipments and the energy of the pre-heating crude trains. The objective function of Equation 6 is subjected to mass balances on each tower such as in Equation 7, to upper and lower limits of feed flow to the towers of Equation 8, to boundaries of operation and quality variables of Equations 9 and 10, to restrictions of products and pumparounds flow according to pump capacity of Equations 11 and 12, to the maximum operation limits of duties of the atmospheric furnaces in Equation 13, to the maximum limit of the fluid temperature h entering the exchangers in Equation 14 and the c fluid coming out the exchangers in Equation 15, to the set of equations representing the networks NET 1 and NET2 of Equations 2, 3, 4 and 5 and, of course, to the process model using Metamodels.

The L and U superscripts indicate the lower and upper level that can take the flow or involved variables. The Operation Variables (V_opt) of the atmospheric tower t, are independent variables inside PRO/II and their amount depends on tower design. Quality variables (Vq_s) included are the EBP of ASTM D-86 for Jet and Diesel. Q_f,t is the furnace duty, f, belonging to t.

The ADT typical model is represented by Equations 16 and 17.

-Total flow of each product stream:

-Property of each product stream:

The function Øs,t  is the Metamodel that calculates flow of product s, of tower t. The function Χ_p,s,t  is the Metamodel that calculated the property p of the products of the tower t. API gravity, temperatures of ASTM D-86 (10, 30, 50, 70, 90%) and Cp are included among the p properties. All Metamodels follow the form of Equation (1).

Besides the typical output variables for an ADT (Equations 16 and 17), Metamodels were constructed to predict the stream temperature leaving the towers, using Equation 18 for both products and pumparounds (alternate sub indexes a/s).  Furthermore, Metamodels were constructed for the pumparound duties using Equation 19 which linked to the networks NET1 and NET2 through Equation 20, depending on the design of each CDU. Equation 21 is the Metamodel of the atmospheric furnaces duty.

The function Η _s/pa,t  is the Metamodel of temperature of flows (s if it is a product and pa if it is a pumparound) exiting tower t. The function μ_pa,t  is the Metamodel of the pumparound duty pa, either UPA, MPA or LPA belonging to the tower t. In Equation (20), the addition of the network duties is conducted on the set of exchangers belonging to the HXpa pumparounds. The function ξ_f,t  is the Metamodel that calculates the atmospheric furnace duty f  of tower t.

RESULTS

The optimization model was formulated in GAMS (Brooke, 1992) containing 3.881 variables, 3.834 equations and 9.833 non-linearity, non zeroes. This model was solved with CONOPT for NLP with a total amount of time between 16 and 32 seconds, depending on the closeness of the initial estimation to the optimum value. The optimization results were compared to the current CDU situation, where operational conditions are determined based on plant engineering heuristics, according to the type of crude to be processed. Considering the operational and quality restrictions as well as the energy exchange within and between CDUs, and the initial estimation, the algorithm could predict a new operational point with better economic results than those obtained with the normal operation. The summary of the main optimum operational conditions of the CDU system for two initial estimations are shown in Table 6.

The product yields are listed in Table 7. The initial estimation 1 (IE-1) is a typical operational point of real plant, and the initial estimation 2 (IE-2) is a converged simulation point in PRO/II. The results of Tables 6-8 were obtained after applying the algorithm described in Figure 2. Therefore, these constitute feasible optimum points (with convergence in PRO/II).

Table 8 shows the system profit (Pr) before and after optimization. The optimum point from EI-1 shows grea-ter profit than the value found by IE-2, with an increase of 4,4%. The optimum results from EI-1 of Tables 6-8 represent an increment in the system profit of about two hundred million dollars (M US$ 200) per year without modifications in installed infrastructure, only changing operational conditions and Crude Oil Total Flow. 

Figure 4 shows the ASTM-D86 distillation curves found in GAMS and validated in PRO/II of some atmospheric products for the CDU system, later optimized by IE-1. These results are evidence of good Metamodels prediction within the optimization model.

CONCLUSIONS

• A new NLP optimization model was proposed for a CDU system that includes the energy restriction of plants and utilizes the metamodel approach to represent the non-linear phenomenon of distillation. The NLP model maximizes the system profit finding the optimum operational conditions for each atmos-pheric tower, calculating products yields and their properties, temperatures and duties of exchangers responsible for crude pre-heating.
• The proposed model included satisfactorily the heat integration of industrial plants as a process restriction without convergence problems by the Metamodels, obtaining feasible optimums in PRO/II spending less time finding the solution.
• The use of Metamodels in this optimization proves the qualities of this modeling since they represented the rigorous distillation phenomenon through a simple function (second degree polynomial), without the convergence and coupling problems characteristic of rigorous simulation models offered by commercial packages. These benefits place them as candidates to be the predictive tools within more complex models such as refinery planning that require the representation of plant reality, thus overcoming the simplified predictions from the LP models used on a routine basis.
• The execution time used by GAMSide to find the solution was between 16 and 32 seconds. This is a low value for a problem modeling plant behavior in a semi-rigorous manner.
• The existence of exchangers in the model of CDU system moves the optimum value to different operational regions from the conditions found if only towers were used, since they add new restrictions to the problem, ensuring safe operation within industrial plants.
• The profit increment of the system without investment and modifications in infrastructure contemplated in this study, demonstrates the importance of optimization models in industrial plants. However, optimization still has the feedstock fixed (constant composition) to the units, an aspect that determines the plant operation. Therefore, the feedstock composition is suggested to be included within operational optimization with heat integration as an optimization variable within an unique model.

ACKNOWLEDGMENTS

The authors express their gratitude to ECOPETROL -ICP, which provided the necessary information for the development of this study and for the construction and tuning of PRO/II models.

REFERENCES

API Technical Data Book (2005). Vol. II: Thermal Properties & Phase Equilibria. 7th edition, Epcon International        [ Links ]

Basak, K. & Abhilash, K. S. (2002). On-Line Optimization of a Crude Distillation Unit with Constraints on Product Properties. Ind. Eng. Chem. Res., 41: 1557-1568        [ Links ]

Brooke, A., Kendrick, D. & Meeraus, A. (1992). GAMS - a user’s guide (release 2.25). San Francisco, CA: The Scientific Press.        [ Links ]

Hartmann, J. C. M. (2001). Determine the optimum crude intake level -a case history. Hydrocarbon Processing, 6: 77-84        [ Links ]

Montgomery, Douglas C. (2001). Design and Analysis of Experiments (5^th  ed.). New York: John Wiley & Sons.        [ Links ]

Palmer, K. & Realff, M. (2002a). Metamodeling Approach to Optimization of Steady State Flowsheet Simulations. Model Generation. Trans IchemE., 80 (A), 760 - 772.        [ Links ]

Palmer, K. & Realff, M. (2002b) Metamodeling Approach to Optimization of Steady State Flowsheet Simulations. Trans IchemE.,  80 (A), 773 - 782.        [ Links ]

Pinto, J.M., Joly, M. & Moro, LFL. (2000). Planning and Scheduling Models for Refinery Operations. Comput. Chem. Eng., 24: 2259-2276        [ Links ]

Wenkai, Li & Chi-Wai., Hui & AnXue, Li. (2005) Integrating CDU, FCC and product blending models into refinery planning. Comput. Chem. Eng., 29 (9), 2010-2028.         [ Links ]

Zhang, J., Zhu, X. X. & Towler, G. P. (2000). A novel mode-ling and decomposition strategy for overall refinery optimization. Comp. Chem. Eng.,  24: 1543-1548.         [ Links ]

Zhang, J., Zhu, X. X. & Towler, G. P. A (2001). Simultaneous Optimization Strategy for Overall Integration in Refinery Planning.  Ind. Eng. Chem. Res.,  40: 2640-2653.        [ Links ]