Integration of Crude-Oil Scheduling and Refinery Planning by Lagrangean Decomposition

In this work, a Mixed-Integer Nonlinear Programming (MINLP) modeling framework for integrating short-term Crude-oil Scheduling (CS) and mid-term Refinery Planning (RP) has been developed and effectively solved by a proposed Lagrangean Decomposition (LD) algorithm. The principles of this integration are based on the fact that both Crude-oil Scheduling and Refinery Planning have their economic net values as their objectives, and that they are physically linked by the Crude Distillation Unit (CDU). A multi-scale approach is proposed in the framework to aggregate continuousand discrete-time formulations in CS and RP, respectively. Compared to hierarchically solving the non-integrated CS and RP, computational results show significant improvement regarding the economic objective values. Moreover, the proposed LD approach requires less CPU time converging to a small (1%-5%) optimality gap when compared to the monolithic approach using state-of-the-art MINLP solvers.


Introduction
In petrochemical or crude-oil refinery industries, two major decision-making problems have been the mid-term Refinery Planning (RP) and the short-term Crude Oil Scheduling (CS). The former problem is typically defined as a midterm downstream unit operating and flowsheet logistic optimization problem for refinery downstream processes, including the crude oil distillation, intermediate separation, reaction processes, and petrochemical product blending. In this problem, it is often assumed that there are several types of crude oil available on stock with certain amounts of volume, so that the optimization can be solved to satisfy the demands of various petrochemicals to maximize the total net sales profit. The latter problem can be defined as a short-term upstream crudeoil inventory management optimization, where the crude transfer volume and tank inventories are optimized in a complex system of vessels, tanks, and crude distillation units (CDUs), in order to minimize the total cost of crude purchasing and/or tank inventory as well as to satisfy the crude demand at CDUs and nonoverlapping processes.
In practice, the decision-making process is commonly performed as follows: 1) solve the RP problem maximizing the total product sales; 2) solve the CS problem satisfying the demand feedstocks of RP and minimizing total cost; 3) determine the total Net Present Value (NPV) by the sales objective in RP subtracted from cost objective in CS. However, from an integrated perspective of the whole refinery, the solutions obtained through this procedure can be suboptimal. To clarify this suboptimality, consider the following scenario of the decision-making processes in a refinery plant. In the RP problem, there is only one distillation unit that has the same capacity generating naphtha and kerosene where the market price of naphtha is higher on the horizon. There are two crude-oil candidates available on stock where Crude A yields more naphtha and Crude B otherwise. Since the purchasing price for Crude A and B are the same, the planner solves the RP problem and decides to produce more naphtha than kerosene by requiring only Crude A without B at feedstocks. Then, in the CS problem, there is only a small amount of Crude A available for charging the distillation unit, where most Crude A is saved in the storage level, but there is plenty of Crude B available. However, the crude transferring rules and the limited number of crude tanks do not allow the scheduler to transfer the required amount of Crude A available for charging in such a short time. Because the feedstocks are fixed by the decision of the planner, the scheduler has to come up with a bad decision of purchasing some immediate costly Crude A. In a refinery-level perspective, the optimal decision is apparently to allow some amount of Crude B until Crude A are available after the crude transferring processes. Indeed, this suboptimality arises because the planner and scheduler make decisions in each of their subsystems without a global view.
Given the recent developments in large-scale optimization algorithms and computing capabilities, there is a large gap between the development of an integrated modeling framework and an integrated solution strategy, which means optimizing simultaneously both the RP and the CS problems.
In this work, these challenging issues are addressed with the following two questions: i) How to perform the integration to connect the crude feedstocks information in terms of volume and time? ii) Once the integration is obtained, what is the most efficient approach to solving the integrated problem to optimal- Figure 1: Simplified flowsheet of a crude oil refinery. Adapted from [18].
ity? To address the first challenge, we consider a real-world refinery problem, stated in Section 3, where the CS is modeled as a continuous-time MINLP in Section 3.2, and RP is modeled as a discrete-time NLP in Section 3.1. In Section 4 we propose an integration including the RP and the CS formulations, as well as their integrated formulation, which contains a set of linking constraints for connecting operating volumes, and a multi-scale modeling approach for connecting operating time. This large-scale full-space formulation is the second challenge we want to address. We propose a Lagrangean Decomposition algorithm to iteratively update the primal/dual objective bounds based on the decomposition of the full-space problem down to scheduling and planning subproblems. Computational experiments are performed in Section 5 on three refinery case studies in terms of different time horizons. In Section 6 the advantage of using the proposed Lagrangean Decomposition is shown from the results, in terms of CPU time and optimality gap. This advantage is significant especially in cases with longer time horizons, as compared to directly using state-of-the-art global MINLP solvers (monolithic approach).

Literature Review
Crude oil refineries are complex systems with multiple operations that rely on the properties of feedstocks and desired products. Typically, a crude oil refinery consists of three segments that make up the production system: a) crude oil management; b) crude-to-fuel transformation (process operations); and c) blending operations, as shown in Figure 1.
The first segment is related to crude oil management and begins with the arrival of crude oil, usually via ships or pipelines. The oil is stored in storage tanks to be further mixed, sent to the feed tanks, and used as feed to the crude oil distillation unit [18]. In the second segment, several unit operations based on physical and chemical transformations such as atmospheric and vacuum distillation, catalytic reforming, fluid catalytic cracking, delayed coker, hydrocracking, hydrotreating, debutanizer, depropanizer, superfractionator, among others, are used to crack the crude oil into smaller and lighter fractions, as well as to im-prove the quality of the final products [6]. In the third segment, known as blend-shops, the final products are blended to meet market specifications regarding their properties, such as specific gravity, sulfur content, octane number for gasoline, cetane number for diesel, among others.

