Dealing with Multiple Tank Outflows and In-Line Blending in

Feb 26, 2019 - The crude oil scheduling problem has been focus of many studies in the past, which is justified by its importance in the oil industry. ...
1 downloads 0 Views 508KB Size
Subscriber access provided by Washington University | Libraries

Process Systems Engineering

Dealing with Multiple Tank Outflows and In-Line Blending in Continuous-time Crude Oil Scheduling Problems Sérgio Mauro da Silva Neiro, Valeria Viana Murata, Bárbara Jahn, Raiana Roland Seixas, Estefane Horn Hollmann, and Cristiane Salgado Pereira Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b03749 • Publication Date (Web): 26 Feb 2019 Downloaded from http://pubs.acs.org on March 4, 2019

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 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 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.

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 44 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

Dealing with Multiple Tank Outflows and In-Line Blending in Continuous-time Crude Oil Scheduling Problems Sergio Mauro da Silva Neiroa,*; Valéria Viana Murataa; Bárbara Jahna; Raiana Roland Seixasa; Estefane Horn Hollmannb; Cristiane Salgado Pereirab a – Uberlândia Federal University, UFU – Av. João Naves de Ávila, 2121. Uberlândia - MG, Brazil 38400-902. b – Petrobras Research Center - CENPES, Cidade Universitária - Rio de Janeiro – RJ, Brazil, 21941-970 * corresponding author

Abstract The crude oil scheduling problem has been focus of many studies in the past, which is justified by its importance in the oil industry. Optimizing crude oil blending is of paramount importance given that it can substantially impact the economic performance of a refinery. In this work operational features of a realworld existing refinery are addressed. Only operations limited to the refinery battery are in scope such as splitting of parcels unloaded from a supplying pipeline segment, tank heels and capacities, brine settling time, multi-quality tracking, Crude Distillation Unit (CDU) straight-run products profile, multiple tank outputs and multiple CDU inputs. Moreover, different policies as to the handling of the refinery tank farm are evaluated. The presented formulation is taking from the Multi-Operation Sequencing (MOS) proposed by Mouret, Grossmann and Pestiaux (Comp & Chem Engng, 35, 1038-1063, 2011). The novelty of the present work relies on how to accurately impose lower and upper bounds on the flowrate of multiple tank outputs given that the formulation is based on a unit-specific time grid. Two approaches are proposed and evaluated. The resulting MINLP models are exhaustively tested with six distinguishing real-world scenarios in which blending can include up to 36 crude oil grades, different tank availability, initial inventories and scheduled parcels to arrive at the refinery. Two solution algorithms that avoid solving the

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 2 of 44

full-scale MINLP problem are used and compared. The computational experiments show that the proposed formulations are able to handle a wide range of problem instances in reasonable computational time, despite of the size dimension.

Keywords Crude oil scheduling, refinery, MINLP, multiple resource outputs, multiple time-grid, continuous time. * corresponding author: [email protected]

1. Introduction A lot of attention has been devoted in finding alternative renewable fuels in recent years in order to diminish the world´s dependency on fossil fuels. The process industry though still heavily relies on nonrenewable energy resources such as petroleum derivatives and coal. The world proven crude oil reserves have not increased significantly in the last years while the demand trend line for crude oil has been consistently increasing since the early 80´s, which has led the world crude oil processing capacity to follow along the demand take off1. As a result, it is imperative that companies operate efficiently in an ever changing competitive global market while environmental regulations become more stringent and the availability of crude oil becomes scarcer and unstable price is observed due to geopolitical situations. Optimizing the whole production supply chain would likely yield the best saving results. Although Enterprise-Wise Optimization (EWO) has been recognized as an important emerging area of research that is able to deal with horizontal and vertical plant operations integration2, the optimization of entire supply chains is still challenging due to multiscale temporal and spatial integration, let alone their complexity and large dimension. Therefore, the traditional approach of addressing major components of the supply chain is still in use these days concerning the oil industry. The same applies for other industry sectors, as

ACS Paragon Plus Environment

Page 3 of 44 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

well. In this work, we focus on the front-end operations of an in-land refinery. However, only the operations comprised in the refinery battery limits are in scope. One of the first contributions concerning the crude oil scheduling problem goes back more than 20 years ago. Lee et al.3 presented a mixed-integer linear programming (MILP) model based on the discrete time representation. The problem addressed in their work deals with a typical in-land refinery comprised of a set of small vessels carrying single oil amounts that are unloaded to dedicated storage tanks available at an on-shore terminal. Discrete amounts of the unmixed crude oils are transferred to charging tanks at the refinery where crudes are blended to satisfy quality restrictions and then fed to crude distillation units (CDU) for processing. The mathematical model was maintained linear by avoiding explicit mixing quality calculations. Material balance for each mix component was stablished individually and composition discrepancy was allowed at the splitting of streams outflowing a tank. As reported by Wenkai et al.4, composition discrepancy can only be eliminated by introducing non-linear quality constraints, which turns the formulation into a mixed-integer nonlinear programming (MINLP) problem and thus more difficult to solve. These authors, besides improving the formulation and including new features to the problem initially addressed by Lee et al.3, also proposed an algorithm that cycles between an NLP and an MILP problem that include nonlinear mixing constraints for eliminating composition discrepancy. However, the algorithm may fail to find a feasible solution even when it is known there exist one7. Jia et al.5 introduced the event-based continuous time representation for the crude oil scheduling problem. The two outstanding benefits of using continuous time representation in comparison to the discrete counterpart are; reduced dimension and better utilization of the time domain, since event points can be positioned anywhere within the scheduling horizon. The discrete time grid, on the other hand, allows operations to take place only at the boundaries of prefixed discrete intervals, so that accuracy is

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 4 of 44

achieved by increasing the number of discrete intervals and thus increasing the problem dimension. Nevertheless, continuous time formulations require the prior postulation of events to be used, which is not trivially determined beforehand. Moreover, continuous time formulations generally make use of big-M constraints that lead to poor relaxations, which in turn affects the computational performance. Some of the just discussed aspects can be observed when comparing the formulation by Lee et al.3 and that proposed by Jia et al.5 . Just like in Lee et al.3, Jia et al.5 proposed an MILP formulation that allows composition discrepancy. However, an instance with introduced bilinear mixing equations was also solved resulting in an MINLP problem in which composition discrepancy is eliminated at the expense of an excessive computational cost. Moro and Pinto6 proposed an MINLP formulation for the crude-oil short-term scheduling concentrated on the refinery tank-farm operations which was based on the uniform-time-grid representation, also known as global time grid. The main objective was to maximize units feed while minimizing the number of tank operations. Their approach relied on subdividing the scheduling horizon in non-uniform time slots whose lengths were determined as part of the optimization decisions. The set of time slots was split in two mutually exclusive subsets; one of the subsets with a predefined number of slots was dedicated to deal with unloading parcels arriving at the refinery through a pipeline, whereas the complement set was used for allocating tanks to CDUs. An attempt was also made to linearize the formulation by discretizing inventory levels, which led to a very large MILP problem that demonstrated to be more difficult do solve than its MINLP version. Reddy et al.7,8 addressed the crude oil short-term scheduling problem involving a coastal refinery. The first contribution presented a formulation based on the discrete time representation that allowed multiple crude transfers in each period, while the latter used a hybrid time representation, discretecontinuous, which shares some resemblance to the time grid framework proposed by Moro and Pinto6 in

ACS Paragon Plus Environment

Page 5 of 44 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

that the discrete intervals were created based on the arrival and departure times of vessels. A number of variable length time slots was allocated to each discrete interval so that operations could be accommodated. According to the definitions from the literature, a coastal refinery differs from an in-land refinery in that the former does not include storage tanks, i.e., only charging tanks are involved, which means that vessels unload directly to charging tanks. Real-world refinery operational features were treated such as including a high-volume single buoy mooring (SBM) pipeline used for unloading multi-parcel very large crude carriers (VLCC) and single-parcel tankers unloading through jetties. Other features include demurrage, settling time, two tanks feeding a CDU and changeover at CDUs. The authors raised the fact that previous reported formulations, including continuous and discrete ones, did not satisfactorily addressed the composition discrepancy in transfer operations. In order to address that issue they proposed to solve their MINLP problem iteratively by a specialized heuristic in which discrepancy is totally eliminated and yet only a series of MILPs are solved. The heuristic did not guarantee global solution. However, one must bear in mind that no nonconvex MINLP problem can probably be solved to optimality unless a global optimization algorithm is employed, which is usually very expensive in terms of solution time and computational resources. This heuristic will be further detailed in the Solution Algorithms section. Li et al.9 expanded the formulation proposed by Reddy et al.7 introducing 15 crude properties that are critical to crude distillation and downstream processing and proposed a robust algorithm based on a backtracking strategy that uses integer cuts to solve the MINLP problem. Moreover, they also developed a partial relaxation strategy that further increase solution speed and improve solution quality.

Furman et al.10 proposed an event-based continuous time formulation that robustly handled the synchronization of time events with material balances. The authors argued that previously reported event based continuous time formulations could potentially generate inconsistent solutions due to the fact that the usual way material balances were laid out only considered the net input and output during some time event and not the timing of the involved transfer operations, leading to possible tank capacities violation.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 6 of 44

Another novelty in their approach was to enable the reduction of the required number of time events by allowing tank input and output flows to occur in the same time event with the only condition that the time spans of the input material transfers did not overlap the time spans of the output ones and that input occurred before output. The authors also developed convex relaxation constraints based on McCormick estimators and a set of valid inequalities designed for pooling problems that could be added to the proposed MINLP formulation as an attempt to improve tightness. However, in both cases the added constraints had the drawback of increasing size dimensions without improving computational time. Karuppiah et al.11 presented an outer-approximation algorithm for finding the global optimum of a nonconvex MINLP problem that was largely taking from Furman et al.10. Saharidis et al.12 presented a novel formulation based on events whose main objective was to minimize tank operations while satisfying operational constraints. The prior known events were the arrival time of vessels carrying crude oil, their departure times which was fixed by contract to be 36 hours after their arrival and the demand of different blend recipes for each CDU along the scheduling horizon (when, which grade and what quantity). The scheduling horizon was then discretized into a set of periods of variable length, forming a common event-based representation, in which period boundaries were enforced to match the known events. Two different recipe preparations were considered; standard that satisfied an exact crude composition mix and flexible, which was a relaxation of the first where only lower and upper quality bounds had be satisfied. Furthermore, blending was admitted to take place in either tanks or manifolds and valid inequalities were included to expedite solution times. It must be emphasized though that the proposed formulation is of limited application since CDU´s demand profile is usually taken to be part of the decisions to be optimized. Mouret et al.13 proposed a continuous-time formulation based on priority-slots named by the authors as Single-Operation Sequencing (SOS) model. The novelty of this contribution was the

ACS Paragon Plus Environment

Page 7 of 44 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

consideration of the schedule as a sequence of transfer operations in which an operation with higher priority should be allocated to a priority-slot before all the others with lower priority, taking into consideration that priority-slots were sequentially ordered. Moreover, if two operations were not allowed to overlap each other they should be allocated to succeeding priority-slots, with the task with higher priority occurring first, while their time spans should not overlap. Variables and constraints were framed based on transfer operations, which were uniquely identified by their origin and destination. That framework allowed the inclusion of symmetry breaking constraints that drastically improved computational times due to a significant reduction of the search space. An outstanding feature of the proposed formulation is that only a single operation could be allocated to each priority-slot, which usually led to the requirement of a high number of pre-postulated priority-slots, although that did not implicate in high computational cost for the used test examples. Two years later, the same authors14 proposed a relaxation to the SOS model, in which multiple operations were allowed to take place in the same priorityslot. The Multi-Operation Sequencing (MOS) model enabled solving the same test problems with a much smaller number of priority-slots. In both cases a two-step MILP-NLP approach was used for solving the resulting nonconvex MINLP problems. This procedure demonstrated to be able to produce near-optimal solutions. More about this procedure is discussed in the Solution Algorithms section. Chen et al.15, conducted a comparative study considering three of the most prominent continuoustime models proposed for addressing the crude oil scheduling problem for in-land refineries. Among the selected models were the Event-Based model by Jia et al.5, the MOS model by Mouret et al.14 and the Unit-Slot (US) model introduced by Hu e Zhu (apud Chen et al.15). The first two models have already been discussed in previous paragraphs whereas the latter can be seen as an Event-Based variant. What is most remarkable about this model is that the time variables are resource-centric. In that case, one cannot tell just by the time variables if an inflow or an outflow occurs at a given resource because time variables

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 8 of 44

