Global Optimal Scheduling of Crude Oil Blending ... - ACS Publications

Article pubs.acs.org/IECR

Global Optimal Scheduling of Crude Oil Blending Operations with RTN Continuous-time and Multiparametric Disaggregation Pedro M. Castro*,†,‡ and Ignacio E. Grossmann‡ †

Centro de Investigaçaõ Operacional, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213-3890, United States

‡

S Supporting Information *

ABSTRACT: This paper addresses the modeling of crude oil operations in refineries assuming that all properties blend linearly. Guidelines are given on how to generate a Resource-Task Network superstructure that implicitly handles the complex logistics, while extending the scope of a well-known continuous-time formulation to variable recipe tasks with multiple input materials. The new single time grid formulation has the advantage of avoiding computationally inefficient big-M constraints, unlike previously proposed unit-specific and priority-slot based models. Through the solution of a set of test problems from the literature, we show that the resulting mixed-integer nonlinear programs can be solved close to global optimality by the commercial solver GloMIQO for the objective of gross margin maximization but not for operating cost minimization. We also show that adopting a two-step MILP-NLP algorithm where the mixed-integer linear relaxation is derived from multiparametric disaggregation can reduce the optimality gap by orders of magnitude.

1. INTRODUCTION Economic and operability benefits associated with improved crude oil blend scheduling can reach multimillion dollars per year.1 Every refinery that charges a mix of crude oils to the atmospheric and vacuum distillation units has an optimization opportunity. From the crude oil purchases perspective, the faster the ability for a refinery to assess whether a particular spot crude oil purchase will result in a feasible operation, the better the refinery can capitalize on short-term market opportunities. From the operational perspective, variation in crude oil mixture quality and flow rate charged to the distillation columns affects downstream production units and is perhaps the single most influential disturbance to a refinery. Overall, crude oil blend scheduling optimization is a relatively inexpensive and timely way to improve performance. Early works on crude oil scheduling2,3 with mathematical programming techniques, addressed the problem with discretetime models. Shah2 proposed a decomposition procedure featuring mixed-integer linear programming (MILP) models for the downstream operation of tanks and crude distillation units as well as for the upstream portside operation. Lee et al.3 proposed a MILP for the whole system that can be viewed as a relaxation of the mixed-integer nonlinear programming (MINLP) model that is required to guarantee that there are no discrepancies in crude compositions inside the tanks with respect to their outlet streams. Additional constraints to prevent widely fluctuating charging rates to distillation units have been suggested by Hamisu et al.4 The turn of the millennium saw increased emphasis on continuous-time models. Jia et al.5 proposed a State-Task Network, unit-specific, nonconvex MINLP formulation. Rather than solving the full MINLP, a two stage MILP-NLP solution procedure was employed, featuring in the first stage a relaxed MILP model without the bilinear blending constraints, followed by the solution of the original MINLP after fixing all binary © 2014 American Chemical Society

variables. Such hierarchical decomposition into logistics and quality subproblems was motivated by the lack of suitable commercially available optimization software and computer power to solve the complete crude oil scheduling problem in reasonable time.1 A more rigorous unit-specific model with respect to the synchronization of time events with material balances is given in Furman et al.6 Moro and Pinto7 opted for a single time grid MINLP continuous-time formulation that was tackled with the augmented penalty version of the outer-approximation method.39 It exhibited a superior performance over the MILP formulation generated after discretizing the inventory levels of the tank farm. Sharing the time representation concept and the inability to provide global optima, Reddy et al.,8 proposed a MILP relaxation along the lines of Lee et al.3 combined with a rolling-horizon algorithm to eliminate the composition discrepancy. The underlying idea is that there is no composition discrepancy for the first slot on each tank. Since the MILP provides the composition in all tanks at the end of the slot, the schedule up to that time can be frozen and the entire procedure repeated for the remaining slots by solving a series of MILPs. Mouret et al.9 considered yet another type of continuoustime formulation, called single operation sequencing model, which represents a schedule as a sequence of priority slots. The advantage is that generic constraints can be used to break symmetries and reduce the size of the search space. The nonlinearities involved by composition constraints were handled by a two stage MILP-NLP procedure identical to Jia et al.,5 which, however, led to small optimality gaps due to the objective function that was selected namely maximizing Received: Revised: Accepted: Published: 15127

May 21, 2014 August 27, 2014 September 1, 2014 September 2, 2014 dx.doi.org/10.1021/ie503002k | Ind. Eng. Chem. Res. 2014, 53, 15127−15145

Industrial & Engineering Chemistry Research

Article

gross margins. Later,10 the MILP part of the model (without the blending constraints) was compared to models relying on different time representation concepts. Using a single time grid to keep track of events taking place was found to be a better approach, with the continuous-time representation scaling better with problem size, and hence outperforming its discrete-time counterpart. An important aspect is that the considered single grid continuous-time model features a set of big-M sequence constraints that are known to decrease the performance of mixedinteger models. Li et al.11 proposed a novel MINLP unit-specific model for marine-access refineries that contrary to inland refineries have just charging tanks to crude distillation units rather than a sequence of storage and charging tanks. Besides this structural limitation, the model incorporates a wide variety of realistic operational features, such as single buoy mooring, multiple jetties, marine vessels containing more than a single crude, brine settling, and the possibility of multiple tanks feeding a column, following the work of Li et al.12 In addition, 15 crude properties are considered that are modeled using linear additive blending indices on either volume or weight basis. A spatial branch-and-bound global optimization algorithm, that on each node uses the MILP-NLP two-stage strategy previously mentioned, is used to solve the MINLP problems. In this way, rather than removing the bilinear blending constraints, a tighter MILP relaxation is created through piecewise McCormick envelopes.13 Overall, the integer feasible solutions obtained for a set of seven test cases from Li et al.12 were guaranteed to be within 2% of global optimality. In this paper, we approach the solution of the crude oil scheduling problem with the Resource-Task Network based single time grid, continuous-time formulation of Castro et al.,14 which can handle both batch and continuous tasks. Contrary to most scheduling problems that feature a fixed recipe with known proportions of input materials to a processing task, the main goal of crude oil scheduling is to identify the optimal mix of low-cost and premium crudes that maximizes profit or minimizes cost while meeting operational constraints. In fact, the proportions of the different available crudes are unknown a priori and can change with time, leading to crude transfer tasks that can be viewed as being of the variable recipe type. More specifically, we show how to (i) extend the generality of the continuous-time formulation to handle variable recipe tasks with multiple input materials (of which crude oil blending is a particular case); (ii) derive a RTN superstructure that implicitly meets the mandatory logistic constraint of no simultaneous inlet and outlet streams from a storage or charging tank, something that is not straightforward whenever multiple inlet/ outlet streams are allowed; (iii) derive a tight relaxation of the bilinear blending constraints using multiparametric disaggregation15,16 (known to be tighter than standard17 McCormick and scaling better with the number of partitions than piecewise13 McCormick envelopes). Combining such elements leads to the main novelty of this work: a two-stage MILP-NLP algorithm that is capable of providing high quality solutions with low optimality gaps for the actual MINLP scheduling problem. As will be shown, this is particularly relevant for the objective of cost minimization. The rest of the paper is structured as follows. In section 2, we give the problem definition and the system configuration as a collection of crudes, marine vessels, storage and charging tanks, distillation units, and existing connections. In section 3, guidelines are provided on how to convert the real system into

the virtual system that is the RTN process model, featuring a variety of tasks and resources. The RTN-based mathematical model is then the subject of section 4, which is divided into linear and nonlinear parts, the latter providing two alternative relaxation approaches for the bilinear constraints that allow computing optimality gaps. The four test problems are provided in section 5, with section 6 describing the computational studies, including comparisons to previous work and commercial global optimization solvers with the objectives of gross margins and total cost. The paper ends with the conclusions in section 7.

2. PROBLEM DEFINITION In crude oil operations in refineries, it is required to schedule the unloading of crude marine vessels into storage tanks as well as the transfer of crude from storage to charging tanks and finally, to the crude oil distillation units (CDUs). Crudes with different properties are involved, and the objective is either to maximize the gross margin in the CDUs, which is accomplished by properly mixing the different crudes at the tanks so as to obtain the most profitable blends, or to minimize operating cost. The crude oil operation systems considered are the ones in Lee et al.3 and Mouret et al.9 Let mv ∈ MV represent a crude marine vessel, st ∈ ST a storage tank, ct ∈ CT a charging tank, cd ∈ CD a crude oil distillation unit, cr ∈ CR a crude oil, and pr ∈ PR a crude oil property (e.g., density). Given are the possible connections between vessels, tanks, and units and the following information. Marine vessels are characterized by arrival time atmv [day] and volume of crude to discharge vinmv,cr [kbbl = 103·bll]. For storage and charging tanks tk ∈ TK = ST ∪ CT, we must know the initial inventory v0tk,cr [kbbl] as well as the lower and upper bounds on storage capacity vtkmin and vtkmax [kbbl]. Tanks involving crude oil mixing must respect the composition limits max for the relevant properties: cmin tk,pr and ctk,pr. Crudes arriving or initially present in storage tanks have known compositions ccr,pr, which are assumed to blend linearly. Each charging tank is assigned a single crude blend that can be viewed as a final product. Final products are subject to minimum and maximum max demands dmin cr and dcr [kbbl]. Transfer flow rates between units max u,u′ ∈ U = MV ∪ ST ∪ CT ∪ CD are between ρmin u,u′ and ρu,u′ [kbbl/day]. In order to ensure that CDUs are operating continuously throughout the scheduling horizon H [day], min ρmin ct,cd > 0; for the remaining units, it is ρu,u′ = 0. Logistic constraints play a very important role, being assumed that (i) a single berth is available at the docking station for vessel unloading; (ii) blending tanks must be in standing-gage operation1 where there can either be flow in or out, not both (a generalization is the mixing-delay restriction); (iii) a charging tank may only feed one CDU at a time; (iv) a CDU can be charged by a single tank at a time. 2.1. Maximizing Gross Margin. Two alternative objectives will be optimized. The first is the maximization of the total gross margin of the distilled crude blends (Mouret et al.9) given the individual crudes margins pcr [$/bbl]. In order to reduce costly CDU switches between crude blends, the maximum number of distillation operations is set to nd. 2.2. Minimizing Operating Cost. The second and more challenging objective deals with the minimization of total operating cost (Lee et al.,3 Jia and Ierapetritou,40 Yadav and Shaik26). Given are sea waiting costs for marine vessels cwsea mv [$/day]; harboring costs for unloading the crude charb mv [$/day]; inventory costs for storage and charging tanks cinv tk [$/kbbl/day]; setup costs for crudes charged to the CDUs cchg [$]. 15128

dx.doi.org/10.1021/ie503002k | Ind. Eng. Chem. Res. 2014, 53, 15127−15145


Article

