A Strategy Based on Convex Relaxation for Solving the Oil Refinery

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

A STRATEGY BASED ON CONVEX RELAXATION FOR SOLVING THE OIL REFINERY OPERATIONS PLANNING PROBLEM Tiago Andrade, Gabriela Pinto Ribas, and Fabrício Oliveira Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.5b01132 • Publication Date (Web): 22 Dec 2015 Downloaded from http://pubs.acs.org on December 26, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A STRATEGY BASED ON CONVEX RELAXATION FOR SOLVING THE OIL REFINERY OPERATIONS PLANNING PROBLEM Tiago Andrade, Gabriela Ribas, Fabricio Oliveira1 Industrial Engineering Department Pontifícia Universidade Católica do Rio de Janeiro – PUC-Rio

Abstract Oil refining is one of the most complex activities in the chemical industry due to its nonlinear nature, which ruins the convexity properties and prevents any guarantees of the global optimality of solutions obtained by local nonlinear programming (NLP) algorithms. Moreover, using global optimization algorithms is often not feasible because they typically require large computational efforts. This paper proposes the use of convex relaxations based on McCormick envelopes for the Refinery Operations Planning Problem (ROPP) that can be used to generate both initial solutions for the ROPP and to estimate optimality gaps for the solutions obtained. The results obtained using data from a real-world refinery suggest that the proposed approach can provide reasonably good solutions for the ROPP, even for cases in which there was no solution available using traditional local NLP algorithms. Furthermore, compared with a global optimization solver, the proposed heuristic is capable of obtaining better solutions in less computational time.

Keywords: Refinery operations planning, Nonconvex optimization, Convex relaxation, McCormick envelopes

1

Corresponding author ([email protected])

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction In recent decades, there has been an intensive increase in global competition in the chemical processing industry (1). To maintain profitability, many new technologies for process operations have been developed, generating a demand for optimized planning (2). In this context, mathematical programming optimization models arose as a useful tool to support refinery-planning decisions. Certain of the pioneering software systems for refinery planning, such as Refinery and Petrochemical Modeling System – RPMS (3) and Process Industry Modeling System – PIMS (4), were based on linear programming (LP) models in their earlier versions. An advantage of this approach is that these models can be efficiently solved. However, the plans generated by these LP models only serve as approximations and, among other limitations, do not allow the consideration of nonlinear mixing properties or more complex process models (2). Although certain of the available tools, such as PIMS, have been updated to address nonlinearities using successive linear programming, most of these tools apply methods that approximate the solutions of nonlinear programming (NLP) models by solving multiple instances of LP models, which can lead to solutions that are far from the global optimum. This study examined a nonconvex nonlinear mathematical model for the Refinery Operations Planning Problem (ROPP) and the difficulties inherent to solve it.The ROPP often considers decisions related to the types and quantities of crude oils to acquire, the intermediate stream levels between the processing units (including the distillation tower), the quantities of final products to be sold, the product specifications (specific product properties, such as viscosity and sulfur concentration), and the storage volumes. Making these decisions in an optimal fashion is a highly complex activity; therefore, a mathematical model becomes an important tool to support the process of building an optimal (or even a feasible) plan, given the decisions of the upper planning levels, the refinery units’ capacities and the product quality level requirements that must be met. In the 1950s, Symonds (5) and Garvin et al. (6) proposed LP models for the ROPP. With the development of solution methods for more complex problems and the constant increases in computational power, more realistic models could be formulated to represent the nonlinearities that are inherent to the chemical processes that occur inside the refinery. Several examples of NLP models that were developed to solve the ROPP can be found in the literature. Moro et al. (2) proposed a Mixed-Integer Nonlinear Programming (MINLP) model for the entire refinery topology. Joly et al. (1) extended that model to consider refinery operations scheduling as well as decisions about the liquefied petroleum gas (LPG) area. Neiro and Pinto (7) proposed an MINLP model for the operational planning of the petroleum supply chains and later expanded it to consider the integrated planning of multiple refineries (8). Alhajri et al. (9) also proposed a model


Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that considers the nonlinear nature of refinery processes, such as blending, by means of an NLP. This model also considered the operating parameters of the crude distillation unit (e.g., cut-point temperatures) as decision variables. Allatas et al. (10) proposed a model for the ROPP that has nonlinear equations associated with the crude distillation unit that are based on the fraction index model instead of the traditional fixed yields. This model also considered the cut-point temperature as a decision variable. Guerra and Le Roux (11) proposed an NLP model for the ROPP using empirical models for the crude distillation units (CDUs) and fluid catalytic cracking (FCC) units. The same authors performed a case study considering small and medium-scale refineries (12). More recently, Kelly et al. (13) proposed a distillation optimization model that considers not only blending but also cut-point temperature as decisions variables using monotonic interpolation. Menezes et al. (14) presented a single period nonlinear model for the ROPP to analyze long-term predictions of the fuel market. Most of the studies that address NLP models recognize that one of the primary difficulties in solving these problems is that the available NLP algorithms can produce infeasible or local optimal solutions. It is well-known that the solution of an NLP problem can be highly sensitive to the starting point provided (i.e., the initial solution provided to the NLP algorithm), and being able to explore the problem structure to determine good starting points can therefore mean the difference between finding a good local solution and not being able to obtain a feasible solution at all. A global optimization method could be applied to overcome the limitations of NLP algorithms, but the solution of NLP problems using global optimization methods requires much more computational effort than using NLP algorithms. In this paper, we present a nonconvex NLP model for the ROPP that is based on the aforementioned literature. Because global optimization algorithms often lead to large computational times and local NLP algorithms may lead to infeasible or local optimal solutions, we propose a strategy to solve the ROPP that aims at improving the quality of the solutions obtained. The strategy relies on a convex relaxation of the original problem that is obtained using McCormick envelopes (15), a randomized multi-start heuristic using off-the-shelf NLP solvers, and the use of starting points that are in the neighborhood of the optimal solution of the convex relaxation. Karuppiah and Grossmann (16) used similar relaxation techniques for the problem of the synthesis of integrated water systems in chemical processes. Floudas and Gounaris (17) provided an extensive review concerning the use of convexification techniques for deriving convex relaxations of NLP problems in the chemical industry. More recently, Kolodziej et al. (18) presented the derivation of a global optimization algorithm for MINLP problems with bilinear terms in their constraints.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