Crude Oil Scheduling
The crude oil scheduling problem, often referred to as part of the upstream operations (crude oil management in Figure 1), includes crude oil unloading from ships or vessels to tanks, transfer from storage to charging tanks to prepare the CDU feed, and feeding operations from the charging tanks to the CDU. The decisions over the scheduling operations are related to continuous CDU charging through the time horizon, with given vessel arrival times, tanks and CDU capacities, crude properties, and non-overlapping constraints of certain operations. The scheduling problem is solved to maximize the gross margin of the crude charging to the CDU and to minimize the cost of inventories. Lee et al. [21] consider an objective of minimizing the cost of tank inventory and vessel demurrage, while Mouret et al. [28] apply another objective of maximizing the gross margin of the amount of crude charged and distilled in the CDU, with the gross margin value estimated by subtracting the purchase value from the potential final product sales in the refinery planning level.
Lee et al. [21] model the crude oil scheduling problem as a discrete-time Mixed-Integer Linear Programming (MILP) formulation using time intervals and the bilinear mixing constraints are avoided by introducing component flows. This model is tested in four cases of different problem sizes, showing advantages in terms of solvability and optimality. Alternatively, a continuous-time MILP formulation has also been developed using continuous event points [16]. Further, continuous-time formulation MINLP models are also developed in [17] and [28] based on event points and priority slots, respectively, and solved with linearization based algorithms. The latter formulation shows advantages in terms of the computational solution time as it requires fewer binary variables. Mouret et al. [28] formulate an MINLP based on priority slots, and use a two-stage solution in which a relaxed MILP is first solved without the bilinear blending constraints followed by the NLP with fixed binary variables from the MILP. Further, Castro and Grossmann [4] developed a single time grid-based continuous-time formulation and made a comparison between event point and priority slots.

Refinery Planning
The RP problem concerns the downstream operations, including the processshops and blend-shops (segments 2 and 3 in Figure 1), and there is an overlap at the CDU level for these problems, including the inventory and time as well. Although many commercial tools such as GRTMPS (Haverly Systems), PETRO-SIM (KBC), PIMS (Aspen Tech), and RPMS (Honeywell Hi-Spec Solutions) [25] have been developed to effectively solve this problem, most of them are simulation-based, do not perform rigorous optimization using equationbased models and/or cannot address multi-period optimization problems [3]. Therefore, optimization-based modeling tools represent a potential benefit when properly used to solve the refinery planning problem.
As it is shown in Figure 1, the process-shop operations begin with the distillation unit, which is one of the most important parts in the refinery planning modeling. This task is usually performed using a fixed-yield linear programming (LP) model, in which the yields and properties of distillates are obtained from the true boiling point (TBP) curve from the crude oil assay profile. The distillates are sent to the other processing units for further processing. The catalytic reforming aims at cracking low octane heavy naphtha into lighter fractions with a higher octane number. In the fluid catalytic cracking, the vacuum gas oil is processed under heat, pressure, and catalyst to produce lighter fractions (mostly gasoline and LPG). The delayed coker thermally processes vacuum residue generating coke and some light fractions. The debutanizer is used to stabilize light naphtha and separate it from light gases such as fuel gas and liquefied petroleum gas. Similarly, the depropanizer separates propane and propylene from butane. The superfractionator (propane-propylene splitter) separates the propane-propylene mix into individual components. Hydrocracking converts heavy gas oil molecules into lighter ones under the presence of hydrogen and catalyst. Hydrotreating is performed to purify the products either before the final blending and specification steps or before units that have requirements regarding their feed properties. The final streams are then blended into the final products to meet specifications and stored in pools for further distribution.

Integration Planning-Scheduling
Both crude oil scheduling and refinery planning problems aim at an economic objective and are linked by the CDU so that there is an incentive to solve an integrated problem, which guarantees global mass flow balances and operating times. Moreover, optimizing the integrated problem will result in improved solutions compared to optimizing each one separately. One of the main challenges in this approach is the time scale difference between these two problems. Ierapetritou and Floudas [16] proposed a rolling horizon integration for planning and scheduling, and successfully solved the problem by a hierarchical decomposition approach. Besides, it is also challenging to solve the large-scale MINLP of the integrated planning and scheduling problem. Lagrangean decomposition has shown its potential advantages in [23]. Further, Mouret et al. [30] integrated a crude oil scheduling model with a continuous-time formulation and a conceptual refinery planning model using an Artificial Neural Network (ANN) for the CDU modeling [12]. In [30], the advantages of considering an integrated approach of the crude oil scheduling and refinery planning problem were shown by achieving a computationally-efficient framework for these problems using Lagrangean decomposition, which will be further described in Section 6. Mouret et al. [29] proposed a new Lagrangean Decomposition approach with a hybrid formulation in the dual problem of updating the Lagrangean multipliers. Oliveira et al. [31] further improved the subgradient formulation in the dual problem by a trust-region approach.

Problem Statement
The problem addressed in this work can be defined in three parts: 1) Crudeoil Scheduling, 2) Refinery Planning, and 3) an integrated framework to link CS and RP. These three parts are discussed in Section 3.1, Section 3.2, and Section 3.3, respectively. Furthermore, the integrated problem is defined and the challenge of modeling the difference of time representation is discussed in Section 3.3.

Refinery Planning problem
The RP problem considered in this work involves the process-shops and blendshops in the crude oil refining process and uses as input the CDU feed information provided by the crude oil scheduling solution. The main decisions are volumetric flow rates between the operating units including CDU, but in a time aggregate level, e.g., volume per week, per two weeks, or longer. The operations include the following information: (i) process units such as atmospheric, flash and vacuum distillation, debutanizer, depropanizer, propane-propylene superfractionator, fluid catalytic cracking, catalytic reforming, hydrotreating, delayed coker, and blenders; (ii) blending pools to store intermediate and/or final products; (iii) inlet and outlet flows between process units and blending pools.
Besides, the following parameters are given for the refinery planning problem: (a) length of the time horizon and number of discrete time steps; (b) feedstocks availability and properties; (c) CDU operational mode(s); (d) capacity of each processing unit as well as holdups for tanks; (e) lower and upper bounds for the flowrates; (f) demands for final products as well as their property specifications. The objective function in the refinery planning is to maximize NPV, with continuous decisions regarding all the flows and properties tracked throughout the process as well as binary decisions regarding operational modes.