do not carry information on the origin or destination. The allocation variables then dictate the flow direction. This strategy thus requires fewer time variables. On the other hand, special timing constraints must be written for coordinating transfer operations between two distinct resources. Li et al.16 proposed a continuous time formulation for addressing a coastal refinery in which all industrial features pointed by Reddy et al.7 and Li et al.9 were also included in their model. Besides, a branch-and-bound global optimization algorithm that makes use of piecewise linear underestimation was devised. In all the previously discussed continuous-time approaches it is assumed that a maximum flowrate is imposed to each transfer operation individually. However, in some instances a tank might have multiple active outflowing operations starting and ending at different times (multiple time-grid) while taking place at the same time slot and still being subjected to a common maximum flowrate due to the fact that only one pump is available (Figure 1). In this case, one must certify that the total output flowrate does not exceed the installed pumping capacity while dealing with the distinct time grids. This operational feature is to be addressed in the present work. Outflowing operation 1 time slot i

Tank

Outflowing operation 2

time slot i Outflowing operation 3 time slot i

time window spanned by multiple output operations

Figure 1 – Illustration of a tank with multiple simultaneous outflows within the same time slot i. On top of having 1 or more flow-out from the same tank, another operational feature to be discussed in the present work is the possibility of having two non-synchronized inflows to the same CDU. As seen in Figure 2a, the multiple time grid representation enables the use of a reduced number of time slots. However, time synchronization among the different transfer operations must be carefully handled.

ACS Paragon Plus Environment

Page 9 of 44 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

Yet, in order for one to take full advantage of the 2 flow-in to the same CDU, no quality constraints should be imposed on the tank outflowing operations. Quality constraints must though be imposed to the resulting CDU inflow operations, requiring that tank outflowing operations be decoupled from CDU inflowing operations (Figure 2b). The tank outputs and CDU input are then interrelated through a fictitious in-line mixer positioned between tanks and CDUs. Our proposal is largely based on the MOS framework proposed by Mouret et al.14. However, major changes are introduced as to parcel unloading and especially to the connections between tanks and CDUs. We also introduce the novel concept of an Enclosing Flexible Time Slot – EFTS (Figure 2a) in which a transfer operation is allowed to be active in a time slot but not spanning the whole enclosing time window for a given time slot. This concept is introduced with the purpose of being able to accurately impose lower and upper flowrate bounds at tanks with multiple simultaneous outputs. A large-scale crude oil scheduling problem representing a real-world inland refinery in Brazil was used for testing the performance of the proposed model. A detailed description of the problem under study, which do not include terminal operations, is given in section 2. The resulting largescale MINLP problems are solved using two different algorithms taken from the literature and compared. The rest of the text is organized as follows. The problem description is given in section 2, some of the highlighting features of the MOS model is presented in section 3, and the proposed approaches as modifications of the MOS model are presented in section 4, solution algorithms are discussed in section 5, whereas a note on the minimum number of time slots is discussed in section 6. Results and discussions are presented in section 7 and conclusions are finally drawn in section 8.

2. Problem description A typical crude oil scheduling problem involving an in-land refinery is treated in the academic literature as being comprised of a set of small vessels each carrying a single crude oil grade, which is unloaded completely to one of the various storage tanks available at an onshore terminal. The storage

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 44

tanks can be either dedicated or hold crude oil mixtures, being the dedicated mode the most common consideration. From the terminal, crude oil parcels composed of either a single crude oil grade or a mix of them, depending on the adoption of tanks being dedicated or not, are transferred to an inland refinery tank farm through a pipeline segment, where final crude oil mixes are blended in charging tanks to satisfy quality specifications before being sent to CDUs for processing. The time and logistics involved in transferring the crude oil between the refinery and the terminal is usually neglected. Transfer from tank 1 to CDU 1

Tank 1

time slot i

Tank 2

Transfer from tank 2 to CDU 1 time slot i Transfer from tank 2 to CDU 2 time slot i

Transfer from tank 1 to CDU 2 time slot i

CDU 1

CDU 2 2

Feeding to CDU 1 time slot i

Feeding to CDU 2 time slot i Enclosing time window for both tanks and CDUs

(a) Tank 1 Tank 2

Transfer from tank 1 to CDU 1 time slot i Mixer

Feeding to CDU 1 time slot i

time slot i Transfer from tank 2 to CDU 1

(b)

CDU 1

Figure 2 – Illustration of a situation with 2 flow-out from the same tank and 2 flow-in to the same CDU within the same time slot (a); Decoupling of tank outflows from CDU inflows (b). The problem addressed in this work is as depicted in Figure 3. Focus is given on the refinery tank farm management because the refinery plant scheduler does not own the decisions involving crude oil supply. That decision is made on a higher management level. Therefore, the refinery scheduler is concerned only with the operations within the refinery battery. It is assumed that a set of crude oil parcels

ACS Paragon Plus Environment

Page 11 of 44 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

composed of different crude oil grades, denominated here as components, are scheduled to arrive at the refinery site at specified instant of times, which are spread along the scheduling horizon. It is also assumed to be known their total volume and composition. A parcel unloading is accomplished at a specified constant flowrate, which can be different for each parcel. Crude oil parcels might be split in multiple subparcels and unloaded to multiple tanks. However, a one-to-one connection between the pipeline and storage tanks is assumed at any time. Moreover, a parcel can only be assigned once for each tank, which means that if a parcel is split the resulting sub-parcels are unloaded to different tanks. What determines the minimal volume of each sub-parcel is the minimum run-length of 3 hours for the unloading operations. Once unloading is started no interruption is allowed even when a parcel is split into multiple sub-parcels, which means that the instant of time at which total unloading is complete is known upfront, i.e., the total unloading time window for each parcel is fixed and determined given that the arrival time, its volume and unloading flowrate are known. Tanks are not dedicated and their total volume can never go below the nominal tank heel due to the pump head. Tanks are not allowed to load and unload simultaneously. Once a tank receive a parcel, a settling period of time of at least 24 hours is imposed. Tanks are allowed to feed 1 or 2 CDUs simultaneously. In that case, an upper flowrate limit is imposed on the sum of both outputs, because there exists only one pump per tank. Distillation units must operate continuously throughout the scheduling horizon. The refinery is comprised of 3 CDUs each with a different capacity. The minimum run-length of each feeding operation is 24 hours so that the unit can reach the new steady state operating point. The minimum run-length time also avoids frequent tank switches to the same CDU. Distillation units can have 2 flow-in, in which case the total feed flowrate must lie between lower and upper processing capacities. As previously mentioned, no quality restrictions are imposed on tank outputs. Quality constraints are though imposed at the inlet of CDUs. Three qualities are considered; specific gravity (API),

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 12 of 44

Total Acid Number (TAN) and Sulphur content (S). The most critical quality for the CDU feed is the TAN, which upper limit is set to be 1.30 mgKOH/g. Two distinct configurations or mode of operations are considered in this study: namely, Full Connectivity – FC (Figure 3a) and Limited Connectivity – LC (Figure 3b). In the former, parcels can be unloaded to any of the active tanks. Moreover, any two tanks (see available connectivity between tanks and CDUs in Figure 3a) can be assigned to feed the same CDU simultaneously. In case of having two tanks assigned to the same CDU, each tank must make up at least 5% of the total CDU feed. In the LC configuration, parcels are classified in two groups, namely: Injection Parcels (IP) and Base Parcels (BP). IPs, identified in orange in Figure 3b, are distinguished from BPs (blue parcels in Figure 3b) for being crudes of better qualities that are typically used in small portions just enough to specify the tanks that present poor quality. Therefore, they must be segregated and used in small quantities. IPs are stored in dedicated tanks identified in this work as Injection Tanks (IT), orange tanks in Figure 3b. On the other hand, BPs can only be unloaded to Base Tanks (BT), blue tanks in Figure 3b. Given what was exposed previously concerning IPs, ITs can only be assigned to feed a CDU if a BT is also assigned and the IT contribution should lie between 5% and 30% of the total CDU feed. The light blue tanks in Figure 3 represent tanks that are currently not holding crude oil. They might be out of service for maintenance or receiving intermediate streams from another refinery. The set of inactive tanks varies along the time which is captured thought the representation of different scenarios, discussed latter in the text. The crude processed at CDUs generate a set of distilled intermediate products such as LPG, naphtha, kerosene, diesel and a tower bottoms (also known as atmospheric residue or reduced crude), which is usually further processed in a Vacuum Distillation Unit (VDU). In this work, we are particularly interested in the diesel cut, for which a maximum sulfur content of 1.70 (% w/w) must be imposed. The

ACS Paragon Plus Environment

Page 13 of 44 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

amount of the diesel cut produced at a CDU is calculated through constant yields provided for each crude oil grade. Inline mixer 1

TA TB

Crude oil parcels

D1

TC

Crude oil parcels

TB

TD

TE

TE

D2

TF

Inline mixer 2

TG

TH TI TK

Inline mixer 3

TL

(a) Full connectivity

Inline mixer 2

TG

D3

TI

D2

TF

Injection parcel

TH

TL

D1

TC

TD

TK

Inline mixer 1

TA

D3 Inline mixer 3

Injection Tank

(b) Limited connectivity

Figure 3 – Schematic representation of the refinery under study. Six different scenarios were used for performance evaluation. Each scenario represents a realworld operational situation at the refinery. Scenarios differ as to the number, volume, composition and time of arrival of parcels, as well as to the initial inventory in storage tanks. Moreover, the set of tanks that are inactive may also vary for each scenario. The crude oil slate available for the refinery is comprised of 36 crude oil types. The problem under study represents a real-world refinery in Brazil and each scenario represents a different operational situation at the refinery. Input data is detailed in the Supporting Information section. The crude schedule produced by the refinery scheduler is very dependent on the initial tank conditions (volume and quality) and on the parcel supply schedule (arrival time, volume, quality and unloading flowrate). The testing scenarios were carefully selected to represent a wide range of initial conditions, supply schedule, as well as scheduling horizons.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 14 of 44

3. The Multi-Operations Sequence (MOS) model The MOS model, introduced by Mouret et al.14 is based on the multiple time grid representation and as such benefits from the fact that each resource has its own time grid, usually leading to the requirement of a lower number of time slots in comparison to other time representations. What is genuine about this formulation is that the model framework is fundamentally conceived based on two key elements: priority time slots and operations performed by units. Within the context of the crude oil scheduling problem, operations designate transfer of materials between resources, as illustrated in Figure 4. An operation establishes a connection between two resources through a material flow. Each connection is unique and thus each operation must univocally be identified. The time window in which a unit performs a task is determined by the time window in which a transfer operation takes place. In that case, both the origin and destination resources involved in the operation must be allocated for the operation to occur. Moreover, if two or more operations cannot occur simultaneously, they are said to form a clique and are considered to be non-overlapping. Cliques are introduced to the model through a non-overlapping matrix, denoted as 𝑁𝑂𝑣,𝑣′, in which an entry in row v and column v’ is 1 if operations v and v’ form a clique. Indeed, only Maximal Cliques are used. Putting it in a very simple way, a Maximal Clique involves the maximum number of non-redundant cliques. The use of Maximal cliques allowed the authors to write timing constraints that produced tighter relaxations (we refer the reader to Mouret et al.13 for a more detailed discussion). Running into symmetric solutions is an inherent characteristic of scheduling problems and the more symmetric solutions exist the more difficult it is to prove global optimality as multiple optimal solutions can be found. In the MOS formulation, a symmetry breaking constraint was introduced. This constraint enforces the allocation of non-overlapping operations to sequenced time slots so that symmetric solutions involving allocation to time slots that are apart from each other are eliminated from the feasible