distillation units CD, that is, REQ = DS ∪ MV ∪ ST ∪ CT ∪ CD. A resource linked to either the tanks or distillation units can be viewed as representing the inlet and outlet pipes associated with the unit. Equipment resources have a maximum availability of one (unary resource), thus ensuring that simultaneous inlet and outlet transfers are forbidden as defined by the logistic constraints. The docking station resource and associated tasks will ensure discharge of a marine vessel only when the contents of its predecessor have been depleted. Material resources RMR include all the different crudes and their possible locations as they move through the system, see Figure 2. Taking crude A as an example, it can be located inside the upper marine vessel A_MV1, inside the top storage tank A_ST1, just after the storage tank A_ST°1 , just before the charging tanks (A_CT1i, A_CT2i) and finally, inside the charging tanks (A_CT1 and A_CT2). Note that the subscripts refer to the numbering of marine vessels and tanks. As for the final crude mixtures X and Y, they are respectively generated after blending A, B and/or C, or A, B, and/or D, and so, they can be viewed as existing just after the charging tanks. Crudes X and Y represent distinct distillation runs that cannot be produced simultaneous since their producing tasks consume the CDU equipment resource (CD). Tasks are primarily divided into fixed recipe and variable recipe tasks IVR, the latter including blending tasks IBL as a subset. In fixed recipe tasks, the standard tasks in STN/RTN scheduling models, the proportion of materials that is consumed/produced with respect to the amount handled by the task is known as a parameter. Fixed recipe tasks can either be continuous IC or storage tasks IS. Storage tasks can be viewed as a special type of batch task, lasting a single time slot regardless of the actual duration of the slot.19 Variable recipe tasks allow for the mixing of different crudes at a proportion to be determined by the optimization model (the selection of a single crude is naturally included). Blending tasks are not as flexible since crude compositions in the inlet stream must match those in the upstream tank. Overall, Figure 2 features all the tasks identified in the system configuration of Figure 1, as well as additional tasks that have the purpose of guaranteeing logistic constraints (i−ii). More specifically, storage tasks (e.g., Hold_in_vessel MV1) ensure that the docking station becomes free only when all the crude

Figure 1. Crude oil operation system for problem 1.

2.3. System Configuration. Figure 1 depicts the refinery configuration for the first problem presented in Lee et al.3 There are two marine vessels: one carrying crude A and another, crude B. Storage tanks are dedicated, the first to crude A and the second to B, and so, there is a one to one correspondence between marine vessels and storage tanks. Since there is a single crude inside a particular storage tank at any given time, the mass balance at the outlet splitter will be linear. In contrast, charging tanks have initially crudes C (upper tank) and D (lower tank) and can receive crudes A and B so as to prepare crude mixtures X (upper tank) and Y (lower tank). Bilinear mass balances will thus be involved. The distinction between tanks that are easy and difficult to model is made by using light and dark blue backgrounds, respectively. Note that in some of the other problems considered in this paper, blending constraints will be required for storage tanks as well, and there will be a case with only one crude in a charging tank. There will also be multiple CDUs and storage vessels with no incoming marine vessel.

3. RESOURCE-TASK NETWORK PROCESS MODEL System configuration and logistic constraints can be captured by a Resource-Task Network18 process model. In this section, we first show how to construct the process model for the problem in Figure 1, while later describing how to ensure that logistic constraints iii and iv, identified in section 2, are met when in the presence of multiple crude distillation units. Overall, a RTN describes the process as a set of resources r ∈ R interacting with a collection of tasks i ∈ I. The first step is to identify the resources, which for now can be divided in two subtypes. Equipment resources REQ include the docking station DS, marine vessels MV, storage ST, charging tanks CT, and

Figure 2. Resource-Task Network process model for problem 1. 15129



Article

column handles Y and Z. The charging tank in the middle can actually handle three crudes A, B, and C, besides the initial and final blends E and Y, but it cannot simultaneously feed the top and bottom columns, which is enforced by making transfer tasks 12 and 13 consume the unary resource CT2; see Figure 7. Even though it is assumed in the following that the crude discharge sequence from marine vessels is fixed to their order of arrival, the given process model allows for such decision to be determined by the optimization. This can be relevant when maximizing gross margin. On the other hand, non-negligible sea waiting costs make different discharge sequences suboptimal when minimizing operating cost. Finally, it should also be noted that while the RTN diagram can become very large with the increase in the number of crude oils and tanks (and hence difficult to show to a scheduler for industrial use), the complete diagram does not need to be generated. In fact, the actual implementation in GAMS features a module that automatically generates the required sets of tasks and resources (and corresponding structural parameters, see section 4.3) given the existing connections between marine vessels and tanks, as well as the location of the different crudes.

has been removed from the vessel. The nonsimultaneous inlet and outlet transfers constraint (ii) is straightforward to model in cases with a single active inlet/outlet stream. However, the existing links between the storage and charging tanks make it possible for multiple active inlet/outlet streams, see Figure 3,

Figure 3. Blending tanks can either be receiving or discharging crude. Storage/charging tanks allow a single/multiple inlet streams and multiple/single outlet streams.

something that is actually observed in the optimal schedules reported by Mouret et al.9 To model this feature, which in the case of inlet streams can be seen as continuous blending,20 auxiliary tasks have been defined that aggregate the outlet streams from storage tanks and the inlet streams to charging tanks. As an example, when active, Transfer_Out_ST 1 continuously produces resource A_ST1° that is consumed at the same processing rate by tasks Transfer 3 and/or Transfer 4 (notice that the maximum processing rate of Transfer_Out_ST1 [1000 kbbl/day] is set so as to match the sum of the maximum rates of the subsequent tasks). Similarly, Transfer_In_CT1 consumes a variable proportion of resources A_CTi1 and B_CTi1, at the same rate as Transfer 3 and Transfer 5. 3.1. Remarks. In Lee et al.,3 logistic constraint (ii) only applies to storage tanks involving crude oil mixing, whereas in Mouret et al.9 it also applies to storage tanks receiving the same crude quality from marine vessels. Figure 2 shows the complete RTN for the latter case, indicating the tasks and resources to remove for the former case. If blending storage tanks can hold various crudes, the Transfer Out tasks may consume multiple resources, thus becoming variable recipe tasks. This is the case of problem 3 in Figure 4. However, there is no need to define

4. MATHEMATICAL FORMULATION (MINLP) The mathematical formulation for optimizing the system of crude oil operations in refineries can be formulated as a mixedinteger nonlinear program (MINLP). Conceptually, it can be viewed as consisting of two parts. The scheduling element is responsible for the mixed-integer linear (MILP) constraints, while the blending constraints and the second objective function are of the nonlinear type (NLP). Compared to the multiperiod blending of refinery end products,21 blending decisions also change from one time period to the next but now the duration of each time period is to be determined by the model rather than being fixed. Furthermore, we do not know a priori how many time periods will be required to achieve the optimum. 4.1. Time Representation. Previous industrial case studies from paper plants22 and pipeline systems23 featuring mostly continuous tasks with variable processing rates have shown a better performance for the continuous-time model. As highlighted in the recent review paper by Harjunkoski et al.,24 continuous-time models have the advantage of requiring fewer slots to represent a schedule, which in this case is particularly important since the complex bilinear terms will feature a time index (Li et al.11 show a significant reduction in the number of bilinear terms compared to a discrete-time formulation). Furthermore, increasing the number of slots may lead to the execution of multiple rather than a single instance of a continuous task to achieve a certain production level, increasing solution degeneracy and hence possibly leading to more changeovers than those effectively required to achieve the same value of the objective function.25 A continuous-time model is thus preferred. The use of a single or multiple time grids remains to be decided. The latter have the advantage of further reducing the number of slots required to find the optimal solution and the drawback of requiring a large number of inefficient big-M constraints to model the interaction between different grids handling the same resource (check sequence constraints in the unit-specific scheduling models for crude oil operations in Jia et al.,5 Furman et al.,6 Li et al.,11 and Yadav and Shaik26). In this case, we opt to use a single time grid featuring t ∈ T event points that span from time zero to the given time horizon H,

Figure 4. Crude oil operation system for problem 3 of Lee et al.3

them as blending tasks since it suffices to enforce the blending constraint for the subsequent transfer tasks, see Figure 5. Notice that to model two main transfer tasks (4 and 5), we already require 25 material location resources. When multiple crude distillation units are involved, each CDU typically handles a subset of the crude blends. Taking the system of problem 2 as an example (Figure 6), the top distillation column can process blends X and Y, while the bottom 15130



Article

Figure 5. Partial RTN process model for problem 3, showing the interaction of storage tank ST1 with charging tanks CT1 and CT2.


see Figure 8. The size of each slot is determined by the optimization. 4.2. Model Variables. The scheduling part of the model requires three sets of binary variables: Ni,t indicates that task i is executed during slot t; Zit,mv assigns the harboring of marine vessel mv to event point t; Z0t,mv assigns the departure of marine vessel mv to event point t. All continuous variables are nonnegative: ξi,t represents the amount of material [kbbl] transferred by task i during slot t; for variable recipe tasks i ∈ IVR, we also need to know the amount consumed of each input resource r, through variables ξ̅r,i,t [kbbl]. Variables Tt hold the absolute time of event point t [day], identifying the start of slot t; excess resource variables Rr,t give the available amount of resource r at event point t; excess variables Rend r,t give the amount just before the end of slot t, that is, immediately before event point t + 1. The two types of excess resource variables have no units for equipment resources r ∈ REQ, while for material resources r ∈ RMR the values are in [kbbl]. In order to compute the number of distillation runs, two additional sets of variables are required: Xcd,t identifies if a sequence independent changeover of crude blend occurs from slot t to t + 1; NDcd gives the number of distillation runs for unit cd ∈ CD. Note that while these variables take respectively binary and integer values, they can be defined as continuous since they appear in constraints with only binary variables, as will be seen later on. The remaining four sets of variables are exclusive to the objective of minimizing operating cost: WSmv and HBmv give respectively the time waiting in sea and harboring for marine vessel mv; DTt represents the duration of time slot t; AItk,t gives the average crude inventory in tank tk during time slot t. 4.3. RTN Structural Parameters. The structural parameters link the RTN graph to the mathematical model. As recently described in Harjunkoski et al.,24 resources can interact with

Figure 7. Partial RTN process model for problem 2, showing the interaction of the charging tanks with the distillation columns.

tasks either discretely or continuously, and five sets of structural parameters are sufficient to model all types of interactions when dealing with fixed recipe tasks and a continuous-time formulation. Parameters μr,i and μ̅r,i are related to the binary extent variables Ni,t and give the discrete interaction at the start or end 15131



Article

upon harboring πiB_MV2,mv2 = 1000, together with the vessel itself: πiB_MV2,mv2 = 1. The equipment resource is then removed from the system upon departure of the marine vessel: πoMV2,mv2 = 1 disabling further crude transfers from it. 4.4. Scheduling Constraints. The RTN continuous-time formulation for the short-term scheduling of the crude oil system is given below. It neglects the blending constraints (given in section 4.5), and so, it can be viewed as assuming that each blending tank, storage or charging, is actually divided into multiple sub tanks, one for each crude type, of the same total capacity as the real tank. Such physical system allows for the compositions in the outlet stream to differ from those inside the tank, which cannot happen in reality. Segregation as a strategy is useful to reduce problem complexity by reducing the number of decisions to be made in terms of where crude oils should be stored.1 4.4.1. Excess Resource Balances. The structural parameters come into action in the excess resource balances, which are multiperiod material balance expressions involving the excess resource variables Rr,t and Rend r,t . In eq 1, the initial resource availability (first term on the right-hand side) is only considered at t = 1. For the remaining event points, the excess value at t is equal to that at the previous event point (t − 1) adjusted by the amounts discretely produced/consumed by all tasks starting or ending at t. The sixth term on the RHS reflects the material input to the system from incoming vessels, which is only accounted for in the case the arrival of marine vessel mv is assigned to event point t through variable Zit,mv. Similarly, the seventh term removes the marine vessel resource whenever the vessel departs at event point t.