This paper has three primary objectives. First, we propose a strategy to solve the ROPP using offthe-shelf NLP solvers that is capable of obtaining improved local optimal solutions. Second, we show that we are able to obtain valid upper bounds to the optimal objective function value of the ROPP, which means that we are capable of calculating optimality gaps for the solutions obtained by the proposed strategy. Last, we apply the proposed strategy by means of a randomized multi-start heuristic to several instances of the ROPP and perform computational experiments using popular NLP solvers, namely, CONOPT, MINOS, and KNITRO, aiming to validate the approach. We also compare the performance of the proposed strategy with the global solver BARON, both in terms of computational time and solution quality. The contributions of this paper are twofold. First, we present a convex relaxation for the ROPP based on McCormick envelopes that can be used to generate the initial points for local NLP solvers. Second, we demonstrate how this relaxation can be used in a randomized multi-start heuristic that relies on the quality of this initial point to generate good local optima for the ROPP using perturbations of the relaxation’s optimal solution. To the best of our knowledge, this report is the first to provide numerical evidence supporting the effectiveness of the combination of these ideas applied to ROPP. This paper is structured as follows. Section 2 describes the ROPP and its main premises. In Section 3, the mathematical model used to represent the ROPP is described. In Section 4, we provide the details of the proposed solution strategy. In Section 5, we describe the computational experiments that we performed and present the results obtained. Last, we present several conclusions and suggestions for future research in Section 6. 2. Problem Description Decision levels in a typical oil supply chain can be divided in multiple ways. Highly often, the decisions are divided into three categories: upstream, midstream, and downstream (19). Upstream segment decisions involve exploration and oil production. Midstream decisions concern oil storage, transportation from the production sites to the refineries to be processed, and the refining activity. The downstream segment includes the transport of the final products to the consumers. These decision levels can also be classified into strategic, tactical, and operational. The differences between these planning levels are primarily related to the time horizon and the level of detail considered in each decision process. Strategic planning usually comprises aggregate long-term decisions, typically comprising several years within the planning horizon. Tactical planning is an intermediate level, considering time horizons typically ranging from few months to a few years, with


Page 4 of 32

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


certain of the decisions being represented at a greater level of detail. Last, operational planning involves decisions at the highest level of detail, and the time horizon considered is usually shorter, typically from one to three weeks (20). The problem considered in this study can be classified as an operational planning model for an oil refinery (ROPP - Refinery Operations Planning Problem), which is part of the midstream segment. The main goal of an oil refinery is to transform crude oil into useful products with higher aggregate value. The refinery topology includes process units, which can be classified as follows: separation, conversion, and treatment units; storage tanks for final or intermediary products or crude oil; and product streams between units. Figure 1 shows a schematic representation of a generic oil refinery, in which the boxes represent units and the arrows represent (intermediate) product flows.

Figure 1 - Refinery schematic flowchart Such processes as crude distillation and vacuum distillation separate the oil into simpler fractions through a physical process. The separation units have different process modes (i.e., different settings of reaction temperature and pressure levels) that can be chosen to change the proportions between the separated fractions. Conversion processes, such as fluid catalytic cracking, convert a fraction or a set of fractions into others, changing its molecular structure through chemical processes. Last, treatment processes, such as hydro treatment, remove contaminants.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The refinery operational planning consists of making decisions concerning the amounts to acquire of each crude oil, the flows between units, the product-processing levels in each process unit, and the inventory levels for intermediate and final products along a given planning horizon. It is also important to determine the blends of the intermediary streams and the process modes of the separation units, as well as selecting which final products will be sold in the market. The blending decision combined with the choice of process mode for each process unit defines the yield and properties of all refined products. The process mode specifies the operational variables, such as the temperature and pressure for the process unit. Thus, the combination of different modes allows the refinery to increase the yield of some desired product. The ROPP also calculates the quantities and properties of all intermediate streams to precisely define the quality of the final products that are subject to severe specifications. A series of nonlinear equations must be modeled to obtain the correct property value for the intermediate and final products. All of these decisions must be made while seeking to maximize the profit obtained from the oil refining.

3. Mathematical Model In this section, we present the mathematical model for the ROPP, which is based on the model originally proposed by Neiro and Pinto (8). However, in our model the property balance constraints were modeled explicitly by using bilinear functions instead of more generic function mass balance constraints are considered independently from inlet stream properties; and the units may have more than one process mode. The process mode that a unit will operate is not considered as discrete decision since it is possible to utilize multiple process modes during the planning horizon. The objective function to be maximized is the profit, which is given by the difference between the revenue, obtained by selling the final products from the shipping points (subset ), and the costs of

acquiring oil from the source points (subset ). The decision variables are the stream flows, the

inventory levels in the storage tanks (subset ) and the stream property values, such as viscosity or

density. These streams can be crude oils, intermediary products, or final products, such as gasoline.

One important decision is the oil blend that will be processed by the distillation tower, a process unit (subset ) of major importance within the refinery and in which different process modes (set )