ACS Paragon Plus Environment

Page 15 of 44 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

space and thus expediting the search for the optimal solution. In this work, the fact that timing variables are attributed to individual operations and not to resources involved with the transfer operation will be exploited in the concept we introduce as Enclosing Flexible Time Slot – EFTS, which is used for addressing the possibility of having multiple outflowing operations from the same resource and yet all added flows subjected to a single upper flowrate. The EFTS concept allows for accurately imposing flowrate limits while still given the possibility of each individual operation having its own time grid.

Storage Tank

Operation v1

Charging Tank

Operation v2 CDU

Figure 4 – Illustration of operations in the MOS model.

4. Mathematical formulation In section 4.1, we present the modified version of the MOS model for addressing the operational features detailed in the problem description section, in which a substantial number of constraints have been modified or added. The modified model is able to comply with both operating configurations: namely, LC and FC. Most of the equations are applied for both configurations. When a specific constraint is necessary for addressing one of the operating configurations it will be explicitly indicated in the text, otherwise it is implied that all constraints are used for both configurations. In section 4.2 the EFTS concept is formally defined and the additional equations for addressing it are introduced.

4.1 – Rigid Time Window Approach - RTW In this first approach a rigid time window is adopted for the multiple output operations from the same tank. The details will be discussed along the text. The modified versions of the MOS model rely on the nomenclature presented at the end of the paper. The model is as follows.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 16 of 44

Parcel Unloading In this work, parcels are not unloaded from vessels. Instead, a set of crude parcels are scheduled to arrive at the refinery through a pipeline segment. The exact time at which a parcel must start unloading is known, 𝑆𝑣. Therefore, any parcel or sub-parcel can only start unloading after its arrival, as given in (1). Parcels are unloaded at constant flowrates, 𝐹𝑅𝑣 = 𝐹𝑅𝑣, so that the exact time at which a parcel must finish 𝐿𝑡0𝑟

unloading is also known, 𝐸𝑣 = 𝑆𝑣 + 𝐹𝑅𝑣. Likewise, each parcel or sub-parcel must finish unloading before its pumping is finished, equation (2). Each parcel may be composed of multiple crudes and when split for being unloaded to multiple tanks the composition of each sub-parcel must be the same as that of the initial parcel, equation (3). 𝑆𝑖𝑣 ≥ 𝑆𝑣

∀ 𝑖,𝑣 ∈ 𝑊𝑢 ― 𝑊𝑖𝑣

(1)

𝐸𝑖𝑣 ≤ 𝐸𝑣

∀ 𝑖,𝑣 ∈ 𝑊𝑢 ― 𝑊𝑖𝑣

(2)

𝐿0𝑟𝑐 𝑉𝑖𝑣𝑐 = 𝑉𝑡𝑖𝑣 𝑡 𝐿0𝑟

∀𝑖,𝑐,𝑟 ∈ 𝑅𝑝,𝑣 ∈ 𝑂𝑟 ― 𝑊𝑖𝑣

(3)

It is important to note that each parcel must be unloaded within its pre-defined time window

[𝑆𝑣,𝐸𝑣 ]. Each time interval is non-overlapping to each other, which means that 𝑆𝑣2 > 𝐸𝑣1, 𝑆𝑣3 > 𝐸𝑣2, 𝑒𝑡𝑐. This allows us to consider the unloading operations to be classified as overlapping operations in the sense that the unloading of different parcels can be allocated to the same time slot, but yet not being overlapping in time. Allowing the unloading operations of each parcel to be overlapping requires a lower number of postulated time slots. However, it must be bore in mind that an unloading operation is non-overlapping to itself and thus, if a parcel is split the unloading of each sub-parcel must be allocated to a different time slot and must not be overlapping in time either. In addition, a parcel cannot be unloaded twice to the same tank, constraints (4), which could happen if a parcel is split. Yet, in each priority-slot, a parcel can only be unloaded to one tank. This is translated as all output operations from the same parcel being non-

ACS Paragon Plus Environment

Page 17 of 44 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

overlapping to each other and is handled by equation (19). Finally, in each time slot, only one parcel can be unloaded to a tank, equation (5).

∑𝑍

𝑖𝑣

≤1

∀𝑣 ∈ 𝑊𝑢 ― 𝑊𝑖𝑣

(4)

≤1

∀𝑖, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

(5)

𝑖

∑𝑍

𝑖𝑣

𝑣 ∈ 𝐼𝑟

Volume balances Constraints (6)-(9) are similar to the equations in the original MOS model. Global volume balances for parcels and tanks are given in (6). Note that 𝐿𝑟 and 𝐿𝑟 for the case of tanks represent the tank heel and its storage capacity respectively. For parcels 𝐿𝑟 = 0 and 𝐿𝑟 = 𝐿𝑡0𝑟. Constraint (7) is similar to (6) except that it is written for each crude component. Because (6) and (7) just take into account i' < i, where i is the actual priority slot, (8) e (9) need to be added to account for the last priority slot. Note that 𝐿𝑟 = 0 in (8) and (9) for parcels so that they are enforced to be totally unloaded.

𝐿𝑟 ≤ 𝐿𝑡𝑖𝑟 = 𝐿𝑡0𝑟 +

∑ ∑𝑉

𝑡 𝑖′𝑣

+

𝑖′ < 𝑖𝑣 ∈ 𝐼𝑟

∑∑𝑉

𝑡 𝑖′𝑣

≤ 𝐿𝑟

(6)

𝑖′ < 𝑖𝑣 ∈ 𝑂𝑟

∀𝑖,𝑟 ∈ 𝑅𝑝 ∪ 𝑅𝑡 ― 0 ≤ 𝐿𝑖𝑟𝑐 = 𝐿0𝑟𝑐 +

∑ ∑𝑉

𝑖′𝑣𝑐

+

𝑖′ < 𝑖𝑣 ∈ 𝐼𝑟

∑∑𝑉

𝑖𝑣𝑐

𝑅𝑖𝑡

≤ 𝐿𝑟

(7)

𝑖′ < 𝑖𝑣 ∈ 𝑂𝑟

∀𝑖,𝑐,𝑟 ∈ 𝑅𝑝 ∪ 𝑅𝑡 ― 𝐿𝑟 ≤ 𝐿𝑡𝑖𝑟 = 𝐿𝑡0𝑟 +

∑∑𝑉 + ∑ ∑ 𝑉 𝑡 𝑖𝑣

𝑖 ∈ 𝑇𝑣 ∈ 𝐼𝑟

𝑡 𝑖𝑣

≤ 𝐿𝑟

(8)

𝑖 ∈ 𝑇𝑣 ∈ 𝑂 𝑟

∀𝑟 ∈ 𝑅𝑝 ∪ 𝑅𝑡 ― 0 ≤ 𝐿𝑖𝑟𝑐 = 𝐿0𝑟𝑐 +

∑∑𝑉 𝑖 ∈ 𝑇𝑣 ∈ 𝐼𝑟

𝑅𝑖𝑡

𝑖𝑣𝑐

+

∑∑𝑉

𝑖𝑣𝑐

𝑅𝑖𝑡

(9)

≤ 𝐿𝑟

𝑖 ∈ 𝑇𝑣 ∈ 𝑂 𝑟

∀𝑐,𝑟 ∈ 𝑅𝑝 ∪ 𝑅𝑡 ― 𝑅𝑖𝑡

As mentioned in the introduction and in the problem description sections, fictitious mixers are placed between tanks and CDUs (see figures 2b and 3). The output of a mixer corresponds to the final blended crude fed to a CDU. Therefore, equations (10) and (11) must be included for establishing global

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 18 of 44

and component volume balances at mixers, respectively. Moreover, the total volume of a transfer operation, (12), or contained in a resource, (13), must equal the sum of the volume of all components. Equation (14) eliminates composition discrepancy of outflowing operations at tanks whereas (3) does the same for unloading operations. Note that (14) is a nonlinear nonconvex expression, whereas (3) is linear, given that parcels composition does not change during unloading. Finally, if an operation is not allocated, no volume should be transferred, constraints (15).



𝑉𝑖𝑣𝑐 =

∀𝑖,𝑟 ∈ 𝑅𝑚

(10)

∑𝑉

𝑖𝑣𝑐

∀𝑖,𝑐,𝑟 ∈ 𝑅𝑚

(11)

𝑣 ∈ 𝑂𝑟

𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣

∑𝑉 = ∑𝐿

𝑉𝑡𝑖𝑣 =

𝑡 𝑖𝑣

𝑣 ∈ 𝑂𝑟

𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣



∑𝑉

𝑉𝑡𝑖𝑣 =

𝑖𝑣𝑐

∀𝑖,𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣

(12)

∀𝑖,𝑟 ∈ 𝑅𝑝 ∪ 𝑅𝑡 ― 𝑅𝑖𝑡

(13)

∀𝑖,𝑐,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ∀𝑖,𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣

(14) (15)

𝑐∈𝐶

𝐿𝑡𝑖𝑟

𝑖𝑟𝑐

𝑐∈𝐶

𝐿𝑡𝑖𝑟𝑉𝑖𝑣𝑐 = 𝑉𝑡𝑖𝑣𝐿𝑖𝑟𝑐 𝑉𝑡𝑣𝑍𝑖𝑣 ≤ 𝑉𝑡𝑖𝑣 ≤ 𝑉𝑡𝑣 𝑍𝑖𝑣

Quality constraints The quality constraints are imposed for the transfer operations between mixers and CDUs (𝑣 ∈ 𝑊𝐷 ). There is only one operation for each mixer-CDU pair. Constraint (16) imposes crude oil blending specifications for properties obeying volumetric mixing rules. Likewise, (17) is the equivalent for properties obeying mass mixing rules. Constraint (18) is applied for setting an upper limit for the sulfur content in the diesel cut produced at CDUs, which follows a mass mixing rule. Note that all these constraints are linear, since the only unknowns are 𝑉𝑡𝑖𝑣 and 𝑉𝑖𝑣𝑐.

ACS Paragon Plus Environment

Page 19 of 44 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

𝑥𝑘.𝑉𝑡𝑖𝑣 ≤

∑𝑥

𝑐𝑘𝑉𝑖𝑣𝑐

≤ 𝑥𝑘.𝑉𝑡𝑖𝑣

(16)

𝑐∈𝐶

∀ 𝑖, 𝑣 ∈ 𝑊𝑑,𝑘 ∈ 𝐾𝑣 𝑥𝑘

∑𝑉

𝑘′ 𝑖𝑣𝑐𝑥𝑐

∑𝑥



𝑐∈𝐶

𝑘′ 𝑐𝑘𝑉𝑖𝑣𝑐𝑥𝑐

𝑐∈𝐶

≤ 𝑥𝑘

∑𝑉

𝑘′ 𝑖𝑣𝑐𝑥𝑐

(17)

𝑐∈𝐶

∀ 𝑖,𝑣 ∈ 𝑊𝑑,𝑘 ∈

∑𝑉

𝑘 𝑘′ 𝑘′′ 𝑖𝑣𝑐𝑥𝑐 𝑥𝑐 𝑥𝑐

𝐷

≤ 𝑥𝑘

𝑐∈𝐶

∑𝑉

𝐾𝑚,𝑘′

= 𝐴𝑃𝐼 𝑔𝑟𝑎𝑣𝑖𝑡𝑦

𝑘′ 𝑘′′ 𝑖𝑣𝑐𝑥𝑐 𝑥𝑐

(18)

𝑐∈𝐶

∀ 𝑖,𝑣 ∈ 𝑊𝑑 𝑘 = 𝑠𝑢𝑙𝑓𝑢𝑟,𝑘′ = 𝐴𝑃𝐼 𝑔𝑟𝑎𝑣𝑖𝑡𝑦, 𝑘′′ = 𝑦𝑖𝑒𝑙𝑑