Figure 8. In the given continuous-time formulation, events report to a single grid.

of the task, respectively. All dashed lines in Figure 2 linking tasks with equipment resources are of this type. More specifically, and taking the first transfer task as an example, we have μST1,Transfer1 = −1 (resource consumed at the beginning of the task) and μ̅ST1,Transfer1 = 1 (resource produced at the end of the task). The task also consumes/produces the first marine vessel (μMV1,Transfer1 = −1; μ̅MV1,Transfer1 = 1). Note that double arrows are used for a particular resource whenever there is both consumption at the start and production at the end. Parameters νr,i and ν̅r,i are also related to discrete interactions but are linked to continuous extent variables ξi,t. For the crude oil system considered, they are used solely by the Hold_in_vessel tasks to ensure that the docking station becomes free only when all the crude has been removed from the corresponding marine vessel. We thus have for the RTN in Figure 2: νA_MV1,Hold_in_vesselMV1 = −1 and ν̅A_MV1,Hold_in_vessel MV1 = 1, meaning that for each bbl handled by the task, 1 bbl is consumed at the start, and 1 bbl is produced at the end. Parameter λr,i is used for the interaction of the fixed recipe continuous transfer tasks with the material resources, which are represented as solid lines. Taking the Transfer 3 task in Figure 2 as an example, we have λA_ST10,Transfer3 = −1 and λA_ST1i,Transfer3 = 1, meaning that for each bbl being transferred between the storage and charging tanks, 1 bbl of crude A is continuously being removed from ST1 while continuously being supplied to CT1 (see also the explanation linked to eq 3). Variable recipes tasks with flexible proportions of output products appear in the generic batch-scheduling problem of Westenberger and Kallrath27 as mentioned in Kallrath.28 Blömer and Günther29 modeled this feature using a discretetime STN scheduling model and highlighted that flexible proportions of input goods may be considered in a similar fashion. Indeed, Castro et al.22 have used flexible inputs to determine the optimal amount of broke to recycle back to the tissue machines of a paper plant, where it is processed together with the corresponding raw material. The set of structural parameters defined for the RTN continuous-time scheduling formulation works well for the simplest case with just two input resources in the flexible recipe. Now, however, there can be up to six input resources in the flexible recipe, see task Transfer in CT2 in Figure 5, and so, an additional index needs to be used. Let λ̅r,r′,i represent the continuous change of resource r due to the consumption of resource r′ by variable recipe task i ∈ IVR. Taking task Transfer 7 in Figure 2 as an example, we have λ̅ r ,r,Transfer7 = −1 and λ̅ X ,r,Transfer7 = 1, ∀r ∈ {A_CT1,B_CT1,C_CT1}. Parameters R0r are used to set the initial inventory v0tk,cr of storage and charging tanks, as well as to define the initial availability of equipment resources other than marine vessels, for example, R0ST2 = 1. In order to model the interaction of marine vessels with the system resources, we rely on parameters πir,mv and πor,mv. Taking problem 1 as an example, the second marine vessel carries νinmv,B = 1000 [kbbl] that is made available to the system

R r , t = R r0|t = 1 + R rend , t − 1|r ∈ RMR + R r , t − 1|r ∈ REQ +

∑ (μr , i Ni , t|t ≠|T| + μr ̅, i Ni , t− 1) + ∑

+

νr , iξi , t

i∈I S t ≠|T |

i∈I

πri, mvZti, mv −

∑ mv ∈ MV t ≠|T |

∑

πro, mvZto, mv

∀ r ∈ R,

t∈T

mv ∈ MV t≠1

(1)

The discrete interaction at the start of storage IS tasks affects excess variable Rr,t, while the interaction at the end affects excess variable Rend r,t , see eq 2. The purpose of the former is to temporarily hide from the system, coupled with eq 3, the amount of crude raw material (RRM ⊂ RMR) that still has not been transferred from the marine vessel so that the docking station remains unavailable for the other vessels (notice that since tasks are not defined for t = |T|, eq 3 forces all marine vessels to be empty at the end of the time horizon). The purpose of the latter is to allow for the capacity constraints of the storage and charging tanks to be enforced, see eq 4, where subset RMR tk holds the material resources that can appear inside tank tk (e.g., A_CT1, B_CT1, and C_CT1 are the material resources associated with charging tank CT1 in Figure 2). The second and third terms on the RHS of eq 2 give respectively the amounts continuously produced/consumed by fixed and variable recipe tasks. For further details concerning the complex interaction of structural parameters with the two types of excess resource variables and the different types of tasks, the reader is referred to Castro et al.,25 where an illustration is given. 15132


Industrial & Engineering Chemistry Research R rend ,t = R r ,t +

∑ λr ,iξi ,t + ∑ ∑ i∈I

+

∑ νr̅ ,iξi ,t i∈I

i∈I

∀r∈R

MR

VR

,

r ́∈ R

MR

According to Figure 8, the first event point corresponds to time zero whereas the last corresponds to the time horizon H (eq 9).

λr̅ , r ,́ iξr̅ ,́ i , t

t ∈ T,

T1 = 0,

t ≠ |T | (2)

S

R r ,t = 0 vtkmin ≤

C

Article

∀r∈R

RM

∪R ,

max R rend , t ≤ vtk

∑

IO

t∈T

∀ tk ∈ TK ,

t ≠ |T |

r ∈ R tkMR

(4)

R r ,t = 0

∀ r ∈ RCD ,

t ∈ T,

t ≠ |T |

(5)

Equation 3 is also useful to model the condition of no intermediate storage for the material resources located just after/before the storage/charging tanks (r ∈ RIO), effectively ensuring that the auxiliary Transfer_Out_ST/Transfer_In_CT tasks in Figures 2 and 5 occur at the same rate as their successors/predecessors and guaranteeing with eq 1 that simultaneous inlet and outlet transfers on tanks are forbidden. Equation 5 certifies that distillation units (r ∈ RCD) process crude in every time slot, becoming available only at the last event point. 4.4.2. Timing Constraints. The timing constraints of the single time grid formulation state that the difference in time between two consecutive event points must be within the minimum and maximum durations of the occurring task. For tasks i ∈ INE not linked to an equipment resource (e.g., Transfer 3−6 in Figure 2), the slot duration must be greater than the amount of material transferred divided by the maximum processing rate (eq 6), where the ρmax are taken from the ρmax i u,u′ data. For the other continuous transfer tasks, previous research 14 has shown that writing the constraint per equipment unit leads to a tighter formulation. Equation 7 applies to the subset RTC ⊂ REQ of the equipment resources that excludes the docking station and marine vessels. Note that the summation on the RHS includes all inlet and outlet tasks of that particular unit (e.g., Transfer 1 and Transfer_Out_ST1 in Figure 2 for unit ST1). The absence of the docking station from the domain is because it is linked to the Hold_in_vessel batch tasks that always adjust to the duration of the time slot. As for the marine vessels, each is linked to a single transfer task that is already taken care in the constraints of the storage tanks. Tt + 1 − Tt ≥

Tt + 1 − Tt ≥

ξi , t ρi

max

∑ μr ̅,i i∈I

∀ i ∈ I NE ,

t ∈ T,

⎡ Zi ⎤ ⎡ Z noi ⎤ t , mv ⎥ t ⎥ ∨̲ ⎢ ∨̲ ⎢ mv ∈ MV ⎢ ⎣Tt ≥ atmv ⎥⎦ ⎢⎣Tt ≥ 0 ⎥⎦

ρi max

∀ r ∈ RTC ,

t ≠ |T |

(10)

∀ t ∈ T,

t≠1

(11)

The single docking station prevents marine vessels from harboring/departing simultaneously and so such occurrences can be linked to a single event point. In eqs 12 and 13 and other equations, the domain of the binary variables does not include the last/first event point for harboring/departure since sufficient time must be given for discharging the vessel contents, which takes at least one time slot. ∨̲ Zti, mv

∀ mv ∈ MV

t∈T t ≠|T |

(12)

∨̲ Zto, mv

∀ mv ∈ MV

t∈T t≠1

(13) 30

The convex hull reformulation of eqs 10 and 11 followed by simple algebraic manipulations31 to eliminate the resulting noo disaggregated as well as the Znoi t and Zt variables, leads to the MILP constraints in eqs 14 and 15. Tt ≥

t ≠ |T |

t ∈ T,

∀ t ∈ T,

⎤ ⎡ Zto, mv ⎥ ⎡ Ztnoo ⎤ ⎢ ⎥ ∨̲ ⎢ max ⎥∨̲ ⎢ in mv ∈ MV ⎢Tt ≥ atmv + ∑ vmv , cr / ρmv , st * ⎥ ⎢ ⎣Tt ≥ 0 ⎥⎦ ⎦ ⎣ cr

Zti, mvatmv

∑

∀ t ∈ T,

t ≠ |T | (14)

mv ∈ MV

(6)

ξi , t

(9)

Material supply from incoming vessels must be assigned to event points through binary variables Zit,mv. If an event point t is linked to the harboring of marine vessel mv, then the absolute time Tt must be greater than the arrival time atmv; otherwise, it is free to vary. This constraint can be written as a disjunction, where each term is exclusive (see eq 10). Notice that marine vessels are allowed to wait before harboring. Similarly, if the vessel departs at event point t, then the time must be greater than the arrival time plus the minimum discharge time, which is calculated dividing the volume of crude to discharge by the maximum processing rate to the downstream storage tank st* (eq 11).

(3)

t ∈ T,

T|T | = H

Tt ≥

t ≠ |T |

∑ mv ∈ MV

in ⎛ ⎞ ∑ vmv , cr Zto, mv⎜⎜atmv + crmax ⎟⎟ ρmv , st * ⎠ ⎝

∀ t ∈ T,

t≠1

(7)

(15)

Equation 8 applies only to the distillation units r ⊂ RCD, for which the inlet flow rate must be greater than a certain minimum value. Compared to the general constraint in Castro et al.14 that applies to semicontinuous flows (either zero or between the bounds), it does not require a big-M term since distillation units are always processing crude, a constraint enforced by eq 5.

Equations 16 and 17 represent the exclusive or in eqs 12 and 13.

Tt + 1 − Tt ≤

∑ μr ̅,i i∈I

ξi , t ρi min

∀ r ∈ RCD ,

t ∈ T,

∑

Zti, mv = 1

∀ mv ∈ MV

t∈T t ≠|T |

(16)

∑ Zto,mv = 1 t∈T t≠1

t ≠ |T |

∀ mv ∈ MV (17)

It is assumed that the unloading sequence of marine vessels follows their order of arrival. Equation 18 thus assigns harboring

(8) 15133



Article

events in an increasing order. t × Zti, mv ≤

∑

Xcd , t ≥ Ni , t + Ni ,́ t + 1 − 1

t × Zti, mv + 1

∑

∀ mv ∈ MV ,

i < i ́,

t∈T t ≠|T |

t∈T t ≠|T |

(18)

i < i ́,

4.4.3. Task Extent Constraints. Every time there is crude being transferred, binary variables Ni,t need to be activated so that the corresponding equipment units are consumed. Naturally, transfer tasks i ∈ INE not consuming equipment units do not need to be a part of eq 19. Parameter vmax is computed from the i max given values of vinmv,cr, vmax tk , and dcr . ∀ i ∈ I \I NE ,

ξi , t ≤ vimax Ni , t

t ∈ T,

t ≠ |T |

∑