Crude Oil Scheduling problem
The crude oil scheduling is presented in Figure 2 which will be used later in the case studies. Crude vessels are represented by r 1 to r 3 , storage tanks as r 4 to r 6 , and charging tanks as r 7 to r 9 . Crude unloading operations are represented by v 1 to v 3 , transferring/blending operations as v 4 to v 10 , and charging operations to the crude distillation unit (CDU) as v 10 to v 14 . The given information and parameters are:   [21] crude to be transferred. Besides, the schedule of operations cannot violate any of the following logic constraints: (i) only one unloading operation is allowed for each vessel; (ii) simultaneous inlet and outlet operations are not allowed for each tank; (iii) each tank can charge no more than one CDU at a time; (iv) each CDU can be charged by no more than one tank at a time; (v) each CDU needs to operate continuously through the time horizon. Regarding the objective function for the crude oil scheduling, it is typically any rational economic value. The objective function considered in this work minimizes the total cost of crude charged to the distillation unit.

Integrated Refinery problem
In this work, the crude oil scheduling and the refinery planning problems are considered to be integrated into a single model that we denote as the full-space problem. Besides the information from both problems, the full-space problem also includes the link between them, which is the mass conservation in the blender unit, which appears both at the end of Crude Oil Scheduling and at the beginning of Refinery Planning. In this section, we discuss how the integration of planning-scheduling is performed regarding different time scales, as well as solution approaches to tackle the full-space problem.
The full-space problem is defined as the combination of 1) multi-period Refinery Planning continuous nonlinear model given in Section 4.1; 2) extended drude-oil Scheduling MINLP model with multi-scale approach in Appendix A; and 3) linking constraints given by Equations 20 and 21, yielding an MINLP optimization model.

Solution Approaches
There are three main solution approaches to tackle the Planning-Scheduling problem. Conventionally, the refinery planning and the crude oil scheduling problems are optimized sequentially in a strategy denoted as the hierarchical approach. In contrast, to exploit the potential benefit of an integrated solution, a single planning and scheduling model is optimized in a monolithic fashion. However, the monolithic problem can be challenging to solve due to its large size, and to integrate planning and scheduling operations at the same level can be practically unnecessary as well as introduce issues of decoupling the important decisions. Thus, in this work, we approach the problem through a decomposition technique, in which the Lagrangean Decomposition method is used to decompose an optimization problem with complicating constraints into more easily solvable subproblems, which are solved separately in an iterative fashion.

Hierarchical Approach
In the hierarchical approach, the Nonlinear Programming (NLP) optimization problem for Refinery Planning is solved first, as the planning objective of the product sales often dominates in the integrated objective. Then, the MINLP problem for Crude Oil Scheduling model is optimized, with the CDU feed flows fixed by using the Refinery Planning solution. This non-integrated approach may return a suboptimal objective for the integrated problem, which is obtained by adding up the planning and scheduling objectives. This suboptimality arises from the fact that when considering the optimization problems independently, the interaction between them through the constraints that involve variables in both subproblems is disregarded. In some cases, when planning solutions are fixed, the scheduling problem is either infeasible or there is a considerable penalty term in its objective. The infeasibility may occur by 1) a mismatch in CDU flows between RP level and CS level due to the time aggregation at which capacities are considered, or by 2) the constraint requiring the flows with respect to meeting the feedstock blends and properties. To avoid these infeasibilities we allow the purchase of crude c in time period t feeding directly to CDU r denoted by f imm t,c,r , and the cost/penalty of using this crude is set to be 3 times more than the nominal crude purchasing values. These penalties can be avoided by considering the integrated problem, considering the constraints simultaneously and obtaining coordinated operations of the RP and the CS.

Monolithic Approach
In the monolithic approach, the full-space MINLP model that integrates the RP and CS models is optimized by using state-of-the-art MINLP solvers, including convex optimization solvers such as DICOPT, Bonmin, and SBB, as well as global solvers such as BARON, ANTIGONE, SCIP, LINDOGLOBAL, and Couenne. A complete comparison and discussion about the solvers for MINLP can be found in [20]. However, this large-scale MINLP model is challenging for the solvers due to a large number of binary variables and logic constraints in the Crude Oil Scheduling as well as the large number of non-convex bilinear terms in the Refinery Planning.

Lagrangean Decomposition Approach
The Lagrangean Decomposition approach used in this work is defined based on [29], which can be considered as a primal-dual algorithm that uses an improved subgradient formulation for the dual problem by exploiting a trust-region approach [31]. The problem is first initialized by defining the initial upper and lower bound. As the linking constraints are the equalities of crude feed, the corresponding Lagrangean multipliers λ t,c are indexed by time period t and crude type c, corresponding to the linking constraints (Equations 20 and 21). They are initialized by using crude purchasing cost V crude c . At each iteration of the algorithm, the Lagrangean dual Z LD is first obtained by solving the decomposed subproblems to optimality. Z LD is the tightest dual bound among Lagrangean relaxation and other continuous relaxations [13]. Because of the nonconvexity of the problem considered here, strong duality is not guaranteed and therefore there might be a gap between the global optimal solution and the Lagrangean dual which this algorithm may not be able to reduce to zero.

Modeling Methodology
The mathematical formulation for the full-space problem is developed in three stages: a) Multi-period discrete-time formulation for Refinery Planning, maximizing the total NPV; b) Priority-slot based continuous-time formulation for Crude-oil Scheduling with multi-period extension, minimizing the total crude purchasing cost; c) Set of linking variables and constraints of crude-oil feedstocks between Crude-oil Scheduling and Refinery Planning, with the objective maximizing the total NPV minus total crude purchasing cost;