Logical Constraints Non-overlapping operations form a clique and thus cannot be allocated to the same priority slot, constraint (19). Constraint (20) is a symmetry breaking constraint in which symmetric solutions are avoided by enforcing allocation of non-overlapping operations to sequencing time slots. Constraint (21) is particularly useful for determining the minimum number of time slots. It enforces the allocation of at least one operation to each time slot. All these constraints retain the same form as the ones in the original MOS model.

∑𝑍

𝑖𝑣

≤ 1 ∀𝑖,𝑊′ ∈ 𝐶𝑙𝑖𝑞𝑢𝑒(𝐺𝑁𝑂)

(19)



(20)

𝑣 ∈ 𝑊′

𝑍𝑖𝑣 ≤

𝑍(𝑖 ― 1)𝑣′ ∀𝑖,𝑖 ≠ 1, 𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣

𝑣′ ∈ 𝑊 ― 𝑊𝑖𝑣 𝑁𝑂𝑣𝑣′ = 1

∑ 𝑣∈𝑊―

(21)

𝑍𝑖𝑣 ≥ 1 ∀𝑖 𝑊𝑖𝑣

For the FC configuration, any two active tanks are allowed to feed the same CDU, so that (22a) is a valid constraint. As for the LC, only one base tank can feed a CDU, (22b), and an additional injection

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 44

tank, (22c). Each tank is allowed to serve two CDUs simultaneously, as given in (23).



𝑍𝑖𝑣 ≤ 2

∀ 𝑖, 𝑟 ∈ 𝑅𝑚|FC configuration only

(22a)

∀ 𝑖, 𝑟 ∈ 𝑅𝑚|LC configuration only

(22b)

∀ 𝑖, 𝑟 ∈ 𝑅𝑚|LC configuration only

(22c)

𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣



𝑍𝑖𝑣 ≤ 1

𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣 ― 𝑊𝑖𝑛𝑗 𝑣



𝑣 ∈ 𝐼𝑟 ∩ 𝑊𝑖𝑛𝑗 𝑣

∑𝑍



𝑍𝑖𝑣 ≤

𝑖𝑣

𝑍𝑖𝑣

𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣 ― 𝑊𝑖𝑛𝑗 𝑣

(23)

≤ 2 ∀ 𝑖, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

𝑣 ∈ 𝑂𝑟

Timing constraints All operations are subject to the following two constraints. The end of an operation is given by the instant of time at which it begins plus its duration, (24). Operations that are not active must have its related timing variables set to zero, (25). These equations are maintained as in the original formulation. Constraint (26) is introduced to satisfy the requirement that unloading operations have a minimum run-length of 3 hours, whereas distillation operations have a minimum run-length of 24 hours. The latter is presumably a very restrictive constraint, which has the purpose of avoiding frequent changes to CDUs load and thus minimizing upset to their operation, in which case off-spec products could be produced. In this work no attempt was made as to identify if the load of a CDU has changed between two consecutive time slots. In case of not having changes, (26) could be relaxed.

𝐸𝑖𝑣 = 𝑆𝑖𝑣 + 𝐷𝑖𝑣 𝐸𝑖𝑣 ≤ H.𝑍𝑖𝑣

∀ 𝑖,𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣 ∀ 𝑖,𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣

(24) (25)

𝐷𝑖𝑣 ≥ 𝐷𝑣𝑍𝑖𝑣

∀ 𝑖,𝑣 ∈ 𝑊 ― 𝑊𝑖𝑣

(26)

ACS Paragon Plus Environment

Page 21 of 44

New variables are also introduced to define the time window within which outflowing operations are active at a tank, as illustrated in Figure 5. Variables 𝑆𝑡𝑖,𝑟, 𝐷𝑡𝑖,𝑟 and 𝐸𝑡𝑖,𝑟 represent the start, duration and end of the enclosing time window. Note that these variables are indexed by resources instead of operations. Thus, each tank has its own enclosing time window. Constraints (27) and (28) are the equivalent to (24) and (25) stated in terms of the tank variables. Note also that theses variables and constraints are applied only for tanks that can feed two distillation units simultaneously (see figure 3).

𝑆𝑡𝑖𝑟 𝐸𝑡𝑖𝑟 𝐷𝑡𝑖𝑟

v3

𝑆𝑡𝑖𝑟 CDU 1

𝑆𝑖𝑣 𝐸𝑖𝑣 𝐷𝑖𝑣 v1

v4

v2

𝑆𝑖𝑣1

𝐸𝑡𝑖𝑟

𝐷𝑡𝑖𝑟 𝐷𝑖𝑣1

CDU 2

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

𝐸𝑖𝑣1 𝑆𝑖𝑣2

𝐷𝑖𝑣2

𝐸𝑖𝑣2

Time Slot i

Figure 5 – Illustration of the time window that encloses outflowing operations from a tank in time slot i. 𝐸𝑡𝑖,𝑟 = 𝑆𝑡𝑖,𝑟 + 𝐷𝑡𝑖,𝑟 𝐸𝑡𝑖,𝑟 ≤ 𝐻

∑𝑍

𝑖,𝑣

∀𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(27) (28)

𝑣 ∈ 𝑂𝑟

Now the timing variables of outflowing operations from tanks must be related to the tank timing variables through (29-32), which enforce all outputs to be synchronized leading to the behavior illustrated in Figure 6a. Note, however, that the solution presented in Figure 6b would be feasible by increasing the number of time slots. Although limiting to some extent, this was though the first approach adopted for being able to accurately satisfying lower and upper flowrate limits with multiple outputs and yet keeping flowrate constraints linear. In the next section, a discussion is conducted in which (29-32) are relaxed and the solution in Figure 5 is feasible using just a single time slot and yet without adding any nonlinear constraint for accurately satisfying lower and upper flowrate limits. 𝑆𝑖𝑣 ≤ 𝑆𝑡𝑖𝑟 + 𝐻(1 ― 𝑍𝑖𝑣)

∀ 𝑖,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

ACS Paragon Plus Environment

(29)

Industrial & Engineering Chemistry Research 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 ― 𝑍𝑖𝑣) 𝐸𝑖𝑣 ≤ 𝐸𝑡𝑖𝑟 + 𝐻(1 ― 𝑍𝑖𝑣) 𝐸𝑖𝑣 ≥ 𝐸𝑡𝑖𝑟 ― 𝐻(1 ― 𝑍𝑖𝑣)

Page 22 of 44

∀ 𝑖, 𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 Tank window

Tank windows in different time slots

v1

v1

(30) (31) (32)

v2

v1 v2

Time Slot i (a)

Time Slot i-1

Time Slot i

v2 Time Slot i+1

(b)

Figure 6 – Solutions found with the frozen time windows. Constraints (33-36) synchronizes start and end times among tanks and mixers, guaranteeing that all operations coming from different tanks have the same duration as that of the CDU. Moreover, once the unloading operations are not considered overlapping, (20) cannot guarantee that unloading be done using time slots in an orderly fashion. This could cause equations (6) and (7) to not properly account for the inventory in tanks given that these equations rely on the fact that time slots are used in increasing order. Therefore, (37) must be added for guaranteeing that multiple loading operations to the same tank be ordered in terms of time slots and in terms of time. Note that (37) can only be written because in each time slot: a) there is only one tank inlet operation active, (5); b) the introduced tank variables account for the multiple outlet operations for the case of multiple outputs (37a), and c) a tank cannot load and unload simultaneously, which is addressed by (19). (37a) applies for tanks with multiple outputs whereas (37b) applies for tanks that connect to a single CDU. Constraints (38-39) are used for determining if there is an outflow operation at the tank, whereas (40-41) determine if the tank has received a load.

𝑆𝑖𝑣 ≤ 𝑆𝑖𝑣 ≥ 𝐸𝑖𝑣 ≤ 𝐸𝑖𝑣 ≥

𝑆𝑖𝑣′ + 𝐻(1 ― 𝑍𝑖𝑣) 𝑆𝑖𝑣′ ― 𝐻(1 ― 𝑍𝑖𝑣) 𝐸𝑖𝑣′ + 𝐻(1 ― 𝑍𝑖𝑣) 𝐸𝑖𝑣′ ― 𝐻(1 ― 𝑍𝑖𝑣)

∀ 𝑖,𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣′ ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚 ∀ 𝑖,𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣′ ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚 ∀ 𝑖,𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣′ ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚 ∀ 𝑖,𝑣 ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣′ ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚

ACS Paragon Plus Environment

(33) (34) (35) (36)

Page 23 of 44 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

)

𝑣 ∈ 𝐼𝑟

∑𝑆



𝑖2𝑣

(



𝐸𝑖1𝑣 + 𝐸𝑡𝑖𝑟 + 𝐵𝑆.𝑌𝑇𝑖1𝑟 +

𝑖1 < 𝑖 < 𝑖2

(

+ 𝑆𝑡𝑖𝑟 + 𝐻 1 ―

𝑣 ∈ 𝐼𝑟

∑𝐷

𝐵𝑆.𝑌𝑇𝑖𝑟 +

𝑣 ∈ 𝐼𝑟

∑𝑍

𝑖2𝑣

)

𝑖𝑣



+

𝐷𝑡𝑖𝑟

(37a)

𝑖1 < 𝑖 < 𝑖2

)

― 𝑍𝑡𝑖𝑟

𝑣 ∈ 𝐼𝑟

∀𝑖1 < 𝑖2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(∑

𝐸𝑖1𝑣 +

𝑣 ∈ 𝐼𝑟

∑ 𝐸 + 𝐵𝑆.𝑌𝑇 ) + ∑ (𝐵𝑆.𝑌𝑇 + ∑ 𝐷 ) + ∑ ∑ 𝐷 ≤ ∑ 𝑆 + ∑ 𝑆 + 𝐻(1 ― ∑ 𝑍 ― ∑ 𝑍 ) 𝑖1𝑣

𝑖1𝑟

𝑖𝑟

𝑣 ∈ 𝑂𝑟

𝑖1 < 𝑖 < 𝑖2

𝑖2𝑣

𝑣 ∈ 𝐼𝑟

𝑖𝑣

𝑣 ∈ 𝐼𝑟

𝑖2𝑣

𝑖2𝑣

𝑣 ∈ 𝑂𝑟

𝑣 ∈ 𝐼𝑟

𝑖𝑣

𝑖1 < 𝑖 < 𝑖2𝑣 ∈ 𝑂𝑟

(37b)

𝑖2𝑣

𝑣 ∈ 𝑂𝑟

∀𝑖1 < 𝑖2,𝑟 ∈ 𝑅𝑂𝑛𝑒 ― 𝑅𝑖𝑡 𝑡 𝑍𝑡𝑖𝑟 ≥ 𝑍𝑖𝑣

∀ 𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

∑𝑍 ≥ ∑𝑍 ― ∑𝑍 ≤ ∑𝑍

𝑍𝑡𝑖𝑟 ≤

(38) (39)

∀ 𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

𝑖𝑣

𝑣 ∈ 𝑂𝑟

𝑌𝑇𝑖𝑟

𝑖𝑣

𝑣 ∈ 𝐼𝑟

𝑌𝑇𝑖𝑟

∀ 𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

(40)

∀ 𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

(41)

𝑖 + 1𝑣

𝑣 ∈ 𝐼𝑟

𝑖𝑣

𝑣 ∈ 𝐼𝑟

As already mentioned, non-overlapping operations cannot be allocated to the same priority slot, nor can their time spans overlap, which is avoided through (42). In our application, (42) includes only unloading and distillation operations. Tank related operations are dealt with (37). Constraint (42) could also include tank operations although regarding non-overlapping unloading operations (37) would still be required. Therefore, we decided to exclude tank operations from (42) and merge those decisions into (37). Finally, CDUs must operate continuously, (43).

∑𝐸

𝑖1𝑣

𝑣 ∈ 𝑊′

+

∑ ∑𝐷 𝑖 ∈ 𝑇 𝑣 ∈ 𝑊′ 𝑖1 < 𝑖 < 𝑖2

𝑖1𝑣



∑𝑆 𝑣′ ∈ 𝑊′

𝑖2𝑣

(

+𝐻 1―

∑𝑍 𝑣 ∈ 𝑊′

)

𝑖2𝑣