∀ i ∈ I VR ,

ξr̅ , i , t

t ∈ T,

t ≠ |T |

∑

≤ ctkmax , pr

R rend ,t ≤

∑ r ∈ R tkMR

∑ ∑ r ∈ R tkMR cr ∈ CR r

R rend ,t

pr ∈ PR ,

t ∈ T,

Xcd , t ∀ cd ∈ CD

(

cr ∈ CR r ∈ R crCT t ∈ T t ≠|T |

∑ i ∈ I VR

∑

λr̅ , r , i =−1

cr ∈ CR r

dcrmax

t ≠ |T |

min

∀ r ∈ RFP +

(22)

cr ∈ CR r

4.4.6. Number of Distillation Runs. A distillation run involving a certain blend of crude can span multiple time slots. In order to determine the number of runs, it is required to relate the changeover variables Xcd,t to the relevant task binary extent variables Ni,t. Given that in all examples there are exactly two transfer tasks linked to each distillation column, we can write the logic propositions in eqs 23 and 24, where ICD cd holds the pair of tasks linked to column cd. These are valid since distillation units are always operating (no idle slots in between). Logic constraints can easily be converted to the inequalities32 in eqs 25 and 26. A unit’s number of distillation runs NDcd is then the number of changeovers plus one (eq 27).

i < i ́, Ni ́

,t

́

,t+1

⇒ Xcd , t

∑

wsea harb (WSmvcmv + HBmv cmv )+

∑

∑ ∑

AItk , t · DTt ·ctkinv

tk ∈ TK t ∈ T t ≠|T |

(NDcd − 1)·c chg

(30)

(i , i )́ ∈

∀ cd ∈ CD ,

t ∈ T \{|T | , |T | − 1}

To relate the sea waiting times with the time values of the event points of the grid, we add a constraint to each disjunction in eq 10. The big-M reformulation30 of eq 31 leads to eq 32, where we have used the tightest31 big-M parameters possible: H − atmv. ⎡ ⎤ ⎡ ⎤ Zti, mv Ztnoi ⎥∨̲ ⎢ ⎥ ∨̲ ⎢ ⎥ ⎢ ⎥⎦ mv ∈ MV ⎢ ≥ ∀ ∈ WS 0 mv MV ≥ − WS T at ⎣ mv mv t mv ⎦ ⎣ ∀ t ∈ T,

t ≠ |T |

(31)

WSmv ≥ Tt − atmv − (H − atmv) (1 − Zti, mv)

IcdCD ,

t ∈ T \{|T | , |T | − 1}

∧ Ni , t + 1 ⇒ Xcd , t

i < i ́,

∀ cd ∈ CD ,

(28)

(29)

cd ∈ CD

Ni , t ∧ Ni

ξi , t pcr )

NDcd ≤ nd

mv ∈ MV

∑

∑ i∈IC λr , i =−1

It should be highlighted that the scheduling model has no bigM constraints when maximizing gross margins. 4.4.8. Objective: Minimize Operating Cost. Equation 30 gives the definition of the total operating cost that includes sea waiting, harboring, inventory, and changeover costs. Notice that the second term is bilinear, involving the product of two continuous variables, the average inventory AItk,t and the slot duration DTt.

4.4.5. Product Demand. The amount available of final product resources RFP (e.g., X and Y in Figure 2) at the end of the time horizon must be within given minimum and maximum demands (eq 22). dcrmin ≤ R r , |T | ≤

ξr̅ , i , t pcr +

cd ∈ CD

(21)

∑

(27)

Specific to this objective function is the constraint that the total number of distillation runs is lower or equal than the maximum specified value (eq 29).

R rend , t ccr , pr

∀ tk ∈ TK BL ,

∑ ∑ ∑

max

(20)

4.4.4. Tank Compositions. The properties of the crude mixture inside a blending tank are calculated by adding for all constituents cr the product of the volume Rend r,t by the property ccr,pr. Recall that we are assuming that all properties blend linearly. Equation 21 then ensures that mixture properties are within given lower and upper limits, where subset CRr holds the crude associated with resource r. r ∈ R tkMR

∑

(26)

4.4.7. Objective: Maximize Gross Margin. Equation 28 defines the objective of maximizing the total gross margin of MR the distilled crude blends, where RCT represents the set cr ⊂ R of material resources inside the charging tanks that are linked to crude blend cr of value pcr (e.g., in Figure 2 these would be resources C_CT1 and C_CT2 for crude C). The two terms reflect the cases of the transfer tasks to the CDUs being variable (more common) or fixed recipe tasks (there is one such task in problem 4).

r ∈ R MR

ctkmin , pr

(i , i )́ ∈ IcdCD ,

∀ cd ∈ CD ,

t ∈ T \{|T | , |T |− 1}

(19)

λr , r , i =−1

(25)

t ∈ T \{|T | , |T | − 1}

NDcd = 1 +

Equation 20 states that the total amount processed by variable recipe tasks is equal to the sum of its components. ξi , t =

t ∈ T \{|T | , |T | − 1}

Xcd , t ≥ Ni ,́ t + Ni , t + 1 − 1

mv ≠ |MV |

(i , i )́ ∈ IcdCD ,

∀ cd ∈ CD ,

∀ mv ∈ MV , (23)

t ∈ T,

t ≠ |T |

(32)

The harboring duration for marine vessel mv is equal to the difference in time between the time values of harboring and departure event points, respectively, t and t′. This constraint can also be written as a disjunction (eq 33), being big-M reformulated into eq 34 (check section 5.1 in Castro and

(i , i ́) ∈ IcdCD , (24) 15134



Article

Grossmann31 for a closely related reformulation). Note that there is no need to use the equality on the left disjunctions of eqs 31 and 33 since both the sea waiting and harboring durations have positive coefficients in the objective function being minimized.

Taking i = {Transfer 4} in Figure 5 as an example, RBL i = {A−ST1,D−AT1} includes the resources corresponding to crudes A and D inside storage tank ST1, whereas resources r′ appearing in eq 37 are those consumed by the blending task, that is, for r = {A_ST1}, r′ ∈ RCR = {A−STO1 }, and for r CR O r = {D_ST1}, Rr = {D−ST1 }. For the charging tanks outlet transfers, resources r and r′ are actually the same (identified through the structural parameter λ̅r,r,i = −1). More specifically, for i = {Transfer 7} in Figure 2, RBL i = {A−CT1,B−CT1,C−CT1}, and so for r = {A_CT1}, r′ ∈ RrCR = {A−CT1}, and so on. It should be noted that there are two choices for modeling the blending constraint. We have opted to consider individual flows and volume fractions as model variables, the alternative being total flows and component compositions.41 This is possible since it is not required to know the compositions explicitly, only that the properties of the crude blends inside the tanks are within given lower and upper bounds, which is ensured through the linear constraint in eq 21, where the compositions of the base crudes ccr,pr are known parameters. The full nonconvex, mixed-integer nonlinear programming (MINLP) formulation for the crude oil scheduling problem thus consists of eqs 1−9, 14−27, and 37, coupled with either eqs 28−29 for gross margin maximization, or eqs 30, 32, and 34−36 for operating cost minimization. 4.6. Relaxation of Bilinear Term in Blending Equation. One way to obtain a good feasible solution to the MINLP problem is to first relax the bilinear terms in eq 37 using one of the methods described next and then to solve the resulting MILP to obtain either an upper bound (UB), when maximizing gross margin, or a lower bound when minimizing operating cost. In the latter case, the bilinear terms in eq 30 need also to be relaxed (see section 4.6.3). The binary variables are then fixed to the values from the MILP, and a constrained NLP version of the original problem is solved to find either a lower (LB) or upper bound. This approach is the main element of the specialized outerapproximation algorithm of Karuppiah et al.33 that allows to compute an optimality gap that can be improved iteratively by adding integer cuts to the upper bounding MILP. While this is, in general, an inefficient approach, for the test problems considered in this paper, the optimality gaps are near zero for the objective of maximizing gross margins and so Mouret et al.9 also relied on a single iteration. The nonlinear blending constraint equivalent to eq 37 was removed from their model to generate the upper bounding problem. In this paper, we follow the single iteration approach since the quality of the MILP relaxation resulting from multiparametric disaggregation15,16 is directly related to the accuracy level of the discretized variables, which is set by the user. Optimality gaps will be calculated for different settings, which could also be changed iteratively in the context of a competitive global optimization algorithm for bilinear up to mixed-integer polynomial programs.16,34,35 The standard McCormick LP17 relaxation will also be tested for comparative purposes. Note that due to the scheduling part of the model in section 4.4, both relaxation approaches lead to MILP problems. 4.6.1. Using McCormick Envelopes (LP). The standard approach to provide a relaxation to eq 37, is to replace the bilinear term ζi,t·Rr,t with a new set of continuous variables Wr,i,t (see eq 38) coupled with four sets of constraints that correspond to the McCormick envelopes.17 These generate the tightest

⎤ ⎡ Zi ∧ Zo ⎤ ⎡ Ztnoio ′,t t , mv t ′ , mv ⎢ ⎥ ⎢ ⎥ ∨̲ ∨̲ mv ∈ MV ⎢ ⎣ HBmv ≥ Tt ′ − Tt ⎥⎦ ⎢⎣ HBmv ≥ 0 ∀ mv ∈ MV ⎥⎦ ∀ t,

t′ ∈ T ,

t ≠ |T | ,

HBmv ≥ Tt ′ − Tt − (H −

t′ > t

(33)

∑ vmvin , cr /ρmvmax, st*) cr

× (2 −

Zti, mv

t ≠ |T | ,

Zto′ , mv)

−

∀ mv ∈ MV ,

t , t′ ∈ T ,

t′ > t

(34)

Tank inventory costs are calculated according to the trapezoid area,3 which requires the computation of the average inventory for each time slot, eq 35. In order to compute the accurate inventory cost, we need to know the duration of each time slot that is determined through eq 36. AItk , t =

∑

R r , t + R rend ,t

∀ tk ∈ TK ,

2

r ∈ R tkMR

t ∈ T,

t ≠ |T | (35)

DTt = Tt + 1 − Tt

∀ t ∈ T,

t ≠ |T |

(36)

As will be seen later on, the bilinear terms involved in the computation of the inventory cost have a major effect in model performance. In order to reduce the complexity, Jia et al.5 approximated the average inventory by adding the inventory levels and dividing by the number of event points in the grid, which translates to eq 37. The approximated operating cost then becomes eq 38. AItk =

∑

(R r ,1 +

r ∈ R tkMR

∑

min

∑

R rend , t )/ |T |

∀ tk ∈ TK

t∈T t ≠|T |

(37)

wsea harb (WSmvcmv + HBmv cmv )

mv ∈ MV

+

∑

AItk ·H ·ctkinv +

tk ∈ TK

∑

(NDcd − 1) ·c chg

cd ∈ CD

(38)

4.5. Matching Compositions Inside Blending Tanks and Their Outlet Streams. The nonlinear part of the model consists of a single set of constraints that ensures that a tank’s outlet stream has the exact same composition of the blend inside the tank. Let 0 ≤ ζi,t ≤ 1 represent the fraction of the volume inside the tank leaving it during slot t through blending task i ∈ IBL. It can be used to relate the amount of resource r ∈ RiBL available inside the tank at the start of the slot t and the amount of corresponding resource r′ ∈ RCR r consumed by the task through eq 37. Notice that the RHS is a bilinear term and that if the task is not active during slot t (Ni,t = 0), no material is being removed from the tank (ξ̅r,i,t = 0 due to eqs 19 and 20), and thus, the volume fraction must equal zero to meet the constraint.