Sets and Indices
The following sets are introduced for the Refinery Planning model.
of CDUs CDU1, CDU2, J L is the subset of other major separation units FLASH, TOWER, VDU, J U is the subset of other units.
• c ∈ C for crude-oil types.
• u ∈ U for all operating flows in the system.
• u i ∈ IN j ⊂ U for inlet flows of unit j.
• u o ∈ OU T j ⊂ U for outlet flows of unit j.
• l ∈ L for all final products.
• p ∈ P = {SG, SUL} for two crude properties tracked, representing Specific Gravity and Sulfur Content, respectively.
• mix ∈ M IX, sep ∈ SEP , for set of mixers and separators.

Parameters
The following parameters are given.
• H is the planning horizon length in each time period t.
• V prod l is the sales value of product l.
is the cost value of blending component b.
• C j is the capacity of unit j.
• ν unit j,u o is the yield of outlet flows u 0 ∈ OU T j of unit j ⊂ J U .
• ν col j,u o ,c is the crude-dependent yields of outlet u o ∈ OU T j of unit j ⊂ J L , which is pre-computed from the micro-cut fractions of crude assay data.
• ν cdu j,k,c is the crude-dependent yield of outlet crude cut k ∈ K of CDU j ⊂ J C , which is pre-computed from the micro-cut fractions of crude assay data.
• ν cdu j,k s ,c is the crude-dependent yield of outlet swing cut k s ∈ K S of CDU j ⊂ J C , which is pre-computed from the micro-cut fractions of crude assay data.
• Q l,p is the upper bound for property p for final product l.

Variables
The following optimization variables are introduced.
• F cdu t,j,c is the volumetric flow rate of crude c to feed CDU j ⊂ J C in time period t.
• F T cdu t,j is the total volumetric flow rate of crude to feed CDU j ⊂ J C in time period t. • F T cdu t,j,k is the total volumetric flow rate of crude cut k in CDU j ⊂ J C in time period t.
• F cdu t,j,k s ,c , F T cdu t,j,k s are the volumetric flow rate in crude type c and in total, from swing cut k s in CDU j ⊂ J C in time period t. and F T cdu,heavy t,j,ks being the total volume over all crudes.
j,u i being the total volume over all crudes.
• F out t,j,u o ,c is the volumetric flow rate of crude c of the output flow u o of unit j ⊂ OU T j in time period t, with F T out t,j,u o being the total volume over all crudes.
is the volumetric flow rate of blending component b used in time period t.
• F prod t,l,c is the volumetric flow rate of crude c in final product l produced in time period t, with F T prod t,l being the total volumetric flow rate of all crudes.

Constraints
Equation (1) restricts the total input volume of unit j under its capacity. As the capacity Cap j of unit j is given on a daily basis, the total capacity can be calculated by multiplying it by the horizon length H.
Equation (2) represents the volumetric flow balance for the total input and output for unit j.
Equations (3) to (5) state the balance for volumetric flow rate variables F and total volumetric flow rate variables F T .
Equations (6) and (7) are used to calculate the production for each unit. The yield data ν col j,u o ,c is crude-dependent and is calculated from the crude assay information, while ν unit j,u o is crude-independent.
Equation (8) and (9) represent the volumetric balance of mixing and separating operations, specifically, one or more out-coming flows u o from units j mixing to one in-coming flow u i towards its unit j for mixing operations, one out-coming flow u o from its unit j separating to one or more in-coming flows u i towards their j for separating, tracked for each crude type c in time period t. Particularly, for simplicity and consistency in the modeling, the case where one out-coming flow proceeds directly towards another one in-coming flow can be treated as either a mixing or a separation operation.
Equations (10) to (12) impose specifications for the final products in order to meet the specified demands. Furthermore, the trilinear terms are replaced by bilinear terms by introducing new variables that represent the product between the sulfur content and specific gravity terms. That strategy reduces the computational effort during the optimization step.

Planning Objective Function
The objective function for the refinery planning is to maximize Net Present Value (NPV), which is the total product sales values subtracted by the total cost of blending components, as shown in Equation 13.

Crude Oil Distillation Unit Modeling
Both rigorous and surrogate representations can be used to predict outputs of distillation processes. Despite the robustness and accuracy in the rigorous predictions, these models demand high computational effort. On the other hand, non-rigorous models are typically simple and have reasonable accuracy, being commonly used in crude oil refineries [22]. Examples of distillation unit models have been addressed by [35,22,1,24,25,19,9], among others. Three different models have been considered in this work: Fixed-Yield (FY), Multiple Fixed-Yield (MFY), and Swing-Cut (SC). The best performance was achieved using the SC model so that its features and equations are presented as follows. For further information about the FY and MFY features as well as their respective models, please see Appendix B.

Swing-Cut
The swing-cuts are essentially internal modeling constructs, not necessarily physically present in the tower. Three swing-cuts sc are created and their adjacent cuts with respective set of micro-cuts CU T sc are shown in Table 1.
The mathematical model for the conventional swing-cut method is given as follows. Equations (B.1) to (B.4) from the fixed yield model are used (see Table 1: Swing Cut design Swing Cut Lighter Cut Heavier Cut Initial Appendix B), which represent a volumetric balance for the CDU, the flow calculation for each distillate and the volume-and mass-based balances, respectively. However, it is important to note that the final cuts in the fixed yield model are now the intermediate cuts in the swing-cut model, which will further be blended with the swing-cuts to create the final cuts. Equation (14) represents the volumetric flow balances for the swing-cuts, whereas Equations (15) and (16) define the total flow variables which are identical to the aggregation of their composition flow variables over all crude types, for the light and heavy swing-cuts, respectively.
When the intermediate cuts are mixed with their respective swing-cuts, e.g., intermediate naphtha cut is mixed with the light fraction of the first swing-cut, Equation (17) calculates the final cut flow while the nonlinear Equations (18) and (19) calculate the volume-and mass-based properties for the final cut.
The conventional swing-cut method creates additional degrees of freedom by allowing the optimization of each swing-cut amount flowing to the lighter and heavier final distillates. Therefore, the search space for the optimization becomes larger and at least an equivalent solution is expected when compared to the (multiple) fixed yield methods.