∀ 𝑖1,𝑖2 ∈ 𝑇, 𝑖1 < 𝑖2,𝑊′ ∈ 𝐶𝑙𝑖𝑞𝑢𝑒(𝐺𝑁𝑂)

ACS Paragon Plus Environment

(42)

Industrial & Engineering Chemistry Research 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 24 of 44

(43)

𝑣 ∈ 𝑊𝑑

𝑖∈𝑇

Flowrate Limits Imposing upper and lower flowrate limits for simultaneous outputs in a multiple time grid model is by far the most challenging feature to be addressed in the present work. Multiple outputs are allowed only for tanks and when multiple outputs are active there should be upper and lower bounds for the sum of all outflows. One might think that having a constraint like (A) would suffice. Indeed, there are two problems with such a constraint; 1) it introduces nonlinear fractional terms, which have the additional problem of division by zero for the unassigned operations and, at the same time, it makes the employment of solution algorithms that rely on solving MILP subproblems (derived from the full scale MINLP), prohibitive, and 2) because of the multiple time grid nature, the simple addition of (A) would not prevent us from having an upper flowrate violation such as the one illustrated in Figure 7, in which the assumed upper bound would be 200 m3/h. Note that the flowrate upper bound is violated in the grayed area where the overlapping occurs. It should though be noticed that the latter situation is prevented with the introduction of (37), as we propose. However, the former point would still not be addressed. Therefore, the addition of (29-36) seems a reasonable appropriate measure for addressing multiple outputs. Unloading and distillation operations flowrate are bounded through (44), whereas (45) bounds output flowrates at tanks. Note that (45) remains linear in this case. 𝐹𝑅𝑟 ≤

𝑉𝑡𝑖𝑣

∑𝐷

𝑣 ∈ 𝑂𝑟

𝑖𝑣

≤ 𝐹𝑅𝑟

∀𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡

ACS Paragon Plus Environment

(A)

Page 25 of 44 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

Overlapping of two sequenced time windows of the same tank VTi1,v1 = 100 m3, Di1,v1 = 1 h

VTi2,v1 = 150 m3, Di2,v1 = 1 h

VTi1,v2 = 100 m3, Di1,v2 = 1 h

t

S i1,Tank1

v1 v2

Operating time window of tank 1 in time slot i1

t

E i1,Tank1

t

S i2,Tank1 Operating time window of tank 1 in time slot i2

t

E i2,Tank1

Figure 7 – Illustration of upper flowrate violation involving multiple tank output in multiple time grid represented models. 𝐹𝑅𝑣𝐷𝑖𝑣 ≤ 𝑉𝑡𝑖𝑣 ≤ 𝐹𝑅𝑣 𝐷𝑖𝑣 ∀𝑖,𝑣 ∈ 𝑊𝑢 ∪ 𝑊𝑑 𝐹𝑅𝑟𝐷𝑡𝑖𝑟 ≤

∑𝑉

𝑡 𝑖𝑣

(44) (45)

≤ 𝐹𝑅𝑟 𝐷𝑡𝑖𝑟 ∀𝑖,𝑟 ∈ 𝑅𝑡

𝑣 ∈ 𝑂𝑟

If a tank is assigned to feed a distillation unit its transferred volume cannot be null (5% is the attributed minimal value), (46). For the case of the LC configuration, the volume transferred from injection tanks cannot be higher than 30% of the total unit load, (47). 𝑉𝑡𝑖𝑣′ ≥ 0.05.𝑉𝑡𝑖𝑣 ― 𝑉𝑡𝑣 (1 ―𝑍𝑖𝑣′) 𝑉𝑡𝑖𝑣′ ≤ 0.30.𝑉𝑡𝑖𝑣 + 𝑉𝑡𝑣 (1 ―𝑍𝑖𝑣′)

∀ 𝑖,𝑖𝑠,𝑣′ ∈ 𝐼𝑟 ― 𝑊𝑖𝑣′, 𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚 𝑖 ∀ 𝑖,𝑖𝑠,𝑣′ ∈ 𝐼𝑟 ∩ 𝑊𝑖𝑛𝑗 𝑣 ― 𝑊𝑣′, 𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚

(46) (47)

4.2 - Enclosing Flexible Time Slot – EFTS In this section, focus is given on modifying the model presented in the previous section for allowing the possibility of having two or more output operations from the same tank each with its own time grid as a true multiple time grid model. Although the RTW model offers the possibility of having a solution like the one shown in Figure 6b by increasing the number of time slots, it is still too restrictive due to (26), which imposes a minimum run-length for each time slot. Thus, we want to be able to have the

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 26 of 44

flexibility shown in Figure 6 and still accurately set upper and lower bounds for the total flowrate. The solution is to sub-divide the tank active time window in three time sub-slots as illustrated in Figure 8. As can be seen by the picture, all variables related to tank outputs are disaggregated into variables representing the operation taking place in each sub-slot. Thus, if the operation takes place in a sub-slot its disaggregated variables are non-null, and they are null otherwise. The length of each sub-slot is defined by the tank enclosing time window, so that 𝐷𝑡𝑖𝑟 = ∑𝑖𝑠∆𝑇𝑖,𝑖𝑠,𝑟 ≥ 𝐷𝑣. Nevertheless, it is possible to have ∆𝑇𝑖,𝑖𝑠,𝑟 ≤ 𝐷𝑣. If no operations are assigned to a sub-slot then its duration must be null, ∆𝑇𝑖,𝑖𝑠,𝑟 = 0. The length of an operation taking place in a sub-slot must be equal to the length of a sub-slot of the tank enclosing window, ∆𝑇𝑖,𝑖𝑠,𝑟 = 𝐷𝑠𝑖,𝑖𝑠,𝑣⟺𝑍𝑠𝑖,𝑖𝑠,𝑣 = 1. Once those conditions are imposed, flowrates can accurately be estimated within each time sub-slot through a linear expression that takes into account the duration of the enclosing time window sub-slot, ∆𝑇𝑖,𝑖𝑠,𝑟, similar to (45). It is not allowed to have all active operations taking place in the first and last sub-slots, so that asynchronous outflows are captured as in Figure 8. Besides, by imposing such a restriction helps preventing symmetric solutions spanning different number of sub-slots when the transfer operations are synchronized. In our problem under study, in which ∑𝑣 ∈ 𝑂 𝑍𝑖𝑣 ≤ 2, 𝑟 ∈ 𝑅𝑡, only one operation can take 𝑟

place in the first and last sub-slots, ∑𝑣 ∈ 𝑂 𝑍𝑠𝑖,𝑖𝑠,𝑣 ≤ 1, 𝑟 ∈ 𝑅𝑡, 𝑖𝑠 ∈ {1,3}. An operation starting in the first 𝑟

sub-slot cannot go onto the third sub-slot and neither can it start on the third sub-slot. Likewise, an operation cannot end on the first sub-slot. The number of active operations in the second sub-slot though must be equal to the number of active tank outputs, ∑𝑣 ∈ 𝑂 𝑍𝑠𝑖,𝑖𝑠,𝑣 = ∑𝑣 ∈ 𝑂 𝑍𝑖𝑣, 𝑟 ∈ 𝑅𝑡, which enforces the 𝑟

𝑟

operations to overlap to some extent. Accordingly, if two operations are assigned to start and end at the same instants of time, like in the RTW model, only the second sub-slot is used plus, 𝐷𝑡𝑖𝑟 = ∆𝑇𝑖,𝑖𝑠2,𝑟 ≥ 𝐷𝑣. Therefore, one can state that the feasible region of the RTW model is contained in the feasible region of the EFTS model, 𝑅𝑇𝑊 ⊂ 𝐸𝐹𝑇𝑆. In other words, all solutions of the RTW model are feasible solutions for

ACS Paragon Plus Environment

Page 27 of 44 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

the EFTS model, but the opposite is not true. The following constraints establish the details just described as well as the relationship between the global and disaggregated variables. Most of the constraints of the RTW model are also applied to the EFTS model. Those which are not applied are indicated in the text. Sub-slot is1 𝑆𝑡𝑖𝑟

𝑍𝑠𝑠𝑖,𝑖𝑠1,𝑣1 𝑆𝑖𝑣1

𝑆𝑠𝑖𝑟

Sub-slot is2

∆𝑇𝑖,𝑖𝑠1,𝑟

∆𝑇𝑖,𝑖𝑠2,𝑟

𝑍𝑠𝑖,𝑖𝑠1,𝑣1 𝐷𝑠𝑖,𝑖𝑠1,𝑣1

𝑍𝑠𝑖,𝑖𝑠2,𝑣1 𝐷𝑠𝑖,𝑖𝑠2,𝑣1 𝑉𝑠𝑡𝑖,𝑖𝑠2,𝑣1 𝑉𝑠𝑖,𝑖𝑠2,𝑣1,𝑐

𝑉𝑠𝑡𝑖,𝑖𝑠1,𝑣1 𝑉𝑠𝑖,𝑖𝑠1,𝑣1,𝑐 𝐷𝑖𝑣1

𝑍𝑠𝑠𝑖,𝑖𝑠2,𝑣2

𝐸𝑠𝑖𝑟

Sub-slot is3 ∆𝑇𝑖,𝑖𝑠3,𝑟

𝐸𝑡𝑖𝑟

𝑍𝑠𝑓𝑖,𝑖𝑠2,𝑣1 𝐸𝑖𝑣1 𝑍𝑠𝑖,𝑖𝑠3,𝑣2 𝐷𝑠𝑖,𝑖𝑠3,𝑣2 𝑉𝑠𝑡𝑖,𝑖𝑠3,𝑣2

𝑍𝑠𝑖,𝑖𝑠2,𝑣2 𝐷𝑠𝑖,𝑖𝑠2,𝑣2 𝑉𝑠𝑡𝑖,𝑖𝑠2,𝑣2 𝑉𝑠𝑖,𝑖𝑠2,𝑣2,𝑐

𝑉𝑠𝑖,𝑖𝑠3,𝑣2,𝑐

𝑆𝑖𝑣2

𝑍𝑠𝑓𝑖,𝑖𝑠3,𝑣2 𝐸𝑖𝑣2

𝐷𝑖𝑣2 Time Slot i, 𝐷𝑡𝑖𝑟

Figure 8 – Enclosing flexible time slots and related variables. Logical Constraints Constraint (48) states that the number of active operations in the second sub-slot is equal to the number of active output operations of a tank, whereas (49) establish a maximum to the number of active operations in the first and third sub-slots based on the number of active output operations. Note that if only one output operation is active, it is assigned to the second sub-slot and the first and third sub-slots will be empty. Moreover, an operation is allowed to span only two time sub-slots, (50). It is noteworthy that constraints presented in this section are applied only to the active tanks that present multiple connections to CDUs, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 . If a tank has only one connection, there is no need for applying the disaggregation.

∑ 𝑍𝑠 𝑣 ∈ 𝑂𝑟

𝑖,𝑖𝑠,𝑣

=

∑𝑍

𝑖𝑣

∀𝑖, 𝑖𝑠 = 2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

𝑣 ∈ 𝑂𝑟

ACS Paragon Plus Environment

(48)

Industrial & Engineering Chemistry Research 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 28 of 44

∀𝑖, 𝑖𝑠 ≠ 2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(49)

∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(50)

𝑣 ∈ 𝑂𝑟

≤2

𝑖𝑠

An operation must start on the first or second sub-slots, (51-53), and finish on the second or third sub-slot, (54-56).

∑ 𝑍𝑠

𝑍𝑠𝑠𝑖,𝑖𝑠.𝑣 ≥ 𝑍𝑠𝑖,𝑖𝑠,𝑣 ―

𝑖,𝑖𝑠′,𝑣

∀𝑖, 𝑖𝑠 < 3,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(51)

𝑖𝑠′ < 𝑖𝑠

𝑍𝑠𝑠𝑖,𝑖𝑠.𝑣

≤ 𝑍𝑠𝑖,𝑖𝑠,𝑣

∑𝑍𝑠

= 𝑍𝑖,𝑣

𝑠 𝑖,𝑖𝑠.𝑣

∀𝑖, 𝑖𝑠 < 3,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(52) (53)