∑

ξr̅ ′ , i , t = ζi , t ·R r , t

∀ i ∈ I BL ,

r ∈ R iBL ,

t ∈ T,

t ≠ |T |

r ′∈ R rCR

(37) 15135



Article

ζi,t associated with power k. The sum of variables R̂ r,i,t,j,k over all digits j must be equal to the original variable Rr,t for every power k, see eq 45. Equation 46 then ensures that a single digit j is active at position k. In case the digit is selected, the value of the disaggregated variable must be lower than the upper bound of the original variable, see eq 47.

linear relaxation for the given bounds on variables ζi,t and Rr,t.

∑

ξr̅ ′ , i , t = Wr , i , t

∀ i ∈ I BL ,

r ∈ R iBL ,

t ∈ T,

t ≠ |T |

r ′∈ R rCR

(38)

In this case, the lower bounds are zero and the upper bounds max are equal to 1 and Rmax gives the r , respectively. Parameter Rr maximum amount of crude that can exist within a blending tank. It can easily be calculated from initial inventory and that of incoming marine vessels and cannot be higher than the tank’s maximum capacity. The McCormick envelopes for this specific case are then given by eqs 39−42. ∀ i ∈ I BL ,

Wr , i , t ≥ 0

r ∈ R iBL ,

t ∈ T,

η

Wr , i , t =

∀ i ∈ I BL ,

t ∈ T,

r ∈ R iBL ,

Wr , i , t ≤

R rmax ·ζi , t

R r ,t =

∀i∈I ,r∈

R iBL ,

∀ i ∈ I BL ,

t ∈ T,

r ∈ R iBL ,

t ∈ T,

k ∈ {ψ , ..., η}

(45)

9

∀ i ∈ I BL ,

∑ Yi ,t ,j , k = 1

t ∈ T,

t ≠ |T |

k ∈ {ψ , ..., η}

j=0

t ∈ T , t ≠ |T |

r ∈ R iBL ,

(44)

∀ i ∈ I BL ,

∑ R̂ r , i , t , j , k

t ≠ |T | ,

r ∈ R iBL ,

(46)

(41)

Wr , i , t ≤ R r , t

t∈T

j=0

(40) BL

∀ i ∈ I BL ,

9

t ≠ |T |

t ≠ |T |

···

k=ψ j=1

(39)

Wr , i , t ≥ R r , t + R rmax(ζi , t − 1)

9

∑ ∑ j·10k ·R̂ r , i , t , j , k + Wr , i , t

0 ≤ R̂ r , i , t , j , k ≤ R rmaxYi , t , j , k

t ≠ |T |

t ∈ T,

(42)

t ≠ |T | ,

∀ i ∈ I BL ,

r ∈ R iBL ,

k ∈ {ψ , ..., η}, {0, ..., 9}

(47)

···

4.6.2. Using Multiparametric Disaggregation (MILP). Multiparametric disaggregation15,16,34,35 involves discretizing one of the variables of each bilinear term over a set of powers k ∈ {ψ,...,η}. In this work, we opt to discretize variables ζi,t, which are fewer in number than variables Rr,t, since a blending task i involves multiple resources r. The upper bound is also constant for all tasks i and slots t (equal to 1) and so parameter η = ⌊log10 1⌋ = 0. Parameter ψ is chosen by the user so as to reach a certain accuracy level for the discretized variables. To achieve a tighter relaxation than the one given by the McCormick envelopes: ψ ≤ −1, with the quality of the relaxation improving as the value of ψ decreases. Binary variables Yi,t,j,k identify if the volume fraction of blending task i at time t features digit j ∈ {0,...,9} in position k of the decimal numerical system. Since there always exists a gap between discretization points for a finite ψ, residual variable

The residual variables W r,i,t actually represent the bilinear ···

term arising from ζ i,t · Rr,t, being relaxed using the McCormick envelopes described in section 4.6.1; see eqs 48−51. Although not strictly necessary, the McCormick envelopes on the original variables (eqs 39−42) have also been added to improve the quality of the relaxation. Equation 38 completes the formulation. ···

∀ i ∈ I BL ,

Wr , i , t ≥ 0

r ∈ R iBL ,

t∈T ···

···

Wr , i , t ≥ (R r , t − R rmax ) × 10ψ + R rmaxζi , t r ∈ R iBL , ···

···

···

∀ i ∈ I BL ,

t∈T ···

Wr , i , t ≤ R rmaxζi , t Wr , i , t ≤ R r , t × 10ψ

0 ≤ ζ i,j ≤ 10ψ is added to obtain a continuous domain; see eq 43 and Figure 9.

(48)

(49)

∀ i ∈ I BL , ∀ i ∈ I BL ,

r ∈ R iBL , r ∈ R iBL ,

t∈T

(50)

t∈T (51)

4.6.3. Remarks. The AItk,t·DTt bilinear term in the operating cost objective function (eq 30) was relaxed in the same manner as ζi,t·Rr,t. For the multiparametric disaggregation relaxation, the option was to discretize the slot duration variables DTt, primarily because their values are in the order of units, closer to those of the other discretized variables (volume fractions ζi,t ∈ [0,1]). Such option reduces the number of binary variables when considering a global accuracy parameter ψ.

Figure 9. Volume fraction variables ζi,t are calculated as a sum of a ···

discrete point (illustration for level ψ = −1) and a residual variable ζ i,t, thus achieving a continuous domain. ···

ζi , t =

η

5. TEST PROBLEMS The four test problems considered are taken from Lee et al.3 with the data being provided as Supporting Information. It should be highlighted that we have multiplied the crude volumes by ten to be consistent with the motivating example.3 In order to be able to generate the same values for the objective function, the cost coefficient parameters were adjusted accordingly. The same system configuration is considered by Jia et al.,5 Ierapetritou et al.,40 Mouret et al.,9 and Yadav and Shaik26 with the exception of P3. While the configuration was not given in Lee et al.,3 it could be inferred from the optimal schedule that

9

∑ ∑ j·10k ·Yi ,t ,j , k + ζi ,t

∀ i ∈ I BL ,

t ∈ T,

t ≠ |T |

k=ψ j=1

(43)

The bilinear variable Wr,i,t is also written as a sum of an ···

approximation and a residual term W r,i,t. In eq 44, R̂ r,i,t,j,k is the disaggregated variable linked to Rr,t and the discrete value j of 15136



Article

The choice of the time representation and the assumptions made also affect the computational performance of the resulting mathematical formulation or optimization algorithm. The main alternatives that have been used in the literature are summarized in Table 2 and compared to the current work. While all listed references present the full set of mixed-integer nonlinear constraints, most focus on the mixed-integer linear programming relaxation of the problem. If there are no composition discrepancies between the inventory of blending tanks and their outlet streams, the solution of the MILP provides the optimal solution to the original MINLP. Otherwise, it is used to guide the search of either the NLP solver (after fixing the binary variables), or of the MINLP, leading to local solutions. Note that previous works have not attempted to solve the MINLP to global optimality.

Figure 10. Crude oil operation system for problem 3 in Mouret et al.9

shows MV2 discharging to ST3. It is also the only way of meeting the composition bounds on that storage tank. Mouret et al.9 do not enforce composition constraints on the storage tanks and use the configuration shown in Figure 10. The other references assume a more complex topology for problem 3 allowing marine vessels to discharge to multiple charging tanks even simultaneously, which is the only way to meet the composition bounds. This feature is not handled in this paper, and so, such problem is not included in comparison with Jia et al.,5 Ierapetritou et al.,40 and Yadav and Shaik.26 Typos in the data were also detected. In terms of system topology, problem complexity has two different sources. From the scheduling point of view, it increases with the number of tasks, since more tasks will typically require more time slots to represent the optimal schedule, leading to a more difficult problem. In this respect, the features in Table 1

6. COMPUTATIONAL STUDIES We now evaluate the performance of different algorithms that can be used to solve the crude oil blending MINLP scheduling problem. These include commercial solvers (i) DICOPT; (ii) BARON; (iii) GloMIQO; and the two-stage MILP-NLP strategies described in section 4.6, relying on (iv) standard McCormick relaxation of the bilinear terms; (v) multiparametric disaggregation (MDT). While all algorithms may fail to find the global optimal solution, DICOPT36 is the most limiting since it works on the assumption that the relaxed nonlinear problem is convex and does not compute an optimality gap. Through spatial branch and bound, BARON37 and GloMIQO38 can theoretically close the gap to zero and prove global optimality, while the same outcome can be obtained by the MDT relaxation for a sufficiently high discretization level16,34 (low ψ value). All mathematical formulations were implemented in GAMS, release 24.1.3, using solvers CPLEX 12.5.1 for the MILP relaxation problems and CONOPT 3.15L for the NLP problems. Solvers DICOPT, BARON 12.3.3, and GloMIQO 2.3 tackled the original MINLP. We have used a single thread and default options up to (a) relative optimality tolerance = 10−6 and (b) maximum computational time = 3600 CPUs. The hardware consisted of a desktop with an Intel i7-3770 (3.40 GHz) processor, with 8 GB of RAM and running Windows 7. In the case of DICOPT, the volume fraction variables ζi,t were initialized to 0.5 to guarantee a feasible solution in all runs. It should be noted at this point that the specified number of time slots |T| for the grid acts as a hidden constraint, meaning that the global optimal solution for |T| is not necessarily the real global optimum solution. In order to overcome this issue, we have considered the standard termination criterion of stopping after no improvement in the objective function following a single increment in |T|. In case of normal termination within the maximum resource limit, we get a best-found solution. For each test problem, we will be showing results for three different |T| settings, starting with the minimum value that ensures feasibility, and ending with lowest value that is able to find the best solution (note that the computational times are for the individual runs and not the accumulated values of all runs leading to |T|). An interesting point is that a small number of time slots contributes to maximizing the performance of the logistics part of the problem that can be viewed as minimizing the number of active movements (i.e., transitions or switchovers),29 something that is not included in the objective function. 6.1. Comparison to Literature Results. In the next subsections we provide detailed comparisons to relevant previous work. Note that the proposed model has been adjusted to

Table 1. Number of Tasks and Resources Per Problem P1 P2 P3 P4

tasks

resources

12 20 23 26

24 38 59 51


predict P4 as the most challenging. Note that according to the system topology and the problem data, the charging tank on the bottom in Figure 11 can only receive crude E, which will be the only component of final blend W. Such charging tank is thus not a blending tank, which explains its light blue background compared to the other charging tanks. From the blending point of view, complexity increases with the number of blending tanks and the number of components for each blend, which is very much related to the number of resources required to model the problem. From this perspective, P3, with its six blending tanks and a total of 59 resources, is the most complex. 15137



Article

Table 2. Main Features of Alternative Scheduling Formulations That Have Tackled the Problem Lee et al.3 time representation # time grids Flow in and out in storage tanks? Composition constraints on storage tanks? objective function inventory cost blending constraints in MILP relaxation Solution of rigorous MINLP attempted?

Yadav and Shaik26

Jia et al.5

Mouret et al.9

this work

discrete single no, if blending yes

continuous multiple yes, always yes

continuous multiple yes, alwaysa yes

continuous multiple no, always no

continuous single no, if blendingb,c yesa

operating cost rigorous (linear term) relaxed

operating cost estimate

operating cost estimate