Crude Oil Scheduling with multi-period extension
The priority-slot based MINLP formulation of Crude Oil Scheduling problem, proposed in [28], is used in this work. Specifically, we consider the Multi-Operation Sequencing (MOS) time representation in [27], which allows several operations assigned to one slot. To be integrated with the Refinery Planning section of the model, we propose a new extension of the priority-slot based formulation, discussed in Section 4.4.1, denoted multi-scale approach. The basic idea is to add the dimension of planning period t to the original scheduling formulation, meaning a short-term scheduling problem is solved in each planning period t. The detailed mathematical formulation can be found in Appendix A.

Integrated full-space problem
The linking constraints are defined as Equations (20) and (21). The left-handside scheduling variables y F represent the amount of crude c charged into CDU r at time period t, while the right-hand-side planning variables F represent the corresponding input flows which are the same flow modeled in scheduling as y F . Thus, these two equations state the equivalency of those two types of variables.
The objective function for the full-space problem is the economic value combining both objective values in crude oil scheduling and refinery planning, specifically, maximizing the NPV in planning subtracted by crude purchase cost in scheduling, as given by Equation (22).

A multi-scale modeling approach for different time representations
The difference in terms of time scales between planning and scheduling remains a major challenge to be addressed. As mentioned in [34], there are not many modeling frameworks that address the integration of "short-term" scheduling and "mid-term" planning where the information between the scheduling-level and planning-level can be connected effectively through their time scales. More generally in [14] and [5], this challenge is considered as one future research direction. Given different time representation in crude oil scheduling (continuous) and refinery planning (discrete), a multi-scale approach is proposed in this work, illustrated by using an example as in Figure 4 (see Su et al. [33]), in which there are three time planning periods, t 1 , t 2 , and t 3 in the planning, as well as  Figure 4 shows the difference between Planning and Scheduling in terms of time representation, which makes it difficult to ensure simultaneous time consistency in the full-space problem. To accomplish that task, the following steps are performed to generate the bottom structure in Figure 4: 1. Add planning period t to scheduling slot i.

Add tank inventory balancing variables and constraints between periods.
All three steps are applied to the scheduling formulation. In step 1, the dimension of planning period t has been added to the scheduling formulation as in Appendix A. In step 2, as the scheduling is extended to multiple planning periods, the initial tank inventories of the future time periods after the first period t 1 need to be the same as the final inventory from the previous period. The following inventory balance constraints are included in the multi-scale scheduling.

Lagrangean Decomposition
The full-space MINLP is given by Equations (30)-(35), where we denote x and f (x) as the planning variables and constraints, y and g(y) as the scheduling variables and constraints. Equation (33) represents the linking constraints, where some variables in both planningx and schedulingŷ need to be equal.
x ∈ X (34) In order to obtain the objective function predicted by the Lagrangean decomposition (Z LD ), the full-space MINLP must be first relaxed through a Lagrangean relaxation, providing an objective function Z LR , by dualizing the linking constraint. In that case the problem becomes the one given by Equations (31), (32), (34) and (35), in addition to Equation (36).
Note that the for a given λ, Z LR (λ) decomposes into two subproblems: one in the space of the planning variables x (Equations (31), (34), and (37)), and the other in the space of the scheduling variables y (Equations (32), (35), and (38)).
where both subproblems are parametrized with respect to the multipliers. Thus, Z LD is obtained by adding up the current two objective values Z LR1 and Z LR2 from the two subproblems. Since the full-space problem is a maximization problem, Z LD will be an upper bound to the original objective function value Z. The Lagrangean dual, Z LD , can be obtained by minimizing, with respect to the multipliers, the sum at Z LR1 (λ) and Z LR2 (λ), as shown in Equation (39).
The bound provided by the Lagrangean dual is a relaxation of the original objective function meaning that Z LD ≥ Z. However, given that there are discrete variables in the models, they result in non-convexities in which strong duality does not hold. Also, the failure to obtain strong duality can be caused by non-convex constraints with bilinear terms from the blending equations. For maximization problems, the upper (dual) bound at least satisfies Z LD > Z, and this difference is defined as the duality gap. We aim at making the dual gap as small as possible so that it can provide a tight relaxation, and therefore generate a better heuristic solution, there is a need to search the space of Lagrangean multipliers λ for the tightest Z LD (λ), i.e., the optimal multipliers. The problem of finding optimal multipliers is denoted as the dual problem and is given by Equations (40) to (43) at a given iteration K.
where LB is a lower bound to the optimal objective function Z, α ∈ [0, 1] is the step size for the multiplier update, and the superscripts denote the iteration for the different quantities evaluated.
Since the lower bounds can be updated by any improved feasible heuristic solution, the full-space problem with fixed binary variables is solved as an NLP to obtain a feasible solution. The stopping criteria are set as the relative gap between the upper and lower bounds being less than , or the current number of iterations exceeding a maximum K max . If converged, the current lower bound is reported as the optimal solution. If not, the algorithm solves the dual problem to update the values of the multipliers and start a new iteration. The dual problem is developed based on [29] and [31], which is a hybrid formulation of cutting planes and sub-gradient step. Notice that the optimal solution of the multipliers λ will be used in the new iteration K + 1 as λ K . The main steps of the Lagrangean Decomposition are shown in Figure 5.