𝑖𝑠

∑ 𝑍𝑠

𝑍𝑠𝑓𝑖,𝑖𝑠.𝑣 ≥ 𝑍𝑠𝑖,𝑖𝑠,𝑣 ―

𝑖,𝑖𝑠′,𝑣

∀𝑖, 𝑖𝑠 > 1,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(54)

𝑖𝑠′ > 𝑖𝑠

𝑍𝑠𝑓𝑖,𝑖𝑠.𝑣 ≤ 𝑍𝑠𝑖,𝑖𝑠,𝑣



𝑍𝑠𝑓𝑖,𝑖𝑠.𝑣

= 𝑍𝑖,𝑣

∀𝑖, 𝑖𝑠 > 1,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ―

𝑅𝑖𝑡



𝑅𝑂𝑛𝑒 𝑡

(55) (56)

𝑖𝑠

Timing constraints The duration of an operation must be equal to the sum of its duration over the multiple sub-slots (57). Likewise, the duration of the enclosing time window is equal to the sum of the durations of each sub-slot, (58). Moreover, if a task is not active in a sub-slot its corresponding duration should be null, (59). In the same way, if no operations are active in a sub-slot, the tank sub-slot duration must be null, (60).

∑𝐷𝑠 = ∑∆𝑇

𝐷𝑖𝑣 =

𝑖,𝑖𝑠,𝑣

∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(57)

𝑖,𝑖𝑠,𝑟

∀𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(58)

∀𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 𝑖 𝑂𝑛𝑒 ∀𝑖,𝑖𝑠,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑡 ― 𝑅𝑡

(59) (60)

𝑖𝑠

𝐷𝑡𝑖𝑟

𝑖𝑠

𝐷𝑠𝑖,𝑖𝑠,𝑣 ≤ 𝐻.𝑍𝑠𝑖,𝑖𝑠,𝑣 ∆𝑇𝑖,𝑖𝑠,𝑟 ≤ 𝐻.𝑍𝑖,𝑣

In the EFTS model, constraints (29) and (32) are replaced by (61) and (62) to indicate that an operation is started in the first sub-slot and finished in the third sub-slot, respectively, while (30) and (31) are still valid constraints regardless of the sub-slot an operation starts or finishes. Furthermore, (63) and

ACS Paragon Plus Environment

Page 29 of 44 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

(64) indicate the instant of time at which an operation starts, whereas (65) and (66) indicate the instant of time at which an operation finishes, respectively, if assigned to start and finish in the second sub-slot. 𝑆𝑖𝑣 ≤ 𝑆𝑡𝑖𝑟 + 𝐻(1 ― 𝑍𝑠𝑠𝑖,𝑖𝑠,𝑣) 𝐸𝑖𝑣 ≥ 𝐸𝑡𝑖𝑟 ― 𝐻(1 ― 𝑍𝑠𝑓𝑖,𝑖𝑠,𝑣)

∀ 𝑖,𝑖𝑠 = 1, 𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 𝑖 𝑂𝑛𝑒 ∀ 𝑖,𝑖𝑠 = 3,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑡 ― 𝑅𝑡

(61) (62)

𝑆𝑖𝑣 ≤ 𝑆𝑠𝑖𝑟 + 𝐻(1 ― 𝑍𝑠𝑠𝑖,𝑖𝑠,𝑣) 𝑆𝑖𝑣 ≥ 𝑆𝑠𝑖𝑟 ― 𝐻(1 ― 𝑍𝑠𝑠𝑖,𝑖𝑠,𝑣)

∀ 𝑖,𝑖𝑠 = 2,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 𝑖 𝑂𝑛𝑒 ∀ 𝑖,𝑖𝑠 = 2,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑡 ― 𝑅𝑡

(63) (64)

𝐸𝑖𝑣 ≤ 𝐸𝑠𝑖𝑟 + 𝐻(1 ― 𝑍𝑠𝑠𝑖,𝑖𝑠,𝑣) 𝐸𝑖𝑣 ≥ 𝐸𝑠𝑖𝑟 ― 𝐻(1 ― 𝑍𝑠𝑠𝑖,𝑖𝑠,𝑣)

∀ 𝑖,𝑖𝑠 = 2,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑖𝑠 = 2,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(65) (66)

Constraints (67-69) are applied for splitting the enclosing time window in three time sub-slots, whereas (70) assures that there will not exist a sub-slot with a null duration if at least an operation is active in that sub-slot. 𝑆𝑠𝑖𝑟 = 𝑆𝑡𝑖𝑟 + ∆𝑇𝑖,𝑖𝑠,𝑟 𝐸𝑠𝑖𝑟 = 𝐸𝑡𝑖𝑟 ― ∆𝑇𝑖,𝑖𝑠,𝑟 𝐸𝑠𝑖𝑟 ― 𝑆𝑠𝑖𝑟 = ∆𝑇𝑖,𝑖𝑠,𝑟 ∆𝑇𝑖,𝑖𝑠,𝑟 ≥ ∆𝑇.𝑍𝑠𝑖,𝑖𝑠,𝑣

∀ 𝑖,𝑖𝑠 = 1, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑖𝑠 = 3, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑖𝑠 = 2, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(67) (68) (69) (70)

If an operation is assigned to a sub-slot, its duration must equal that of the enclosing time window corresponding to the same sub-slot, (71-72). Although not mandatory, (73) is also introduced as a valid constraint that ends up expediting solution time. 𝐷𝑠𝑖,𝑖𝑠,𝑣 ≤ ∆𝑇𝑖𝑟 + 𝐻(1 ― 𝑍𝑠𝑖,𝑖𝑠,𝑣) 𝐷𝑠𝑖,𝑖𝑠,𝑣 ≥ ∆𝑇𝑖𝑟 ― 𝐻(1 ― 𝑍𝑠𝑖,𝑖𝑠,𝑣)

∑𝐷

𝑖𝑣

≤ 𝐷𝑡𝑖𝑟

∀ 𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀ 𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀𝑖,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(71) (72) (73)

𝑣 ∈ 𝑂𝑟

Once the global variables are disaggregated for tank output operations, some of the decisions involving CDU feeding must also be disaggregated. The decisions that we are particularly interested are the inlet flowrate and the qualities in each sub-slot. For determining flowrate, sub-slot durations of the

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 30 of 44

CDU inflowing operation must coincide with those of the tanks assigned for feeding the unit, which is established through (74 -77). 𝐷𝑚𝑖,𝑖𝑠,𝑣 ≤ 𝐷𝑠𝑖,𝑖𝑠,𝑣′ + 𝐻(1 ― 𝑍𝑠𝑖,𝑖𝑠,𝑣′) 𝐷𝑚𝑖,𝑖𝑠,𝑣 ≥ 𝐷𝑠𝑖,𝑖𝑠,𝑣′ ― 𝐻(1 ― 𝑍𝑠𝑖,𝑖𝑠,𝑣′)



𝐷𝑚𝑖,𝑖𝑠,𝑣 ≤ 𝐻

𝑍𝑠𝑖,𝑖𝑠,𝑣′

∀ 𝑖,𝑖𝑠,𝑣′ ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣 ∈ 𝑂𝑟, 𝑟 ∈ 𝑅𝑚 ∀ 𝑖,𝑖𝑠,𝑣′ ∈ 𝐼𝑟 ― 𝑊𝑖𝑣, 𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚 ∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚

(74) (75) (76)

∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚

(77)

𝑣′ ∈ 𝐼𝑟 ― 𝑊𝑖𝑣

𝐷𝑖𝑣 =

∑𝐷𝑚

𝑖,𝑖𝑠,𝑣

𝑖𝑠

Volume Balances and Flowrate The total volume of an operation is disaggregated into the volumes transferred in each sub-slot, (78), and the total volume in each sub-slot is disaggregated into the components, (79). If an operation is not active in a sub-slot, its associated volume must be null, (80). 𝑉𝑡𝑖𝑣 =

∑𝑉𝑠 = ∑𝑉

𝑡 𝑖,𝑖𝑠,𝑣

∀𝑖,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(78)

∀𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(79)

∀𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(80)

𝑖𝑠

𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣 𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣



𝑐 𝑉𝑡𝑣

𝑖,𝑖𝑠,𝑣,𝑐

𝑍𝑠𝑖,𝑖𝑠,𝑣

The upper bound for the flowrate in the first and third sub-slots is imposed through (81), whereas (82) is used for the second sub-slot. The lower bound for the flowrate in the first and third sub-slots is imposed through (83), whereas (84) and (85) are used for the second sub-slot. The maximum processing capacity of each CDU is guaranteed in each sub-slot through (86) 𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣 ≤ 𝐹𝑅𝑣.𝐷𝑠𝑖,𝑖𝑠,𝑣

∑ 𝑉𝑠

𝑡 𝑖,𝑖𝑠,𝑣

≤ 𝐹𝑅𝑟.∆𝑇𝑖,𝑖𝑠,𝑟

∀𝑖,𝑖𝑠 ≠ 2,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡 ∀𝑖,𝑖𝑠 = 2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(81) (82)

𝑣 ∈ 𝑂𝑟

𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣 ≥ 𝐹𝑅𝑣.𝐷𝑠𝑖,𝑖𝑠,𝑣

∀𝑖,𝑖𝑠 ≠ 2,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(83)

𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣 ≥ 1/8.𝐹𝑅𝑣.𝐷𝑠𝑖,𝑖𝑠,𝑣

∀𝑖,𝑖𝑠 = 2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(84)

∑ 𝑉𝑠

∀𝑖,𝑖𝑠 = 2,𝑟 ∈ 𝑅𝑡 ― 𝑅𝑖𝑡 ― 𝑅𝑂𝑛𝑒 𝑡

(85)

𝑡 𝑖,𝑖𝑠,𝑣

≥ 𝐹𝑅𝑟.∆𝑇𝑖,𝑖𝑠,𝑟

𝑣 ∈ 𝑂𝑟

ACS Paragon Plus Environment

Page 31 of 44 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



𝑉𝑠𝑡𝑖,𝑖𝑠,𝑣′ ≤ 𝐹𝑅𝑟.𝐷𝑚𝑖,𝑖𝑠,𝑣

∀𝑖,𝑖𝑠,𝑣 ∈ 𝑂𝑟,𝑟 ∈ 𝑅𝑚

(86)

𝑣′ ∈ 𝐼𝑟 ― 𝑊𝑖𝑣

Additional constraints Besides the previous constraints presented in this section, constraints (14), (16), (17), (18), (46) and (47) must also be written for each sub-slot and will not be presented here for sake of readability and to avoid a too lengthy text.

4.3 - Objective Function In both models RTW and EFTS, the objective function is to maximize the gross margin given in (87) as proposed by Mouret et al.13, where GMc is the gross margin of crude c. max

∑ ∑ ∑𝐺𝑀 .𝑉 𝑐

(87)

𝑖,𝑣,𝑐

𝑖 𝑣 ∈ 𝑊𝑑 𝑐

5. Solution algorithms All full scale formulations presented in section 4 result in nonconvex MINLP problems, which are usually hard to solve. In this work, two solution algorithms were used for solving the obtained models that avoid solving the full scale problems, namely: the Two-Step algorithm proposed by Mouret et al.13, and the Rolling-Horizon Solution Algorithm (RHS) proposed by Reddy et al.8. Both algorithms are detailed next. 5.1 The Two-Step Algorithm The Two-Step algorithm is summarized in Figure 9. It consists of solving an MILP problem in the first step using a Branch-and-Bound based algorithm (e.g., CPLEX). The linear problem is derived from the full scale MINLP by dropping the nonlinear constraints (14). In the second step, the integer variables

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 32 of 44

found in the first step are fixed in the original MINLP problem, resulting in an NLP, which is then solved using a local NLP solver (e.g., CONOPT). FIRST STEP GAMS – CPLEX12 Max CPU = 3600 s GAP = 0% Solution Pool Parameters solnpoolintensity = 4 solnpoolpop = 2 solnpoolcapacity = 100 populatelim = 100 solnpoolgap = 50 solnpoolreplace = 2

Solve MILP collecting n MILP feasible solutions