can be used. The planning horizon is uniformly divided into discrete time periods, and all decisions must be made for each of these time periods. Regarding the constraints, there are several of different natures. There are market constraints, limiting how much oil can be bought and how much of each product can be sold. Other classes of constraints are the volume balances in the units and the property balances, which are typically


Page 6 of 32

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonlinear because of the mixture of different streams. Several constraints have their origin in the

units’ capacities for processing or storage. In addition, certain properties (set ) must be met in the final products or even in the intermediary ones for multiple reasons, such as federal regulations or the proper functioning of the refinery units and chemical processes. The objective function (1) defines the profit of the refinery operation, which is given by the difference between the revenue and the costs. The revenue comes from selling the final products and the costs of acquiring oil. Maximize OF = ,

,,

− , ,,

∈ ∈ ∈, ∈

(1)

∈ ∈ ∈, ∈

Constraints (2) to (5) represent material balances. Equations (2) and (3) are definitions of auxiliary

variables ,,

and ,,

, which represent the total stream volumes that are entering and leaving

the units, respectively. Constraint (4) represents the mass balance in the storage tanks, stating that the

flows that enter the unit (∑ ∈, ,,

) plus the volume in storage at the beginning of time period , must be equal to the stream flow that leaves the unit (∑ ∈

,

,, ) plus the volume in

storage at the end of period , . Equation (5) describes the mass balance of the blending units !,

which is similar to Equation (4), but without the storage alternative. ,,

= ,,

=

#, #, ,,

,, , # , #

$#, #, ,,%∈&

$,, ,0,0%∈&

, = 21., + ,

+ ,,

− ,,

∈,

0∈,

,,

= ,,

∈,

∈,

