Multiple Optima in Gasoline Blend Planning

Article pubs.acs.org/IECR Multiple Optima in Gasoline Blend Planning Shefali Kulkarni-Thaker and Vladimir Mahalec* School of Computational Engineering and Science, McMaster University, Hamilton, Ontario, Canada L8S 4L7 ABSTRACT: Gasoline is produced by blending several different components in ratios such that the blended mixture meets the required quality specifications. The blender produces different batches of gasoline by switching operation from one grade of gasoline to another. Blend planning horizon usually spans 10 to 14 days. Blend plan optimization minimizes the total blend costs by solving a multiperiod problem, where demands need to be satisfied in each period and some inventory is carried into the future time periods to meet the demands. Since blend component production is determined by a longer range refinery production plan, inventory carrying costs are not included in the objective function. It is shown that nonlinear programming (NLP) as well as mixed integer nonlinear programming (MINLP) solvers lead to different blend recipes and different blend volume patterns for the same total cost. The new algorithm described in this work systematically searches for multiple optimum solutions; this opens the way for blend planners to select from different blend plans based on additional considerations (e.g., blend more of regular gasoline earlier in the planning horizon thereby creating an opportunity to meet more demand for it in early periods) instead of having to use only one solution that varies with the choice of the solver. Inherent structure of the proposed algorithm makes it well suited for implementation on parallel CPU machines. The total amount of gasoline to be blended over the planning horizon is determined from the supply commitments to the customers. Since shipments of different grades are uneven from day to day, the blender operates in multibatch mode: for example, produce a batch of midgrade gasoline, then a batch of premium gasoline, etc. Switching from one gasoline grade to another requires calibration of analytical instruments; that is, there is a setup time that reduces the total blend-capacity. From the blender, the produced gasoline is sent to the product storage tanks. Product shipments (in batches) are from the storage tanks to pipelines or trucks or ships. Simplistic approach to gasoline blending would be to add together all expected product shipments (for each grade separately) and compute how much additional gasoline needs to be produced for each grade and also compute the corresponding optimal blend recipes for each grade. However, such an approach most of the time does not produce feasible solutions, since product shipments vary over the planning horizon and there are limits on the storage of blend components and the storage of the gasoline products. Hence, gasoline blending must consider when specific products need to be delivered along the blending time horizon, production rates, and properties of blend components, as well as the blender capacity and the inventory constraints. Work process currently prevalent in the industry to optimize gasoline blending consists of the following: 1. INTRODUCTION Gasoline blending produces several different grades of gasoline by blending various intermediate refinery streams. Costs of intermediate streams (blend components) depend on the cost of operating the process units that produce them. Each grade of gasoline has different specifications that vary with season and the geographical location of the target market. Gasoline quality specifications are inequalities, that is, a specific property either has to be greater than or equal to the specification (e.g., octane number greater than or equal 91) or less than or equal to the specification (e.g., Reid vapor pressure less than 11). Since components for blending gasoline are complex mixtures of many different chemical compounds, it is possible to meet the same product specifications by using different blend recipes while using the same components. These different blend recipes will result in gasoline batches (mixtures of components) that are closer or further away from the constraints. An optimal percentage of components in a blend (optimal blend recipe) is computed in such a manner that the product quality specifications are satisfied, while minimizing the cost of the blend. Gasoline costs are largely driven by the cost of high octane components. Hence, if a batch of gasoline is blended even slightly above the minimum required octane, there is a significant loss of profit. An average refinery can lose millions of dollars per year1 if it blends consistently 0.1 octanes above the minimum required octane (e.g., 91.1 vs 91.0). Such economic importance of gasoline blending has made it a subject of research for the last several decades. Simplified gasoline blending system is shown in Figure 1. The components for gasoline blending are stored in separate storage tanks; that is, there is a separate tank for each component. Materials from the component tanks are blended in ratios that correspond to the specified blend recipe. Blending is carried out in a blend header, a device where gasoline components are mixed. In order to simplify the model, and without losing accuracy of the results, one can think of blending occurring in the tank that contains the respective grade of gasoline. © 2013 American Chemical Society • Blend planningminimize total cost of all blends over a multiperiod horizon. ◦ Blend planning time horizon is divided into periods (discrete time representation). Received: Revised: Accepted: Published: 10707 August 29, 2012 June 20, 2013 June 21, 2013 June 21, 2013 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Figure 1. Simplified gasoline blending system. consisted of 1104 continuous variables, 150 binary variables, and 640 nonlinear terms; it was solved in 5274 CPU seconds with an optimality gap of 0.5%. When blend properties are computed by (quality*volume) or (quality*mass) constraints, and the blend component properties are known, the resulting blending model for a single period is linear. However, the corresponding multiperiod model is nonlinear, since the finished product volumes and the qualities are not known from the second period onward. Possible simplification is to assume that the properties of the starting inventory of each grade have properties that are exactly on the quality constraints. This means that the subsequent blends do not have to adjust the qualities of the already blended gasoline in its storage tank (no “heel compensation”). Under these conditions, one can formulate a multiperiod LP or MILP model (Mendez et al.2). Linear constraints to describe blend properties are commonly used (e.g., Li and Karimi5 or Mendez et al.2). If product tank heel properties need to be considered in blend planning, even (quality*quantity) constraints the model is nonlinear, since it contains the bilinear terms that are unknown in second and subsequent time periods. Current industrial practice employs by MINLP or NLP algorithms6 to solve such problems. As illustrated in paragraph 3, using such approach, one can arrive at different optimal solutions, depending on the algorithm is used. This is clearly caused by the existence of the multiple optima. In order to eliminate multiple optima, one can add to the objective function terms that include some other criteria, for example, minimize deviations of the blend recipes from one period to another or from some favorite recipe (Mendez et al.2). However, such formulation means that the objective function is a combination of the costs and the penalties for deviations from the average blend recipes; that is, the blend recipe corresponding to the solution might not be the same as the real economic optimum of the unconstrained recipe problem. An alternative way to avoid multiple solutions is to include the cost of carrying product inventory in the objective function. Since the refinery production plan determines the amounts of blend components to be produced and the amounts gasoline to be blended, the inclusion of inventory carrying costs is not justified, because the refinery has to carry either the blend components or the finished product in the quantities specified by the production ◦ Period boundaries are set as fixed points on the planning horizon (either by the calendar or by the beginning of product liftings). ◦ Within each period, one or more grades are blended. ◦ An algorithm computes how much of each grade to produce in each period and the blend recipe for each batch. The model is solved via NLP (if no constraints on batch sizes are considered) or MINLP. Alternatively, if product properties (or their transformations, so-called blending indices) can be computed via (quality*volume) constraints, then an MILP formulation (Mendez et al.2) can be employed. • Blend schedulingsequence the batches computed in the blend plan so that there is a minimal number of switches from blending one grade of gasoline to another. In addition, minimize switching of swing tanks from storing one grade of gasoline to another. It is solved via MILP or heuristic algorithms. Blend planning model typically deals with a planning horizon of 7 to 14 days. Since the product inventory is typically larger than the amount shipped in a single time period, some inventory is often carried from one period to the next, as needed to meet future demands. Glisman and Gruhn3 proposed a two level algorithm corresponding to the two tasks (planning task, scheduling task) in gasoline blending. At the top level, they employed multiperiod NLP to compute the blend plan. Period boundaries in at the top level are set based on calendar time unit (e.g., days). The lower level of the algorithm uses fixed recipes from the top level in a MILP with very short periods (2 h) to produce a schedule. Adjustments to the blend recipes are carried out if the scheduling problem is infeasible. Environmental Protection Agency model for reformulated gasoline introduces bilinear and polynomial terms, as well as the terms raised to a fractional power. Misener et al.4 presented a MINLP model that integrates these constraints into the pooling problem. In the model, binary variables are used to represent the quality breakpoints defined by the EPA model. To solve this highly nonlinear mixed-integer model to global optimality, they introduced a specialized algorithm. One of their case studies 10708 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Figure 2. Gasoline blending model structure. This work focuses on computation of multiple optimal solutions to multiperiod blend planning models, in order to provide the blend planner with options to select among the available resulting inventory profiles of blend components and finished products. The remainder of this paper is organized as follows: Section 2 describes a gasoline blending model which contains bilinear terms and has multiple minimal cost solutions. Section 3 demonstrates that solutions for gasoline blend planning obtained by several different NLP and MINLP solvers are significantly different from each other. This shows that the current practice of adhering to the blend plan computed by an NLP or MINLP solver ties the hands of the blend planner to a rather arbitrary set of blend recipes and inventory profiles. A two-level optimization algorithm is presented in Section 4. At the upper level, we compute the volume of each blend for each period via Differential Evolution (DE). Once the volume of each batch is determined, quality constraints become linear and the multiperiod planning problem can be solved as an LP. Section 5 presents the results and compares them to solutions via NLP and MINLP solvers. Evolutionary search via DE computes a number of optimal solutions. Even though we cannot prove that DE finds the global optimum, it is likely that these solutions are globally optimal, since the same value of the objective function is determined after large number of generations and starting from different initial populations. plan. While addition of penalty terms or carrying costs to the objective function enables us to avoid having to deal with the existence of multiple optima, it does not allow the decision maker (blend planner) to select those optimal solutions that may be preferred due to some additional operational considerations. Existence of multiple solutions of multiperiod gasoline of blend planning models can be intuitively deduced by the following reasoning. Let us assume that the production rates of blend components, blender capacity, and the product liftings are such that one can decide, for instance, to (a) blend all regular grade gasoline during the first 3 days, then do nothing for 2 days, then blend premium gasoline for three days, then blend midgrade gasoline for three days, (b) blend midgrade gasoline for three days, then blend premium gasoline for three days, then do nothing for 2 days, then blend regular gasoline for 3 days, (c) etc. Clearly, (a) and (b) have the same number of switches and would compute the same blend recipes for respective grades. The total blended volume for each grade would be equal to the demand for that grade of gasoline. In the presence of inventory constraints, some of these choices (how much to blend and when) may not be feasible, since the inventory constrains limit the range of choices. Nevertheless, there still exist multiple solutions; there is just a smaller set of them. This is illustrated by the fact that, for a given blend planning problem, different solvers arrive at different blend plans with the same value of the objective function. An alternative to discrete time period approach to gasoline blending is to employ a continuous time representation: individual blends are represented by slots that are placed at such points along the time horizon as required to produce optimal operation. Li and Karimi5 presented such MILP formulation with a schedule adjustment algorithm for the recipe determination and scheduling of gasoline blending operations. Their results show that some realistic-scale problems (2 to 3 blenders, 4 to 6 products, 9 components, 9 properties, 11 product tanks, 10 to 45 orders, and a planning horizon of 8 days) can be solved, but others were unsolvable to optimality within the allocated CPU times (10800 s to 108000 s, depending of the problem). 2. GASOLINE BLEND PLANNING Figure 2 presents representation of gasoline blending in the blend planning model. Each grade of gasoline is produced by a separate instance of the blender. Only one of the grades can be blended at any time, since, in reality, there is only one blender switching from one grade to another. Changeover times between product runs on the blender are product and sequence independent. It is assumed that the demand and supply of every gasoline grade and raw material are known. The opening inventories in all the component and blending tanks are known. It is also assumed that the blend components have consistent composition during each time period. In order to simplify the notation, we will assume that the component quality is constant 10709 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article products is equal to the orders. Each order involves only one product, and each order is completed during the planning horizon. The product inventory tank may receive a new blend while gasoline is being withdrawn from the tank. A component inventory tank may receive additional amounts of the component while the material is being withdrawn from the tank and sent to the blender. We will assume in the example presented that the component qualities are constant for the whole planning horizon. Since blend components are produced by upstream units that are operated to meet target qualities of the blend components, this assumption is valid so long as the process units produce blend components to their target qualities. If this is not the case, then the point in time when the quality changes is taken as a start of a new period; that is, component quality is constant during any time period. Quality constraints are either linear or can be made linear via blend index transformation (Mendez et al.2). It is assumed that a perfect mixing occurs in the blender. If the initial inventory (heel) of some product tank is offspec, it will be used as blend component to produce on-spec gasoline. Blend planning for a gasoline blending systems as shown in Figure 1 is described by 1. a planning horizon divided into fixed-duration time periods 1,2,...,N. 2. a set of required shipments for each product and for each period along the planning horizon (demand profile). If there are multiple orders for the same grade of gasoline in a given time period, these orders are lumped into a total demand for that grade of gasoline in that period. 3. a set of blend components, their properties, initial inventories, costs, and flow rates for every period along the horizon. 4. a set of products (i.e., gasoline grades) with prescribed minimum and maximum quality specifications, their initial inventories and corresponding initial quality. 5. minimum and maximum inventories (for every time period) for each blend component and for each grade of gasoline. If there are multiple tanks available for storage of the same material, they are treated as one aggregate inventory capacity for that material. 6. the blender capacity loss due to switching the blender from one grade to another. 7. if a specific grade of gasoline is to be blended, then the amount blended must be greater than or equal to some threshold amount. We need to compute: 1. How much of each grade of gasoline to make and how much of each blend component to use for each gasoline blend (blend recipes) for each period Table 1. Explanation of Symbols symbol G g I i K k T t f k,g,t Ck Cr dt Fmax Fmax NP Qopen,i,g,t Qclose,i,g,t Qi,k,t Qi,g,min Qi,g,max Qi,g,t Sg,t ug,t Vblend,g,t Vclose,g,t Vin,k,t Vopen,g,t Vopen,k,t Vclose,k,t Vout,g,t Vout,k,g,t Vmin Vmax description set of gasoline grades index in set G set of qualities index in set I set of component tanks index in set K set of time periods index in set T fraction of component k in blend of gasoline grade g, period t unit cost of component k, [USD] crossover probability threshold in differential evolution length of the time period [time] mutation scaling factor in differential evolution maximum blender flow rate [volume/time] number of population members in differential evolution opening quality, i, for gasoline grade g in period t [quality units] closing quality, i, for gasoline grade g in period t [quality units] quality i of component k in period t [quality units] minimum and maximum specifications on the quality i gasoline g [quality units] calculated quality, i, for gasoline grade g in period t [quality units] binary switching variable (=1 if grade g in period t is blended) trial vector generated by DE. Used to calculate the blend volume Vblend,g,t blend volume of gasoline grade g in period t [volume] closing inventory for gasoline grade g in period t [volume] supply for component k in period t [volume] opening volume in time t for gasoline grade g [volume] opening inventory of component k in period t [volume] closing inventory of component k in period t [volume] demand for gasoline grade g in period t [volume] volume contributed by component k to blend gasoline g, period t [volume] minimum and maximum tank capacities [volume] throughout the blend planning horizon. If that was not the case, the only change in the equation presented below would be the values of the qualities in each period. Problem Statement. The refinery production plan determines the amounts of gasoline blend components and their qualities, which are required to meet the total gasoline demand throughout the time horizon covered by the production plan. The blend planning horizon is typically 1 to 2 weeks long with periods being one day or half a day in duration. Product orders (demands) need to be met; that is, lifting (shipments) of the Table 2. Component Properties, Initial Inventories, and Costs component properties, initial inventories, and costs property ALK BUT HCL HCN LCN LNAP REF aromatics [%] benzene [%] MON [octane] RON [octane] olefins [%] RVI [psi] SPG sulfur [%] inventory [BBl] unit cost [USD/BBL] 0.00 0.00 93.7 95.0 0.00 5.15 0.703 0.0 500 $29.2 0.00 0.00 90.0 93.8 0.00 138.00 0.584 0.0 400 $11.5 0.00 0.00 79.8 82.3 0.00 22.33 0.695 0.0 160 $20.0 25.00 0.50 75.8 86.7 14.00 2.38 0.791 0.490 350 $22.0 18.00 1.00 81.6 93.2 27.00 13.88 0.744 0.080 700 $25.0 2.97 0.59 81.6 67.8 0.00 19.90 0.677 0.010 200 $19.7 74.9 7.50 90.8 103.0 0.00 3.62 0.818 0.000 500 $24.5 10710 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Table 3. Product Property Specifications, Costs, Initial Inventories, and Their Properties product inventories, qualities, and constraints initial inventory regular midgrade premium property regular midgrade premium min max min max min max aromatics [%] benzene [%] MON [octane] RON [octane] olefins [%] RVI [psi] SPG sulfur [%] inventory [BBl] 21.061 3.821 85.0 92.0 13.00 5.46 0.750 0.008 7000 21.061 3.821 86.0 94.0 16.00 13.00 0.800 0.005 6000 14.569 3.032 87.0 97.0 10.00 15.60 0.800 0.005 6000 0.00 0.00 82.00 92.00 0.00 0.00 0.70 0.00 0 60.00 5.90 100.00 100.00 24.20 15.60 0.81 0.01 50000 0.00 0.00 84.00 94.00 0.00 0.00 0.70 0.00 0 60.00 5.90 100.00 100.00 24.20 15.60 0.81 0.01 50000 0.00 0.00 86.00 96.00 0.00 0.00 0.70 0.00 0 60.00 5.90 100.00 100.00 24.20 15.60 0.81 0.01 50000 introduces nonlinearity, since neither the quality nor the volume in the second and subsequent periods are known. Total volume blended of gasoline grade g in period t is Vblend,g,t: Table 4. Component Supply and Product Demand for Each Period component supply for each period component 1 2 3 ALK BUT HCL HCN LCN LN REF 5000 4000 5600 3000 3000 3000 5000 grade 1 2 3 4 regular midgrade premium 9000 6000 7000 5000 7000 9000 3000 8000 5000 9000 8000 4000 4000 3000 8000 3000 6000 7000 4000 7000 4300 3200 3000 8000 9000 6000 product demand for each period K 4 ∑ Vin,k ,g ,t − Vblend,g ,t = 0; 5000 3200 4000 7000 5000 6300 3000 G G Vin, k , t − Vclose, k , t + Vopen, k , t − The closing inventories in the component and blending tanks at the end of time period, t, must be within the volumetric capacity of the tank. ∀ k, t (5a) Vmin ≤ Vclose, g , t ≤ Vmax ; ∀ g, t (5b) Q i ,min ≤ Q i , g , t ≤ Q i ,max ; ∀ i, g , t (6) The closing inventory of the gasoline tanks is given by the volumetric balance equations: Vblend, g , t − Vclose, g , t + Vopen, g , t − Vout, g , t = 0; ∀ g, t (7) Maximum production capacity is limited by the maximum blender flow rate Fmax: G ∑ Vout,g ,t ≤ Fmax dt ; g=1 (1) ∀t (8) If the volume of the blend Vblend,g,t is specified, since Vout,g,t is equal to the specified product lifting (order) and since Vopen,g,t is known, then the closing volume Vclose,g,t can be computed from eq 7 [Vopen,g,t is specified for t = 1; Vclose,g,t for t = 1 becomes Vopen,g,t for t = 2, etc.]. The algorithm developed in this work (section 4) employs the fact that eq 2 becomes linear if Vblend,g,t is known for every period t (unknown variable is Qi,g,t). Equations 1−7 describe continuous aspects of blend planning. They can be solved by an NLP solver. MINLP formulation includes an additional constraint specifying the switching time required to prepare between blending two different grades of Vopen, g , tQ open, i , g , t − Vout, g , tQ i , g , t − Vclose, g , tQ i , g , t K k=1 Vmin ≤ Vclose, k , t ≤ Vmax ; To ensure that each gasoline grade maintains its physical composition calculated in eq 2 the values of the properties are governed by the following constraints: The physical composition of every gasoline grade is calculated by quality balance equations around the blending tank, eq 2. ∑ Vin,k ,g ,tQ i ,k ,t = 0; ∀ k, t (4) min ∑ ∑ ∑ CkVin, k , g , t + ∑ Vout,k ,g ,t = 0; g=1 T k=1 g =1 t=1 (3) The closing inventory of the blend component tanks is given by the volumetric balance equations. 2. Inventory profiles of the components and of finished gasoline. Note that the multiperiod model implies that, within each time period, the products required during that period are first produced in the required quantities (if they are not available in the storage tanks) and then the products are shipped (“lifted”) in the required quantities. The multiperiod blend planning model contains the unknown closing inventory of each gasoline grade and the unknown associated qualities for the second and subsequent periods. Product terms of these two unknown variables are the nonlinear terms in the model. 2.1. Blend Planning Model. The nomenclature used in the blend planning equations is shown in the Table 1. The objective of blend planning is to minimize the cost of blend recipes over the entire time horizon. K ∀ g, t k=1 ∀ i, g , t (2) In the first time period (t = 1), the initial quality Qopen,i,g,1 and the initial volume Vopen,g,1 are given. The term Vclose,g,tQi,g,t 10711 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Table 5. Blend Volume (% of demand) Patterns from Three NLP Algorithms regular: demand = 19 000 BBl midgrade: demand = 23 000 BBl premium: demand = 19 000 BBl period KNITRO IPOPT CONOPT KNITRO IPOPT CONOPT KNITRO IPOPT CONOPT 1 2 3 4 13.5 23.5 18.8 44.2 29.1 22.8 27.6 20.5 10.5 35.7 6.5 47.4 1.9 28.6 34.9 34.6 18.9 29.2 29.0 22.9 67.8 4.7 27.5 0.0 5.6 47.1 28.5 18.8 24.4 36.7 21.6 17.2 5.3 47.4 26.3 21.1 Figure 3. Blend volume profiles for regular gasoline for three NLP solvers. K gasoline. During the time switching, the blender is idle; that is, some production capacity is lost. This introduces a binary variable that is zero if a given grade is not blended, and it is 1 if a particular grade is blended. Another constraint is the minimum amount to be blended: if a specific grade of gasoline is to be blended, then the amount to be blended has to be greater than or equal to a specified threshold. This introduces another binary variable. 2.2. Lower Bound on the Optimum Solution. Imagine that one can wait until the end of the planning horizon before delivering any amount of the products. In other words, accumulate all components in their storage and then blend the total required amount for each gasoline grade. At that point, it is possible to compute the lowest cost recipes for all gasoline grades that have to be blended over the entire planning horizon. Constraints for this model are • total available amount of each component • total amount of each grade of gasoline that needs to be shipped • available gasoline inventories at the start of the planning horizon • minimum allowed inventories at the end of the planning horizon (for each gasoline grade and each component). The equations for the corresponding model are eqs 9−13 K ∑ Vin,k ,g − Vblend,g = 0; (10) Vblend, g + Vopen, g − Vout, g − Vclose, g = 0; ∀g (11) K Vopen, gQ open, i , g + ∑ Vin,k ,gQ i , k − Vout,gQ i , g − Vclose,gQ i ,g k=1 = 0; ∀ i, g Q i ,min ≤ Q i , g ≤ Q i ,max ; (12) ∀ i, g (13) The solution to the model represents the lower bound on the optimum solution of the multiperiod planning model, since it is computed in the absence of the inventory constraints and it minimizes quality deviations from the constraints. 2.3. Gasoline Blending Example. The example used in this article is based on data provided by AspenTech and Honeywell Canada. There are seven components that are used to produce regular, midgrade, and premium gasoline. Table 2 provides data on the blending components (properties, initial inventories). It is assumed that the properties of the blending components remain constant throughout the planning horizon. Table 3 provides data for the starting inventories (“opening inventories”) for each the three grades of gasoline. Since maximum inventory limits have not been limiting in these experiments, for the sake of brevity, they are not presented. It is assumed that the cost of the starting G min ∑ ∑ Vin, k , gCk k=1 g =1 ∀g k=1 (9) 10712 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Table 6. Optimal MINLP Solution Calculated by α-ECP blend recipes and blend volumes calculated by α-ECP (MINLP) period ALK [%] BUT [%] HCL [%] 1 2 3 4 0.00 0.00 5.11 0.00 29.15 1.87 2.13 1.86 13.57 50.50 48.09 50.58 1 2 3 4 36.82 0.00 0.00 0.00 6.13 3.26 3.23 2.91 24.02 39.81 39.63 36.30 1 2 3 4 0.00 72.30 37.64 70.96 5.78 7.99 6.43 7.83 8.14 0.47 14.60 0.00 HCN [%] LCN [%] Regular 0.01 56.19 1.97 0.42 1.70 1.70 1.98 0.36 total blended volume = 19 000.00 Midgrade 3.02 0.26 1.76 1.74 1.66 2.34 0.02 12.39 total blended volume = 23 000.00 Premium 5.56 1.04 1.93 0.70 1.94 0.60 1.42 3.78 total blended volume = 19 000.00 LNAP [%] RFT [%] blend vol. [BBI] 0.00 0.03 0.00 0.00 1.07 45.21 41.28 45.22 2000.00 8085.68 0.24 8914.08 0.00 0.01 0.00 0.00 29.75 53.43 53.14 48.37 11 950.42 1823.46 8837.22 388.90 16.32 0.00 0.00 0.00 63.16 16.60 38.80 16.01 3045.79 6954.21 8162.55 837.45 LNAP [%] RFT [%] blend vol. [BBI] 0.79 0.78 0.97 0.71 21.16 37.95 36.99 34.96 3339.89 5038.58 5104.75 5516.78 0.33 1.36 0.90 1.29 34.85 43.89 42.56 39.17 10 781.95 3342.09 5386.92 3489.05 1.22 0.58 0.85 1.69 20.41 47.04 46.80 45.71 2660.96 8179.30 5488.45 2671.30 Table 7. Optimal MINLP Solution Calculated by Bonmin blend recipes and blend volumes calculated by Bonmin (MINLP) period ALK [%] BUT [%] HCL [%] 1 2 3 4 18.79 9.43 9.80 16.66 20.08 2.27 2.25 2.70 29.09 44.17 43.12 42.54 1 2 3 4 28.49 17.36 18.42 23.53 5.90 4.25 4.22 4.47 26.46 30.88 31.01 27.98 1 2 3 4 45.53 25.09 21.80 24.19 5.71 5.92 5.63 5.84 15.02 18.61 17.71 15.21 HCN [%] LCN [%] Regular 5.56 4.53 1.37 4.04 1.07 5.80 1.95 0.48 total blended volume = 19 000.00 Midgrade 3.01 0.96 1.96 0.30 1.85 1.04 1.71 1.84 total blended volume = 23 000.00 Premium 5.54 6.57 1.89 0.87 1.01 6.19 0.96 6.39 total blended volume = 19 000.00 Multiperiod solution has been computed by several NLP and MINLP solvers. NLP solvers that have been able to find the optimum are CONOPT, IPOPT, and KNITRO. While the total cost is the same for all solutions, each of these solvers leads to a solution that differs by the amount of each gasoline grade blended in different time periods (see Table 5) and different blend recipes. For the sake of brevity, the blend recipes are not shown, since similar results are presented for MINLP solutions. Table 5 shows what percentage of the total product demand is blended in each period. KNITRO and CONOPT compute blend volumes (blend profile) that vary widely from period to period. IPOPT computes solutions that blend similar amounts in each period, for a given grade. Figure 3 presents the blend profiles for regular gasoline. The MINLP model includes minimum bed size threshold constraints and the capacity lost due to blend switching (we chose not to include blend switching costs in the objective inventory of each gasoline grade is zero. This ensures that the blend plan first uses already existing inventory and then blends additional amount of each grade of gasoline. Data are presented in English system of units, since it is prevalent in the refining industry. Table 4 presents component supply and product demand for each period. 3. BLEND PLAN OPTIMIZATION VIA MINLP AND NLP The same gasoline blending problem has been solved by different solvers, in order to verify whether the different solvers arrive at different blend recipes while the objective function (total cost of all blends) remains the same. We first solved the aggregate blending problem, given by eqs 9−13. The optimal objective function value is 1.411906 × 106; this is a lower bound on the optimum feasible solution. The multiperiod period model, which includes all attributes of the blending system, has optimal solution of 1.428969 × 106. 10713 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Figure 4. Two level algorithm to compute globally optimal gasoline blending solutions. function in order to compare MINLP results with the NLP solutions). The problem has been solved by α-ECP and by Bonmin MINLP solvers (Tables 6 and 7, respectively). Again, the solvers present significantly differing blend recipes and the volumes blended in each period. Given that we arrive at different solutions, solely due to a different choice of a solver, a natural question to ask is whether a blend planner potentially has many solutions from which to choose. Such a question is not answered by the commercial software,6 which presents only a single solution. The algorithm presented in the next section systematically searches for these multiple optimal solutions. 4. TWO LEVEL ALGORITHM TO IDENTIFY MULTIPLE OPTIMUM In order to find multiple optimal solutions to the blend planning, we need to find different values of the blended amounts for every time period, for each grade. As explained in Section 2.1, if the blend volume is specified for each gasoline grade in every period, then one can compute Vclose,g,t from eq 7. Consequently, eq 2 becomes linear, and the entire set of eqs 1−7 can be solved as an LP problem. Hence, a two level algorithm can be constructed: • The top level searches for blend volumes Vblend,g,t of each grade and each period, such that the resulting multiperiod model will have the optimal (lowest) cost of all blends. • The lower level solves the LP model described by eq 2−7 to find the blend recipes that lead to the lowest total cost of all blends. If a specific blending problem is infeasible due to qualities or available volumes of blend components, or due to inventory constraints, the LP at the lower level will be infeasible. If blend properties are to be computed via nonlinear constraints then the top level of the algorithm would remain Figure 5. Differential evolution algorithm. the same, while the lower level would solve an NLP. Note that in such a case the lower level would not have multiple solution, since the blend volumes would be specified by the upper level of the algorithm. In this work, we have selected Differential Evolution7−9,14 (DE) as the algorithm at the top level. Overall structure of this two level algorithm is shown in Figure 4. 4.1. Differential Evolution. Traditional optimization techniques compute single optimum. In contrast, meta-heuristic algorithms simultaneously pursue many possible solutions by creating successive generations of populations, where each member of a population is one solution. The goal is to generate a population where all members are optimal solutions. These algorithms begin by evaluating all members of an initial population, retaining the best members of the current generation, and diversifying the population by mating and mutation of the member of the current population. Some of the algorithms in this class are Evolution Strategies,10,11 Genetic Algorithms,12,13 Differential Evolution, and others. 10714 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Numerical experiments have shown that DE outperforms most of the time genetic algorithms or particle swarm optimization algorithms. These results were the motivation for choosing DE in this work. The DE algorithm is summarized in Figure 5. DE begins with initializing a population pool of Np members. It then randomly selects three mutually exclusive vectors that are not equal to the current target vector being evaluated. These vectors are used to form a trial vector. The trial vector is selected as the member of the new generation if its fitness is no worse than the target vector. One variation of DE, known as “Scheme DE1”, introduced by Storn and Price7 is as follows: Initialization of the Population Members: Optimization variables are represented as a vector of real numbers. Each element of a population vector is randomly generated between specified upper, bupper, and lower, blower, bounds. Every element j of each population member p is randomly initialized as Differential Evolution Algorithm. Differential Evolution (DE) is an optimization algorithm that was originally proposed by Price and Storn in 1995.7 An important property of DE is the ability to “move” across regions where objective function is flat. D{j,p,initial} = rand j(0, 1)(bupper − b lower) + b lower (14) where randj(0,1) is uniformly distributed random number between [0,1]. Figure 6. Global optimum solutions and their neighbors within 1%. Figure 7. Blend volumes and closing inventories for five of the optimal solutions. 10715 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article generating the mutant vector and the type of crossover used. A variation is specified as DE/base/n/crossover where base indicates how the base vector is chosen, n represents the number of difference vectors that participate in mutation, and crossover represents the type of crossover used. Many variations of DE are researched; the most popular is DE/rand/1/bin, where the base vector is randomly chosen, one difference vector participates to produce the mutant vector, and uniform binomial (bin) crossover is used. Population Members and Fitness Function in Gasoline Blend Planning. Each population member represents one solution. Hence, in the work presented here, each member contains the blend volumes for each gasoline grade for each period. Blending three grades of gasoline over the period of 13 days is described by a vector with 39 elements. Each element was scaled to be between 0 and 1, showing what fraction of the maximum possible production capacity Fmax is to be blended for a given grade in a specific time period. If this value was less than a certain threshold, α, then the grade was not blended. The threshold α = (Fmin/Fmax), where Fmin is the minimum acceptable volume of the blend. All parameter values for all initial members of the population were initialized randomly. If the value of the DE element, ug,t is less than α, then the g grade of gasoline is not blended; otherwise, it is. If this grade of gasoline is to be blended, the switching (calibration) time is calculated as well the amount of gasoline g to blend, Vblend,g,t. Mutation: The scaled difference of two randomly chosen vectors is added to a third randomly chosen vector to produce a mutant vector. Three mutually exclusive vector indices: first, second, and base are chosen, with none of them equal to target vector index, p. The difference of the first and the second is scaled by a factor, F. This weighted difference vector is added to the base vector to produce a mutant vector, v. vi , g = D base, g + F *(Dfirst, g − Dsecond, g ) (15) Crossover: The user specifies a crossover probability, Cr, which is a threshold used to decide whether or not recombination of some elements will be carried out between two vectors. The mutant and the target vectors are recombined to produce a trial vector, u. For instance, uniform crossover examines every pair of corresponding elements of the mutant and the target vectors; the element of the trial vector u is chosen from either mutant or target vector, as follows: uj , p , g = vj , p , g = Dj , p , g rand j(0, 1) < Cr if otherwise (16) Selection: If the trial vector has fitness function value that is equal or less than the value of fitness function of the target vector, it is selected to become a member of the next generation population in order to maintain diversity in the population, and the target vector is discarded. This choice enables the DE explore flat surfaces. Dp , g + 1 = ui , g = Dj , p , g ug , t = 0 otherwise ug , t < α (18) Vblend, g , t = Fmax dtug , t fitness(ui , g ) ≤ Dp , g if if (19) Blend Unit Recalibration. Since a single blender is used to blend all grades, every time a different grade of gasoline is blended, the blender needs to be recalibrated for the next blend. The algorithm accounts for the recalibration time as the time when no gasoline can be blended. As a result, the production capacity of the unit is reduced in that period by the amount Vlost,t (lost capacity): (17) Once the new generation population is created, it is evaluated and tested for convergence. This process is repeated until the stopping criterion is satisfied. Stopping criteria is either the maximum number of generations, or some measure of how close are the fitness function values of different members of a population. The variations in DE are introduced by the choice of base vector, number of weighted difference vectors that participate in G ∑ [Vlost,g ,tSg ,t + Vblend, g ,t ] ≤ Fmax dt (19A) g=1 Table 8. Premium Gasoline: Blend Recipes for Five Samples of Globally Optimal Solutions premium gasoline percent blend recipes [%] solution period ALK BUT HCL HCN LCN RFT LNAP A 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 61.16 18.66 0 0 23.88 8.71 0 0 82.76 4.24 0 0 83.85 11.29 57.41 52.04 77.17 34.8 74.04 22.25 5.78 6.1 6.83 4.74 5.41 4.94 9.26 4.74 8.01 4.93 6.83 4.35 7.78 5.25 7.34 7.1 7.75 5.92 8.09 5.35 0 15.16 0 30.1 23.45 24.54 0 30.1 0 28.38 0 25.97 0 25.51 6.76 8.94 2.34 11.82 0 16.92 0 1.95 1.69 2.04 3.86 1.05 37.72 2.04 4.15 2.04 1.69 0 6.4 2.04 2.04 2.04 4.19 0 2.04 0 32.01 0 0 0 0 6.1 0 0 0 0 0 12.5 0 0 0 0 0 12.5 0 12.5 1.06 53.9 74.16 63.11 43.41 54.66 53.02 63.11 5.07 60.41 74.16 57.19 0 55.9 26.45 29.88 8.55 34.96 15.83 42.98 0 4.23 17.32 0 0 0 0 0 0 0 17.32 0 1.96 0 0 0 0 0 0 0 B C D E 10716 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Table 9. Midgrade Gasoline: Blend Recipes for Five Sample of Globally Optimal Solutions midgrade gasoline percent blend recipes [%] solution period ALK BUT HCL HCN LCN RFT LNAP A 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 40.82 39.71 0 78.2 34.77 39.71 0 35.43 0 74.85 0 39.88 43.56 0 0 14.56 0 0 0 0 7.67 5.1 3.31 6.85 9.3 5.1 3.31 4.91 5.31 6.3 3.31 4.71 8.26 3.31 3.31 3.97 5.66 5.19 3.31 3.31 17.71 24.25 40.4 8.61 5.17 24.25 40.4 25.99 38.13 5.84 40.4 20.05 19.3 40.4 40.4 34.48 37.72 13.24 40.4 40.4 3.06 2.04 2.04 2.04 0 2.04 2.04 2.04 3.8 0 2.04 0 4.66 2.04 2.04 2.04 4.11 1.72 2.04 2.04 8.81 0 0 0 43.81 0 0 0 0 12.5 0 12.5 0 0 0 0 0 0 0 0 21.92 28.89 54.25 4.3 6.95 28.89 54.25 31.62 52.77 0.52 54.25 22.85 24.22 54.25 54.25 44.95 52.5 64.22 54.25 54.25 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15.62 0 0 B C D E Table 10. Regular Gasoline: Blend Recipes for Five Samples of Globally Optimal Solutions regular gasoline percent blend recipes [%] solution period ALK BUT HCL HCN LCN RFT LNAP A 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 13.66 0 30.27 0 38.88 0 27.39 0 1.47 38.76 30.79 0 0 0 0 42.11 0 28.21 32.8 0 9.6 5.39 3.24 1.87 11.79 5.39 3.11 1.87 15.68 3.63 3.96 1.48 7.97 5.39 1.87 3.78 12.42 3.15 3.36 1.87 40.61 0 38.39 50.7 19.13 0 39.56 50.7 41.34 34.94 27.85 46.57 32.92 0 50.7 33.58 41.42 39.23 37.36 50.7 3.79 1.45 2.04 2.04 3.73 1.45 2.04 2.04 5.42 2.04 1.86 0 0.59 1.45 2.04 2.04 3.5 2.04 2.04 2.04 0 0 0 0 0 0 0 0 0 0 0.37 12.5 17.48 0 0 0 7.26 0 0 0 32.34 63.99 26.05 45.38 20.12 63.99 27.89 45.38 36.1 20.63 29.29 39.46 36.37 63.99 45.38 18.49 35.39 27.37 24.43 45.38 0 29.17 0 0 6.36 29.17 0 0 0 0 5.88 0 4.67 29.17 0 0 0 0 0 0 B C D E local exploration.9 The proposed initial values8,9 of Cr and F are F ε [0.5, 1], Cr ε [0.8, 1], and NP = 10*nDim, where nDim is the dimensionality of the problem (number of elements in each vector). The values of Cr and F are limited to [0, 1].8 When the value of NP is increased, values of Cr and F should also be increased.8,9 The value of Vlost,g,t is reset to 0 for each period, and it becomes nonzero only when we blend a different grade of gasoline. Mutation and Crossover. DE involves three main parameters: Cr (crossover factor), NP (number of population members), and F (scaling factor). DE is affected by all these three parameters, and hence, the performance of it varies. There is always a tradeoff between the rate of convergence and robustness of an algorithm. Cr determines the number of new components that can be introduced in the next generation. The higher Cr is the faster the convergence is, but the convergence toward the optimum becomes less robust.8,9 Cr acts as a fine-tuning parameter. Scaling factor F determines the scope of exploration of the space of decision vriables. In general, large F favors global exploration, while small F favors 5. DIFFERENTIAL EVOLUTION RESULTS Since every member of DE population represents one set of blend volumes for each period, each generation of the population contains many solutions. Optimal solutions and the solutions that are within 1% of the optimum are shown in Figure 6. Within the first 20 generations, the algorithm finds five optimal solutions, with 10717 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article Figure 8. Exponential crossover rate of convergence. Figure 9. Binomial crossover rate of convergence. value (1 428 969 USD) as it was found by NLP and MINLP algorithms. Each of these solutions blends different amounts of each gasoline grade in each period, leading also to different closing inventories in every time period (Figure 7). Gasoline Blend Recipes. Blend recipes for each optimal solution (computed via LP) differ one from another, since additional optimal solutions being found as 120 generations are generated. Blend Volumes and Closing Inventories. Five optimum solutions (labeled A, B, C, D, and E) from Figure 6 were randomly chosen to illustrate the results. These solutions correspond to population members with same objective function 10718 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719 Industrial & Engineering Chemistry Research Article different amounts are blended in each time period and there are inventory constraints. Tables 8−10 show blend recipes for the five sample solutions. 5.1. Impact of Mutation and Crossover on Convergence Rate. Generations with smaller population sizes (10) were studied to understand the effect of crossover, mutation, and strategy on the algorithm’s convergence. With exponential crossover, the entire population converged in less than 25 generations when the mutation probability was fixed at 0.8 and the crossover probability was varied in the range [0.7, 0.95]. The rate of convergence for each population member using an exponential crossover is shown in Figure 8. As can be seen, the population members converged slightly faster toward optimum solutions as crossover probability was increased from 0.7 to 0.95. Similarly, the population members converged at a faster rate with increase in the crossover probability with binomial crossover (Figure 9). However, the two strategies differed in that the population members converged in lesser number of iterations (≤25) with exponential crossover than for the binomial crossover (≤50). ACKNOWLEDGMENTS This work has been supported by funding from NSERC Canada and from Ontario Center for Excellence. We particularly acknowledge the help from Jeff Kelly of Industrial Algorithms (formerly of Honeywell Canada), who has been a valuable sounding board for this work. ■ ■ 6. CONCLUSIONS We have illustrated existence of multiple optimum solutions for multiperiod blend planning problems by using different NLP and MINLP solvers. Among three NLP solvers (CONOPT, IPOPT, KNITRO), IPOPT arrives at solutions that have the least variation in blend recipe from one period to another. MINLP solvers (α-ECP and Bonmin) also arrive at significantly different solutions; their blend recipes also vary widely from one period to another. These results illustrate that the current practice of formulating multiperiod problem and solving it by some specific solver, simply produces one of many possible options for a blend plan. Common practice to avoid multiplicity of optimal solutions is to augment the economic objective function with some penalty terms, for example, minimize deviation in blend recipes from one period to another or minimize deviations from a preferred blend recipe. Such approaches leave unexplored alternative blend plans that may contain inventory profiles that could be advantages from an operational viewpoint. In order to compute systematically multiple optimum solutions and present the refinery blend planner with multiplicity of choices, we have employed Differential Evolution algorithm combined with a multiperiod LP. Differential Evolution algorithm selects the blend volume for each gasoline grade and each time period. With blend volumes specified, the gasoline blend planning problem becomes linear; hence, it can be solved by an LP. Multiplicity of optimal bend plans provide the refinery planner with options to select blend plans that meet additional considerations, for example, selecting a blend plan that has preferred profile of inventory for some grade or a blend plan that uses lower amounts of alkylate in the earlier periods. Even though the algorithm requires extensive computations, its inherent characteristics make it well suited for parallel computations on multicore CPU machines. Rapidly decreasing price of multicore CPU machines will promote the use of algorithms such as the one presented here. ■ REFERENCES (1) Singh, A.; Forbes, J. F.; Vermeer, P. J.; Woo, S. S. Model-based realtime optimization of automotive gasoline blending operations. J. Process Control 2000, 10, 43. (2) Mendez, C. A.; Grossmann, I. E.; Harjunkoski, I.; Kabore, P. A simultaneous optimization approach for off-line blending and scheduling of oil-refinery operations. Comput. Chem. Eng. 2006, 30, 614. (3) Glismann, K.; Gruhn, G. Short-term scheduling and recipe optimization of blending processes. Comput. Chem. Eng. 2001, 25, 627. (4) Misener, R.; Gounaris, C. E.; Floudas, C. Mathematical modeling and global optimization of large-scale extended pooling problems with EPA complex emissions constraints. Comput. Chem. Eng. 2010, 34, 1432. (5) Li, J.; Karimi, I. A. Scheduling gasoline blending operations from recipe determination to shipping using unit slots. Ind. Eng. Chem. Res. 2011, 50, 9156. (6) Aspen Refinery Multi-Blend Optimizer; Aspen Technology; Burlington, MA, 2009. (7) Storn, R.; Price, K. Differential EvolutionA Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces, Report TR-95-012; International Computer Science Institute; Berkeley, CA, 1995. (8) Storn, R.; Price, K. Differential evolutionA simple and efficient heuristic for global optimization over continuous spaces. J. Global Optimization 1997, 11 (4), 341. (9) Price, K. V.; Storn, R. M.; Lampinen, J. A. Differential Evolution: A Practical Approach to Global Optimization; Springer; Berlin Heidelberg, 2005. (10) Rechenberg, I. Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution; FrommannHolzboog; Stuttgart, 1971. (11) Schwefel, H. P. Kybernetische Evolution als Strategie der experimentellen Forschung in der Stromungs technik. Ph.D. dissertation, Technical University of Berlin; Hermann Fottinger Institute for Fluid Dynamics; Berlin, 1965. (12) Goldberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison-Wesley Publishing; Reading, 1989. (13) Holland, J. H. Outline for a logical theory of adaptive systems. J. Assoc. Comput. Mach. 1962, 9, 297. (14) Das, S.; Konar, A.; Chakraborty, U. K. Two improved differential evolution schemes for faster global search. GECCO’05, Washington, DC, June 25−29, 2005. AUTHOR INFORMATION Corresponding Author *E-mail: mahalec@mcmaster.ca. Notes The authors declare no competing financial interest. 10719 dx.doi.org/10.1021/ie3011963 | Ind. Eng. Chem. Res. 2013, 52, 10707−10719