gross margin not applicable

both rigorousc (bilinear term)

ignored

ignored

ignored

no

yes, fix binaries and solve NLP

yes

yes, fix binaries and solve NLP

relaxed with either McCormick or Multiparametric Disaggregation yes, to global optimality

The model of Yadav and Shaik26 actually presents the constraints for the different possibilities of mixing and material flow in a tank but assumes in the examples that flow in and out is always allowed. bCurrent model is configurable to match assumptions of Mouret et al.9 cCurrent model is configurable to match assumptions of Jia et al.5 and Yadav and Shaik.26 a

Table 3. Best Found Solutions for Gross Margin Maximization Showing Minimum Number of Slots, Optimality Gap, and Total Computational Time of MILP-NLP Strategy priority slots (Mouret et al.9)

RTN, single time grid (this work)

MINLP without blending constraints

blending constraints relaxed using multiparametric disaggregation

problem

# slots

gross margin [k$]

gap [%]

CPUs

# slots

gross margin [k$]

gap [%]

CPUs

P1 P2 P3 P4

5 6 9 4

7975 10117.5 8544.9 13254.8

0 0 2.2 0

0.35 0.93 5.96 0.47

6 9 7 6

7982.5 10246.1 8544.9 13254.8

0 0 0.37 0

1.02 55.6 3600 8.15

transfer from each marine vessel, whereas the optimal schedules for P1 and P2 given in Figures 12 and 13, feature respectively one and three discharge interruptions. In Figure 12, the filling task of storage tank ST1 (Transfer 1) starts at time zero and receives 50 kbbl of crude A from the first marine vessel at the maximum rate of 500 kbbl/day up to T2 = 0.1 days. ST1 then sends 300 kbbl of A to CT1 at the same rate up to T3 = 0.7 days, before receiving the remaining 950 kbbl of A. During that time, CT1 is also receiving 150 kbbl of B from ST2, following 50 kbbl of B from ST2. This mixing of crudes in a destination tank, one after the other, can be classified as batch blending.20 The resulting blend of 30% A, 20% B and 50% C is charged to the distillation column up to T5 = 4.396 days. The crude from CT2 in the first two slots is pure D and is transferred to the CDU at the minimum rate of 50 kbbl/day. In slots five and six, the crude composition from CT2 is 0.7% A, 52.8% B, and 46.5% D. The incoming crude from the second marine vessel is only discharged in the last slot, starting at T6 = 5.965 days, and so, it does not have time to contribute to the profit. At the end of the time horizon, ST1 holds 943 kbbl of A, ST2 is full, while CT1 and CT2 are almost empty. Interrupting the discharge of marine vessels is important in the context of a system that does not allow simultaneous inlet and outlet streams from a tank. In the case of P2 in Figure 13, emptying a vessel takes a minimum of 2 days, a time during which no crude can be sent from the storage tank to its downstream charging tanks. The no interruption constraint9 makes the system less flexible to generate optimal blends leading to a decrease in profit. As can be seen in the optimal schedule, ST1 receives 100.5 kbbl of A at the maximum rate, so that in the second slot, 300.5 kbbl can be sent to CT1, emptying ST1 in the process. CT1 also receives 399.5 kbbl of B during slots 1−2, to make 1000 kbbl of crude with composition 30.05% A,

match the conditions previously reported (to the best of our knowledge). 6.1.1. Priority Slots Continuous-Time Approach of Mouret et al.9 We start by comparing the results to the other work that has tackled the gross margin maximization objective function using a MILP-NLP solution strategy. The priority slots based, continuous-time formulation of Mouret et al.9 is freely available at the MINLP Cyber-infrastructure http://minlp.org/library/ problem/index.php?i=117&lib=MINLP and so the reported computational times in Table 3 are for the same hardware and software. Recall that their MILP relaxation ignored blending constraints whereas in this work, the blending constraints are relaxed using the multiparametric disaggregation technique16,35 (MDT). In Table 3, we show the best solutions obtained by the two approaches, listing the suitable number of time slots and the computational time for that specific run (bear in mind that the solutions from Mouret et al.9 were also found after solving the problem for various number of slots). Notice that global optimality was proven for all but P3. For this problem, we were able to reduce the gap between the lower and upper bounds by roughly 1 order of magnitude to 0.37% and new best solutions were obtained for P1 and P2 (0.1 and 1.3% increase in profit). The computational times required are nevertheless much longer, but this is due to the additional binary variables and constraints of MDT since with the McCormick relaxation the scheduling model can be solved in just a few seconds as well (see Table 6). The fact that the priority slots model by Mouret et al.9 is unable to find the same solutions in P1 and P2 (despite the zero gap) is due to cardinality constraints on the maximum number of slots for each transfer task from the marine vessels to the storage tanks, which are not part of the problem definition but lead to a simplification of the model. Indeed, the model9 in the MINLP Cyber-infrastructure shows a maximum of a single 15138



Article

Figure 12. Global optimal solution for problem 1 and gross margin maximization ($7,982,500).

39.95% B, and 30.0%, which feeds column 1 until the end of day 10. The second parcel from the first marine vessel occurs during slot 3 (T3 = 0.802, T4 = 0.962 days) and transfers 79.763 kbbl of A at near maximum rate to ST1, which is then sent to CT2, creating 959.90 kbbl of crude with composition 8.31% of A, 10.47% of B, 33.31% of C, and 47.91% of E. Tank CT2 feeds column 1 with pure E in the first couple of slots, switching to column 2 at T5 = 1.441 days and becoming empty at T9 = 8.144 days. Finally, part of crude C from the third marine vessel still makes it to column 2 through CT3, contributing to create a blend of 2.37% B, 74.84% C, and 22.79% F. The final remark about this solution is that it will not be possible in the subsequent scheduling period to keep the condition of always operating distillation columns, since the two possible charging tanks to column 1 reach the end of the time horizon with no crude inside. 6.1.2. Discrete-Time Approach of Lee et al.3 Lee et al.3 minimized the operating cost of the system considering an accurate inventory term. The objective function also includes harboring costs, and it is allowed for nonblending storage tanks to receive crude from marine vessels while simultaneously transferring it to the downstream charging tanks. It is thus no surprise that the interruption of marine vessels discharge does not appear in the optimal solutions. For the same number of time slots, the proposed continuoustime model can be seen as a rigorous relaxation of the discretetime model of Lee et al.,3 since the location of events is selected by the optimization rather than being fixed a priori. As can be seen by the results in Table 4, this is an important limitation of the discrete-time formulation that returns considerably worse solutions even when considering twice as many time slots.

Figure 13. Global optimal solution for problem 2 and gross margin maximization ($10,246,075).

Table 4. Best-Found Solutions for Operating Cost Minimization and Accurate Inventory Costs Lee et al.3

this work

discrete

continuous

problem

# slots

cost [k$]a

# slots

cost [k$]

reduction [%]

P1 P2 P3 P4

8 10 12 15

217.667 352.55 296.56 420.99

7 7 7 7

210.538 320.496 287.000 365.088

3.3 9.1 3.2 13.3

a

Reported values are for the MILP relaxation and hence may feature composition discrepancies at blending tanks that can possibly lead to further increments in cost when corrected.

We can thus conclude that the continuous-time representation is by far a better choice for this continuous process, which is consistent with the results in Castro et al.23 for pipeline systems. The most important cost reductions are achieved for P2 and P4. In our best found schedule for P2 in Figure 14, sea waiting times are equal to 1.4, 3, and 2 days, whereas in Lee et al.,3 MV1 wait times are longer (2 days), leading to a higher sea waiting cost ($35,000 vs $32,000). Harboring and changeover costs are the same in both cases ($48,000 and $150,000), the latter corresponding to two oil changes in CD1 and one in CD2. The main difference is thus the inventory cost, particularly due to 15139



Article

Figure 15. Best-found solution for problem 4 and operating cost minimization ($365,088).

Figure 14. Best-found solution for problem 2 and operating cost minimization ($320,496).

relaxation ignoring blending constraints (also given in ref 5). In their later work,40 only the costs for the MILP relaxation are reported, which are higher than those from current work. Yadav and Shaik26 allow for the possibility of incoming crudes from marine vessels bypassing storage tanks with the same crude (not allowed in this work). This seems to be the case in P1 and P2 and the reason why better solutions are obtained. The other possible explanation is given next. The approximate inventory cost function in eq 37 features in the denominator the number of event points in the grid, which is not a good option considering that the optimal number of events is changed in the search for the global optimal solution. More specifically, it was observed experimentally that after a certain number of slots, the optimization will “cheat” by duplicating event points corresponding to the points of minimum inventory so that they contribute more to the objective function and the total cost is reduced. This can be spotted easily by checking if multiple consecutive event points (e.g., t and t + 1) with the same absolute time value (e.g., Tt = Tt+1) exist. The number of time slots reported in Table 5 for our current work are before this happened. Using more slots leads to a further decrease in cost, $213,000 in P1 for |T| = 7 (6 slots) and $304,694 in P2 for |T| = 13. These values are already lower than those reported by Yadav and Shaik.26 The drawback of the single time grid formulation is that it requires a higher number of events to represent a schedule compared to the unit-specific approaches. The advantage is that

CT2 that is now empty most of the time. The other aspect worth highlighting is that time slots 3−5 do not start at the beginning of the days, which is not possible for the specified discretization level in Lee et al.3 (1 day). We also show the best-found solution for P4 in Figure 15, where the 67.6% of the reduction in cost is due to quicker harboring, with all vessels discharging in 1.2 days, compared to the 3 days in the seminal paper. 6.1.3. Unit-Specific Continuous-Time Approaches. The unit-specific continuous-time approaches of Ierapetritou and co-workers,5,40 and Yadav and Shaik,26 have solved the problem considering an approximate representation of inventory cost (eq 37) that has the advantage of being linear, compared to the nonlinear formula of the accurate cost in eq 35. This choice has a major impact in the complexity of the resulting problems, which are much simpler to solve. As can be seen in the last column of Table 5, GloMIQO can solve to global optimality in roughly 1 min the MINLPs resulting from our single grid formulation (previous unit-specific works do not report on attempts to solve the rigorous MINLP with a global solver). Note that simultaneous flow in and flow out is allowed in storage tanks to provide the same basis for comparison. Better solutions are obtained compared to Ierapetritou and coworkers5,40 in P1 and P2 but not in P4. The results in Jia et al.5 reported in Table 5 are for the MINLP model that oddly correspond to lower costs than the values from the MILP 15140



Article

Table 5. Best-Found Solutions for Operating Cost Minimization and Approximate Inventory Costs Jia et al.5

Jia and Ierapetritou40

Yadav and Shaik26

this work

multiple time grids

multiple time grids

multiple time grids

single time grid

problem

# slots

cost [k$]

CPUsa

P1 P2 P4

3 3 3

225.00 325.80 341.10

0.96 1383 21912

# slots

cost [k$]b

CPUsa

# slots

cost [k$]

CPUsa,b

# slots

cost [k$]

CPUsc

0.28 4.89 7.87

3 4

217 305.8

0.59 51.84

3

247.0 413.48 387.15

5 8 9

217.667 320.358 362.022

0.76 26.6 79.4

a

Computational times are for different hardware and software than current work. bReported values are for the MILP relaxation. cResults for global optimization solver GloMIQO.

Table 6. Computational Statistics for Different Solution Approaches, Commercial Solvers, and Gross Margin Maximizationa |T| = (5,6,5,7)b approach P1