∀( ∈ ∪ * ∪ , ∀+ ∈ , ∀, ∈ -., , ∀/ ∈ * ∀( ∈ ∪ * ∪ , ∀+ ∈ , ∀, ∈ -1, ∀/ ∈ *

(2)

(3)

∀( ∈ , ∀+ ∈ , ∀, ∈ -., , ∀/ ∈ *

(4)

∀( ∈ !, ∀+ ∈ , ∀/ ∈ *

(5)

Constraint (6) defines the material balance in process unit ( ∈ . In these units, every inlet stream

generates multiples outlet streams, and the proportion is given by the yield parameter 4,, , 0 .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀/ ∈ *

,,

= ,, 0 4,, , #

# ∈,

Page 8 of 32

(6)

The nonlinear constraint used to calculate the property value of all intermediate streams and final products are defined in Equations (7) to (10), where the multiplication of the variables flow and

property value is needed. Equation (7) calculates the level of property of outlet stream , in process unit ( ∈ . Equation (7) is similar to Equation (6), with a property variable ,, ,5 on the left side and a parameter with a known property 1,, ,

# ,5 on the right side. ,, ,5 ,,

= ,, 0 4,, , 0 1,, ,

# ,5

# ∈,

(7)

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1,, 0 ∀/ ∈ *

Equation (8) computes the property value for outlet stream , in tank unit ( ∈ . The outlet property

,, ,5 depends on the property of the inlet stream ,, ,5 and the property level of the amount

stored at the end of the previous time period. ,, ,5

621.,

+

,

+

# ∈,

,, 0 7

(8)

= 1.,,5 21., + , , + ,, 0,5 ,, 0

# ∈,

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1,, , ∀/ ∈ *

Constraint (9) computes outlet stream property ,, ,5 using the inlet stream property ,,

# ,5 in

the blending tanks ( ∈ !.

,, ,5 ,,

,,

# = ′ ,5 ,, ′

# ∈,

(9)

# ∈,

∀( ∈ !, ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1,, ∩ .,, # , ∀/ ∈ *

Equation (10) computes the inlet flow property of the tank and process units using the outlet flow

properties that leave other units and enter into unit ( ∈ .


Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


,, ,5 ,,

=

$#, # , ,,%∈&

#, #, ,, #, #, ,5

(10)

∀( ∈ ∪ *, ∀+ ∈ , ∀, ∈ -., , ∀ ∈ .,, , ∀/ ∈ *

Constraints (11) to (18) define the known bounds for the decision variables. As stated in (19), all flow variables, including the auxiliary variables and the inventory level variables, are nonnegative. The property variables are defined as real numbers. ,,

≤ :1,,

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀/ ∈ *

∀( ∈ , ∀+ ∈ , ∀, ∈ -., , ∀/ ∈ *

,< , ;!,,

≤ ,, ≤ ;!,,

∀( ∈ , ∀+ ∈ , ∀/ ∈ *

,< , 21, ≤ , ≤ 21,

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀/ ∈ *

,< , 1,, ,5 ≤ ,, ,5 ≤ 1,, ,5 ,< , .,, ,5 ≤ ,, ,5 ≤ .,, ,5

,< , :,, , # , # ≤ ,, ,# , # ≤ :,, , # # ,< , :.,,

≤ ,,

≤ :.,,

,< , :1,,

≤ ,,

≤ :1,,

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀/ ∈ * ∀(, +, ,, ( 0 , + 0 ∈ =

∀( ∈ , ∀+ ∈ , ∀, ∈ -., , ∀/ ∈ *

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀/ ∈ *

? ,, , # , # , ,, , ,, , , ∈ ℜ , ,,5, , ,, ,5 ∈ ℜ

(11) (12) (13) (14) (15) (16) (17) (18) (19)

4. Proposed Strategy The model stated in the previous section is an NLP model that is nonconvex due to the bilinear terms present in Constraints (7) to (10). This class of problem is well-known for presenting challenges, especially in regard to applying those models to real-world problems. First, these problems typically present multiple local optima, i.e., a solution that is the best in a given neighborhood is usually not the best solution available over the complete feasible region. Second, it is often difficult to even obtain a feasible solution using current off-the-shelf solvers for this class of problem. In addition, the algorithms to solve general nonconvex problems have non-polynomial complexity, implying a large computational time for typical real-world (i.e., large) problems. Nevertheless, the aforementioned difficulties can be eased if one can provide a good initial solution as an input for the selected NLP algorithm. In this sense, we propose a strategy that relies on a convex relaxation of the feasible space to provide such an initial solution. The relaxation is performed using McCormick envelopes (15), which can be described as a collection of supporting hyperplanes for the nonconvex feasible region. The relaxation



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

obtained using the McCormick envelopes is an LP model, and therefore convex, for which there are known efficient algorithms to solve it.

4.1 Convex relaxation for the ROPP Constraints (7) to (10) in the model presented in Section 2 are responsible for the nonlinearity (and nonconvexity) of the problem. These constraints enforce the product property balances, which are defined by means of bilinear terms. Note that Constraints (7) to (10) are nonlinear due to the bilinear terms; therefore, McCormick envelopes can be formulated for that set of constraints to obtain a convex relaxation for the proposed model. The convex relaxation for model (1-19) using McCormick envelopes can be obtained as follows. For more detailed information regarding how to obtain the envelopes, please, refer to the appendix A. Equation (20) is the relaxed version of Constraint (7), in which the product between the outlet

and outlet flow ,,

is replaced by the new variable ,, ,5 . property ,, ,5 ,, ,5 = ,,

# 4,, , # 1,, , # ,5

# ∈,

(20)

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1(,+,,′ ∀/ ∈ *

Constraint (21) presents the relaxed version of equation (8). The product between the variables

, ,, ,5 and , was replaced by the variable ,, ,5 . Time index / refers to storage variable #

, , and / 0 refers to property variable ,, ,5 . The product of ,, ,5 and ,,

was replaced

, while the product of inlet property ,,5 and inlet flow , was replaced by with ,, ,5 ,,5 .


Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


, ,, ,5 21., + ,, ,5 + ,,

# , ,5

(21)

# ∈,

, = 1.,,5 21., + , + ,,

# ,5

# ∈,

∀( ∈ , ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1,, , ∀/ ∈ *

Constraint (22) states the relaxed version of equation (9). The product of ,, ,5 and the inlet flow ,,

# is replaced by the variable ,, # ,A,B . Following the same idea, the product of the inlet

property ,,

. # ,5 and the inlet flow ,, # was replaced by ,, ′ ,5 ,,

,,

# , ,5 = # ,5

0∈,

# ∈,

(22)

∀( ∈ !, ∀+ ∈ , ∀, ∈ -1, , ∀ ∈ 1,, ∩ ., , ∀/ ∈ *

Constraint (23) represents the relaxation for constraint (10). The product of variables ,, ,5 and

,,

was replaced with ,, ,5 , while the product of flow variable #, #, ,, and outlet property

#, ′, ,5 was replaced with variable #, #, ,,,5 . ,, ,5 =

$# , # , ,,%∈&

#, #, ,,,5

(23)

∀( ∈ ∪ *, ∀+ ∈ , ∀, ∈ -., , ∀ ∈ .,, , ∀/ ∈ *

Once all of the reformulations are performed, we can then state the convex relaxation of model (119) as follows: (1-6) (11-19) (20-23) ,< ,< ,< ,< ,, ,5 ≥ ,,

1,, ,5 + :1,,

,, ,5 − :1,,

1,, ,5 , , , , ,, ,5 ≥ ,,

1,, ,5 + :1,,

,, ,5 − :1,,

1,, ,5


(24) (25)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

, ,< ,< , ,, ,5 ≤ ,,

1,, ,5 + :1,,

,, ,5 − :1,,

1,, ,5 ,< , , , ,, ,5 ≤ ,,

1,, ,5 + :1,,

,, ,5 − :1,,

1,, ,5 , ,< ,< ,< ,< ,, ,5 ≥ , 1,, ,5 + 21, ,, ,5 − 21, 1,, ,5 #

#

#

#

#

#

#

#

#

#

, ,< , , ,< ,, ,5 ≤ , 1,, ,5 + 21, ,, ,5 − 21, 1,, ,5 #

#

(30)

#

#

(31)

, ,< ,< ,< ,< ,,

# , ,5 ≥ ,, 0 1,, ,5 + :.,, 0 ,, ,5 − :.,, 0 1,, ,5 #

#

#

#

, , , , , ,,

# , ,5 ≥ ,, 0 1,, ,5 + :.,, 0 ,, ,5 − :.,, 0 1,, ,5 #

#

#

#

, , ,< ,< , ,,

# , ,5 ≤ ,, 0 1,, ,5 + :.,, 0 ,, ,5 − :.,, 0 1,, ,5 #

#

#

#

, ,< , , ,< ,,

# , ,5 ≤ ,, 0 1,, ,5 + :.,, 0 ,, ,5 − :.,, 0 1,, ,5 #

#

#

(27)

(29)

, , ,< ,< , ,, ,5 ≤ , 1,, ,5 + 21, ,, ,5 − 21, 1,, ,5 #

(26)

(28)

, , , , , ,, ,5 ≥ , 1,, ,5 + 21, ,, ,5 − 21, 1,, ,5 #

Page 12 of 32

#

,< ,< ,< ,< ,, ,5 ≥ ,,

.,, ,5 + :.,,

,, ,5 − :.,,

.,, ,5 , , , , ,, ,5 ≥ ,,

.,, ,5 + :.,,

,, ,5 − :.,,

.,, ,5 , ,< ,< , ,, ,5 ≤ ,,

., ,,5 + :.,,

,, ,5 − :.,,

.,, ,5

,< , , ,< ,, ,5 ≤ , , .,, ,5 + :.,,

,, ,5 − :.,,

.,, ,5

,< ,< ,< ,< ,, , ,, ,5 − :,, , # , # ,5 ≥ ,, , # , # 1,, ,5 + : # ,0 1,, ,5 ,, ,#, # , , , , ,, , ,, ,5 − :,, , # , # ,5 ≥ ,, , # , # 1,, ,5 + : # ,0 1,, ,5 ,, ,#, #

,< ,< , , ,, , ,, ,5 − :,, , # , # ,5 ≤ ,, , # , # 1,, ,5 + : # , # 1,, ,5 ,, ,#, # , ,< ,< , ,, , ,, ,5 − :,, , # , # ,5 ≤ ,, , # , # 1,, ,5 + : # ,0 1,, ,5 ,, ,#, #

(32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43)

The domain of Constraints (24) to (43) are omitted for the sake of clarity. The proposed relaxation for model (1-19) using McCormick envelopes will be hereinafter referred to as the relaxed model, while model (1-19) itself, containing the original bilinear terms, will be referred to as the original model.


Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.2 Solution heuristic To obtain solutions for the original problem, we propose the following heuristic. First, we solve the

∗ relaxed model, from which we obtain a valid upper bound D$EF$G% % for the original problem and an

∗ ∗ optimal solution EF$G% . Next, we solve the original model using EF$G% as an initial point. To

increase the chance of finding good locally optimal solutions, we rely on a multi-start strategy, ∗ solving the problem in parallel with different starting points, all derived from EF$G% .

To obtain additional starting points, we generate small random perturbations in the optimal solution

∗ of the relaxed model EF$G% , seeking to obtain different points that lie within a ball centered at ∗ ∗ EF$G% . The idea behind the multi-start centered at EF$G% is to search for multiple local solutions and

to return the best of them. Notice that this local search is performed in the neighborhood of

∗ EF$G% because if the function does not diverge much in this neighborhood (which is rather true in

our case, given that the objective function is uniformly continuous), then the local maximum ∗ objective function value might be near the relaxed problem maximum objective function D$EF$G% %,

which is known to be a global maximum estimator.

In principle, we can repeat this procedure as many times as desired, as in a greedy search-based algorithm. In that algorithm, we keep an incumbent solution that is defined as the best solution E ∗

found during the iterative process. Note that this procedure has no convergence guarantee. However, early computational experiments have shown that, for practical purposes, the solution can be improved with a small number of repetitions before stabilization (typically 10 to 30, depending on the instance). Figure 2 presents a flowchart representing the proposed solution heuristic.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2 - Heuristic flowchart Notice that using the proposed heuristic, we are not only able to obtain reasonably good solutions, but also an optimality gap can be assessed for the current solution. The solution of the relaxed model ∗ ∗ EF$G% provides us with a valid upper bound (D$EF$G% %) for the original NLP model that can be

combined with any local optimum obtained from the NLP to calculate a global optimality gap

(please refer to Propositions 1 and 2 in the Appendix section for a detailed explanation). In practical situations, this optimality gap information can be as important as finding the solution itself because it helps to elucidate the quality of the current solution by allowing us to estimate how far the solution value is from the global maximum in the worst-case scenario.

5. Numerical experiments The numerical experiments performed are divided in two parts. First, we present a simplified version of the problem to illustrate the use of the proposed heuristic in a reproducible example. Next, we evaluated the performance of the proposed heuristic considering six distinct real-world instances. All numerical experiments were performed on an Intel i7 (2.8 GHz) with 12 GB RAM. The optimization model and the proposed heuristic were implemented in GAMS 23.7. The versions of the solvers used were CPLEX 12.3 (21), MINOS 5.51 (22), CONOPT 3 (23), KNITRO 7.0 (24) and BARON 9.3.1 (25).


Page 14 of 32

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.1. Example 1 To illustrate the use of the proposed heuristic, we devised an example instance of a simplified refinery. Our objective with this is to provide a reproducible example that illustrates the benefits of the proposed heuristic and how an optimality gap (i.e., valid upper and lower bounds) can be calculated for the problem using the solution obtained. In this example, a hypothetical refinery has three oils available for purchasing (namely O1, O2 and O3), two products (namely S1 and S2), and two properties (namely P1 and P2) that must be controlled for each product. The hardware is composed by a distillation tower with a single process mode and a single time period is considered. All flows are limited in 2 units. Figure 3 schematically illustrates this hypothetical refinery. Oil yields, properties, and prices are presented in Table 1, while final product prices and property ranges are presented in Table 2.

Figure 3 – Example 1 - refining scheme Yields (%) UPs properties Price S1 S2 P1 P2 O1 70 30 1 1 1 O2 60 40 3 2 0.2 O3 50 50 2 3 0.3 Table 1- Yields, UPs tabled properties and oil prices P1 Lower Bound

S1 S2

P2 Upper Bound

Lower Bound

Upper Bound

1.3 1.7 1.2 1.5 1.6 1.4 Table 2- Properties bounds and final products prices


1.9 1.8


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

For the sake of reproducibility, a single starting solution (the one obtained directly from the relaxed model) is provided to the original NLP model. The optimal solution for the relaxed problem is to acquire 2 units of oil O1, 0.75 units of oils O2, and 0.25 units of oil O3, selling 1.875 units of product S1 and 1.125 units of product S2, with a final profit of 1.85 units. As previously stated, given that the globally optimal solution for the original nonlinear problem belongs to the feasible region of the relaxation, the relaxed problem optimal value 1.85 is therefore an upper bound for the original problem’s optimal value. The property values P1 and P2 are 1.602 and 1.6 for S1, and 1.804 and 1.8 for S2. Notice that, even though the property constraints are relaxed (it is straightforward to check that the property P1 is not balanced for neither S1 or S2), they still limit the solution, as additional O2 and O3 could be acquired otherwise (and thus a larger profit could be obtained with the production of additional S1 and S2). Next, the original NLP model was solved using the relaxed model’s optimal solution as a warm start. In this case, the same solution is obtained in terms of the flows (i.e., amounts of O1, O2, O3, S1, and S2) and therefore, the same optimal objective function value is obtained (1.85 units). However, the properties values so the original balance constraints were met in this case. Provided that this solution is feasible for the original NLP model, it is thus a valid lower bound for the original problem. Since we have both an upper and a lower bound for the problem, we can calculate an optimality gap for the obtained solution. Because of the simplified nature, the optimality gap obtained for this problem is of 0% in this case, meaning that the solution obtained from the NLP problem is actually globally optimal. Notice, however, that this conclusion could only be reached because the upper bound obtained from the relaxation was available as the solution obtained by the NLP solver is proved to be only locally optimal. 5.2. Example 2 In this example, we present computational results for the real-world refinery that motivated this study. The input data considered is based on a Brazilian refinery, which is located in the Northeast region and was designed to produce asphalt, lubricants, and fuel oil. This refinery typically processes heavy and asphaltic oils that are specific for this purpose. The refinery also produces diesel, natural gas, liquefied petroleum gas, gasoline, maritime oil, fuel oil, and insulating oil for electrical transformers. The refinery’s production capacity is 1300 m3/day. Figure 4 presents the refining scheme for the refinery considered in this case study, showing all units considered and possible flows between them. This refinery has the following process units: vacuum distillation, hydrogen


Page 16 of 32

Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


generation, hydro-treatment, and a natural gas plant. Table 3 lists the code of the process units (Code), their names (Unit), type of refining process (Process Type), and how many process modes they have (Process Modes). The remaining units not presented in Table 3 have a single process mode.

Code VAC HG HDT NGP

Unit Process Type Process Modes Vacuum distillation Separation 4 Hydrogen Generation Conversion 2 Hydro-treatment Treatment 3 Natural gas plant Conversion 1 Table 3 - Refinery's process units

The units shown in Figure 4 that are not identified in Table 3 correspond to raw material supply points (Source), storage tanks of the intermediates and the final mixture (gray circles), and demand point for final products (Shipping Point). In this refinery, four final products are produced, namely two maritime and two fuel oils, each of which having its viscosity measured at 50 ºC and at 60 ºC as controlled properties. All the other data concerning operational characteristics of the refinery under study is omitted due to confidentiality reasons. However, masked data for research purposes can be made available upon request to the authors. From the data available, we derived six distinct instances that represent snapshots of the operational reality faced by the refinery during the period of one and a half years. The instances differ in terms of processing unit capacities, market parameters, initial storage and its properties, time-horizon length (two or three months, discretized into one-month discrete periods), and oil availability, which are different from each other in each instance. Although the one-month time period may seem too long for operational planning, it was a requirement of the company that supplied the data because this model is supposed to be used considering a midterm view of the refining process (notice however that the model could be straightforwardly adapted to consider shorter planning horizons by simply scaling the input data accordingly). The units capacities is measured in m³/month, and its value varies depending on the total number of days in a given month (28, 29, 30 or 31 days) if there is some maintenance scheduled during a given month or if there is a process rate reduction forecasted due to supply shortage or the lack of available storage capacity at marine terminals. Table 4 presents the model size of both the original and relaxed models. As seen from this table, the number of variables and constraints are in the thousands, which presents a challenge for the available NLP algorithms. Moreover, the relatively large number of nonlinear constraints makes the problem even harder to solve.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Original Model Instance

Variables

Constraints

1 2 3 4 5 6

1159 1015 1114 1207 1746 1513

1097 963 1059 1180 1647 1444

Page 18 of 32

Relaxed Model Nonlinear Constraints 336 280 336 294 504 420

Variables

Constraints

1889 1625 1844 1861 2850 2437

4017 3403 3979 3796 6063 5140

Table 4 - Model Size


Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Figure 4 – Example 2 - refining scheme



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

To test the proposed heuristic for obtaining a solution to the problem at hand, we must solve a linear programming problem and then repeatedly solve a nonlinear programming problem with different initial points. For this purpose, we have tested several popular offthe-shelf solvers. The relaxed model was solved using CPLEX (21), while the nonlinear problems were solved using MINOS (22), CONOPT (23), and KNITRO (24). MINOS (22) uses Reduced Gradient and an Augmented Lagrangian algorithm. CONOPT (23) uses Generalized Reduced Gradient. KNITRO (24) uses the Barrier Method, also known as the Interior Points Method. These three solvers used for the nonlinear problem are local NLP algorithms, meaning that they return solutions that are not guaranteed to be either a local or a global optimum. Because the solvers use different methods, they often return different solutions and present different computational performances in finding locally optimal solutions. However, it is not a trivial task to decide in advance which method will perform better on a certain problem. Therefore, to assess the method that would perform best in our case, we tested these three solvers. For test purposes, the original model was solved directly using the three nonlinear solvers by considering their standard settings. Table 5 summarizes the results we obtained when solving the original model using the standard (zero) initial point and the local solvers MINOS, KNITRO, and CONOPT.

Instance 1 2 3 4 5 6

MINOS KNITRO Elapsed Solution Elapsed Solution Value Time (s) Value ($IJK ) Time (s) ($IJK ) 2.76 0.45 9.86 Infeasible 2.97 1.01 10.03 Infeasible 2.76 1.31 10.07 Infeasible 2.47 0.49 9.98 Infeasible 3.07 0.98 10.21 Infeasible 3.39 0.91 11.12 Infeasible Table 5 – Results obtained with local solvers

CONOPT Elapsed Solution Value Time (s) ($IJK ) 2.64 Infeasible 2.64 Infeasible 2.69 Infeasible 2.44 Infeasible 2.86 Infeasible 2.62 Infeasible

As seen in Table 5, MINOS was the only solver capable of directly solving the problem considering its standard settings. Note that both KNITRO and CONOPT were not able to find feasible solutions in any of the six instances considered, which might be explained by the nature of the problem and how difficult it is to obtain a good initial solution for it.


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Furthermore, using local solvers to search for optimal solutions in a nonconvex problem does not provide any information regarding the global quality of the solution (i.e., optimality gaps cannot be obtained given that no valid lower bound is provided by these solvers). One possible alternative is to use global optimization solvers, such as BARON, that are capable of finding global solutions within a prespecified tolerance gap. A point in the feasible region is said to be a global solution if it is guaranteed to be the best solution available in the same feasible region. For more information about these methods, please refer to Nocedal and Wright (26) and Biegler (27). To use BARON as a benchmark we set the same stop criteria that were used for the proposed heuristic, namely 10 minutes (600 seconds), 50 iterations, and an optimality gap of 10%. In order to verify whether the criteria set were overly restraining for the solver, we performed additional experiments only considering a time limit of 1 hour (3600 seconds). The results obtained using BARON in both experiments are presented in Tables 6 and 7, respectively. Instance Elapsed time (s) UB ($104) LB ($104) Opt. gap (%) 555.36 1.14 1.03 10.00 1 600.00 1.50 1.26 19.16 2 5.17 1.99 1.89 5.15 3 212.69 1.42 1.38 3.02 4 600.00 1.64 1.40 17.58 5 600.00 2.04 1.59 28.47 6 Table 6 - BARON results with all stop criteria Instance UB ($104) LB ($104) Opt. gap (%) 1.14 1.06 7.02 1 1.50 1.38 8.69 2 1.99 1.99 0.00 3 1.42 1.38 2.94 4 1.64 1.53 7.01 5 2.04 1.63 24.97 6 Table 7 - BARON results after 1 hour without iteration and gap limits In the first tests performed, BARON was able to obtain a feasible solution for all instances. In three of the instances, the solver was not able to reach the desired optimality gap of 10%



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

in less than 10 minutes. As can be seen in Table 7, the additional experiments have shown that the solver was not capable of improving the upper bounds obtained in the first experiments, being thus the observed optimality gap improvement caused by the improvement of the lower bound. The results obtained with the proposed heuristic are presented in Tables 8, 9, and 10, in which we present the results using MINOS, KNITRO, and CONOPT, respectively, to solve

the instances. In all experiments, L (i.e., the number of distinct initial solutions in the neighborhood of the relaxed problem solution generated at each iteration) was set to 50,

and the alternative initial points were generated by perturbing each variable with a random increase or decrease of 10% using a uniform distribution. Stop criteria used were the same as in the experiments using BARON (i.e., 10 minutes (600 seconds) and optimality gap of 10%), with an additional 50 iteration limit. Instance 1 2 3 4 5 6

Elapsed Time (s) UB ($IJK ) LB ( $IJK ) Opt. Gap (%) 25.21 1.14 0.98 16.01 23.53 1.50 1.34 11.72 3.66 1.99 1.98 0.77 23.78 1.42 1.22 16.16 32.12 1.64 1.46 12.73 27.89 2.04 1.65 24.10 Table 8 - Proposed heuristic with MINOS


Elapsed Time (s) UB ($IJK ) LB ($IJK ) 442.86 1.14 1.07 319.29 1.50 1.39 24.25 1.99 1.99 44.00 1.42 1.38 489.53 1.64 1.60 578.52 2.04 1.77 Table 9 - Proposed heuristic with KNITRO


Elapsed Time (s) UB ($IJK ) Opt. Gap (%) LB ($IJK ) 13.68 1.14 1.04 9.93 13.39 1.50 1.38 8.88 3.13 1.99 1.98 0.77 4.87 1.42 1.38 2.99 20.76 1.64 1.54 6.90 20.50 2.04 1.72 18.44 Table 10 - Proposed heuristic with CONOPT


Opt. Gap (%) 6.94 8.05 0.00 2.81 2.70 15.54

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As seen in Tables 8, 9, and 10, the proposed heuristic was capable of providing feasible solutions, (and thus information regarding optimality gaps), for all instances considered, which was not the case when we attempted to solve the problem using KNITRO or CONOPT. In addition, in all instances we observed an improvement in the objective function value when we compare the solution with the solutions obtained using MINOS as solver, some being particularly remarkable, as in Instances 1 and 4. In cases in which the optimality gap reported is higher than 10%, the procedure stopped because of the iteration limit set to 50. Another advantage of using the proposed heuristic is the possibility of establishing optimality gaps for the obtained solution, which consist valuable information that would not be available with the use of standard NLP solvers. As previously stated, the optimal objective function value of the relaxed problem is a valid upper bound for the original NLP problem, which allows us to use it to calculate an optimality gap for the locally optimal solution obtained from the NLP solver and thus verify the quality of the solution obtained. It is important to highlight that, despite that possibility, in the proposed heuristic only the lower bound is updated iteratively, while the upper bound remains the same throughout the execution of the heuristic. Regarding the differences between the alternative solvers, the heuristic achieved better gaps using KNITRO, although the optimality gaps obtained with CONOPT were comparable and the elapsed time was smaller in this case. In all cases, the proposed heuristic provided the same upper bound obtained by the global solver BARON. However, the heuristic was capable of finding better feasible solutions and, thus, better lower bounds for all cases using KNITRO and MINOS, and for cases 1, 3, 5 and 6 using MINOS. Even when BARON was let running for one hour (additional experiments), the heuristic using KNITRO still outperformed or obtained the same solution as BARON; CONOPT presented better performance for instances 2, 4, 5, and 6; and MINOS was better in instances 6. Moreover, using CONOPT or MINOS as the local NLP solver within the proposed heuristic presented smaller computational times than BARON, while using KNITRO presented a competitive solving time. 5. Conclusions



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In this paper, we presented a strategy to solve the ROPP, which, because of the nonlinear nature of the processes involved, presents several challenges in terms of obtaining reasonably good plans. The proposed strategy consists of the solution of a convex relaxation of the original problem and subsequent use of its solution as a starting point for the local NLP algorithms. The proposed strategy was able to provide feasible solutions for the problem considered with data from a real-world refinery. In addition, we demonstrated that the solutions obtained from the relaxation can be used to calculate optimality gaps for the solutions obtained, which, in practice, consists of valuable information that is not provided by traditional off-the-shelf local solvers. The numerical results suggest that this strategy can be used to develop a heuristic that is capable of generating good starting solutions for local solvers, solving cases that were originally unsolvable using standard local solvers. In addition, the proposed heuristic appears to be competitive with the current state of the art global solvers, such as BARON, for this particular problem, identifying better optimality gaps in less computational time. Future research could be performed to both extend the mathematical model and the proposed strategy. Both the mathematical model and the strategy could be extended to other refineries and their particularities, such as properties that depend on the stream mass instead of volume, which will generate extra constraints with trilinear terms. Another direction of investigation is to use alternative forms of relaxation, such as piecewise McCormick envelopes or multiparametric disaggregation with a slack variable, which would generate mixed-integer relaxations instead of purely linear ones.

Acknowledgment We gratefully acknowledge the support of the Tecgraf Institute (PUC-Rio), in particular of Prof. Silvio Hamacher, for the discussions and insightful suggestions. This research has been financially supported by the CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brazil) Grant 455013/2014-4.


Page 24 of 32

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nomenclature Sets Set of process units (u, u') Set of operational modes (c, c') Set of streams (s, s') Set of properties (p, p') Time periods {n | n = 1,..., T} Operational modes c performed at unit u Properties p of inlet streams s at unit u Properties p of outlet streams s at unit u Inlet streams s at unit u Outlet streams s at unit u Tanks of raw materials Delivery units for final products Tank units (storage and blending) Storage tanks Blending tanks Process units (separation and conversion) Flows (u,c,s,u',c'), where pair (u,c) is the origin and pair (u',c') is the destination Variables Flow rate of stream s between (u,c) and (u’,c’) Inventory level at unit u Property p of inlet stream s at unit u Property p of outlet stream s at unit u Inlet flow rate of stream s at unit u Outlet flow rate of stream s at unit u Parameters Product price Cost of additional raw material Process unit yield Initial stock at unit u Minimum storage capacity at unit u Maximum storage capacity at unit u Minimum demand Maximum demand Maximum offer of additional raw material Property value p of outlet stream s at process unit u as a function of inlet stream s’ Property value p of initial stock Minimum feed flow rate at unit u


* ⊂ .,, ⊂ 1,, ⊂ -., ⊂ -1. ⊂ ⊂ ⊂ * ⊂ ⊂ ! ⊂ ⊂

=,, ,#, # ⊂ ,, , # # , ,,A,5 ,, ,5 ,,

,,

,

,

4,, , # 21., ,< 21, , 21, ,< ;!,,

, ;!,,

:1,,

1,, ,

# ,5 1.,,5 ,< :.,,


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Maximum feed flow rate at unit u Lower bound of outlet property p at unit u Upper bound of outlet property p at unit u Lower bound of inlet property p at unit u Upper bound of inlet property p at unit u Minimum stream s flow rate leaving unit u Maximum stream s flow rate leaving unit u

Minimum stream s flow rate between unit u and unit (′

Maximum stream s flow rate between unit u and unit (′


Page 26 of 32

, :.,,

,< 1,, ,5 , 1,, ,5 ,< .,, ,5 , .,, ,5 ,< :1,,

, :1,,

,< :,, , ′ , ′ , :,, , # , #

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Appendix A - Convex relaxation using McCormick Envelopes Consider the following compact representation of a nonconvex NLP problem: : ! P D$E, Q%

(44)

T< ≤ E ≤ T

(46)

R $E, Q% = 0

(45)

4< ≤ Q ≤ 4

(47)

where R $E, Q% = EQ + RU $E, Q%; and D$E, Q% and RU $E, Q% are affine functions. First, let us

define two functions, ℎ $E% = E − T < and ℎU $Q% = Q − 4 < . Those functions are

nonnegative by construction, and it is easily shown that the product of two nonnegative functions is also a nonnegative function, that is,

$ℎ ℎU %$E, Q% = $E − T < %$Q − 4 < % = EQ − E4 < − T < Q + T < 4 < ≥ 0. In this case, we have that

EQ ≥ E4 < + T < Q − T < 4 < .

The remaining supporting hyperplanes can be obtained in a similar way by defining the

other possible combinations of the variable bounds as functions ℎ $E% and ℎU $Q%, resulting in the following constraints:

EQ ≥ E4 < + T < Q − T < 4

A Strategy Based on Convex Relaxation for Solving the Oil Refinery

Recommend Documents