SECOND STEP GAMS – CONOPT Max CPU = 3600 s

Store m ≤ n feasible solutions

Solve m NLP By fixing integer variables in original MINLP problem

Report best NLP (if any)

Figure 9 – The two step MIP-nNLP algorithm. By dropping constraints (14) in the first step, composition discrepancy is allowed between the mixture in the transfer operation and in the remaining content of the tank, which is a splitting point. Composition discrepancy is always possible in the problem addressed in the present work because all tanks are subject to holding a minimum volume equivalent to the tank heel at all times. Upper and lower property constraints could be added to both the tank content and the mixture under transfer to mitigate the discrepancy. However, once again, in our problem, that measure cannot be employed because the final crude mix fed to a CDU is not necessarily blended at the tanks. Therefore, it is expected to always have a gap between the MILP and the NLP solutions. Moreover, the chances for running into infeasible NLP solutions after fixing the integer variables found in the first step are not minimal depending on the scenario under consideration. As a remediation, we increase the possibilities of finding a feasible solution by collecting 100 (or less, if solution time tolerance is reached before 100 solutions are found) using the Solution Pool feature available in the solver CPLEX12. After solving all NLPs, only the best NLP solution is reported.

ACS Paragon Plus Environment

Page 33 of 44 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

5.1 Rolling-Horizon Solution Algorithm (RHS) This algorithm was devised by Reddy et al.8 and it is based on the fact that a tank cannot load and unload simultaneously, bearing in mind that a tank output is the only point where splitting takes place in the problem under study. Although a parcel arriving at the refinery might be composed of multiple components, their composition do not change along the unloading operation even if a parcel is split, equation (3). So, if one takes the first time slot as an instance, (14) becomes linear for each tank that has an active output operation because 𝐿𝑡𝑖1𝑟 = 𝐿𝑡0𝑟 and 𝐿𝑖1𝑟𝑐 = 𝐿0𝑟𝑐. Taking that into consideration, the algorithm starts by solving the full scale MINLP problem with (14) applied only to the first time slot, so that composition discrepancy is eliminated for tanks feeding CDUs in that time slot. For all other time slots composition discrepancy is allowed by not applying (14). After solving the resulting MILP problem, all decisions with respect to the first time slot are fixed. Tanks that might have receive a parcel in the first time slot must have their composition updated with the amount of each crude component received. In the next iteration, (14) is applied only for the second time slot, in which case 𝐿𝑡𝑖2𝑟 = 𝐿𝑡𝑖1𝑟 and 𝐿𝑖2𝑟𝑐 = 𝐿𝑖1𝑟𝑐, and making (14) linear again since 𝐿𝑡𝑖1𝑟 and 𝐿𝑖1𝑟𝑐 have been fixed in the previous iteration. The algorithm proceeds until the last time slot is reached. As summarized in Figure 10, the full scale problem is solved from the first to the last time slot using only MILP problems. The method presented in the picture is slightly different from the original one in that at each iteration only the variables corresponding the current time slot are fixed. Just like the two-step algorithm, the RHS does not guarantee optimality. It should be noted that because variables are fixed, smaller MILP problems are solved in each iteration since the scheduling horizon is progressively shortened. Moreover, we set the relative gap to 0% for the first iteration and 1% for all other iterations.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 34 of 44

Initialize the algorithm with t=1

Apply (14) for i = t tanking

𝐿𝑡𝑖,𝑟 = 𝐿𝑡𝑖 ― 𝑟 𝐿𝑖,𝑟,𝑐 = 𝐿𝑖 ― 1,𝑟,𝑐 GAMS – CPLEX12 Max CPU = 3600 s GAP = 0% 1st iteration 1% for the others

Solve MILP Fix all variables for current i

No

feasible?

Report: infeasible

Yes No

i = |I|?

Yes

Report last solution

Figure 10- Rolling-Horizon Solution Algorithm (RHS).

6. Note on the minimum number of time slots It is a common practice to adopt the Additive Approach for determining the minimum number of time slots for solving production scheduling problems. When the two-step algorithm is employed to solve a nonconvex problem with a local NLP solver though, there is no guarantee that the NLP solution will be improved as the number of time slots is increased. Therefore, in this work all problems are solved for the range of time slots varying from 3 to 7 for a better chance of finding the best solution using local solvers.

7. Computational Results All problem instances were implemented in the GAMS system distribution 24.3 and run on an Intel® Core™ i7 2860QM 2.5 GHz 8GB RAM. CPLEX was the MILP solver and CONOPT the NLP solver. The maximum CPU time was set to 3600 s and relative gap set to 0, whichever happens first. Statistical results are presented in Tables S22-S33 and discussion is presented next.

ACS Paragon Plus Environment

Page 35 of 44 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

In tables S22-S25, results are presented for the problems solved with the two-step algorithm. In these tables, the Nº Changeover is an estimation of the total number of changeovers at all CDUs calculated with the optimal solution, Nº sub-parcels is the total number of sub-parcels created by the model for unloading the initial scheduled parcels and is equal to the sum of all active unloading integer variables. These numbers are intended to capture operational aspects of the solution to some extent. The larger these numbers the more handling are demanded for executing the solution in real life. It must be borne in mind though that none of these parameters are influenced by the current objective function (87). Attempts were made to include such decisions in the objective function as penalties. However, estimating the real cost penalty associated exclusively with those parameters is not a trivial task. Therefore, we decided to monitor instead of affecting them. The number of changeover is though impacted by the minimum run-length of 24 hours once a tank is assigned to feed a CDU. A minimal number of changeovers is highly desirable given that it directly impacts the amount of out-of-spec products produced during tank switching at CDUs. Neither of these parameters can be evaluated absolutely. They are directly dependent on the initial inventory, the number of parcels scheduled to arrive at the refinery and on the quality of the initial inventory in tanks and parcels. Taken table S22, in which results for the RTW model considering the FC configuration is presented with the two-step algorithm, it can be seen that not all instances returned an NLP feasible solution even considering the pool of 100 solutions. In fact, for scenario 6 solved with 3 time slots not even a relaxed MILP solution was found. Scenario 1 is by far the most difficult problem to solve. Only two instances were solved. Scenario 1 with 3 slots returned a single NLP solution, whereas the same scenario with 6 time slots returned 5 NLP solutions. Note that, because of (22), the solution found with 3 time slots is not contained in the feasible space of the same problem solved with 6 time slots. Scenario 2 is the second most difficult problem to be solved as the number of NLP solutions is notably reduced in comparison to

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 36 of 44

the other scenarios. For all other scenarios, the number of NLP solutions found depend on the number of time slots used, but in most cases all collected MILP solutions resulted in feasible NLP solutions. The MILP problem demonstrated to be way more time consuming than the corresponding NLP problems. Despite the number of nonlinear constraints, NLP problems were solved in just few CPU seconds. The largest gap between the MILP found in the first step of the two-step algorithm and the MILP solution corresponding to the best NLP solution was 0.63%, which suggests that instead of using solnpoolgap = 50, we could have used a much smaller value. In other words, the MILP solution that produced the best NLP solution was frequently very close to the MILP solution without using the Solution Pool feature of CPLEX. In addition, a smaller number of MILP solutions could have being collected as well, which would impact computational time favorably. The high computational expense solving the MILP problems is due the null relative gap. The computational time could be improved by allowing to exist a nonzero relative gap. Indeed, all instances were also solved with a 1% relative gap (see tables S26S29). Solution time was reduced significantly whereas the objective function was reduced for some instances. However, there were also instances for which the solution found was better in terms of the final objective function, which can only be explained by the nonconvex nature of the problem and the use of a local NLP solver. For the results in table S22, the largest gap between the MILP and best NLP solutions was 4.6%, whereas the largest difference between the best and worst NLP solutions was 9.5%. Comparing results in table S22 with the results in table S23, the latter also corresponding to solutions produced with the RTW model but considering the LC operation solved with the two-step algorithm, it can be readily seen that feasibility is tremendously affected by limiting connectivity between tanks and CDUs (segregating crudes). In addition, the objective function is reduced by as much as 9%. As a general conclusion it could be stated that the FC configuration outperforms the LC configuration.

ACS Paragon Plus Environment

Page 37 of 44 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

Figures S34 and S35 show the Gantt chart for scenario 2 using 5 time slots for the RTW model and the FC and LC configurations, respectively, both solved with two-step algorithm. In these figures, only tanks that had active operations in or out are shown. Red bars represent unloading operations whereas green bars represent loading operations. The top section corresponds to parcel unloading, the intermediate section corresponds to tank operations and the bottom section corresponds to CDU operations. Every parcel unloading operation has a matching tank loading operation. Likewise, every tank unloading operation has a matching CDU loading operation. For sake of simplicity, mixer operations are not shown. On top of bars, depending on whether it is a loading or unloading operation, the destination or origin is presented together with the time slot in which the operation takes place. Numbers underneath bars refer to the volume/flowrate performed in that operation. Scenario 2 has tanks E, G and K inactive and tank L as an injection tank. In figure S34 one can see that tank L is the only tank feeding distillation columns D1 and D3 in time slot 3. That solution is not feasible for the LC configuration given that an injection tank cannot make up more than 30% of a column feed. Alternatively, the solution found for the LC configuration is shown in Figure S35, where one can see that tank L is not used at all. The configuration in which the refinery is operated causes solutions to be different not only in terms of the objective function but mainly in terms of allocation. Figure S36 shows the inventory profile of some of the tanks for solutions presented in Figures S34 and S35. The grey area represents the working volume of the tank, being the bottom the tank heel and the top the tank capacity. It should be noted that in the LC configuration, tank L receives no parcels since there is no injection parcel scheduled for scenario 2 and that tank is never used along the scheduling horizon. On the other hand, in the FC configuration, tank L receives a sub-parcel in time slot 2 filling the tank to a level over its half capacity and, after its settling period, the whole working volume is used to feed distillation columns D1 and D3. The fact that the available working volume of a tank is almost always totally consumed once a

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 38 of 44

tank is allocated to feed a CDU is an interesting behavior that complies with common operational practice for, in that case, less handling at the tank farm is required. Yet, another aspect to be noticed is that when tanks are switched, there can be a significant difference in quality of the final crude blend fed to a CDU, as shown in Figure S37. Smoothening transitions when tanks are switched is desirable for that would mitigate CDU upset. This is the subject of another work that we have under development whose results will be soon submitted for publishing. Figure S38 shows the solution for scenario 2 using 5 time slots when the EFTL is used for the LC configuration and solved with the two-step algorithm. In this figure, tank operations are shown by subslots for those tanks that have multiple outputs. Likewise, CDU loading is also shown by sub-slots for those cases in which they can be loaded by tanks that have multiple outputs so that it can be shown that CDU capacities are not violated in each sub-slot either. Tanks A, B and H actively make use of the flexible enclosing window. Tank B, for instance, starts feeding distillation column D1 3.6 hours before it starts feeding unit D3 and it ceases feeding D1 0.5 hours before it is switched in D3. The flexibility in terms of the time grid of individual operations though, did not necessarily always produced an improved objective function, as can be seen by comparing the results in Tables S23 and S24. Also, the computational effort is in general increased, as a great number of additional constraints and variables are introduced. One of the main motivations for developing the EFTL was to enable the possibility of increasing feasibility by given more flexibility to the problem. However, only one additional instance was solved with the EFTL that could not be solved with the RTW model for the LC configuration. Despite the fact that a lower number of feasible NLP solutions were obtained for some instances and the higher burden in computational requirement, we still believe that the proposed EFTL approach is more compliant with actual common practice operations given that the alignment of multiple tanks starting or ending at exactly

ACS Paragon Plus Environment

Page 39 of 44 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