Case Studies
Three case studies for different horizon lengths, 12 days, 24 days, and 48 days are proposed to represent from short-term to mid-term horizon of the integration.
The layout of the refinery configuration is based on Figure 6. Three feedstocks, three storage tanks and three charging tanks are the upstream operations for the scheduling whereas several unit operations such as blenders for gasoline (GBLENDER), diesel (DBLENDER) and fuel oil (FOBLENDER), hydrotreaters for light cracked naphtha (LCNHT), coker light naphtha (CLNHT), and diesel (DHT), fluidized catalytic cracking unit (FCCU), delayed coker unit (DCU), stabilizer unit (STABILIZER), catalytic reforming unit (CRU), and delayed coker (DCU), are the downstream operations for planning. The distillation unit is composed of five towers: a flash distillation tower (FLASH), two distillation units (CDU and CDU2), a vacuum distillation unit (VDU) and a debutanizer unit (DEBUT).
The case data for Scheduling are shown in Table 2 and the data for Planning        The model statistics can be found in Table 7. In Table 8 we show the supply data at the scheduling level. The 12-day case will consider the scenario only in Period 1. 24-day case will consider Periods 1 and 2, while the 48-day case considers all four periods in the table.

Computational Results
The computational experiments are performed in GAMS 28.1.0, in a Dell Pow-erEdge T410 Ubuntu Server 16.04.2 LTS with an Intel® Xeon® CPU (24 cores) 2.67 GHz and a DIMM DDR3 Synchronous 1333 MHz (128 GB) RAM. The solvers used in this comparison are ANTIGONE 1.1 [26], BARON 18.5.8[32], and SCIP 6.0 [10] as global MINLP solvers. The solvers CONOPT 3.17 [7] and CPLEX 12.8 [15] were used when solving NLP and MILP problems within the solver DICOPT 2.0 [11,2]. In the implementation of the Lagrangean decomposition algorithm, we solved the corresponding MILP and NLP solvers using CPLEX and CONOPT, respectively.

Overall Solution
The overall solutions of the three cases, 12-days, 24-days, and 48-days, are presented in Table 9. The objective functions for both Monolithic and LD approaches are larger than for the Hierarchical approach (75%, 60%, and 6%, for the 12-days, 24-days, and 48-days cases, respectively), corroborating their advantage.
The CPU time of Lagrangean Decomposition refers to the total time required to obtain the bounds for all iterations. In all cases, the Hierarchical approach resulted in the smallest CPU time, since it is a non-integrated approach that does not consider the overall integrated CS and RP objective function. However, this method returned an infeasible solution for the 48-days case. Furthermore, even when a feasible solution was found, its objective function was much lower compared to the Monolithic and Lagrangean Decomposition approaches as seen in Table 9, as they aim to solve the integrated full-space problem and hence, there is a wider search space in the optimization.
To have a fair comparison between the two integrated approaches, Monolithic and Lagrangean Decomposition, the stopping criteria is set to be the optimality gap lower than = 1%. Lagrangean Decomposition can also be terminated when the iteration number exceeds K max = 100. Both approaches yield the same solution in the 12-day case whereas in the 24-day and 48-day cases the Lagrangean Decomposition converges to a slightly better near-optimal solution and much faster than the Monolithic Approach, as it iteratively solves smaller subproblems which is more efficient than considering the full problem at once. The detailed iterative results are presented in Section 6.3.
The Gantt charts for the solutions of these three cases solved by Lagrangean Decomposition are presented in Figures 7, 8, and 9. Notice that the total number of iterations ranges between 33 and 78. The x-axis represents the full refinery planning-scheduling time length (days), while the y-axis represents the crude-oil transfer operations marked in Figure 2.The two pairs of CDU charging operations, (v 11 , v 12 ) towards CDU1 and (v 13 , v 14 ) CDU2, are continuous. The non-overlapping scheduling constraints are satisfied. Moreover, the operation schedules in the first 12 days are different, which means that the inclusion of a longer planning horizon affects the early decisions in the optimal schedule, even if the scenario parameters for the first periods are identical for all three cases.
As can be seen in Figure 10, the production volumes are not constant through time. This shows that the integrated optimization is reacting to the variation in the crude oil availability from the scheduling part while maximizing the overall economic objective.

Hierarchical and Monolithic approaches
For the Hierarchical approach, the sequence of solving the two sub-problems CS and RP is a key decision. Both methods are tested in three cases and the computational results are presented in Table 10. The two options are stated as Scheduling-Planning and Planning-Scheduling, which explicitly represent the     Table 10 show that Scheduling-Planning is the best option in all cases. This reflects the fact that the Scheduling constraints are tighter and their penalty cost for using extra crude in the objective is higher than those in Planning. The full objective value returned by the Planning-Scheduling option is a poor negative value in the 12-day case, and in the 24-day and 48-day cases, infeasible solutions are obtained meaning that in terms of the solving sequence Scheduling is predominant compared to Planning, based on the three cases presented here.
For the Monolithic approach, the different alternative solutions come from the state-of-the-art MINLP solvers that are used. The computational results returned by a set of MINLP solvers are presented in Table 11. In the 12day case, all selected global MINLP solvers return the optimal solution with an objective value of 2.704 million dollars, whereas the Outer-Approximationbased convex solver DICOPT also returns the same objective value. Finding an optimal solution for the 24-day and 48-day cases is computationally much more expensive, as the optimality gap cannot be closed within the given time limit of 6 hours using the global solvers. In the 24-day case, the solver DICOPT was unable to find a solution to the problem given that the local NLP subsolver CONOPT failed at optimizing the continuous relaxation of the MINLP. When the termination condition is relaxed to be less than a = 5% relative gap, nearoptimal solutions are obtained within the time limit but still require a certain computational effort.