P2

P3

P4

McCormick MDT (ψ = −1) GloMIQO BARON DICOPT McCormick MDT (ψ = −1) GloMIQO BARON DICOPT McCormick MDT (ψ = −1) MDT (ψ = −2) GloMIQO BARON DICOPT McCormick MDT (ψ = −1) GloMIQO BARON DICOPT

discrete variablesd

total variablesd

total equationsd

56 152 43

396 812 372

488 944 392

110 350 110

88 616 1056 88

168 528 168

828 1998 758

1081 3473 5313 941

1331 3203 1217

1119 2439 839

1495 4147 5731 935

1780 3922 1324

|T| = (6,7,8,−)c

margin [k $]

gap [%]

7975e

0.00e

7750 9000

no sol. 9000 8250

13254.76

0.00

0.00

0.00

no sol. 13215.56

CPUs 0.37 0.36 0.31 37.5 0.45 0.36 1.33 0.75 3600 1.33 0.39 0.99 3.80 0.55 2573 0.68 0.85 5.26 1.94 3600 1.94

margin [k $]

gap [%]

7982.323

0.00

9461.05 9639.475 no sol. 9524.5 8250 8369.601

8250 7922.857

3.19 0.00

2.37 0.005 0.002 0.001 8.67

|T| = (7,10,8,−)c CPUs 0.35 0.48 0.63 150 0.73 1.04 18.4 596 3600 1.98 0.61 7.15 268 3600 3600 0.86

margin [k $]

gap [%]

7982.5

0.00

0.18 7975 10246.08

no sol. 9000 8540 8544.891

no sol. 8250

0.00

2.29 0.56 0.37 0.70

CPUs 0.49 0.89 0.76 3600 0.44 4.27 46.0 16.5 3600 31.7 1.59 43.8 3600 3600 3600 1.32

a Optimal solution in bold. bNumber of event points used for P1 (5), P2 (6), P3 (5), and P4 (7), respectively. cResults for P4 for runs with |T| > 7 are not listed since the gross margin did not improve and the optimality gap remained zero. dWhile the problem size increases with the number of slots, the number of variables and constraints is only given for the lowest setting in order to save space. eTo facilitate interpretation of the results, values for gross margins and optimality gaps are merged for algorithms leading to the same outcome and are listed on the first row of the group (e.g., DICOPT is the only one not proving optimality for P1 at |T| = 5).

the constraints required to model inventory in tanks are much simpler, not requiring big-M terms. A trade-off is thus involved, which in this particular problem seems to favor the single time grid approach. This may not be the case in other refinery applications and it is certainly not the case for scheduling in general.24 6.2. Comparison to Different Optimization Approaches. Tables 6 and 7 provide detailed computational statistics for the objectives of gross margin maximization and total cost minimization, respectively. The comparisons are written for the two proposed MILP-NLP solution approaches, as well as for the solution of the original MINLP with commercial solvers. The most important result is that the objective function plays a major role. In fact, while three of the four gross maximization problems can be solved to global optimality, none of the operating cost minimization problems is even close to be solved to provable global optimality when considering the largest number of time slots (last three columns in Tables 6 and 7). This is primarily due to the quality of the MILP relaxation resulting either from McCormick envelopes that are also part of commercial solvers GloMIQO and BARON, or multiparametric

disaggregation (MDT). Overall, we conclude that the simplest MILP-NLP relaxation using McCormick envelopes (see first row for each problem) is acceptable for gross margin maximization, only failing to find the global optimum for P3 (2.29% gap), while operating cost minimization demands a global optimization algorithm, since all found solutions were highly suboptimal and the returned optimality gaps are as high as 112%. Using the two-stage MILP-NLP approach with the MDT relaxation for operating cost minimization (Table 7) leads to quite good optimality gaps for the minimum number of time slots (≤0.071%). A single increase in the number of slots is still manageable for P1−P3 (≤0.03%) but is already too much for P4 (25.8%). Increasing the number of event points to eight widens the gap to 7.8, 16.2, 48.1, and 55.2%, but such gaps are still below the ones reported by the solver GloMIQO, with the exception of P3 (48.1 vs 45.8%). Overall, the optimality gaps obtained by the proposed MDT approach can be up to 4 orders of magnitude smaller than the values from GloMIQO (P1 for |T| = 6). The two-stage MDT approach was able to find slightly better solutions than GloMIQO in P1 and P2, but returned suboptimal solutions in P3 and P4. This is to be expected in the 15141



Article

Table 7. Computational Statistics for Operating Cost Minimization with Accurate Inventory Costs |T| = (4,5,5,6) approach P1

P2

P3

P4

a

McCormick MDT (ψ = −1) MDT (ψ = −2)a MDT (ψ = −4)a GloMIQO BARON DICOPT McCormick MDT (ψ = −1) MDT (ψ = −2)a MDT (ψ = −4)a GloMIQO BARON DICOPT McCormick MDT (ψ = −1) MDT (ψ = −2)a GloMIQO BARON DICOPT McCormick MDT (ψ = −1) MDT (ψ = −2)a GloMIQO BARON DICOPT

discrete variables

total variables

total equations

36 165 255 435 36

309 921 1311 2091 279

434 1085 1424 2102 314

76 356 556 956 76

88 704 1184 88

110 520 820 110

678 2258 3258 5258 598

1139 4175 6295 975

1104 3929 5679 959

1015 2783 3683 5483 695

1676 5040 6892 1020

1692 4947 6572 1112

|T| = (6,6,6,7)

|T| = (8,8,8,8)

cost [$]

gap [%]

CPUs

cost [$]

gap [%]

CPUs

cost [$]

gap [%]

CPUs

245357 239000

69.4 0.46 0.046 0.0005 0.0001 0.0001

0.25 2.2 6.9 20.9 113 47.1 0.41 0.37 25.4 50.8 1185 3600 3600 0.67 0.59 42.6 555 3600 3600 0.74 0.90 1145 2420 3600 3600 3.15

230500 211588

74.6 0.30 0.029 0.0004 4.6 208

0.48 67.7 148 481 3600 3600 1.06 0.70 331 1549

233707 210538

77.1 7.8

0.46 3600

210838 230121 214000 360208 320496

45.3 421

3600 3600 3.52 1.38 3600

320504 no sol. 329400 292168 288867

37.0

382600 361800

374488 361800 338400 331700

354155 351267 395170

78.8 0.46 0.045 0.0006 9.6 126 50.4 0.26 0.026 0.031 215 88.9 0.71 0.071 11.9

no sol. 395170

229144 214000 348300 337950

71.6 0.32 0.03 16.3

no sol. 339659 310400 305800

no sol. 348200 388800 374102

59.2 0.22 0.022 9.4

104 25.8 33.4

no sol. 401498

81.9 16.2

3600 3600 1.31 1.03 497 3215 3600 3600 2.13 1.56 3600

287000 no sol. 300000 379020 366006

45.8

3600 3600 3.60

365088 no sol. 393486

112

62.3 48.1

99.3 55.2

3600 3600 6.03 5.38 3600 3600 3600 2.35 1.58 3600 3600 3600 9.96

No point in solving the problem for higher accuracy settings, if optimality cannot be proven in less than 1 h for ψ = −1.

priority-slot continuous-time formulation from the literature, the model allows interrupting crude discharge from marine vessels, leading to higher gross margins in two of the four test problems solved. The new MINLP formulation has shown to be computationally very efficient since GloMIQO could solve to global optimality in just a few seconds three of the gross maximization and all of the cost minimization problems when considering approximate inventory costs. This was no longer possible for accurate inventory costs but better solutions than those in the literature could still be found. It has also been shown that relaxing the bilinear terms appearing in blending constraints and inventory cost calculations using multiparametric disaggregation can lead to orders of magnitude lower gaps than GloMIQO. However, because of the use of a single starting point for the solution of the original MINLP, solution quality from the proposed MILP-NLP decomposition strategy was sometimes inferior. The main focus of future work will thus be on developing a spatial branch-and-bound global optimization algorithm featuring MILP relaxations derived from multiparametric disaggregation, so as to generate multiple high quality starting points for the solution of the original MINLP and thus increase the likelihood of finding the optimal solution in cases where the optimality gap cannot be closed in reasonable time. Other interesting research directions concern the handling of nonlinear blending rules for some of the crude properties, expanding the system boundary to include downstream refinery operations, and developing solution approaches for large-scale problems that allow quickly obtaining good feasible solutions.

presence of a weak MILP relaxation, considering that a single starting point is used to find a feasible solution to the original MINLP. In contrast, global solvers GloMIQO and BARON generate multiple starting points naturally, as part of the spatial branch-and-bound procedure, being thus more likely to find the optimum. GloMIQO was the best performer in terms of solution quality, which was not totally unexpected given that we are dealing with a recent solver that was developed for the solution of blending problems. In contrast, the poor performance of BARON came as a surprise. Notice in Tables 6 and 7 that BARON is incapable of finding a feasible solution in 12 out of 22 problems, proving optimality in just three of the remaining eight problems (all dealing with gross margin maximization). DICOPT, which assumes that the resulting MINLPs become convex once the binary variables are relaxed to continuous, can still find five best solutions for the lowest settings of |T| but becomes less effective with the increase in the number of slots. Still, it is a very fast solver, a feature shared with the McCormick MILP-NLP algorithm.

7. CONCLUSIONS This paper has shown that a resource-task network continuoustime formulation can effectively capture the most important scheduling constraints of the crude oil scheduling problem. The highlights include the handling of tasks with variable proportions of input materials, essential for blending, and the extensive use of material location resources to deal with the logistic condition that forbids simultaneous inlet and outlet transfers on storage and charging tanks. When compared to a 15142



■

Article

ASSOCIATED CONTENT

ST/st = storage tank T/t = Event points (time slots) of the single time grid TK/tk = storage and charging tanks TKBL = blending tanks U/u = system units excluding docking station

S Supporting Information *

Tables with all the necessary data for problems P1−P4. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

Parameters

AUTHOR INFORMATION

atmv = arrival time of marine vessel mv [day] ccr,pr = composition of raw-material crude cr in property pr cchg = cost involved for each change in crude to a distillation column [$] charb mv = harboring costs for unloading the crude from marine vessel mv [$/day] cinv tk = inventory cost for tank tk [$/kbbl/day] cmax tk,pr = maximum composition in blending tank tk for property pr cmin tk,pr = minimum composition in blending tank tk for property pr cwsea mv = sea waiting cost for marine vessel mv [$/day] dmax cr = maximum demand of final product crude cr [kbbl] dmin cr = minimum demand of final product crude cr [kbbl] H = time horizon [day] nd = maximum number of distillation operations pcr = gross margin of crude cr [$/bbl] R0r = initial availability of resource r Rmax mv,cr = maximum excess value for of resource r vinmv,cr = incoming volume from marine vessel mv of crude cr [kbbl] v0tk,cr = initial volume inside tank tk of crude cr [kbbl] vmax = maximum amount of material processed by task i i [kbbl] vmax tk = upper bound on total volume inside tank tk [kbbl] vmin = lower bound on total volume inside tank tk [kbbl] k η = position of most significant digit for discretized variables λr,i = continuous interaction of resource r during the execution of task i λ̅r,r′,i = continuous interaction of resource r due to resource r′ during the execution of task i μr,i = discrete interaction of resource r at the start of task i acting on variables Ni,t μ̅r,i = discrete interaction of resource r at the end of task i acting on variables Ni,t vr,i = discrete interaction of resource r at the start of task i acting on variables ζi,t vr̅ ,i = discrete interaction of resource r at the end of task i acting on variables ζi,t πir,mv = amount of resource r entering the system due to harboring of marine vessel mv πor,mv = amount of resource r leaving the system due to departure of marine vessel mv ρmax = maximum processing rate of task i [kbbl/day] i ρmin = minimum processing rate of task i [kbbl/day] i ρmax u,u = maximum transfer flow rate between units u and u′ [kbbl/day] ρmin u,u = minimum transfer flow rate between units u and u′ [kbbl/day] ψ = position of least significant digit for discretized variables