the same time (RTW approach) would be very restrictive for implementing in tank farms deficient in automation. An increase in feasibility was only achieved by means of dropping the symmetry breaking constraints (20). As expected, the solution time increased as a consequence of avoiding exploring symmetric solutions in the search tree, as can be seen by comparing solution times in Tables S24 and S25. On the other hand, feasible solutions for scenario 6 were obtained. Figure S39 shows the solution found for scenario 6 with 5 time slots using the EFTL model for the LC configuration without (20) solved with the two-step algorithm. The encircled operations illustrate an instance for which constraint (20) would be violated, since both operations are non-overlapping and are allocated to time slots that are not sequenced. The reasoning for dropping constraint (20) was supported by an additional task we undertook in which we have implemented the RTW concept for the Unit Slot model15 (formulation provided in the Supporting Information S40 together with results in tables S41). This model does not make use of symmetry breaking constraints. The only imposition is that time variables must be monotonically increasing, even for nonallocated operations. In that case, it can be seen that solutions for scenario 6 were obtained. This led us to believe that constraint (20) seamed to exclude feasible solutions for that scenario. One of the reasons why the symmetry breaking constraint might have not worked for our problem might be attributed to the fact that we have made the unloading operations overlapping. The theory supporting that constraint relies on the fact that one has non-overlapping operations throughout the system13. However, no effort was made to prove that fact. As for scenario 1, solutions were only obtained for the RTW model solved for the LC configuration when the upper bound on TA is relaxed. The quality constraints can be transformed in soft constraints adding slack variables, which is penalized in the objective function.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 40 of 44

All problems were also solved with the Rolling-Horizon Solution Algorithm and the results are shown in tables S30-S33. When results in those tables are compared with results in Tables S22-S25, a general conclusion can be withdraw; the Rolling-Horizon Solution Algorithm always produce less favored objective solutions as it would be expected, since it offers less flexibility after all variables have been fixed for the initial time slots, i. e., this algorithm presents a myopic behavior. Therefore, after solving for the first time slots the solution for the last time slots becomes more restrictive at each iteration and so compromising the overall solution quality. Solutions obtained for the different proposed approaches follow similar pattern already discussed when the two-step algorithm is used, as the points already discussed are intrinsic of modeling approach and not on the employed solution strategy.

8. Conclusions In this work, the crude oil scheduling problem was addressed with focus on operational features of a real-world refinery from Brazil. Two models were derived from the MOS formulation as proposals to address the fact that tanks have multiple outputs subject to common lower and upper flowrates, namely: the RTW and the EFTS. Moreover, two configurations as to the connectivity between parcels and tanks and between tanks and CDUs were studied (FC and LC), being the latter the one in greater accordance with operational practices in the studied refinery. The resulting models were tested against a number of distinct scenarios presenting different levels of difficulty to be solved. In all cases, the LC configuration demonstrated to be more challenging and producing inferior economical solutions. Even though, the RTW formulation seemed to be robust enough to solve a great number of problem instances, especially when the symmetry breaking constraints were dropped and the two-step algorithm was used. As for the EFTS model, although it enables finding solutions that seem to be more achievable from an operational perspective, it still requires improvements.

ACS Paragon Plus Environment

Page 41 of 44 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

Supporting Information The supporting information for this article is available and includes the common data for all scenarios such as tank capacities and output flowrates (S1), crude oil adopted marginal values (S2), and distillation unit flowrate limits (S3). Qualities for each crude oil grade and qualities of the diesel cut produced by each crude oil grade cannot be presented for confidentiality reasons. In S4-S21 are given input data specific for each scenario. Tables S22-S33 present computational results for the two presented models (RTW and EFTW) solved with the two solution algorithms (Two-Step and RHS). Figures S34-S39 present Gantt charts, inventory profiles and quality of resulting mixture fed into CDUs. In S40 the RTW version applied under the ideas of the Unit Slot model is presented. Finally, S41 shows solutions found with the model presented in S40. This information is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement The authors gratefully acknowledge the financial support from Petrobras under the grant 0055.0089001.14.9.

Nomenclature Sets i is v r c k 𝐼𝑟 𝑂𝑟 𝑅𝑝 𝑅𝑖𝑛𝑗 𝑝 𝑅𝑡 𝑅𝑖𝑡 𝑅𝑖𝑛𝑗 𝑡 𝑅𝑂𝑛𝑒 𝑡 𝑅𝑚

set of time slots (priority-slots) set of time sub-slots 𝐼𝑆 = {𝑖𝑠1, 𝑖𝑠2,𝑖𝑠3} set of transfer operations set of resources set of crude oil grades set of properties (API gravity, sulfur content, TAN, yield) set of inlet operations at resource r set of outlet operations at resource r set of resources r that represent crude oil parcels set of injection parcels set of resources r that represent tanks set of tanks that are inactive set of injection tanks set of tanks that connect to a single CDU set of resources r that represent virtual mixers

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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

𝑅𝑑 𝑊𝑢 𝑊𝑖𝑣 𝑊𝑖𝑛𝑗 𝑣 𝑊𝑡 𝑊𝑑 𝐾𝑚 𝐾𝑣

set of resources r that represent CDUs set of unloading operations set of inactive operations set of operations related to either an injection parcel or an injection tank set of inflow and outflow operations related to tanks set of distillation operations set of properties calculated on a mass basis set of properties calculated on a volume basis

Parameters 𝑆𝑣, 𝑆𝑣 𝐸𝑣, 𝐸𝑣 𝐷𝑣 𝑉𝑡𝑣, 𝑉𝑡𝑣 𝐹𝑅𝑣, 𝐹𝑅𝑣 𝐹𝑅𝑟, 𝐹𝑅𝑟 𝐿𝑡𝑟, 𝐿𝑡𝑟 𝐿𝑡0𝑟 𝐿0𝑟𝑐 BS 𝑥𝑘, 𝑥𝑘 𝐷

𝑥𝑘 𝑦𝑟𝑐

lower/upper bound on start time of operation v lower/upper bound on end time of operation v lower bound on the duration of operation v lower/upper bound on total volume transferred with operation v lower/upper flowrate bound for operation v lower/upper flowrate bound for resource r lower/upper limits for the volume held in resource r initial total inventory in tank r initial inventory for crude c in tank r brine settling time lower/upper bound on property k of crude mixes upper bound on property k of the diesel cut produced at CDUs composition of crude c contained in resource r

Continuous variables 𝑆𝑖𝑣 start time of operation v allocated to priority-slot i 𝐷𝑖𝑣 duration of operation v allocated to priority-slot i 𝐸𝑖𝑣 end time of operation v allocated to priority-slot i start time of enclosing time window of tank r in priority-slot i 𝑆𝑡𝑖𝑟 𝑡 duration of enclosing time window of tank r in priority-slot i 𝐷𝑖𝑟 𝑡 end time of enclosing time window of tank r in priority-slot i 𝐸𝑖𝑟 𝑡 start time of second sub-slot of enclosing time window of tank r in priority-slot i 𝑆𝑠𝑖𝑟 𝑡 end time of second sub-slot of enclosing time window of tank r in priority-slot i 𝐸𝑠𝑖𝑟 𝐷𝑠𝑖,𝑖𝑠,𝑣 duration of operation v in sub-slot is and priority-slot i 𝐷𝑚𝑖,𝑖𝑠,𝑣 duration of operation v in sub-slot is and priority-slot i, where v is an output of a mixer ∆𝑇𝑖,𝑖𝑠,𝑟 duration of sub-slot is in priority-slot i for tank r 𝑡 total volume of crude mix transferred through operation v in priority-slot i 𝑉𝑖𝑣 𝑉𝑖𝑣𝑐 volume of crude c transferred through operation v in priority-slot i 𝑡 total volume of crude mix held in resource r at the beginning of priority-slot i 𝐿𝑖𝑟 𝑡 total volume of crude mix transferred through operation v in time sub-slot is and 𝑉𝑠𝑖,𝑖𝑠,𝑣 priority-slot i 𝑉𝑠𝑖,𝑖𝑠,𝑣,𝑐 volume of crude c transferred through operation v in time sub-slot is and in priorityslot i

ACS Paragon Plus Environment

Page 42 of 44

Page 43 of 44 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

volume of crude c held in resource r in priority-slot i (0,1) continuous variable that assumes 1 if a tank has receive crude mix in priorityslot i 1 if operation v is assigned to start in sub-slot is and priority-slot i 1 if operation v is assigned to finish in sub-slot is and priority-slot i (0,1) continuous variable that assumes 1 if resource r has at least one active output in priority-slot i

Binary variable 𝑍𝑖𝑣 1 is operation v is assigned to priority-slot i 𝑍𝑠𝑖,𝑖𝑠,𝑣 1 if operation v is assigned to sub-slot is and priority-slot i

References (1) OPEC Montly Oil Market Report, 11 May 2017 available at http://www.opec.org/opec_web/en/publications/4054.htm (2) Grossmann, I.E., Advances in mathematical programming models for enterprise-wide optimization, Computers and Chemical Engineering 2012, 47, 2– 18 (3) Lee, H.; Pinto, J. M.; Grossmann, I. E.; Park, S., Mixed-Integer Linear Programming Model for Refinery Short-Term Scheduling of Crude Oil Unloading with Inventory Management, Ind. Eng. Chem. Res. 1996, 35, 1630-1641 (4) Wenkai, L.; Hui, C.W.; Hua, B.; Tong, Z., Scheduling crude oil unloading, storage, and processing, Ind. Eng. Chem. Res. 2002, 41 (26), 6723-6734. (5) Jia, Z; Ierapetritou, M.; Kelly, J. D. Refinery Short-Term Scheduling Using Continuous Time Formulation: Crude-Oil Operations 2003, Ind. Eng. Chem. Res., 42, 3085-3097 (6) Moro, L. F. L.; Pinto, J. M. Mixed-Integer Programming Approach for Short-Term Crude Oil Scheduling, Ind. Eng. Chem. Res. 2004, 43, 85-94 (7) Reddy, P. C. P.; Karimi, I. A.; Srinivasan, R. Novel Solution Approach for Optimizing Crude Oil Operations, AICHE Journal 2004, 50 (6), 1177-1197. (8) Reddy, P. C. P.; Karimi, I. A.; Srinivasan, R. A new continuous-time formulation for scheduling crude oil operations, Chemical Engeneering Science 2004, 59, 1325 – 1341. (9) Li, J.; Karimi, I. A.; Srinivasan, R. Improving the robustness and efficiency of crude scheduling algorithms, AIChE Journal, 2007, 53 (10), 2659-2680.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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

(10) Furman, K. C.; Jia, Z.; Ierapetritou, M. G. A Robust Event-Based Continuous Time Formulation for Tank Transfer Scheduling, Ind. Eng. Chem. Res. 2007, 46, 9126-9136. (11) Karuppiah, R.; Furman, K. C.; Grossmann, I. E. Global optimization for scheduling refinery crude oil operations, Computers and Chemical Engineering 2008, 32, 2745–2766. (12) Saharidis, G. K. D.; Minoux, Dallery, Y. Scheduling of loading and unloading of crude oil in a refinery using event-based discrete time formulation, Computers and Chemical Engineering 2009, 33, 1413–1426. (13) Mouret, S.; Grossmann, I. E.; Pestiaux, P. A Novel Priority-Slot Based Continuous-Time Formulation for Crude-Oil Scheduling Problems, Ind. Eng. Chem. Res. 2009, 48, 8515–8528. (14) Mouret, S.; Grossmann, I. E.; Pestiaux, P. Time representations and mathematical models for process scheduling problems, Computers and Chemical Engineering 2011, 35, 1038–1063. (15) Chen, X.; Grossmann, I. E.; Zheng, L. A comparative study of continuous-time models for scheduling of crude oil operations in inland refineries, Computers and Chemical Engineering 2012 44, 141– 167. (16) Li, J.; Misener, R.; Floudas, C. A. Continuous-time modeling and global optimization approach for scheduling of crude oil operations. AICHE Journal 2012, 58 (1), 205-226.

For Table of Contents Only (Tiff – 300 DPI) Crude oil parcels

CD1 T A TB TC CD2 TD TE Inline TF Mixer TG H TI T CD3 TK Refinery L

Multiple tank outflow

v

CD1

1

TA

v

2

CD3

𝑣1

𝑣2

Unsynchronized operations within the same time slot

For Table of Contents Only (Word Art)

ACS Paragon Plus Environment

Page 44 of 44