Lagrangean Decomposition properties
To obtain a better understanding of Lagrangean Decomposition, detailed primal/dual bounds against Lagrangean iteration numbers are presented in Figures  11, 12, and 13. The Lagrangean Dual objective Z LD is obtained by solving the dualized sub-problems and Upper Bounds are updated to be Z LD if it improves the current Upper Bound. As the full problem is a maximization problem, the Lower Bound is obtained by solving the full problem with the binary variables fixed, yielding an NLP. All three cases show that the lower bounds updated at  Figure 11: Computational results of bounds vs. Iteration Number in 12-day case around the fourth Lagrangean iteration are the optimal solution, although it can not be guaranteed based on the bounds obtained at that point. In the following iterations, the upper bounds are efficiently improved by updating Z LD , and the relative gap is closing to = 1% gap, which terminates the algorithm. Note that the Z LD might eventually fail to improve the current upper bound since the quality of Z LD is strongly dependent on the Lagrangean Multipliers λ K t,c that are used to formulate the Lagrangean sub-problems in the iteration K. The multipliers λ t,c are returned from solving the Dual Problem in iteration K and there is no guarantee that they always update in a direction that improves Z LD , given the non-convexity of the MINLP. The evolution of the multipliers with respect to the iterations are shown in Appendix C ( Figures C.1, C.2, and C.3). In all cases, after changing abruptly in the first iterations, the multipliers stabilize and converge in the final iterations. This behavior is desired since the first exploration steps allow the solution algorithm to explore the different feasible solutions while the last stable steps are there to find the optimality guarantees.
To assess the efficiency of the approach considered in this manuscript, we compare the time evolution of the lower and upper bounds predicted by the Lagragean decomposition with the ones obtained by the global MINLP solver BARON applied to the monolithic model. The bounds profiles for the cases with time horizon of 12, 24, and 48 days are presented in Figures 14, 15,    In all three cases, the Lagrangean decomposition approach was able to obtain a feasible solution and a valid relaxation (upper) bound in less than 100 seconds. Another observation is that the best solution found was always found by the Lagrangean decomposition approach. Notice that the relaxation gap for the Lagrangean decomposition is steadily improving while using BARON provides good quality starting bounds but these do not seem to improve considerably afterward. This is particularly noticeable in the 48-day case, where the upper bound predicted by BARON is virtually the same as the best solution found by the Lagrangean decomposition. Although BARON is not able to find such a good quality solution, using the validity of its bound we can conclude that the solution found by the Lagrangean decomposition, in this case, is the global optimal solution. The Lagrangean decomposition can converge to the desired optimality gap = 1% for the 12day case, and = 5% for the remaining cases, in significantly less time than BARON.

Conclusion
In this paper, both Crude-oil Scheduling and Refinery Planning are formulated as optimization problems and are integrated into a single full-space MINLP. The integrated formulation is constructed based on the overlap between the CDU charging operations in the scheduling level and the CDU feeding decisions  at the planning level. An MINLP priority-slot based continuous-time Crudeoil Scheduling formulation proposed by [28], and a multi-period discrete-time NLP Refinery Planning formulation are used. A multi-scale framework has been proposed to integrate different time representations, for which the Hierarchical, Monolithic and Lagrangean decomposition solution approaches were tested. Computational results were obtained on 12-day, 24-day, and 48-day cases regarding the time horizon length, where a 12-day case represents a short-term refinery planning and scheduling scenario, while 24-and 48-day represent the mid-term scenarios. In all cases, the Monolithic and the Lagrangean Decomposition approaches return better objective values compared to those returned by the Hierarchical approach. Significant economic and operational gains could be had by integrating the two-levels of decision making. Moreover, the proposed Lagrangean Decomposition can yield near-optimal solutions much faster compared to the Monolithic approach using global MINLP solvers.

A.2 Parameters
• H is the scheduling horizon

A.3 Variables
• Z tiv ∈ {0, 1} t ∈ T, i ∈ I, v ∈ W Z tiv = 1 if operation v is assigned to priority-slot i in time period t, Z tiv = 0 otherwise.
• S tiv ≥ 0, D tiv ≥ 0, E tiv ≥ 0 t ∈ T, i ∈ I, v ∈ W S tiv , D tiv , E tiv are the start time, duration, end time of operation v, respectively, if it is assigned to priority-slot i in time period t, S tiv = 0, D tiv = 0, E tiv = 0 otherwise.
tiv is the total volume of crude transferred during operation v if it is assigned to priority-slot i in time period t, V t tiv = 0 otherwise. V tivc is the volume of crude c transferred during operation v if it is assigned to priority-slot i in time period t, V tivc = 0 otherwise.
• L t tir and L tirc t ∈ T, i ∈ I, r ∈ R, c ∈ C L t tir is the total accumulated level of crude in tank r ∈ R S ∪ R C before the operation assigned to priority-slot i. L tirc is the accumulated level of crude c in tank r ∈ R S ∪ R C before the operation assigned to priority-slot i.
• y F trc t ∈ T, i ∈ I, r ∈ R C , c ∈ C y F trc is the total accumulated level of crude c ∈ C charged into CDU r ∈ R C in time period t.

Supplementary material
Appendix B Crude-Oil Distillation Unit Modeling To calculate CDU yields using non-rigorous models, the temperature distribution from the crude oil TBP curve may be used. The TBP curve, also known as crude oil assay, represents how the crude oil yields and properties (such as specific gravity and sulfur content) vary with the distillation temperate [19]. Due to operational limitations regarding the reflux rate and the number of stages, there is an overlap in the TBP boiling ranges of adjacent fractions in distillation columns [22], as can be seen in the TBP distillation curve presented in Figure  B.1. Besides, an example of crude oil assay data in which the temperature range is discretized in 89 small sections of 10 ºC each (micro cuts) is presented in Figure B.2. This data can be used to calculate how much material would be produced in each micro-cut range as well as its respective properties.

