An Approximate Method for the Optimization of Long-Horizon Tank Blending and Scheduling Operations
Benjamin Beach, Robert Hildebrand, Kimberly Ellis, Baptiste Lebreton
AAn Approximate Method for the Optimization of Long-Horizon TankBlending and Scheduling Operations
Benjamin Beach ∗ , Robert Hildebrand † , Kimberly Ellis † , and Baptiste Lebreton ‡ June 16, 2020
Abstract
We address a challenging tank blending and scheduling problem regarding operations for a chemicalplant. We model the problem as a nonconvex MIQCP, then approximate this model as a MILP usinga discretization-based approach. We combine a rolling horizon approach with the discretization ofindividual chemical property specifications to deal with long scheduling horizons, time-varying qualityspecifications, and multiple suppliers with discrete arrival times. This approach has been evaluatedusing industry-representative data sets from the specialty chemical industry. We demonstrate that theproposed approach supports fast planning cycles.
For speciality chemical manufacturing, companies purchase and blend supply streams from a variety ofsuppliers. Because the supply streams are derivatives of refinery operations, the chemical composition(specifications) of the streams may vary significantly by supplier as well as by season. Supply streams mayhave both desirable properties and less desirable properties that prevent their exclusive use during theproduction process. Furthermore, some streams might have limited available quantities or intermittentavailability. Thus, supply streams of differing composition are often combined into intermediate storagetanks. The resulting compounds are then blended to meet demand for production feed (see Figure 1).For our motivating application, supply streams typically arrive in barges from suppliers and may remain atthe inbound docks for several days before unloading into one or multiple storage tanks. The time windowsfor the arrival and departure of barges are committed several months in advance. As the supply bargesare unloaded into the tanks, the volume and chemical composition of the compounds in the storage tanksare changed. The compounds from one or multiple storage tanks are then blended to create the desired ∗ Grado Department of Industrial and Systems Engineering, Virginia Tech. Corresponding author, e-mail: [email protected] † Grado Department of Industrial and Systems Engineering, Virginia Tech ‡ Operations Research Team Leader, Eastman Chemical Company, Kingsport, TN, email: [email protected] a r X i v : . [ m a t h . O C ] J un INTRODUCTION quantity and specification of production feed. The production feed is planned in product campaigns(runs), typically spanning 3-14 days per product. Approximately 3 to 12 months of product campaignsare planned in advance due to the wide variety of spec requirements over time for very different non-concurrent production runs, the non-periodic nature of the supply acquisition process, and the limitedcapacity of the storage tanks.Given the volume and characteristics of the compounds on the inbound barges, the time windows for theavailability of barges, and the desired production feed (both quantity and specifications), the companymust determine an assignment of supply barges to tanks, an unloading plan from barges to tanks, anda blending plan from tanks to production feed. Although similar problems are often inherently multi-objective, the most important objective by far is to satisfy production requirements and to unload supply.Hence, the objective is to determine an operational schedule for the given barge arrivals and productionplan. Rapid solutions are required for long scheduling horizons to support effective planning decisions.Thus, we propose an efficient optimization-based approach to provide solutions for this assignment, mix-ing, and blending problem within ten minutes for a planning horizon of up to a year.Barge | S | Barge s v s , f s,q Barge 1 Tank 1Tank kv k,t , f k,q,t Tank | K | Production Feed y in s,k,t v prod t , f prod q,t Production y out k,t Figure 1: We schedule the the front-end portion of the problem, before the production stage (as indi-cated by the large box). This includes the scheduling of when and where to unload each barge and thedetermination of flow plans from tanks to production feed.Within the context of supply chain management, this problem encompasses short-term decisions as high-lighted in the Supply Chain Planning Matrix (Figure 2, adapted from [24]).2
PROBLEM DESCRIPTION
Figure 2: Context of this work in terms of supply chain planning (adapted from [24]).In this work, we employ a discrete-time mixed-integer linear program (MILP) approximate formulationof the core nonconvex mixed-integer quadratically-constrained program (MIQCP), formulated via directtracking of specs, using single-day scheduling periods with scheduling horizons of 90 to 368 days. The‘binary’ version of this discretization combines aspects of multi-parametric disaggregation (MDT) [60]with normalized MDT (NMDT, [4]), while avoiding the explicit tracking of spec qualities. The demandspec and spec ratio requirements are tightened to help ensure that the approximate representation of thespecs does not lead to violations of the demand composition requirements. To satisfy the long-horizonscheduling goals within a reasonable computation time, we use a rolling-horizon solution approach. Afteroptimization of the model approximation, the scheduling solution is simulated to obtain the correct tankinventory and demand feed specs.In the following sections, we provide a formal problem description, highlight the related research, formulatea mixed integer non-linear representation of the problem, propose approaches for solving the model, andcompare the approaches with computational studies based on industry-representative data sets.
For a given a set S of supply nodes (barges), assume each supply node s contains a given volume v s of raw material, with a level for chemical property q of f s,q , where q ∈ Q . The levels of these chemicalproperties are often specified as a percentage of the volume for a material component of interest. Thepercentages may not necessarily sum to 100% of the material, because some components may overlap andothers may not be of interest. Each supply node s is available during a set T s of consecutive time periods,from t min s to t max s . With the availability of all supply nodes known, the set of supply nodes available each3 PROBLEM DESCRIPTION time period t is denoted as S t . The number of barges that are unloaded each day is limited to ul max . If abarge is unloaded in a given day, the minimum percentage of inventory that can be unloaded is ulp min s . Ifa supply barge has less than ulp min s remaining, then the inventory is not unloaded and the barge departs.A set K of intermediate storage tanks is available to hold raw material or a blend of raw materials. Eachtank k has an initial inventory volume v k, , minimum inventory level v min k , and maximum inventory level v max k . The initial quality level for chemical property q is denoted by f k,q, . Each tank can be suppliedby the supply nodes in K s . If the mixture in a tank k is used to fulfill production feed, the minimumpercentage of feed that can originate from the tank is tkp min k . The demand for production feed is knownfor each time period t ∈ T . This demand is characterized by a volume d t , bounds f min t,q and f max t,q onthe quality level of the chemical properties, and bounds r min t,q ,q and r max t,q ,q on certain ratios of chemicalproperties. Production feed is planned in campaigns or runs of consecutive days. During a run, f min t,q andf max t,q remain constant, and we require the daily flow plan for this feed to remain constant.In order to meet production feed d t , each time period, the volume of material transferred needs to bedetermined: • volume to unload from supply node s to tank k during time t , y in s,k,t and • volume to unload from tank k to meet production feed in time t , y out k,t .In addition, the following assignment variables on these flows are needed to address various operationalconstraints: • if supply barge s is unloaded during time t , then γ s,t = 1. • If tank k provides mixture to meet production feed in time t , then σ k,t = 1.To model our objective of meeting demand and unloading requirements, we maximize a weighted sum ofthe amount of demand met and product unloaded.Additional operational constraints include the following. The number of times a barge can be unloadedis limited to two, and all unload operations for the barge must be performed within a range of ult max time periods (e.g., seven days). Flow from tanks to meet production feeds must be continuous for eachcampaign (or demand run).Conservation of flow must be maintained, such that the sum of the initial tank inventory level and thevolume of supply that flows into the tank less the volume that flows out for production feed equalsthe final tank inventory level. The resulting levels of chemical properties from this flow must meet thespecified composition requirements. Due the relatively small sizes of the storage tanks in this work, weallow storage tanks to both receive supply and provide material for production feed simultaneously. Ifthis occurs in a given time period, for the sake of simplicity, we model that complete mixing of inflowsoccurs before any outflows within that time period. The problem and constraints are formally specifiedin Section 4. 4 LITERATURE REVIEW
The blending and scheduling problems, also known as multi-period pooling problems, studied in theliterature vary widely in network size, scope, objective, and detail. Such problems, similar to the standingpooling problem, are typically characterized as bilinear nonconvex MIQCP’s due to the balance of specs(or properties) in mixing and splitting of product. This is referred to as ‘linear blending’, and is oftenan approximation of what occurs in practice (as e.g. the mixing process may not finish completelybefore material is split into multiple streams). However, similar problems can also involve linear [33]representation for spec balance, or a more highly nonlinear [65, 58, 42, 13] representations for balance ofsome specs (‘nonlinear blending’, as in the extended pooling problem), depending on the nature of theoperational network or the specifications to track.The underlying resource networks can be small and separated into discrete layers, and can representfront-end acquisition, as in e.g. [45, 50, 11, 9, 6] and the current work; back-end operations and sales,as in e.g. [13]; or the full process from acquisition to sales, as in [45, 64]. The networks can also spanlarger, more complex supply chains, as in [52, 2]. The scope of the problem can include scheduling only,as in [44, 41, 30, 6] and the current work, but can also extend to include more high-level acquisition [50],production planning decisions [38, 37, 61], or more detailed equipment control decisions [15]. Some works,such as [22, 23], integrate all three levels of planning, scheduling, and control in a single model.Common examples of objective functions used in similar problems include maximization of profits [53,54, 34, 29, 9, 3, 36, 2, 64] or gross margins [45, 50, 7, 16, 6], or minimization of operating costs [45, 33,35, 32, 7, 6]. The objective function used in [14] in particular is similar to that in this work, but includesadditional penalties for performing mixing operations, missing the demand composition requirements,and for not respecting the maintenance requirements. [13] maximizes the yields of final products. Manyproblems have only a single type of supply source, which can be characterized either by discrete arrivals[45, 50, 11, 7, 6] or continuous supply streams e.g. [33, 32, 13, 3, 16, 36]. Some works require multipletypes of supply sources. For example, in [53, 54, 34, 9, 64], there are both small barges with a singletype of inventory and larger barges, referred to as very large crude carriers (VLCCs), that carry multiplepackages with varied compositions. Some works explicitly include vessel scheduling decisions to be madeat unloading points, such as docks [53, 34]. Most works assume known supply and demand quantities orbounds, but some, such as [48, 38, 13, 19, 12, 15, 49], explicitly take into account stochastic supply anddemand.Muti-period pooling problems classically come in the forms of P, Q, or PQ formulations. The Q formula-tion is source-based, the P formulation is spec-based, and the PQ is a hybrid of the two. That is, the Qformulation tracks compositions based on origin fractions or volumes, while the P formulation tracks com-positions explicitly based on the concentrations of the physical qualities to track. See [19] for a discussionof these models. Some examples of spec-based formulations can be found in [13, 11, 3], while examples ofsource-based formulations can be found in [53, 45, 16, 9]. Comparing to the traditional pooling problem, aspec-based formulation mirrors the traditional P-formulation while a source-based formulation mirrors theQ or PQ-formulations. This choice of representation affects both the required number of bilinear termsand the MILP tightness of the formulation. Overall, the number of bilinear terms for a source-based for-5 .1 Comparison to Similar Research 3 LITERATURE REVIEW mulation scales with the number of mixing or splitting nodes in the network and the number of distinctsource compositions, while the number of bilinear terms for a spec-based formulation scales instead withthe number of specs that require tracking. As such, this choice of representation can have a significanteffect on both the tightness of the MIP and the number of required bilinear terms. While source-basedformulations can be tighter [36], they are only viable in situations where the number of distinct sources isrelatively small. In this work, as many as 35 distinct source compositions can be represented within thescheduling horizon, rendering source-based formulations computationally expensive.
We highlight and compare three research streams that are similar to this current work. Note that manydifferences between this work and others stem from the specialty chemical nature of our problem, combinedwith the properties of the production site for the application of interest.Castro [6] studies a similar problem motivated by crude oil operations in refineries. This problem isalso addressed in [31, 46, 45, 7], among others. The problem involves barge-based supply where thebarges carry a single type of crude, arrive within given time intervals, and unload according to the arrivalsequence. The barges unload to storage tanks. The material in the storage tanks are mixed into chargingtanks that supply crude oil into charging tanks that supply crude oil distillation units (CDUs).The similarities of the core problem include a similar network structure, a barge-based supply stream, andthe capacity of each storage tank ( ∼ ≤
15 days, ≤ ∼
30 barges) in the current work necessitates a heuristic such as a rolling horizon scheme. Atthe same time, the larger network in [6] reduces the length of the scheduling horizon over which theproblem can be solved within a reasonable time. A significant observation of [6] is that, for the problemconsidered, the choice of objective function significantly impacted the relative efficacy of the discrete-timeand continuous-time formulations. Another key difference is that the production schedule for each CDUis not fixed in [6], but is chosen based on gross margins per crude or on operating costs. On the otherhand, a CDU can be sourced from only a single charging tank at a time.In Castro’s work, the source-based choice is convenient because there are few composition types to track.In the current work, the source-based formulation is prohibitive due to the large amount of differentpossible composition types, and a spec-based formulation is used instead. As mentioned in [6], whenviable, source-based formulations yield tighter MILP relaxations [36], and can yield better computational6 .2 Related Solution Methodologies 3 LITERATURE REVIEW performance [3].Castro [3] addresses a similar multi-period pooling problem with composition constraints on the flowsto demand stream (tested on both per-flow and a mixture bases). Similar to the solution approachused in [6], Castro employs an iterated MILP-NLP approach, where the MILP stage is represented usingbase-10 NMDT. However, the utilized objective functions and formulations differ. In [3], a discrete-timeformulation is used, and both spec-based and source-based formulations are presented, and comparedwith commercial global solvers. The objective function is profit maximization as opposed to gross marginmaximization or cost minimization. The problem instance tested has two source streams, up to fourblending/storage tanks for which tank-to-tank transfer operations are allowed, and two demand streams,and are solved with a time horizon of up to four periods.Oddsdottir [50] addresses a similar network, with the same structure (barges to storage tanks to demand(CDUs)). The main differences in this respect involve barge availability (less restrictive time windows),the number of storage tanks (6), and the number of demand streams (2). Another common feature isthe requirement of a much longer scheduling horizon (3 months) and in the flexibility in barge arrivals.However, the goal in [50] is related to procurement planning instead of operational scheduling. As such,the choice of which barges to order and when is a problem decision, and a courser time-grid is used (threedays per period), resulting in a total of 30 time periods to model. Additionally, in both works, multiplebarges might be unloaded within a single time period (up to a limit), possibly contributing to the choiceof a discrete-time formulation. Further, a given barge may unload to multiple storage tanks in a giventime period. As with [6], however, the limited number of crude types in the system in [50] allow for areasonably compact source-based formulation, which would be unwieldy in the current work.
Rolling planning aside, the solution approach selected in [50] is comparable to the approach in [6], butwith some key differences. As in [6], a source-based formulation is used and composition constraintsare formulated with respect to crude volumes instead of crude volume fractions. However, instead ofdiscretizing, [50] applies a MILP approximation to the problem in which the constraints enforcing correctproportions for the outflow composition are reformulated in a manner that allows composition discrepancyin the outflows. A nonconvex MIQCP stage is then applied to resolve the composition discrepancy, yieldingnear-optimal feasible solutions in a much shorter time frame than the full nonconvex MIQCP: Running inGAMS with BARON, over 20 datasets, the average solution gap for this PRONODIS method (comparedto the full nonconvex MIQCP) was 1.1%, with an average computation time of 20s (compared to over13000s for the full nonconvex MIQCP).A wide variety of formulations and approaches have been developed to tackle these difficult schedulingproblems, with varying time-scheduling formulations, spec representations, and algorithms to deal withthe problem nonlinearity, which usually incorporate some form of MILP approximation to the basicnonconvex MIQCP. 7 .2 Related Solution Methodologies 3 LITERATURE REVIEW
A number of continuous-time and discrete-time scheduling formulations have been presented. The choiceof time-formulation can have a significant impact on problem performance, as showcased in works suchas [54] and [46, 45]. The preferred scheduling formulation can depend on the type of objective functionused, as observed in [6], as well as on the types of operational constraints to be considered. Compared tocontinuous-time formulations, discrete-time formulations frequently yield simpler representations, with atrade-off of more integer variables. Even so, problems with certain types of operational constraints, suchas variable flow rates and constraints on the number of simultaneous operations (such as barge unloads),can be significantly easier to model and solve in discrete-time (as seen in e.g. [6]). A discrete-timeformulation is also used in our work.A number of solution methodologies have emerged to deal with the handling of the nonlinear natureof these scheduling problems, particularly over long scheduling horizons. Such methods include single-step [63, 45, 50, 7, 3, 6] and iterated MILP-NLP [27, 26, 9, 14, 64] approaches, which have emerged assome of the most efficient methods to ensure small optimality gaps for relatively short scheduling horizons.Other methods include MILP-(smaller MINLP) iterations [50, 36] and iterative MILP-refinement methods[53, 53, 41, 34, 30].For longer scheduling horizons, rolling planning methods are often used to limit the number of integervariables. Such models typically yield good feasible solutions far faster than the full model, with far betterscaling of computation times with the scheduling horizon. In general, a rolling horizon model operateson a discrete time grid by representing in detail only the ‘present’ segment { t h , . . . , t h + H roll − } of thescheduling horizon { , . . . , H } at each iteration, then tanking a step forwards in time of size t step ≤ t h forthe next iteration, freezing the decisions made in the interval { t h , . . . , t h + t step − } and proceeding untilthe full horizon is reached.In some rolling planning models, such as those in [57, 13, 50, 59, 37] and heuristic H2 in [61], all aspectsof the ‘past’ decisions in { , . . . , t h − } are frozen; in others, such as [9] and heuristics H1, H3, and H4in [61], only certain variables are frozen (most commonly the binary/integer decisions). The treatment ofthe ‘future’ portion of the scheduling horizon { t + H roll , . . . , H } varies considerably between models: Somemodels exclude the ‘future’ from the current model entirely, as in [57, 13, 59, 39, 9], whereas some relaxcomplicating constraints [38, 61] and/or treat integer variables as continuous [37, 61]. Some works, suchas [13, 50, 10, 56], deal with uncertain events by updating uncertain data or events at each (or certain)rolling planning steps(s). For bi-level integrated planning and scheduling models, it is common to solve thescheduling subproblem only for the ‘present’ segment, while including (a subset of) the future segments inthe planning subproblem, as is done in [35, 66, 57, 56]. Some works use an alternate framework based onsupply arrivals to determine the time steps for each rolling iteration, as in [53, 54, 34, 9]; [34] in particularallows for some backtracking to correct spec quality discrepancies in the current plan.Finally, many ideas for MILP approximations to deal with bilinear terms have emerged in recent years,especially for the pooling problem. Some methods are either binary-free, or add a number of binary vari-ables that scales linearly with the desired level of precision. These include various linearization techniques[25, 44], piecewise-linear approximation [51, 16], outer-approximation methods [64], and McCormick orpiecewise McCormick envelopes [40, 1, 47, 5, 3, 14]. [17] in particular gives a numerical comparison for fif-teen different piecewise McCormick-related schemes, split into three classes of schema. For other methods,8 MATHEMATICAL MODEL often utilizing digit-based representations for a single variable, the number of binaries scales logarithmi-cally with the desired precision. Such methods include multi-parametric disaggregation [60, 28, 7, 8, 3, 4, 6]and other techniques [62, 43, 20, 30, 21]. Approximation via generalized disjunctive programming (GDP)is another, far more expressive, approach that has been used to tackle pooling-type problems. Castroand Grossmann et. al. have derived a MILP-based formulation for variants of multi-parametric disag-gregation as a disjunctive program [18, 28, 7, 3, 4], and have also applied it separately [55, 6, 36]. Inthis work, we choose to directly discretize the inventory specs of each tank using a formulation based onnormalized multi-parametric disaggregation (NMDT), as set forth in [4], yielding a relatively simple andhighly efficient model in conjunction with a rolling planning approach.We were unable to find any other work in the literature which combines a version of the multi-parametricdisaggregation technique with a version of rolling planning; this combination is a key novelty of thiswork. The primary results of this work include the speed of the method achieved while yielding near-optimal solutions without compromising the qualities of the specs in the demand feeds, combined withthe incorporation of constraints on the ratios of specs in the demand feeds.
We describe a nonconvex MIQCP model in this section that introduces variables to handle the constraints.Data in the problem is written with an underscore, continuous variables are written with lower-case letters,binary variables are written with Greek letters, and sets are written with capital letters.
Run Production run, consisting of a number of consecutive days for which the production feedmust be constant in terms of volume.Feed Flow from a storage tank to the downstream production processes.Mixing The mixing of material from a number of different sources to a single destination.Splitting The splitting of material from a single source to a number of different destinations.Spec Concentration of a certain type of chemical property.9 .2 Sets and Parameters 4 MATHEMATICAL MODEL
Sets S : Set of supply nodes (barges with raw material), indexed by s . S k : Set of supply barges allowed to unload to tank k ∈ K . S t : Set of supply barges available during time period t ∈ T . K : Set of tanks, indexed by k . K s : Set of tanks to which supply barge s ∈ S is allowed to unload. Q : Set of chemical properties, indexed by q . Q rat : Pairs of specification for which ratio bounds are to be enforced. R : Set of non-overlapping demand runs (consecutive days for which the production feed is constant). T : Set of time periods (days), such that T := { , , . . . , H } . T r : Time periods corresponding to run r where runs do not overlap, such that T r = { t ir , t ir + 1 , . . . , t fr } ⊆ T . T s : Time periods during which supply s ∈ S can be unloaded, such that T s = { t min s , t min s + 1 , . . . , t max s } . T d : Time periods where there is some demand, such that T d = (cid:70) r ∈ R T r .Parametersc inb s : Penalty cost per volume of not unloading supply from supply node s .c dmd t : Penalty cost per volume of not meeting production feed volume for time period t .d t : Demand to meet in time step t ∈ T . Must be constant throughout each run.Thus, d t = d t (cid:48) for all t, t (cid:48) ∈ T r for r ∈ R . Also, d t = 0 for t / ∈ T d .f s,q : Value of chemical characteristic q ∈ Q in supply source s ∈ S .f k,q, : Initial value of chemical characteristic q ∈ Q in tank k ∈ K .f min q,t , f max q,t : Bounds on the value of chemical characteristic q ∈ Q for the demand to be filledin time step t ∈ T . H : Scheduling time horizon (implicit in the definition of T ).r min q ,q ,t , r max q ,q ,t : Bounds on the ratio of chemical characteristics f q and f q ,( q , q ) ∈ Q rat for the demandto be filled in time step t ∈ T .tkp min k : Minimum percentage of feed each day that can originate fromtank k (if used).ul max : Maximum number of supply barges that can be unloaded each day.ulp min s : Minimum percentage of inventory from supply barge s .that can be unloaded each day that unloading occurs for s .ult max : Maximum permissible gap between the first and last unloadsfor each supply barge.v s : Volume of supply available for unloading from supply node s ∈ S .v k, : Initial inventory volume of tank k ∈ K .[v min k , v max k ] : Bounds on permissible inventory for each tank k ∈ K. Variables (All variables are non-negative) 10 .3 Constraints 4 MATHEMATICAL MODEL
Assignment: γ s,t : = (cid:40) s is unloaded in time step t s ∈ S, t ∈ Tσ k,t : = (cid:40) k provides material for production feed in time step t k ∈ K, t ∈ T Transfer: y in s,k,t : Volume of raw material to transfer from supply s ∈ S totank k ∈ K in time step t ∈ Ty out k,t : Volume of material to feed to demand from node k ∈ K in time step t ∈ T Resource: v end k,t : Volume of inventory in tank k ∈ K at the end of time step t ∈ Tv mid k,t : Volume of inventory in tank k ∈ K after material is supplied from supply barges,but before feed operations, in time step t ∈ T Demand: v prod t : Total volume fed to production in time step t ∈ T Tank specs: f k,q,t : Value of chemical characteristic q ∈ Q in tank k ∈ K at the end of time step t ∈ T Final specs: f prod q,t : Value of chemical characteristic q ∈ Q at production at the end of time step t ∈ T Unload time: t first s : Earliest unload date for supply barge st last s : Latest unload date for supply barge s Slack: v unused s : Volume of material that is not unloaded from supply barge s ∈ Sd unmet t : Volume of production demand missed in time step t ∈ T Objective: Our objective is to minimize the amount of material that is not unloaded from barges ( v unused s )and the demand that is unmet ( d unmet t ). To accomplish this, we maximize the sum of the value of theraw material unloaded from the barges and the value of the production feed demand that is met. Theamount of material unloaded for a barge s is exactly v s − v unused s , and the amount of demand met at time t is d t − d unmet t . Hence, we havemax (cid:88) s ∈ S c inb s · (v s − v unused s ) + (cid:88) t ∈ T c dmd t · (d t − d unmet t ) (1) The following constraints limit the solution space of the assignment, blending, and mixing problem. Inthe following discussion, note that variables y in , y out correspond to flows, while the variables v correspondto stationary volumes at certain times including volumes on barges v s , volumes in tanks v k,t , and volumesrequired for production feed v t . Initial Conditions : The volume of material in the tanks and the chemical characteristics of the materialare initialized. Specifically, (2a) enforces that the initial volumes in each storage tank are correctlyrepresented, and (2b) enforces that the initial chemical characteristics in each storage tank are specified. v end k, = v k, k ∈ K (2a) f k,q, = f k,q, k ∈ K, q ∈ Q (2b)11 .3 Constraints 4 MATHEMATICAL MODEL Flow : For these flow constraints, we assume that inbound supply material is unloaded and blended inthe tanks before any material is drawn from the demand tanks to meet production feed on any givenday. The model assumes that in the first half of a time period, barges unload to tanks (3a), adding tothe volume previously in the tank and creating an intermediate volume v mid k,t . In the second half of thetime period, material is unloaded from the tanks (3b) to support production (3c), while within requiredvolume bounds at production (3d) and (3e). (cid:88) s ∈ S k y in s,k,t + v end k,t − = v mid k,t k ∈ K, t ∈ T (3a) v mid k,t − y out k,t = v end k,t k ∈ K, t ∈ T (3b) (cid:88) k ∈ K y out k,t = v prod t t ∈ T (3c) v end k,t ≥ v min k k ∈ K, t ∈ T (3d) v mid k,t ≤ v max k k ∈ K, t ∈ T (3e) Chemical Characteristics : Bi-linear The following constraints capture the blending and mixing of chem-ical characteristics. Specifically, (4a) ensures the characteristics of the material in a tank at the end ofthe previous time period are combined with the material from supply barges to create a mix. Expression(4b) determines the characteristics of the material provided from tanks to meet production feed. f k,q,t · v mid k,t = (cid:88) s ∈ S k f s,q · y in s,k,t + f k,q,t − · v end k,t − k ∈ K, q ∈ Q, t ∈ T (4a) f prod q,t · v prod t = (cid:88) k ∈ K f k,q,t · y out k,t q ∈ Q, t ∈ T (4b)Thus, while supply unloads and production feeds can occur in the same time period (day), they are notconsidered to occur simultaneously. Demand : The following constraint set captures the production feed volume that is not met each timeperiod, d unmet t . v prod t = d t − d unmet t t ∈ T (5) Feed Characteristics : Both the chemical characteristics (6a) and ratios of characteristics (6b) need to bewithin the required bounds. Note that in practice for (6b) the denominator would be multiplied through,as numerical issues can arise when the denominator is small.f min q,t ≤ f prod q,t ≤ f max q,t t ∈ T, q ∈ Q (6a)r min q ,q ,t ≤ f prod q ,t f prod q ,t ≤ r max q ,q ,t t ∈ T, ( q , q ) ∈ Q rat (6b) Total Inflow : Ideally, the volume of supply on each barge is fully unloaded. Expression (7a) captures thevolume of supply that is not unloaded ( v unused s ). In addition, (7b) ensures no material is unloaded from12 .3 Constraints 4 MATHEMATICAL MODEL a supply barge when the barge is not available. (cid:88) k ∈ K s (cid:88) t ∈ T s y in s,k,t = v s − v unused s s ∈ S (7a) y in s,k,t = 0 s ∈ S, k ∈ K s , t ∈ T \ T s (7b) Continuous Feed During Runs : The volume of supply provided by each tank to meet production feedneeds to remain constant during a production run. y out k,t = y out k,t − r ∈ R, k ∈ K, t ∈ T r \{ t ir } (8) Tank Feed Share Constraints : If production feed is supplied by tank k during time period t (such that σ k,t = 1), then the volume of production feed ( y out k,t ) must be at least tkp min k of the demand d t at time t ,as enforced by (9a). If σ k,t = 0, then (9b) ensures no feed is provided by tank k for time t . y out k,t ≥ (tkp min k · d t ) · σ k,t k ∈ K, t ∈ T (9a) y out k,t ≤ d t · σ k,t k ∈ K, t ∈ T (9b) Supply Limitations : Although it is preferred that barges are completely unloaded in one day, due to tankstorage capacity, this may need to be done over several days. Hence, we constrain the number of unloadsto be bounded by bul max , which are 2 in all of our instances (10a). The number of times any barge isunloaded within a time period is also limited (10c). Expression (10b) serves as a tightening constraintto ensure that barges are not unloaded when they are not available. If a supply barge is not unloadedduring a time period (such that γ s,t is zero), then no volume is unloaded from the barge (10d). (cid:88) t ∈ T s γ s,t ≤ bul max s ∈ S (10a) γ s,t = 0 s ∈ S, t ∈ T \ T s (10b) (cid:88) s ∈ S γ s,t ≤ ul max t ∈ T (10c) y in s,k,t ≤ v s · γ s,t s ∈ S, t ∈ T, k ∈ K s (10d) Unload Inventory Percentage Constraints : (11) enforces that the percent of inventory unloaded per supplybarge remains above required minimum ulp min s each day, if unloading is to occur that day. (cid:88) k ∈ S k y in s,k,t ≥ (v s · ulp min s ) · γ s,t s ∈ S, t ∈ T (11) Unload Time Gap Constraints : In an effort to minimize the cost associated with barges remaining in thearea for an extended time, these constraints limit the time duration for a supply barge to remain in theunloading area. To accomplish this, constraints are included to limit the interval [ t first s , t last s ]. Expression(12a) ensures that t first s is no later than the first period that any inventory is unloaded for barge s .Expression (12b) enforces that t last s is no earlier than the last period that any inventory is unloaded for13 .4 Reformulation without spec variables f prod q,t barge s . Finally, (12c) enforces that the difference between the first and last unloading times for eachbarge is no larger than the allotted time gap. Although t first s may assume a value earlier than the actualearliest unload time and t last s may assume a value after the last unload time, if the interval is shorter thanthe required length, then the actual interval will also be shorter. Thus, this set of constraints is sufficient. t first s ≤ t · γ s,t + H (1 − γ s,t ) s ∈ S, t ∈ T (12a) t last s ≥ t · γ s,t − H (1 − γ s,t ) s ∈ S, t ∈ T (12b) t last s − t first s ≤ ult max s ∈ S (12c)All of the constraints described above are linear except for the chemical characteristic constraints describedby (4a),(4b), which contain bilinear terms. MILP approximations for these constraints are formulated inSection 5.2. f prod q,t For the implementation of this model, we adjust the model to eliminate terms consisting of only a single f prod variable from the model, restricting the model to product terms f · v or f prod · v , referred to asspec-volume terms. In this way, after discretization, we no longer have need to explicitly keep track ofthese specs at all, relying solely on these spec-volume product terms. In addition, we eliminate the f prod q,t and v prod t variables.The initial constraint (2b) is replaced with f k,q, · v end k, = f k,q, · v k, k ∈ K, q ∈ Q (2b-prod)By multiplying (3b) by f k,q,t we obtain f k,q,t · v mid k,t − f k,q,t · y out k,t = f k,q,t · v end k,t k ∈ K, t ∈ T (3b-prod)Note that the v mid k,t variables could also be eliminated via the constraint (3b). However, as this eliminationchoice can impact the quality of the MILP approximation, we defer this discussion to a later section.We reformulate the constraints (6a),(6b) by multiplying through by volume (and any denominators) toobtain the following:f min q,t · v prod t ≤ f prod q,t · v prod t ≤ f max q,t · v prod t t ∈ T, q ∈ Q (6a-prod)r min q ,q ,t · f prod q ,t · v prod t ≤ f prod q ,t · v prod t ≤ r max q ,q ,t · f prod q ,t · v prod t t ∈ T, ( q , q ) ∈ Q rat (6b-prod)Finally, we substitute the definitions of v t and f k,q,t · v t in (3c) and (4b) into constraints (5), (6a-prod),(6b-prod)14 .5 Model MINLP-SPLIT: Volume-based reformulation 4 MATHEMATICAL MODEL to obtain: (cid:88) k ∈ K y out k,t = d t − d unmet t t ∈ T (5-sub)f min q,t · (cid:88) k ∈ K y out k,t ≤ (cid:88) k ∈ K f k,q,t · y out k,t ≤ f max q,t · (cid:88) k ∈ K y out k,t t ∈ T, q ∈ Q (14a)r min q ,q ,t · (cid:88) k ∈ K f k,q ,t · y out k,t ≤ (cid:88) k ∈ K f k,q ,t · y out k,t ≤ r max q ,q ,t · (cid:88) k ∈ K f k,q ,t · y out k,t t ∈ T, ( q , q ) ∈ Q rat (14b)To summarize, the MINLP-Mixing reformulation is the initial model, but with (2b), (3b), (3c), (4b), (5),(6a), and (6b) replaced by (2b-prod), (3b-prod), (5-sub), (14a),(14b). The model presented above can be reformulated to keep track of volumes of chemical characteristics insteadof percentages f k,q,t , removing f k,q,t from the model. This reformulation results in a linear representationfor the mixing constraints in (14a),(14b),(4a) and modified flow balance constraint in (3b-prod). However,there is now a nonlinear representation for the splitting process, in which some volume is split into multipleparts (e.g. future inventory and demand feed). In the original model, the definition and usage of f k,q,t ensures that the outflow characteristics are correct. However, after reformulating without f k,q,t , this isnot longer the case, so we must ensure that the outflow characteristics are correct.To enact the reformulation, we make the following substitutions: f k,q,t · v end k,t = v f , end k,q,t k ∈ K, q ∈ Q, t ∈ T (15a) f k,q,t · v mid k,t = v f , mid k,q,t k ∈ K, q ∈ Q, t ∈ T (15b) f k,q,t · y out k,t = y f , out k,q,t k ∈ K, q ∈ Q, t ∈ T (15c)Thus, (3b-prod), (4a), (14a) and (14b) are written as v f , end k,q,t = v f , mid k,q,t − y f , out k,q,t k ∈ K, t ∈ T (16a) v f , mid k,q,t = (cid:88) s ∈ S k f s,q · y in s,k,t + v f , end k,q,t − k ∈ K, q ∈ Q, t ∈ T (4a-vol)f min q,t · (cid:88) k ∈ K y out k,t ≤ (cid:88) k ∈ K y f , out k,q,t ≤ f max q,t · (cid:88) k ∈ K y out k,t t ∈ T, q ∈ Q (16b)r min q ,q ,t · (cid:88) k ∈ K y f , out k,q ,t ≤ (cid:88) k ∈ K y f , out k,q ,t ≤ r max q ,q ,t · (cid:88) k ∈ K y f , out k,q ,t t ∈ T, ( q , q ) ∈ Q rat (16c)In addition, the initial condition (2b-prod) is replaced with v f , end k,q, = f k,q, · v k, k ∈ K, q ∈ Q (17)Now, to enforce the requirement that the chemical characteristics sent to production feed match matchthe chemical characteristics in the tank; that is: f k,q,t = v f , mid k,q,t v mid k,t = y f , out k,q,t y out k,t (18a)15 MILP APPROXIMATION METHODS we add the following nonlinear constraint v f , mid k,q,t · y out k,t = y f , out k,q,t · v mid k,t k ∈ K, q ∈ Q, t ∈ T (19a)To summarize, the MINLP-Splitting model is the initial model, but with (2b),(3b),(3c),(4a),(4b),(5),(6a),(6b)replaced by (17), (16a), (5-sub), (16b), (16c), and (19a). We consider approaches of changing the nonlinear portions of the model in Section 4.4 to an MILP tomore efficiently solve the problem and more easily interact with a rolling planning approach. The corechange is that we will replace all products y in · f and y out · f with approximations. The new models, givensufficiently accurate approximations, will produce feasible flows y in and y out that completely determinethe actions of the system.For a spec variable f ∈ [ l, u ] we describe it by a discete part ˙ f ∈ { l, l + (cid:15), l +2 (cid:15), . . . , l + ζ(cid:15) } and a continuouspart ∆ f ∈ [0 , (cid:15) ], where we choose (cid:15) such that u − (cid:15) < ζ(cid:15) ≤ u . The discrete part ˙ f will be modeled withbinary variables in a way known as Normalized Multi-Parametric Disaggregation Techniques (NMDT).For the products with ∆ f , we will consider two options - restricting the value to a midpoint or using theMcCormick relaxation. That is, the two main ideas considered are the models: x · ˜ f = x · ˙ f (cid:124)(cid:123)(cid:122)(cid:125) NMDT + x · ∆ f (cid:124) (cid:123)(cid:122) (cid:125) ∆ f ← ˆ (cid:15) (Center) x · ˜ f = x · ˙ f (cid:124)(cid:123)(cid:122)(cid:125) NMDT + x · ∆ f (cid:124) (cid:123)(cid:122) (cid:125) McCormick (McCormick)When using (Center), we relax (4a) to maintain flexibility on inflow volumes. We do this by replacingthe left-hand side with its upper and lower bounds with respect to ∆ f . These choices, along with variousrelaxations of redundant constraints are described in detail in this section.We will first describe common linearizations that we will consider, then describe our two competingmodels, and then have a discussion of rolling planning approaches. We define here two types of linearizations that we will consider. First, the McCormick envelope is a3-dimensional polytope M ( x, β ) that is the smallest convex set containing the equation x β = x · β withupper and lower bounds on β and x . We will assume 0 ≤ β ≤ min ≤ x ≤ x max . The McCormickenvelope is M ( x, β ) = (cid:40) ( x, β, x β ) ∈ [x min , x max ] × [0 , × R : x min · β ≤ x β ≤ x max · βx − x max · (1 − β ) ≤ x β ≤ x − x min · (1 − β ) (cid:41) . (20)16 .2 Discretization 5 MILP APPROXIMATION METHODS More generally, suppose we have β , . . . , β m ∈ { , } such that (cid:80) mj =0 β j = 1 hence (cid:80) mj =0 x · β j = x .Defining x β j = x · β j , we can then model the convex hull of all such solutions as M m ( x, β ) = ( x, β, x β ) ∈ [x min , x max ] × [0 , m +1 × R m +1 : m (cid:88) j =0 β j = 1 , m (cid:88) j =0 x β j = x x min · β j ≤ x β j ≤ x max · β j x − x max · (1 − β j ) ≤ x β i ≤ x − x min · (1 − β i ) . (21)In the special case where m = 1, the set M ( x, β ) is a bijective lifting of M ( x, β ). We discuss two methods of linearizing or relaxing bilinear products: Normalized Multiparametric De-scretization Technique (NMDT) [4] and McCormick envelopes. The discrete variable ˙ f , contained in theinterval [ l, u ], is converted into binary variables as˙ f = λ + ε n (cid:88) i =1 m (cid:88) j =0 ( j λ i ) · α ij , m (cid:88) j =0 α ij = 1 for all i = 1 , . . . , n, (22)where α ij ∈ { , } for all i, j and λ i are appropriately chosen scalars. Given a desired level of precision0 < ˆ ε ≤ u − l , the choices of λ i (and resulting choices for ε , m , and n ) will affect the number of binaryvariables α ij that are needed to reach the desired precision.1. Multiparametric Disaggregation Technique with base b provided l ≥ λ = 0 , λ i = b i − , m = b − → ε = b (cid:98) log b (ˆ ε ) (cid:99) , n = (cid:6) log b ( uε ) (cid:7) . (MDT)2. Normalized Multiparametric Disaggregation Technique, with base b : λ = l, λ i = b i − , m = b − n = (cid:6) log b (cid:0) u − l ˆ ε (cid:1)(cid:7) , ε = ( u − l ) b − n . (NMDT)3. Monolithic Disaggregation: λ = l, n = 1 , λ = 1; m = (cid:6) u − l ˆ ε (cid:7) , ε = ( u − l ) /m. (Mono)(MDT) discretizes ˙ f independent of the lower bound l , while (NMDT) and (Mono) use the lower boundand can be interpreted as discretizing ˙ f after normalizing the bounds [ l, u ] to [0 , ε . Furthermore, the exponential approach of (NMDT) is vastlypreferred over the unary approach of (Mono), allowing only ( O (log(1 /ε )) as opposed to O (1 /ε )).In the originating work for (NMDT) [4], the same number of digits n (and hence, precision ε ) is appliedto all discretized variables (to maintain the same precision in the normalized space), while in the version17 .2 Discretization 5 MILP APPROXIMATION METHODS presented here, n is computed based on the desired precision ˆ ε for each individual discretized variable (tomaintain similar precision in the original space). D ( x, f, α ) = ( x, f, ˙ f , ∆ f,α, x f , x α ) : f = ˙ f + ∆ f ˙ f = λ + ε n (cid:88) i =1 m (cid:88) j =1 ( j λ i ) · α ijm (cid:88) j =0 α ij = 1 i = 1 , . . . , n x f = x · f = λ x + x · ∆ f + ε n (cid:88) i =1 m (cid:88) j =0 ( j λ i ) · x α ij x α ij = x · α ij i = 1 , . . . , n, j = 0 , . . . , m x min ≤ x ≤ x max α ij ∈ { , } i = 1 , . . . , n, j = 0 , . . . , m . (23)Most of the nonlinear products are now linearized and we drop the variables ˙ f and f . In this way, f isnow only recorded implicitly (approximately) through product terms x α ij , x f , and x ∆ f . This produces theset˜ D ( x, f, α ) = ( x, α, x f , x ∆ f , x α ) : x f = λ · x + x ∆ f + ε n (cid:88) i =1 m (cid:88) j =0 ( j λ i ) · x α ij ( x, α i , x α i ) ∈ M m ( x, α i ) α ij ∈ { , } i = 1 , . . . , n, j = 0 , . . . , m , (24)where λ , n , and ε are computed from a specified ˆ ε , x min and x max .The remaining nonlinear product x · ∆ f , approximated by x ∆ f , will be dealt with in two different ways.The binarization in [21] is similar to (NMDT) with b = 2 after projecting out the α i variables using theequation (cid:80) j =0 α ij = 1, i.e., α i + α i = 1. b = 2 reduces variables. We suggest the best model uses b = 2 since it uses asymptotically fewer binary variables than any otherchoice b > ε → + ). Proposition 1.
Let z ˆ ε ( b ) := ( b − · n = ( b − (cid:24) log b (cid:18) u − l ˆ ε (cid:19)(cid:25) be the number of binary variables usedto represent ˙ f , and let b , b be two different bases. Then z (cid:48) ( b , b ) := lim ˆ ε → + z ˆ ε ( b ) z ˆ ε ( b ) = ( b −
1) ln( b )( b −
1) ln( b ) . .2 Discretization 5 MILP APPROXIMATION METHODS Proof.
Asymptotically, as ˆ ε → + , we can remove the ceiling function in the definition of z ˆ ε . Hence z (cid:48) ( b , b ) = lim ˆ ε → + ( b −
1) log b (cid:18) u − l ˆ ε (cid:19) ( b −
1) log b (cid:18) u − l ˆ ε (cid:19) = ( b − b −
1) log b ( b ) = ( b −
1) ln( b )( b −
1) ln( b ) . Because z ˆ ε ( b ) is an increasing function of b , z (cid:48) ( b, >
1, for all b >
2. Using the proposition, we compute z (cid:48) ( b,
2) for various values of b in 3. In [4], Castro uses b = 10. As the table shows, this binarization willuse nearly three times as many binary variables. b z (cid:48) ( b,
2) 1.26186 1.5 1.72271 1.93426 2.13724 2.33333 2.52372 2.70927Table 3: Asymptotic number of binary variables required for base b , as a multiple of the number requiredfor base 2.In summary, we suggest to use NMDT with a base 2 because it will require the fewest number of binaryvariables to achieve a desired accuracy. Constraint (3b) (and subsequently (3b-prod)) couples the variables v mid k,t , y out k,t , and v end k,t . For convenience,let us write x = v mid k,t , x = y out k,t , x = v end k,t , and f = f k,q,t . That is, x = x + x . (25)As a consequence of (25), by multiplying through by f, ∆ f , and α i , we obtain three sets of valid equations x f = x f + x f , (26a) x ∆ f = x ∆ f + x ∆ f , (26b) x α i = x α i + x α i i = 1 , . . . , n. (26c)Consider the approximation of f and the resulting product form of the equations f = l + ∆ f + ε n (cid:88) i =1 i − · α i (27a) x f κ = l · x κ + x ∆ f κ + ε n (cid:88) i =1 i − · x α i κ κ = 1 , , x α i κ = x κ · α i i = 1 , . . . , n, κ = 1 , , .3 MILP Models 5 MILP APPROXIMATION METHODS Of key importance here for model reduction and a stronger formulation are the observations: (1) Bilinearterms can be eliminated using the above equations (2) Applying McCormick relaxations to these bilinearterms before variable elimination produces a tighter formulation.First, although other model choices are possible, we project out the variables associated with κ = 1.Conveniently, equations (27b) and (27c), for κ = 1, are both implied by all other equations. In fact,variables x f , x ∆ f and x α i by substitution, can be projected out of the model using equations (25), (26b),and (26c), then the equations (27b) and (27c) become tautologies. This observation provides convenientreduction choices in the model in terms of both variables and constraints.Second, applying the projection/linear relaxation in the right order can improve the model. In particular,let F represent the feasible set given by the above equations, let M ( x, α κ ) be the McCormick relaxationin the space of all variables (similarly, M ( x, α κ ) in the space of variables associated with κ = 2 , be the projecting into the space of variables with κ = 2 ,
3, and denote the removal of bilinear equations(27c).Then we have the strict inclusionrelax proj F ∩ (cid:92) κ =1 , , M ( x, α κ ) (cid:40) relax proj ( F ) ∩ (cid:92) κ =2 , M ( x, α κ ) . (28)That is, using the inequalities from the McCormick relaxation from the κ = 1 variables improves theformulation. We will now outline the two models we propose to consider and then discuss how to integrate them withrolling planning.
We create a mixed integer linear approximation to the main problem that will only approximately trackconstraints on specifications. To check feasibility of the resulting flows y in and y out , these solutions must beplugged into the main model to recover spec values and ensure the constraints are met. The approximationerror will be guided by how fine the discrete approximation of the spec variables is. For the model, aprevision ˆ ε q is chosen such that, given an accurate discretization, we should have | f k,q,t − ˙ f k,q,t | ≤ ˆ ε . Thechoice of ˆ ε q , as discussed earlier, defines the choices of (cid:15) q , n q and m q which are implied parameters inthe sets D ( x, f k,q,t , α ). Thus, for each spec q , there are (potentially) unique (cid:15), n and m that are used forthe set D ( x, f k,q,t , α ). We write all products in terms of single variables and we relax (4a) to increase thenumber of feasible solutions to the approximate model.20 .3 MILP Models 5 MILP APPROXIMATION METHODS We thus model (3b-prod),(4a),(14a),(14b) as v f , mid k,q,t − y f , out k,q,t = v f , end k,q,t k ∈ K, t ∈ T, (29a) v f , mid k,q,t − ε q · v mid k,q,t ≤ (cid:88) s ∈ S k f s,q · y in s,k,t + v f , end k,q,t − k ∈ K, q ∈ Q, t ∈ T, (29b) v f , mid k,q,t + ε q · v mid k,q,t ≥ (cid:88) s ∈ S k f s,q · y in s,k,t + v f , end k,q,t − k ∈ K, q ∈ Q, t ∈ T, (29c)f min q,t · (cid:88) k ∈ K y out k,t ≤ (cid:88) k ∈ K y f , out k,q,t ≤ f max q,t · (cid:88) k ∈ K y out k,t t ∈ T, q ∈ Q, (29d)r min q ,q ,t · (cid:88) k ∈ K y f , out k,q ,t ≤ (cid:88) k ∈ K y f , out k,q ,t ≤ r max q ,q ,t · (cid:88) k ∈ K y f , out k,q ,t t ∈ T, ( q , q ) ∈ Q rat , (29e)The remaining constraints are generated using ˜ D with λ = f min k,q and ∆ f k,q,t = ε q and b = 2,( v mid k,q,t , α k,q,t , v f , mid k,q,t , v ∆ f , mid k,q,t , v α, mid k,q,t ) ∈ ˜ D ( v mid k,q,t , f k,q,t , α k,q,t ) k ∈ K, q ∈ Q, t ∈ T, (30a)( v end k,q,t , α k,q,t , v f , end k,q,t , v ∆ f , end k,q,t , v α, end k,q,t ) ∈ ˜ D ( v end k,q,t , f k,q,t , α k,q,t ) k ∈ K, q ∈ Q, t ∈ T, (30b)( y out k,q,t , α k,q,t , y f , out k,q,t , y ∆ f , out k,q,t , y α, out k,q,t ) ∈ ˜ D ( y out k,q,t , f k,q,t , α k,q,t ) k ∈ K, q ∈ Q, t ∈ T. (30c)Note that each α k,q,t is a vector of n k,q,t binary variables, which may vary in number depending on theaccuracy used to approximate f k,q,t . Defining both sets of constraints in (30a), (30b) creates some redundant inequalities. In particular, afteromitting the k, q, t indices for convenience, the inequalities v α, mid i,j ≥ v min k · α i,j i ∈ I k,q,t , j ∈ { , } (31a) v α, end i,j ≤ v max k · α i,j i ∈ I k,q,t , j ∈ { , } (31b)are redundant due to (3b) and the non-negativity of all volume variables. Because the proposed model is an approximation of the original problem, the model might produce aschedule with flows that result in infeasible specs - violating either the spec bounds or spec ratio bounds.However, in most cases, the magnitude of these violations we observed to be within ˆ ε/ ε/ min t,q , f max t,q ]is replaced with [f min t,q + ˆ ε q / , f max t,q − ˆ ε q / .3 MILP Models 5 MILP APPROXIMATION METHODS As for the spec ratios, we used differentials of the ratio to obtain an upper approximation for the magnitudeof the buffer to use: If the spec ratio to control is f q f q , where f q and f q are spec qualities to control, thebuffer to use is given via the differential:∆ (cid:18) f q f q (cid:19) = 1 f min2 ∆ f q + f max1 (cid:0) f min2 (cid:1) ∆ f q , (32)where ∆ f q and ∆ f q are the computed buffers ˆ ε/ f q and f q , respectively. Thus, the interval[r min q ,q ,t , r max q ,q ,t ] enforced for (cid:16) f q f q (cid:17) is replaced by [r min q ,q ,t +∆ (cid:16) f q f q (cid:17) , r max q ,q ,t − ∆ (cid:16) f q f q (cid:17) ]. These buffers provedeffective in largely resolving the issue of incorrect feed specs and spec ratios for the problems tested. Thebuffers could be reduced if they would result in possible ranges for specs or spec ratios that are smallerthan the range provided by a single discretization bin (based on initial tank and supply inventories andthe parameter ˆ ε q specified by the user). However, for the purposes of this work, we will assume that ˆ ε q issmall enough that the resulting ranges for the demand specs are sufficiently large (with width ≥ ˆ ε aftertightening) that feasibility issues do not become problematic for the proposed discretization scheme.These buffers result in a replacement of the parameters f max t,q , f min t,q , r max t,q , and r min t,q in Section 4 with their‘buffered’ counterparts, tightened further by the possible feed specs based on initial tank inventory andinbound supply specs. In a similar way, the final model for the original NMDT method (with McCormick envelopes around thecontinuous product terms) is given as follows: v f , mid k,q,t − y f , out k,q,t = v f , end k,q,t k ∈ K, t ∈ T (33a) v f , mid k,q,t = (cid:88) s ∈ S k f s,q · y in s,k,t + v f , end k,q,t − k ∈ K, q ∈ Q, t ∈ T ((4a)-relax)f min q,t · (cid:88) k ∈ K y out k,t ≤ (cid:88) k ∈ K y f , out k,q,t ≤ f max q,t · (cid:88) k ∈ K y out k,t q ∈ Q, t ∈ T (33b) (cid:88) k ∈ K y f , out k,q ,t · r min q ,q ,t ≤ (cid:88) k ∈ K y f , out k,q ,t ≤ (cid:88) k ∈ K y f , out k,q ,t · r max q ,q ,t t ∈ T, ( q , q ) ∈ Q rat (33c)Where, for all k ∈ K, q ∈ Q, t ∈ T ,( v mid k,q,t , ∆ f k,q,t , α k,q,t , v f , mid k,q,t , v α, mid k,q,t ) ∈ ˜ D ( v mid k,q,t , f k,q,t , α k,q,t ) (34a)( v end k,q,t , ∆ f k,q,t , α k,q,t , v f , end k,q,t , v α, end k,q,t ) ∈ ˜ D ( v end k,q,t , f k,q,t , α k,q,t ) (34b)( y out k,q,t , ∆ f k,q,t , α k,q,t , y f , out k,q,t , y α, out k,q,t ) ∈ ˜ D ( y out k,q,t , f k,q,t , α k,q,t ) (34c)∆ v f , mid k,q,t ∈ [0 , ε · v max k,q,t ] (34d)∆ v f , end k,q,t ∈ [0 , ε · v max k,q,t ] (34e)∆ y f , out k,q,t ∈ [0 , ε · d k,q,t ] (34f)22 .4 Rolling Horizon Approach 5 MILP APPROXIMATION METHODS Finally, the continuous product terms ∆ v f , mid k,q,t , ∆ v f , end k,q,t , and ∆ y f , out k,q,t are modeled via their McCormickenvelopes, for all k ∈ K, q ∈ Q, t ∈ T , ∆ f k,q,t ∈ [0 , ε ] (35a)(∆ f k,q,t /ε, v mid k,q,t , ∆ v f , mid k,q,t ) ∈ M (∆ f k,q,t , v mid k,q,t ) (35b)(∆ f k,q,t /ε, v end k,q,t , ∆ v f , end k,q,t ) ∈ M (∆ f k,q,t , v end k,q,t ) (35c)(∆ f k,q,t /ε, y out k,q,t , ∆ y f , out k,q,t ) ∈ M (∆ f k,q,t , y out k,q,t ) (35d) From a perspective of numerical error propagation, the McCormick-based approach results in an exactly-represented tank-inflow blending process (constraint (4a)) for our problem: since only two ‘spec · volume’terms contain variable specs, the volumes of the specs are coupled correctly. On the other hand, thesplitting processes in (6a) are inexact, since the McCormick envelopes around the ∆ f product termseffectively allow the specs leaving a tank to differ slightly from the specs propagated to the next timeperiod. However, due to the tightness of the McCormick envelopes, once a post-blending value ˜ f is chosenfor a tank, the represented specs will remain in the interval [ ˜ f , ˜ f + ε ] until the next blending operationin the tank, allowing a maximum possible spec error in the interval [ ε/ , ε ] to accrue between blendingoperations.On the other hand, the Center approximation presented here incurs an error of up to ε/ As mentioned in Section 1, we need a blending and scheduling optimizer that can efficiently optimizeover a scheduling horizon of at least 3 months and up to 12 months. For the instances of interest, thetarget time for this optimization is typically preferred to be no more than five minutes (with a ten minutemaximum). However, even the linear discrete approximation above is extremely slow even for a schedulinghorizon of 120 days, requiring several hours of computation time in Gurobi 8.1, even for extremely coursespec discretization (with a total of roughly six binary variables per tank per day for only three tanks).As such, we employ a rolling planning approach.We present three major alternatives for rolling planning. These alternatives are characterized treatment ofthe ‘past’, ‘present’, and ‘future’, coupled with the choice of time periods (‘present’) to use for each rollingplanning step. The ‘past’ is the portion of the planning horizon before the current rolling horizon window,the ‘present’ is the current horizon window, and the ‘future’ is the portion of the planning horizon afterthe current rolling horizon window. For the considered schemes, the ‘future’ is split into two partitions:23 .4 Rolling Horizon Approach 5 MILP APPROXIMATION METHODS the ‘near future’ and the ‘far future’. We consider two alternatives for the treatment of the ‘past’ and‘future’, and two major schemes for the choice of time periods.We will compare two schedules for rolling planning - fixed-length periods and run-based periods. Inboth schemes, periods are chosen to be disjoint since previous performed test with overlapping periodsperformed slower.Our fixed-length periods are simply equal length periods of some length ∆ t , except the last period maybe truncated to not exceed the horizon length H . See Figure 3. Figure 3: This figure shows specified demand runs, the barge time windows, and the periods that we usein the rolling horizon calculation for the fixed-length periods scheme (and considering a horizon H = 30days).Although this period schedule is consistent, it subdivides runs. This makes modeling consistent volumefeed throughout a run (8) complicated in the rolling planning process and can lead to poor choices earlyin the rolling planning. To alleviate this issue, we propose a ’run-based’ schedule.The run-based periods either • Continue in a run to the end of the run or, • Encompass all periods and gaps (days between runs) that are completely contained within a timechange of ∆ t from the beginning of period.This choice makes it so that periods and gaps are never subdivided. This is depicted in Figure 4, wherea sample partition resulting from this procedure, using ∆ t = 7, is illustrated. Formally, our choice ofperiods is described by Algorithm 1. 24 .4 Rolling Horizon Approach 5 MILP APPROXIMATION METHODS Algorithm 1
Rolling Planning with Run Based Sub-periods.
Input:
Desired approximate length ∆ t for a single rolling planning sub-period Output:
Sub-periods P i for rolling planning. t ← i ← (cid:46) Initialize beginning of current rolling planning step and index. while t i < H do t end ← first ending time period of a run that is ≥ t i t change ← latest ending time period of a run or gap between runs that is ≤ t i + ∆ tt i +1 ← max( t end , t change ) + 1 P i ← [ t i , t i +1 ) i ← i + 1 Figure 4: This figure shows specified demand runs, the barge time windows, and the periods that we usein the rolling horizon calculation for the ‘run-based periods’ scheme (and considering H = 30 days).Next, we will discuss the major schemes for the treatment of the past, present, near future, and futuresegments of the scheduling horizon. These schemes include a full-horizon scheme and a partial-horizonscheme. Full-Horizon Scheme
Given the current rolling horizon period P i = [ t i , t i + 1) and a near futurehorizon H NF , let t NF i = min( H, t i + H NF − H NF = 90.1. In the past, we freeze the binary decisions made, but allow all continuous variables to continue tovary. This allows the model some additional ability to compensate for poor early decisions, which(as we will see) helps to improve the optimality gap substantially.2. In the present, the full model is used.3. In the near future, the binary variables related to demand decisions and discretized tank compo-sitions are relaxed to the interval [0 , , .4 Rolling Horizon Approach 5 MILP APPROXIMATION METHODS The decision not to relax supply decision binaries in the ‘near future’ is due to the inclusion of the ‘limitedunloading window’ constraints for barges. These constraints rendered the inclusion of small objectivefunction terms based on supply unloading decisions computationally prohibitive, due to poor branchingperformance of the solver. Before these were included, small objective function terms had been addedto slightly favor earlier unloading times (all else equal), encouraging barges to be unloaded more quickly(and allowing for a model that more consistently avoided small unnecessary supply misses). After theinclusion of these constraints, the solver performance with these small objective terms became extremelypoor, requiring at least an order of magnitude more computational time to solve to the same optimalitygap in preliminary tests via Gurobi 8.1. When these terms were removed, the model began to exhibit verylarge extraneous supply misses, due in part to the poor linear relaxation of the ‘limited unloading window’constraints, and in part to relaxations of the constraints limiting the number of unloads per barge. Thesepoor relaxations caused the rolling planning model to favor unloading during the relaxed-binary periods, adecision which the model could not recover from. However, by including the supply-based binaries in the‘near future’ (defined to end 90 days after the start of the ‘present’ period), preliminary testing showedthat these unnecessary misses became small enough that their cost was considered to be acceptable, whilethe inclusion of these extra binary variables had only a moderate impact on computational cost (close toa 50% increase).For all rolling horizon methods, we use 0.5% as the MIP gap. For the Full-Horizon scheme, thoughexperimental results we learned that increasing this, on average, does not significantly reduce computationtime, but may have a large effect on resulting objective values. In particular, we want to optimize as muchas we can in each step since, due to the way we relax variables, the objective value in each subsequentstep can only decrease. In particular, since we are targeting the best objective value possible, we attemptto make optimal decisions at every step.
Partial-Horizon Scheme
1. In the past, we fix all decisions made, and use them to re-compute accurate values for spec com-positions to use for the next rolling planning period. Possible [ l, u ] for the specs in each tank arere-computed, and the number of remaining unloads (and remaining volumes) for barges that havebeen unloaded in the past are reduced accordingly.2. In the present, the full model is enforced.3. In the near future, the binary variables related to demand decisions and discretized tank compo-sitions are relaxed to the interval [0 , .5 Post-processing and Error Mitigation 5 MILP APPROXIMATION METHODS Full Horizon SchemeVars Past Present Near Far[0 , t i ) P i (cid:2) t i +1 , t NF i (cid:3) (cid:0) t NF i , H (cid:3) y, v active active active active y , v active active active active d active active active active t s active active active active γ fixed active active relaxed σ fixed active relaxed relaxed α fixed active relaxed relaxed Partial Horizon SchemeVars Past Present Near Far[0 , t i ) P i (cid:2) t i +1 , t NF i (cid:3) (cid:0) t NF i , H (cid:3) y, v fixed active active omitted y , v fixed active active omitted d fixed active active omitted t s fixed active active omitted γ fixed active active omitted σ fixed active relaxed omitted α fixed active relaxed omittedTable 4: Description of how variables are either fixed, fully active, relaxed to be continuous, or omittedin the two versions of the rolling planning horizon schemes.In the partial-horizon scheme, we note that it is possible that only a small part the arrival window, evenjust a single period, for a single barge overlaps with the ‘present’ and ‘near future’. To compensate forthis, we pro-rate the inventory to be unloaded from such barges: If only 10% of the volume on a bargeoverlaps with the currently-considered period, then we enforce that only 10% of the volume from thebarge is to be unloaded in the current rolling planning step. In this way, as the model rolls forwards,more of the barge’s unloading window enters the visible horizon window, allowing for improved decisionsregarding the volume on the barge. After solving the approximation model, we check feasibility of the scheduled flows in terms of requiredspec quality bounds. To do so, using the flow variables y in , y out , we plug their computed values into theoriginal model from Section 4 to uniquely determine the spec values f k,q,t .After the optimization of the discretized model with rolling planning, we compute specs f k,q,t from simulatethe scheduling solution to reveal the true specs throughout the scheduling horizon, and report the truetank and feed volumes and compositions throughout the horizon in addition to the basic scheduling plan.During preliminary computational testing, after simulating the scheduling solution to reveal the truespecs throughout the scheduling horizon, we found that, without tightening the spec requirements, theMILP discretization techniques for the model often resulted in solutions for which the actual spec andspec ratios for the demand feeds were not within the specified ranges. However, this was largely remediedafter tightening the feasible region as described in Section 5.3.327 NUMERICAL RESULTS
Tank 1 Tank 2 Tank 3Total Volume 1179 metric tons 2948 metric tons 1225 metric tonsMinimum Volume 158 metric tons 272 metric tons 136 metric tonsTable 5: The capacities of the three tanks we consider. Minimums are about 10% of total volumes. Notethat barges hold between 1200 - 1400 metric tons. Hence a barge is comparable to our small tanks, andour large tank can hold slightly more than 2 barges.Barge Type Type 1 Type 2 Type 3 Type 4Volume 1240 1360 1182 1360c inb
To evaluate the performance of the proposed methods, we constructed a realistic problem with a 368-day scheduling horizon from real operational data with a 119-day scheduling horizon, extended to thedesired horizon in a periodic fashion. The resulting problem had a total of 33 supply arrivals over thescheduling horizon, three available storage tanks, two specs ( S and S ) to track, and one spec ratio tocontrol. Note that the range of possible spec qualities was relatively narrow for some tanks, resulting inrelatively few binary variables used to approximate the specs for each tank. See (5) for example data onthe storage tanks. Other relevant parameters include ul max = 2, ult max = 7, tkp min k = 10% for each tank,and ulp min = 10% for each supply barge. For the objective coefficients, we globally use c dmd t = 3000.c inb t depends on the barge type, and is either $800 or $1000 for each type, as specified in Table 6. Figure6 showcases the unloading and arrival windows for a sample instance with H = 150, while Figures 5,7,showcase sample solution schedules with H = 40. Note that the high variability in bounds for the S /S To perform our computational tests, we obtained three years of demand data for the same system, andgenerated randomized supply data for those three years. We generated ten different such randomizations,and performed each computational test with five different starting dates (separated by half a year),obtaining a total of fifty test problems for each scheduling horizon.28 .1 Comparison with McCormick 6 NUMERICAL RESULTS
Figure 5: Volume profiles with some bounds describing allowable ranges of spec and spec ratios.29 .1 Comparison with McCormick 6 NUMERICAL RESULTS F i g u r e : E x a m p l e o f t i m e w i nd o w s f o r b a r g e a rr i v a l s a ndd a t e s o f r un s o n a150 d a y t i m e h o r i z o n . F i g u r e : E x a m p l e o f s o l u t i o n f o r s c h e du l e s o f l oa d i n g t o t a n k s a ndun l oa d i n g ( f ee d ) f r o m t a n k s t o p r o du c t i o n f o r a40 d a y t i m e h o r i z o n . .1 Comparison with McCormick 6 NUMERICAL RESULTS In this section, we compare the effectiveness of two different methods for handling the product terms v · ∆ f :McCormick and the proposed Center method, each using a base-2 binarization of the specs. The tests wererepeated using several different values of H , and with several values for ˆ ε (chosen to be the same for allspecs). For ˆ ε = 1, we used H ∈ { , , , , , } , and for ˆ ε = 0 .
25, we used H ∈ { , , , , } (since McCormick became computationally prohibitive over longer scheduling horizons). A total timelimit of 600 seconds was enforced for all computations in this section. The average-case results are shownin Figure 8. In this section, we investigate performance in terms of computation time and percentage ofvalue lost (%loss) in the objective, with respect to solution values v unused s and d unmet t , computed via val target = (cid:88) s ∈ S c inb s · v s + (cid:88) t ∈ T c dmd t · d t (36a) val missed = (cid:88) s ∈ S c inb s · v unused s + (cid:88) t ∈ T c dmd t · d unmet t (36b)%loss = val missed val target · ε = 1 yielded bythe Center method is consistently quite a bit worse that that yielded by McCormick. As showcased inFigure 9(c), the majority of this difference seems to be caused by significantly worse solutions in roughlya quarter of the test cases. However, this disadvantage disappears at higher fidelity; with ˆ ε = 0 .
25, theoptimality results are very close (well within Gurobi’s default optimality gap of 10 − ), with the Centermethod taking a slight edge over McCormick. At the same time, the Center method proved to yieldsimilar feasibility results as McCormick, with each becoming infeasible for no more than one out of the50 trials in all test cases.Figure 9 gives another perspective on the performance advantage of the proposed Center method, basedon the number of instances solved by a given time for the longest scheduling horizons tested: for ˆ ε = 1and H = 90, the advantage is overwhelming: almost all of the 50 test instances have converged for theCenter method before any have converged for McCormick. The results for ε = 0 . H = 60 also showa strong advantage for the proposed method, though this advantage is somewhat less pronounced due tothe shorter scheduling horizon.These results suggest that, for longer scheduling horizons, the new Center idea fares far better computa-tionally than does the McCormick-style discretization in NMDT for the same level of precision.31 .1 Comparison with McCormick 6 NUMERICAL RESULTS (a) (b)(c) (d)Figure 8: Comparison of computational performance of the Center approximation vs. the McCormickrelaxation. Averaged over 50 similar instances with different randomized supply parameters. (a) and (c):Computational time using step sizes of ˆ ε = 1 and 0 .
25. (b) and (d): % Loss computed as fraction ofupper bound on objective value not obtained (see (36c)), where model accuracy is set to ˆ ε = 1 and 0 . .1 Comparison with McCormick 6 NUMERICAL RESULTS (a) (b)(c)Figure 9: Time-completion (a,b) and % value missed (c) profiles for the comparison of computationalperformance for McCormick vs. the proposed approximate method, using the largest values of H testedfor (a,c) ˆ ε = 1, and (b) ˆ ε = 0 . ε = 1, we find little difference between the performance of the two methods, asshowcased in figure Figure 10. As a significant performance difference between the methods had alreadymanifested by such small scheduling horizons in the cases with rolling planning, we find that, for theseproblem instances, the primary advantage with the Approx method lies in its interaction with rollingplanning. 33 .2 Comparison to Nonconvex MIQCP in Gurobi 9.0 6 NUMERICAL RESULTS (a) (b)Figure 10: Time-completion profiles for the comparison of computational performance for NMDT vs. theproposed approximate method, using ε = 1 tested for (a) H = 30 and (b) H = 40. In this section, we compare the computational performance of the methods presented, using ε = 1 andwithout rolling planning, to solving the models MINLP-Mix and MINLP-Split defined in Section ?? withGurobi 9.0, using scheduling horizons of H = 20, H = 30, and H = 40. The MINLP-Mix model timedout (5-minute limit) for H = 30 on all instances but one, so we omit it from our comparisons beyondthe first plot with H = 20. The results are shown in Figure 11. Note that the MINLP-Mix results areomitted from (b), as many instances also didn’t finish for H = 20.The performance of MINLP-Mix is far inferior to that of MINLP-Split: For H = 20, ∼
44 instancesfinished for MINLP-Split within 50s, while only 10 had finished by that point for MINLP-Mix, with only19 finishing within 5 minutes. For H = 30, no more than a single instance finish within 5 minutes forMINLP-Mix.It is interesting to note that, in quite a few instances, Gurobi is able to solve the MINLP-Split modelvery quickly ( ∼
26 for H = 20, ∼
21 for H = 30, and ∼
15 for H = 40. The results show that moreinstances get stuck with long solution times than with the NMDT-based options, with the slowest 40-60%of solutions slower than the NMDT-based methods. However, the solutions it returns are precise in termsof spec feasibility, while the level of solution precision ε = 1 for the discrete options is quite coarse. Itseems likely that, using to ε = 0 . ε = 0 .
25, the Gurobi-based solver may overtake the NMDT-basedvariants in both solution time and quality for these methods. As a result, one might imagine that arolling planning scheme based on this MINLP-split model, e.g. relaxing the bilinear constraint (19a) forthe ‘future’ portion of the horizon, could be quite fruitful. We leave this exploration as future work.The %loss results show that, even at such course fidelity, the approximation scheme without rollingplanning consistently yields near-optimal results. Indeed, the relative difference between the solutions isconsistently within the MIP-gap, 0.5%, used for the optimization.34 .2 Comparison to Nonconvex MIQCP in Gurobi 9.0 6 NUMERICAL RESULTS
Computation Time % Loss(a) (b)(c) (d)(e) (f)Figure 11: Time-completion (a,c,e) and percent-loss (b,d,f) profiles for the comparison of computationalperformance for Gurobi 9.0.1 vs. the proposed NMDT-based methods, using ε = 1 tested for (a,b) H = 20, (c,d) H = 30, and (e,f) H = 40. 35 .3 Performance Enhancement 6 NUMERICAL RESULTS We explore how to enhance the performance of the method via two simple changes to the Center modelthat preserve the mixed-integer feasible solutions. The first of these is to add the coupling constraintsdescribed in Section 5.2.2 v α, mid i,j = v α, end i,j + y α, out i,j i ∈ I k,q,t , j ∈ { , } (37a)where I k,q,t is defined as in Section 5.3.1. We then remove constraints (31a) and (31b), since they arenow made redundant in the associated LP by (37a) combined with the nonnegativity of v α, mid p,i, , v α, end p,i, ,and y α, out p,i, .The second of these is to relax the MILP approximation of the v α, mid p,i,j terms. This corresponds thenremoving constraints (31a), (31b), and the constraints v α, end i,j ≥ v min k · α i,j i ∈ I k,q,t , j ∈ { , } , (38a) y α, out i,j ≤ d t · α i,j i ∈ I k,q,t , j ∈ { , } , (38b)which are part of the representation of (30a) and (30c). The viability of the removal of (38b) is shown inSection 5.2.2: since any one of the three associated variables v α, mid p,i, , v α, end p,i, , and y α, out p,i, can be eliminatedcompletely without compromising the MILP, a valid model is still obtained when relaxing any number ofconstraints related to the definition of y α, out p,i, . Constraint (38a) is always safe to remove, as it is impliedby the fact that only one v α, mid i,j variable can nonzero for each i , combined with the fact that v α, mid i,j sum to v mid k,t , which is bounded below by v min k . Note that either (31b) or (38b) must be included if the couplingconstraints (37a) are not included. In this case, we choose to remove (38b).The results, displayed in Figure 12, show that while both modeling changes make a large difference incomputational cost, the simple relaxation of the MILP yields a larger improvement in computationaltime. This improvement is increased further when the coupling constraints are added, yielding a methodin which the majority of instances finish within the desired 10 minutes when using a scheduling horizonof a full year with ˆ ε = 1. 36 .4 Comparison of Rolling Planning Ideas 6 NUMERICAL RESULTS Figure 12: Time-completion profiles using different combinations of relaxing (constraints) and adding(coupling constraint), with H = 365 and a 7-1-1 fixed-step-size rolling planning scheme. In this section, we compare the computational performance of different rolling planning regimes. Werestrict ourselves here to the full-horizon scheme, and compare the performance obtained using twodifferent ideas. For the first, we use fixed-length periods with H int = 7. For the second, we use run-based periods using ˜ H int = 4 to generate the periods, with n H = ∆ n = 2 as our stepping parameters.We compare the methods using a time limit of 1800 s using ˆ ε = 1 and scheduling horizons of H ∈{ , , , , , } (in days).The results, as shown in Figure 13, indicate that the regime with a fixed step size yields a more consistent,and often lower, optimality gap on average while incurring a higher computational cost. This outcomeis emphasized in Figure 14: when H = 365, the 7-1-1 fixed-step scheme results in nearly all test casesfinishing within about a 5-10 minute time window, while incurring high optimality gaps in a handfulof instances. On the other hand, for the 4-2-2 run-dependent-step scheme, most instances finish withinabout a 10-20 minute time window, but optimality gap for the worst instances is far more controlled.37 .4 Comparison of Rolling Planning Ideas 6 NUMERICAL RESULTS (a) (b)Figure 13: Comparison of fixed-step-length and run-dependent-step rolling planning methods. Averagedover 50 similar instances with different randomized supply parameters, using spec step size ˆ ε = 1. (a):Computational cost, (b) Optimality gap.(a) (b)Figure 14: Time-completion profiles of fixed-step-length and run-dependent-step rolling planning methods,using H = 365 days and ˆ ε = 1. Run using both the relaxation and coupling things. (a) Computationalcost, (b) Optimality gap. 38 EFERENCES
In this work, we have developed an approximation method for the tank mixing and scheduling problemthat combines rolling horizon planning with two different discretization schemes, yielding reasonably high-quality solutions very quickly over long scheduling horizons, while with high probability ensuring thatdemand specs remain within the required ranges. Due to this performance, the model could be especiallyuseful for obtaining rough, fast plan feasibility results during the planning of supply acquisitions andproduct production runs in an industry setting (counting solutions with no (or only small) supply anddemand inventory misses as ‘feasible’).Comparing the in-house “Center” discretization scheme to the original “McCormick”-based NMDT scheme[4], we find that the in-house method yields much faster performance when combined with a rolling hori-zon scheme, while sometimes yielding worse solutions. However, without a rolling horizon scheme, themethodologies yield comparable performance, suggesting that the “Center” scheme may yield relaxationsmore compatible with a rolling horizon scheme. Further, when comparing rolling horizon schemes, wefind that a scheme using run-based rolling horizon periods yields better solutions than a fixed-durationscheme, at some cost to performance. To extend this work, we plan to explicitly incorporate acquisitionplanning of barges, and possibly planning of production runs.
References [1] Maria Lorena Bergamini, Pio Aguirre, and Ignacio Grossmann. Logic-based outer approximation forglobally optimal synthesis of process networks.
Computers & Chemical Engineering , 29(9):1914–1933,aug 2005.[2] Pedro Castillo Castillo, Pedro M. Castro, and Vladimir Mahalec. Global optimization algorithmfor large-scale refinery planning models with bilinear terms.
Industrial & Engineering ChemistryResearch , 56(2):530–548, January 2017.[3] Pedro M. Castro. New MINLP formulation for the multiperiod pooling problem.
AIChE Journal ,61(11):3728–3738, sep 2015.[4] Pedro M. Castro. Normalized multiparametric disaggregation: an efficient relaxation for mixed-integer bilinear problems.
Journal of Global Optimization , 64(4):765–784, jul 2015.[5] Pedro M. Castro. Tightening piecewise McCormick relaxations for bilinear problems.
Computers &Chemical Engineering , 72:300–311, jan 2015.[6] Pedro M. Castro. Source-based discrete and continuous-time formulations for the crude oil poolingproblem.
Computers & Chemical Engineering , 93:382–401, oct 2016.39
EFERENCES REFERENCES [7] Pedro M. Castro and Ignacio E. Grossmann. Global optimal scheduling of crude oil blending op-erations with RTN continuous-time and multiparametric disaggregation.
Industrial & EngineeringChemistry Research , 53(39):15127–15145, sep 2014.[8] Pedro M. Castro and Ignacio E. Grossmann. Optimality-based bound contraction with multipara-metric disaggregation for the global optimization of mixed-integer bilinear problems.
Journal ofGlobal Optimization , 59(2-3):277–306, mar 2014.[9] Jaime Cerd´a, Pedro C. Pautasso, and Diego C. Cafaro. Efficient approach for scheduling crude oiloperations in marine-access refineries.
Industrial & Engineering Chemistry Research , 54(33):8219–8238, aug 2015.[10] Billy R. Champion and Steven A. Gabriel. A multistage stochastic energy model with endogenousprobabilities and a rolling horizon.
Energy and Buildings , 135:338–349, jan 2017.[11] Xuan Chen, Simin Huang, Dejing Chen, Zhihai Zhang, Li Zheng, Ignacio Grossmann, and ShihuiChen. Hierarchical decomposition approach for crude oil scheduling: A SINOPEC case.
Interfaces ,44(3):269–285, jun 2014.[12] Yunfei Chu, Fengqi You, John M. Wassick, and Anshul Agarwal. Integrated planning and schedulingunder production uncertainties: Bi-level model formulation and hybrid solution method.
Computers& Chemical Engineering , 72:255–272, January 2015.[13] Cao Cuiwen, Gu Xingsheng, and Xin Zhong. A data-driven rolling-horizon online scheduling modelfor diesel production of a real-world refinery.
AIChE Journal , 59(4):1160–1174, sep 2012.[14] Leonardo Salsano de Assis, Eduardo Camponogara, Bernardo Zimberg, Enrique Ferreira, and Igna-cio E. Grossmann. A piecewise McCormick relaxation-based strategy for scheduling operations in acrude oil terminal.
Computers & Chemical Engineering , 106:309–321, nov 2017.[15] Lisia S. Dias and Marianthi G. Ierapetritou. Integration of scheduling and control under uncertainties:Review and challenges.
Chemical Engineering Research and Design , 116:98–113, December 2016.[16] Xiaoyong Gao, Yongheng Jiang, Tao Chen, and Dexian Huang. Optimizing scheduling of refineryoperations based on piecewise linear models.
Computers & Chemical Engineering , 75:105–119, apr2015.[17] Chrysanthos E. Gounaris, Ruth Misener, and Christodoulos A. Floudas. Computational comparisonof piecewise-linear relaxations for pooling problems.
Industrial & Engineering Chemistry Research ,48(12):5742–5766, June 2009.[18] Ignacio E. Grossmann and Juan P. Ruiz. Generalized disjunctive programming: A framework forformulation and alternative algorithms for MINLP optimization. In
Mixed Integer Nonlinear Pro-gramming , pages 93–115. Springer New York, November 2011.[19] Vijay Gupta and Ignacio E. Grossmann. Multistage stochastic programming approach for offshoreoilfield infrastructure planning under production sharing agreements and endogenous uncertainties.
Journal of Petroleum Science and Engineering , 124:180–197, December 2014.40
EFERENCES REFERENCES [20] Akshay Gupte, Shabbir Ahmed, Myun Seok Cheon, and Santanu Dey. Solving mixed integer bilinearproblems using MILP formulations.
SIAM Journal on Optimization , 23(2):721–744, January 2013.[21] Akshay Gupte, Shabbir Ahmed, Santanu S. Dey, and Myun Seok Cheon. Relaxations and discretiza-tions for the pooling problem.
Journal of Global Optimization , 67(3):631–669, apr 2016.[22] Miguel Angel Guti´errez-Lim´on, Antonio Flores-Tlacuahuac, and Ignacio E. Grossmann. MINLPformulation for simultaneous planning, scheduling, and control of short-period single-unit processingsystems.
Industrial & Engineering Chemistry Research , 53(38):14679–14694, September 2014.[23] Miguel Angel Guti´errez-Lim´on, Antonio Flores-Tlacuahuac, and Ignacio E. Grossmann. A reactiveoptimization strategy for the simultaneous planning, scheduling and control of short-period contin-uous reactors.
Computers & Chemical Engineering , 84:507–515, January 2016.[24] Michael Wagner J. Rohde, Herbert Meyr. Die supply chain planning matrix.
PPS Management ,5(1), jan 2000.[25] M. Joly and J.M. Pinto. Mixed-integer programming techniques for the scheduling of fuel oil andasphalt production.
Chemical Engineering Research and Design , 81(4):427–447, apr 2003.[26] Ramkumar Karuppiah, Kevin C. Furman, and Ignacio E. Grossmann. Global optimization forscheduling refinery crude oil operations.
Computers & Chemical Engineering , 32(11):2745–2766,nov 2008.[27] J. D. Kelly and J. L. Mann. Crude oil blend scheduling optimization: an application with multimilliondollar benefits - part 1.
Hydrocarbon Processing, International Edition , 82(6):47–53, jun 2003.[28] Scott Kolodziej, Pedro M. Castro, and Ignacio E. Grossmann. Global optimization of bilinear pro-grams with a multiparametric disaggregation technique.
Journal of Global Optimization , 57(4):1039–1063, jan 2013.[29] Scott P. Kolodziej, Ignacio E. Grossmann, Kevin C. Furman, and Nicolas W. Sawaya. A novelglobal optimization approach to the multiperiod blending problem. In
Computer Aided ChemicalEngineering , pages 1492–1496. Elsevier, 2012.[30] Scott P. Kolodziej, Ignacio E. Grossmann, Kevin C. Furman, and Nicolas W. Sawaya. Adiscretization-based approach for the optimization of the multiperiod blend scheduling problem.
Computers & Chemical Engineering , 53:122–142, jun 2013.[31] Heeman Lee, Jose M. Pinto, Ignacio E. Grossmann, and Sunwon Park. Mixed-integer linear program-ming model for refinery short-term scheduling of crude oil unloading with inventory management.
Industrial & Engineering Chemistry Research , 35(5):1630–1641, January 1996.[32] Jie Li and I. A. Karimi. Scheduling gasoline blending operations from recipe determination to shippingusing unit slots.
Industrial & Engineering Chemistry Research , 50(15):9156–9174, August 2011.[33] Jie Li, I. A. Karimi, and Rajagopalan Srinivasan. Recipe determination and scheduling of gasolineblending operations.
AIChE Journal , pages NA–NA, 2009.41
EFERENCES REFERENCES [34] Jie Li, Wenkai Li, I. A. Karimi, and Rajagopalan Srinivasan. Improving the robustness and efficiencyof crude scheduling algorithms.
AIChE Journal , 53(10):2659–2680, 2007.[35] Zukui Li and Marianthi G. Ierapetritou. Rolling horizon based planning and scheduling integrationwith production capacity consideration.
Chemical Engineering Science , 65(22):5887–5900, nov 2010.[36] Irene Lotero, Francisco Trespalacios, Ignacio E. Grossmann, Dimitri J. Papageorgiou, and Myun-Seok Cheon. An MILP-MINLP decomposition method for the global optimization of a source basedmodel of the multiperiod blending problem.
Computers & Chemical Engineering , 87:13–35, apr 2016.[37] Shan Lu, Hongye Su, Yue Wang, Lei Xie, and Quanling Zhang. Multi-product multi-stage productionplanning with lead time on a rolling horizon basis.
IFAC-PapersOnLine , 48(8):1162–1167, 2015.[38] Chunpeng Luo and Gang Rong. A strategy for the integration of production planning and schedulingin refineries under uncertainty.
Chinese Journal of Chemical Engineering , 17(1):113–127, feb 2009.[39] Julien F. Marquant, Ralph Evins, and Jan Carmeliet. Reducing computation time with a rollinghorizon approach applied to a MILP formulation of multiple urban energy hub system.
ProcediaComputer Science , 51:2137–2146, 2015.[40] Garth P. McCormick. Computability of global solutions to factorable nonconvex programs: Part i— convex underestimating problems.
Mathematical Programming , 10(1):147–175, dec 1976.[41] Carlos A. M´endez, Ignacio E. Grossmann, Iiro Harjunkoski, and Pousga Kabor´e. A simultaneousoptimization approach for off-line blending and scheduling of oil-refinery operations.
Computers &Chemical Engineering , 30(4):614–634, feb 2006.[42] Ruth Misener and Christodoulos A Floudas. Advances for the pooling problem: Modeling, globaloptimization, and computational studies.
Applied and Computational Mathematics , 8(1):3–22, 2009.[43] Ruth Misener, Jeffrey P. Thompson, and Christodoulos A. Floudas. APOGEE: Global optimizationof standard, generalized, and extended pooling problems via linear and logarithmic partitioningschemes.
Computers & Chemical Engineering , 35(5):876–892, May 2011.[44] Lincoln F. L. Moro and Jos´e M. Pinto. Mixed-integer programming approach for short-term crudeoil scheduling.
Industrial & Engineering Chemistry Research , 43(1):85–94, jan 2004.[45] Sylvain Mouret.
Optimal Scheduling of Refinery Crude-Oil Operations . PhD thesis, Carnegie MellonUniversity, mar 2010.[46] Sylvain Mouret, Ignacio E. Grossmann, and Pierre Pestiaux. A novel priority-slot based continuous-time formulation for crude-oil scheduling problems.
Industrial & Engineering Chemistry Research ,48(18):8515–8528, sep 2009.[47] Sylvain Mouret, Ignacio E. Grossmann, and Pierre Pestiaux. Tightening the linear relaxation ofa mixed integer nonlinear program using constraint programming. In
Integration of AI and ORTechniques in Constraint Programming for Combinatorial Optimization Problems , pages 208–222.Springer Berlin Heidelberg, 2009. 42
EFERENCES REFERENCES [48] Sergio Neiro and Jose Pinto. Lagrangean decomposition applied to multiperiod planning of petroleumrefineries under uncertainty.
Latin American Applied Research , 36:213–220, 10 2006.[49] Chao Ning and Fengqi You. A data-driven multistage adaptive robust optimization framework forplanning and scheduling under uncertainty.
AIChE Journal , 63(10):4343–4369, May 2017.[50] Thordis Anna Oddsdottir, Martin Grunow, and Renzo Akkerman. Procurement planning in oilrefining industries considering blending operations.
Computers & Chemical Engineering , 58:1–13,nov 2013.[51] Viet Pham, Carl Laird, and Mahmoud El-Halwagi. Convex hull discretization approach to the globaloptimization of pooling problems.
Industrial & Engineering Chemistry Research , 48(4):1973–1979,February 2009.[52] J.M. Pinto and L.F. L. Moro. A planning model for petroleum refineries.
Brazilian Journal ofChemical Engineering , 17(4-7):575–586, December 2000.[53] P. Chandra Prakash Reddy, I. A. Karimi, and R. Srinivasan. Novel solution approach for optimizingcrude oil operations.
AIChE Journal , 50(6):1177–1197, 2004.[54] P.Chandra Prakash Reddy, I.A. Karimi, and R. Srinivasan. A new continuous-time formulation forscheduling crude oil operations.
Chemical Engineering Science , 59(6):1325–1341, mar 2004.[55] Juan P. Ruiz and Ignacio E. Grossmann. Strengthening of lower bounds in the global optimiza-tion of bilinear and concave generalized disjunctive programs.
Computers & Chemical Engineering ,34(6):914–930, June 2010.[56] Javier Silvente, Georgios M. Kopanos, Vivek Dua, and Lazaros G. Papageorgiou. A rolling horizonapproach for optimal management of microgrids under stochastic uncertainty.
Chemical EngineeringResearch and Design , 131:293–317, mar 2018.[57] Javier Silvente, Georgios M. Kopanos, Efstratios N. Pistikopoulos, and Antonio Espu˜na. A rollinghorizon optimization framework for the simultaneous energy supply and demand planning in micro-grids.
Applied Energy , 155:485–501, oct 2015.[58] A. Singh, J.F. Forbes, P.J. Vermeer, and S.S. Woo. Model-based real-time optimization of automotivegasoline blending operations.
Journal of Process Control , 10(1):43–58, February 2000.[59] Raik Stolletz and Emilio Zamorano de Acha. A rolling planning horizon heuristic for schedulingagents with different qualifications.
SSRN Electronic Journal , 2014.[60] Jo˜ao P. Teles, Pedro M. Castro, and Henrique A. Matos. Multi-parametric disaggregation tech-nique for global optimization of polynomial programming problems.
Journal of Global Optimization ,55(2):227–251, nov 2011.[61] S. Torkaman, S.M.T Fatemi Ghomi, and B. Karimi. Multi-stage multi-product multi-period produc-tion planning with sequence-dependent setups in closed-loop supply chain.
Computers & IndustrialEngineering , 113:602–613, nov 2017. 43
EFERENCES REFERENCES [62] Juan Pablo Vielma and George L. Nemhauser. Modeling disjunctive constraints with a logarithmicnumber of binary variables and constraints.
Mathematical Programming , 128(1-2):49–72, July 2009.[63] Li Wenkai, Chi-Wai Hui, Ben Hua, and Zhongxuan Tong. Scheduling crude oil unloading, storage,and processing.
Industrial & Engineering Chemistry Research , 41(26):6723–6734, dec 2002.[64] Jialin Xu, Shujing Zhang, Jian Zhang, Sujing Wang, and Qiang Xu. Simultaneous scheduling offront-end crude transfer and refinery processing.
Computers & Chemical Engineering , 96:212–236,jan 2017.[65] Yale Zhang, Dayadeep Monder, and J. Fraser Forbes. Real-time optimization under parametricuncertainty: a probability constrained approach.
Journal of Process Control , 12(3):373–389, April2002.[66] Edwin Zondervan, Michiel Kaland, Martijn A.H. van Elzakker, Jan C. Fransoo, and Jan Meuldijk.Integrating planning and scheduling in an oil refinery with a rolling horizon approach. In