Corresponding Author

*Email: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Pedro Castro gratefully acknowledges financial support from the Luso-American Foundation under the 2013 Portugal−U.S. Research Networks Program, from Fundaçaõ para a Ciência e Tecnologia (FCT) through the Investigador FCT 2013 program, and from FEDER (Programa Operacional Factores de Competitividade-COMPETE) and FCT through project FCOMP-01-0124-FEDER-020764. Ignacio Grossmann acknowledges funding from the Center of Advanced Process Decision-making at Carnegie Mellon.

■

NOMENCLATURE

Sets/Indices

CD/cd = crude oil distillation unit CR/cr = crude oil CRr = crude associated with material resource r CT/ct = charging tank I/i = tasks IBL = variable recipe blending tasks IC = fixed recipe continuous tasks IS = fixed recipe storage tasks ICD cd = pair of transfer tasks to distillation column cd INE = transfer tasks not consuming an equipment resource IVD cr = variable recipe transfer tasks to distillation columns involving final product crude cr IVR = variable recipe tasks j = digits for decimal numeric representation system ∈ {0,...,9} k = positions in decimal representation of discretized variables ∈ {ψ,...,η} MV/mv = crude marine vessels PR/pr = crude oil property R/r = resources RBL = inside tank material resources located immediately i upstream of blending task i RCD = resources corresponding to distillation units RCT cr = material resource inside charging tanks corresponding to crude cr REQ = equipment resources RFP = material resources corresponding to final crude blends RIO = in or out material resources, located just after/before storage/charging tanks RCR r = out material resource (storage tanks) or inside tank resource (charging tanks) corresponding to the same crude of inside tank resource r RMR = material resources RMR tk = material resources than can appear inside tank tk RRM = raw-material resources RTC = equipment resources that appear in timing constraints

Variables

Ni,t = binary variable indicating if task i is executed during slot t (starts at event point t) Yi,t,j,k = binary variable indicating if the volume fraction of task i at slot t features digit j at position k 15143



Article

Zimv,t = binary variable indicating harboring of marine vessel mv at event point t i Zno = binary variable indicating that no vessel harbors at t event point t io Zno = binary variable indicating that no vessel arrives at t′,t event point t or departs at t′ o Zno = binary variable indicating that no vessel departs at t event point t Zomv,t = binary variable indicating departure of marine vessel mv at event point t AItk = average inventory of crude in tank tk during slot t [kbbl] AI tk = approximation of average inventory of crude in tank tk over the time horizon [kbbl] DTt = duration of time slot t [day] HBmv = duration of harboring for marine vessel mv [day] NDcd = number of distillation runs for column cd Rr,t = excess amount of resource r at event point t Rend r,t = excess amount of resource r immediately before the end of interval t R̂ r,i,t,j,k = disaggregated variable linked to Rr,t and discrete value j of Fi,t associated with power k Tt = absolute time of event point t [day] Wr,i,t = bilinear term variable involving resource r consumed by variable recipe task i during slot t

(9) 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. (10) Mouret, S.; Grossmann, I. E.; Pestiaux, P. Time representations and mathematical models for process scheduling problems. Comput. Chem. Eng. 2011, 35, 1038. (11) Li, J.; Misener, R.; Floudas, C. A. Continuous-time modeling and global optimization approach for scheduling of crude oil operations. AIChE J. 2012, 58, 205. (12) Li, J.; Li, W.; Karimi, I. A.; Srinivasan, R. Improving the robustness and efficiency of crude scheduling algorithms. AIChE J. 2007, 53, 2659. (13) Karuppiah, R.; Grossmann, I. E. Global optimization for the synthesis of integrated water systems in chemical processes. Comput. Chem. Eng. 2006, 30, 650. (14) Castro, P. M.; Barbosa-Póvoa, A. P.; Matos, H. A.; Novais, A. Q. Simple continuous-time formulation for short-term scheduling of batch and continuous processes. Ind. Eng. Chem. Res. 2004, 43, 105. (15) Teles, J. P.; Castro, P. M.; Matos, H. A. Multiparametric disaggregation technique for global optimization of polynomial programming problems. J. Global Optimization 2013, 55, 227. (16) Kolodziej, S.; Castro, P. M.; Grossmann, I. E. Global optimization of bilinear programs with a multiparametric disaggregation technique. J. Global Optimization 2013, 57, 1039. (17) McCormick, G. P. Computability of global solutions to factorable nonconvex programs. Part I. Convex underestimating problems. Mathematical Programming 1976, 10, 147. (18) Pantelides, C. C. Unified frameworks for the optimal process planning and scheduling. In Proceedings of the Second Conference on Foundations of Computer Aided Operations; Cache Publications: New York, 1994; pp 253. (19) Castro, P. M.; Barbosa-Póvoa, A. P.; Novais, A. Q. Simultaneous design and scheduling of multipurpose plants using resource task network based continuous-time formulations. Ind. Eng. Chem. Res. 2005, 44, 343. (20) Kelly, J. D.; Mann, J. L. Crude oil blend scheduling optimization: An application with multimillion dollar benefitsPart 2. Hydrocarbon Processing 2003, 82 (7), 72. (21) Kolodziej, S. P.; Grossmann, I. E.; Furman, K. C.; Sawaya, N. W. A discretization-based approach for the optimization of the multiperiod blend scheduling problem. Comput. Chem. Eng. 2013, 53, 122. (22) Castro, P. M.; Westerlund, J.; Forssell, S. Scheduling of a continuous plant with recycling of byproducts: A case study from a tissue paper mill. Comput. Chem. Eng. 2009, 33, 347. (23) Castro, P. M. Optimal scheduling of pipeline systems with a resource-task network continuous-time formulation. Ind. Eng. Chem. Res. 2010, 49, 11491. (24) Harjunkoski, I.; Maravelias, C.; Bongers, P.; Castro, P. M.; Engell, S.; Grossmann, I.; Hooker, J.; Méndez, C.; Sand, G.; Wassick, J. Scope for industrial applications of production scheduling models and solution methods. Comput. Chem. Eng. 2014, 62, 161. (25) Castro, P. M.; Harjunkoski, I.; Grossmann, I. E. New continuous-time scheduling formulation for continuous plants under variable electricity cost. Ind. Eng. Chem. Res. 2009, 48, 6701. (26) Yadav, S.; Shaik, M. A. Short-term scheduling of refinery crude oil operations. Ind. Eng. Chem. Res. 2012, 51, 9287. (27) Westenberger, H.; Kallrath, L. Formulation of a Job Shop Problem in Process Industry, Internal Report; Bayer AG, Leverkusen and BASF AG: Ludwigshafen, 1995. (28) Kallrath, J. Planning and scheduling in the process industry. OR Spectrum 2002, 24, 219. (29) Blömer, F.; Günther, H. Scheduling of a multi-product batch process in the chemical industry. Computers Ind. 1998, 36, 245. (30) Balas, E. Disjunctive programming and a hierarchy of relaxations for discrete optimization problems. SIAM J. Algebraic Discrete Methods 1985, 6 (3), 466−486. (31) Castro, P. M.; Grossmann, I. E. Generalized disjunctive programming as a systematic modeling framework to derive scheduling formulations. Ind. Eng. Chem. Res. 2012, 51, 5781.

···

W r,i,t = residual variable of bilinear term involving resource r, variable recipe task i and slot t WSmv = time marine vessel mv waits in sea before harboring [day] Xcd,t = continuous variable identifying if a crude blend changeover occurs for unit cd at the end of slot t ζi,t = volume fraction associated with the execution of variable recipe task i during slot t ···

ζ i,t = residual value of volume fraction linked to task t and slot t ξi,t = total amount of material handled by task i during slot t [kbbl] ξ̅r,i,t = amount of resource r consumed by task i during slot t [kbbl]

■

REFERENCES

(1) Kelly, J. D.; Mann, J. L. Crude oil blend scheduling optimization: An application with multimillion dollar benefitsPart 1. Hydrocarbon Processing 2003, 82 (6), 47. (2) Shah, N. Mathematical programming techniques for crude oil scheduling. Comput. Chem. Eng. 1996, 20, S1227. (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. (4) Hamisu, A. A.; Kabantiok, S.; Wang, M. Refinery scheduling of crude oil unloading with tank inventory management. Comput. Chem. Eng. 2013, 55, 134. (5) Jia, Z.; Ierapetritou, M.; Kelly, J. D. Refinery short-term scheduling using continuous time formulation: Crude-oil operations. Ind. Eng. Chem. Res. 2003, 42, 3085. (6) Furman, K.; Jia, Z.; Ierapetritou, M. G. A robust event-based continuous time formulation for tank transfer scheduling. Ind. Eng. Chem. Res. 2007, 46, 9126. (7) Moro, J. F. L.; Pinto, J. M. Mixed-integer programming approach for short-term crude oil scheduling. Ind. Eng. Chem. Res. 2004, 43, 85. (8) Reddy, P. C. P.; Karimi, I. A.; Srinivasan, R. A new continuoustime formulation for scheduling crude oil operations. Chem. Eng. Sci. 2004, 59, 1325. 15144



Article

(32) Raman, R.; Grossmann, I. E. Modeling and computational techniques for logic based integer programming. Comput. Chem. Eng. 1994, 18, 563. (33) Karuppiah, R.; Furman, K. C.; Grossmann, I. E. Global optimization for scheduling refinery crude oil operations. Comput. Chem. Eng. 2008, 32, 2745. (34) Castro, P. M.; Teles, J. P. Comparison of global optimization algorithms for the design of water-using networks. Comput. Chem. Eng. 2013, 52, 249. (35) Teles, J. P.; Castro, P. M.; Matos, H. A. Univariate parameterization for global optimization of mixed-integer polynomial problems. Eur. J. Oper. Res. 2013, 229, 613. (36) Kocis, G. R.; Grossmann, I. E. Computational experience with DICOPT solving MINLP problems in process systems engineering. Comput. Chem. Eng. 1989, 13, 307. (37) Tawarmalani, M.; Sahinidis, N. V. A polyhedral branch-and-cut approach to global optimization. Mathematical Programming 2005, 103 (2), 225. (38) Misener, R.; Floudas, C. A. GloMIQO: Global mixed-integer quadratic optimizer. J. Global Optimization 2013, 57, 3. (39) Viswanathan, J.; Grossmann, I. E. A combined penalty-function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14, 769. (40) Jia, Z.; Ierapetritou, M. G. Efficient short-term scheduling of refinery operations based on a continuous time formulation. Comput. Chem. Eng. 2004, 28, 1001. (41) Quesada, I.; Grossmann, I. E. Global optimization of bilinear process networks with multicomponent flows. Comput. Chem. Eng. 1995, 19, 1219.

15145


Global Optimal Scheduling of Crude Oil Blending ... - ACS Publications

Recommend Documents