B.1 Fixed-Yield modeling
When the fixed-yield formulation is addressed, the crude oil assay data presented in Figure B.2 is used to calculate the yields and properties for the final cuts. Hence, the distillates are assigned to a certain range of micro-cuts, as shown in The yield for cut k from crude c can be calculated as a summation of microcuts cp which are assigned to cut k in Equation B.2. The set CU T k can be determined from Table B.1 and represents the final cut flows (fuel gas (FG), liquefied petroleum gas (LPG), liquid and heavy naphtha (LN and HN), kerosene (K), light and heavy diesel (LD and HD), and atmospheric residue (ATR)).
Volume-and mass-based mixing rules are used to calculate the properties of cut k from crude c. Equation B.3 represents the volumetric rules for SG, and Equation B.4 represents the mass rules for SUL. The implementation of these constraints includes a product of the property with the flow given that they are weight-based measurements.
Therefore, both yields and properties (volume-and mass-based) for the final cuts can be calculated from the crude oil assay data using Equations B.1 to B.4. Regarding the operational modes, M d01 refers to a standard mode (same yields as in the FY model) whereas M d02 refers to a mode which prefers lighter crude (higher yields for lighter cuts) and M d03 refers to a mode which prefers heavier crude (higher yields for heavier cuts). In Figure B.3, the yields for each cut in the three modes are given for crude C01. Yields for other crude are calculated in the same way. For the study case presented here, the fixed yield coefficients are presented in Table B.2.

B.3 Distillation unit modeling and optimization
To compare the performance of these CDU models, single-period optimizations are performed, which are expected to showcase the advantages and disadvantages of the methods previously described. In general, the following assumptions are considered in the modeling:

Md01
Md02 Md03 2) There are seven crude oil feedstocks and their assay data was generated using the process simulator PetroSIM; 3) Each crude availability is 40 Mbbl/day; 4) The CDU processing capacity is 100 Mbbl/day; 5) Two properties are tracked, specific gravity (SG), as a representative of volume-based properties, and sulfur content (SUL), as a representative of massbased properties; 6) The objective function maximizes the total net present value of the eight distillates (to each one has been assigned a potential value).
Three sets of potential values for the distillates are used: Normal-Producing (NP), Light-Crude-Dominated (LCD), and Heavy-Crude-Dominated (HCD), as shown in Table B.3. These three cases aim to be representative of real-world scenarios, such as seasonal changes of the crude marketing value. The models are implemented in the General Algebraic Modeling System (GAMS) version 25.1.3, using a standard machine with macOS High Sierra, 2.2 GHz 6-core Intel Core i7, and 16 GB Memory. Commercial state-of-the-art nonlinear solver BARON 18.5.8 has been used to provide global solutions. The models have around 100 variables and 100 constraints, and the optimizations are performed within 1 second returned by BARON.
The objective value for the three CDU models in each scenario is compared in Figure B.4. The x-axis refers to the NP, LCD, and HCD scenarios while the y-axis refers to the objective value of NPV. In general, the bar chart shows that the NPV is strongly determined by different cases, which refer to different  producing seasons. Given the nature of the selected 7 crude candidates from the assay data, the refinery will make more profit in the case of the Heavy-Crude-Dominated market. Besides, to compare three CDU models, we can see that the FY model returns the lowest value in all cases. MFY and SC return better values than FY, and SC is slightly better than MFY. If we independently analyze each case. In NP, the difference between MFY and SC is very close, while on LCD, the difference is more considerable. In Figure B.5, we give the crude cut solution profile in the NP case. The optimal solutions of crude cut volumes vary in three CDU models. The SC model decided to lose some production volume of HN and HD so that the CDU could contribute more to producing LD, which is more valuable in this case of NP. That's how SC could return the best objective NPV. As seen in Figure B.4, MFY, and SC perform better than FY, because the freedom of CDU operating condition is allowed in MFY and SC. Since a binary decision of mode selection is considered in MFY, and a continuous flow rate swinging is considered in SC, it is difficult to see how well MFY and SC will perform in different cases. As SC has a better performance in all three cases, this model will be embedded in our full-space refinery problem.  T1Crude3CDU1  T1Crude5CDU1  T1Crude7CDU1  T2Crude2CDU1  T2Crude4CDU1  T2Crude6CDU1   T1Crude1CDU2  T1Crude3CDU2  T1Crude5CDU2  T1Crude7CDU2  T2Crude2CDU2  T2Crude4CDU2  T2Crude6CDU2   T1Crude2CDU1  T1Crude4CDU1  T1Crude6CDU1  T2Crude1CDU1  T2Crude3CDU1  T2Crude5CDU1  T2Crude7CDU1   T1Crude2CDU2  T1Crude4CDU2  T1Crude6CDU2  T2Crude1CDU2  T2Crude3CDU2 T2Crude5CDU2 T2Crude7CDU2  T1Crude5CDU1  T2Crude2CDU1  T2Crude6CDU1  T3Crude3CDU1  T3Crude7CDU1  T4Crude4CDU1   T1Crude1CDU2  T1Crude5CDU2  T2Crude2CDU2  T2Crude6CDU2  T3Crude3CDU2  T3Crude7CDU2  T4Crude4CDU2   T1Crude2CDU1  T1Crude6CDU1  T2Crude3CDU1  T2Crude7CDU1  T3Crude4CDU1  T4Crude1CDU1  T4Crude5CDU1   T1Crude2CDU2  T1Crude6CDU2  T2Crude3CDU2  T2Crude7CDU2  T3Crude4CDU2  T4Crude1CDU2  T4Crude5CDU2   T1Crude3CDU1  T1Crude7CDU1  T2Crude4CDU1  T3Crude1CDU1  T3Crude5CDU1  T4Crude2CDU1  T4Crude6CDU1   T1Crude3CDU2  T1Crude7CDU2  T2Crude4CDU2  T3Crude1CDU2  T3Crude5CDU2  T4Crude2CDU2  T4Crude6CDU2   T1Crude4CDU1  T2Crude1CDU1  T2Crude5CDU1  T3Crude2CDU1  T3Crude6CDU1  T4Crude3CDU1  T4Crude7CDU1   T1Crude4CDU2  T2Crude1CDU2  T2Crude5CDU2  T3Crude2CDU2  T3Crude6CDU2 T4Crude3CDU2 T4Crude7CDU2