Alternative Formulations and Improvements with Valid Inequalities for

Article pubs.acs.org/IECR

Alternative Formulations and Improvements with Valid Inequalities for the Refinery Production Scheduling Problem Involving Operational Transitions Lu Zhang,† Yongheng Jiang,*,† and Dexian Huang†,‡ †

Department of Automation, Tsinghua University, Beijing 100084, China Tsinghua National Laboratory for Information Science and Technology, Tsinghua University, Beijing 100084, China

‡

S Supporting Information *

ABSTRACT: Operation mode switching of production units would result in long transitions with fluctuant product yields, production costs, and key product properties. To formulate the remarkable process dynamics, a mixed-integer linear programming (MILP) model involving operational transitions of mode switching was proposed for the refinery production scheduling problem in Shi et al. (Ind. Eng. Chem. Res. 2014, 53 (19), 8155−8170), which can describe transitional behaviors and provide implementable schedules. However, the model is very computationally expensive because it involves a large number of discrete and continuous variables, and the constraints about operational transitions are numerous and complex. In this paper, we study a scheduling problem similar to that in Shi et al. (Ind. Eng. Chem. Res. 2014, 53 (19), 8155−8170) and aim at improving the computational efficiency of MILP models by alternative formulations and valid inequalities. First, we redefine some sets, variables and introduce new variables to propose the basic formulation, with which the nonlinear items are avoided and the variable definitions are simplified. Second, by reformulating constraints to explicitly describe the Fixed Charge Network Flow characteristic of operation mode switching, three Fixed-Charge-Network-Flow-based reformulations are proposed. Third, valid inequalities are developed for lot-sizing relaxations derived from the reformulations. Computational results show that the basic formulation is much smaller-sized and its computational performance is significantly better. The Fixed-Charge-Network-Flowbased reformulations together with the valid inequalities can further reduce by up to more than 95% the computational time while tightening the linear relaxation of the MILP model by more than 45%.

1. INTRODUCTION Refinery production is a complex process which converts crude oils into components through distillation, cracking, refining, and other reactions, then blends components into final products in blender units, and finally stores final products for external demands. There are two types of operation modes defined for different units. For the primary, secondary, and most hydro-upgrading processing units, with the unit-level model predictive control and online optimization technologies,2 the operating conditions are adjusted to guarantee that the unit produces some kind of outflow to the maximum while all the outflows are qualified, and we define an operation mode corresponding with the kind of outflow reaching its maximal production. Since production units can produce several kinds of products simultaneously, a finite number of operation modes are defined. The only determinant of the operation mode is what kind of outflow reaches its maximal production. Therefore, a certain mode asks for different operating conditions for different inflow properties, and mode switching means the operating condition is adjusted accordingly. However, for hydrogen desulfurization units (HDS) and etherification units (ETH), the production target is just to ensure the quality specifications of products. Different inflow properties would result in different yields of the outflows while operating conditions can guarantee the product quality. Therefore, several operation modes are defined according to the inflow properties, and their modes change with the qualities © 2015 American Chemical Society

of their inflows. Because of the dynamic characteristic of the oil refining process, the transitional process from one steady operation mode to another has to be considered for all the production units. Take the atmospheric distillation units (ATM) as an example, when the reflux ratio is adjusted to switch modes, the amount of reflux would vary and the vapor− liquid equilibrium at all trays are disrupted; the concentrations of vaporous materials and liquid materials change, and the outflows would fluctuate; before ATM enters the destined steady mode, the process with obviously fluctuant product yields, production costs, and key product properties would last a long period of time, which is called the operational transition.1 Therefore, the state of the production unit is either a steady operation mode or an operational transition. Refinery scheduling determines which operating states run on the production units, what are the production rate and the consumption rate for each material, and what are the deliveries and the inventory levels of products at any time. The objective of scheduling is to minimize the total cost or to maximize the total profit while satisfying the restrictions of material balance, production capacity, etc. Refinery is a typical energy intensive industry which consumes crude oils as raw materials and utilizes Received: Revised: Accepted: Published: 7871

February 10, 2015 May 19, 2015 July 15, 2015 July 15, 2015 DOI: 10.1021/acs.iecr.5b00600 Ind. Eng. Chem. Res. 2015, 54, 7871−7889

Article

Industrial & Engineering Chemistry Research

linear programming relaxation,29 polynomial algorithm, as well as extend formulation30 for the multi-item capacitated lot-sizing problem. For the refinery production scheduling problem, the FCNF characteristic of operation mode switching is observed. That is, by reformulating mode switching constraints with standard FCNF structure of the flow conservation constraints, the operating states of production units can also be explicitly described as “the virtual flows” through directed graphs. Agra et al.31 studied a single-product maritime inventory routing problem and stated that a FCNF reformulation can lead to a much tighter linear programming relaxation than an original formulation which involved the implicit network construction. In this paper, we develop FCNF-based reformulations as alternative formulations for the refinery production scheduling problem and analyze their effects on lessening computational efforts. Tightening linear relaxations of MILP models via valid inequalities is an effective approach to improve the solution of MILP models. There have been considerable research efforts on constructing valid inequalities for lot-sizing problems in discrete manufacturing, and numerous results have been achieved over the last decades. Pochet and Wolsey32 pointed out that it is a popular method to identify valid inequalities or facet-defining valid inequalities for relaxations of MILP formulations, and two typical relaxations called low-level and high-level relaxations appear in lot-sizing models widely. Lowlevel relaxations consist of the bounds on the variables and a single constraint derived from MILP models, like the n-mixing set33 and the network-dual set.34 Whereas high-level relaxations are canonical production models such as the single-item singlelevel,35 the multi-item single-level36 and the multi-item multilevel lot-sizing model,37 which can be viewed as submodels for lot-sizing problems. Recently, a growing quantity of active researchers have committed themselves to developing valid inequalities involved in process industry planning38 and chemical production scheduling.39,40 Saharidis et al.41 presented an event-based discrete-time formulation for the scheduling of crude oil loading and unloading and combined it with multiple sets of valid inequalities which are implied by operational rules or represent the production information like loading and unloading data. To solve a comprehensive integrated continuous-time optimization model for the scheduling problem involving production units and final products blending, Shah and Ierapetritou42 proposed several valid inequalities to enforce the connections between material flows and utilization of equipment, which can improve the computational performance significantly. In this paper, based on valid inequalities for lot-sizing problems, the tailored valid inequalities for discrete-time MILP models are constructed for the refinery production scheduling problem. In order to solve real-world refinery production scheduling problems involving operational transitions more efficiently, we present alternative formulations and valid inequalities and test their effectiveness in this paper. For a scheduling problem similar to that in Shi et al.,1 we first develop a discrete-time MILP model named the basic formulation with redefined sets, variables, and new variables, which contributes to avoiding nonlinear items and simplifying variable definitions. Second, three FCNF-based reformulations are presented, which involve the constraints with standard FCNF structure to describe the FCNF characteristic of operation mode switching. Third, valid inequalities of lot-sizing relaxations are constructed for the FCNF-based reformulations. Computational tests show that the

fuels for production operations. Considering the great challenges from the crude oil supply, market competition, and environmental protection, operators have viewed schedule as an important approach to obtain better economic benefits and social benefits for refineries. An optimized schedule can help cut down operational costs, improve economic benefits, increase the utilization of production capacities, and reduce the negative effects on the environment. Over the years, various methodologies have been developed to address the scheduling problems and a great many research results have been reported. Most researchers focused on the scheduling of crude oil operations,3−5 essential processing units,6,7 or overall refinery production.8,9 The systems under study were modeled by mathematical programming models employing different time representations10−13 or hybrid Petri nets representing operation decisions from a control perspective.14,15 Also, iterative,16 heuristic,17 hierarchical,18 and inventory pinch19 optimization algorithms were proposed to solve the scheduling problems in reasonable times. However, each operating state of production units in the existing papers referred to the fixed input and operating condition, which was too accurate to be realized exactly in practice. Also, the inevitable operational transitions of mode switching were seldom discussed to reflect the process dynamics in refinery production. With the advanced control and online optimization technologies presented by our research group,2 Shi et al.1 developed an overall refinery production scheduling model (denoted RPSMiT) involving realizable operation modes and operational transitions of mode switching, which can provide implementable schedules. Nevertheless, when longer scheduling horizons are considered for real-world refineries, it is a challenging issue to solve RPSMiT involving a great many binary and continuous variables in reasonable times. This is particularly true when the constraints about operational transitions are included so that the complexity of RPSMiT is increased further. A review20 on advanced solution methods of MIP models for chemical production scheduling mentioned that both alternative formulations and valid inequalities can enhance the computational efficiency of MIP models. A number of references in the open literature21,22 also showed that the improvements of mathematical formulations are crucially important to effectively solve scheduling problems. In this paper, we focus on a formulation with optimized sets and variables, reformulations based on the Fixed Charge Network Flow characteristic, and valid inequalities which can strengthen the reformulations further, so that real-world refinery production scheduling problems involving operational transitions can be solved more efficiently. As an important category of MIP problems, the Fixed Charge Network Flow (FCNF) problem has been a research hotspot in the past decades. The practical applications of FCNF exist in network design,23 transportation,24 and supply chain25 areas extensively. The production planning problem is also a typical class of FCNF due to the flow conservation constraints, which are essential in the standard network flow formulations of lotsizing problems. The flow conservation constraints describe production and inventories as feasible flows along arcs and represent time intervals as nodes in networks. On the basis of them, almost all the research efforts have been made to develop optimization algorithms or solution techniques for FCNF lotsizing problems, like multicommodity reformulation,26 valid inequalities27 and branch-and-cut algorithm28 for the singlecommodity uncapacitated lot-sizing problem and strengthened 7872

DOI: 10.1021/acs.iecr.5b00600 Ind. Eng. Chem. Res. 2015, 54, 7871−7889

Article


Figure 1. Flow sheet of the refinery production system.

RD1, and RD2). The gasoline blender unit and the diesel blender unit blend these components into gasolines and diesels, respectively. A total of five kinds of gasoline products (JIV93, JIV97, GIII90, GIII93, and GIII97) and three kinds of diesel products (GIII0, GIIIM10, and GIV0) are produced to meet the external demands. There are five kinds of blend components (HG, EG, LD, RD1, and RD2) whose key properties vary with the operating states of production units, and these blend components from outlet ports of production units are directly transferred to blender units by oil pipelines, whereas there is a dedicated storage tank for each of the other blend components and final products. All the facilities are interconnected by oil pipelines. The following information are available for the refinery production scheduling problem: a finite number of operation modes and operational transitions with fixed transition times, indicators of operating states such as yields of outflows, technical requirements on operations such as maximum and minimum proportions of components in blending, equipment capacity limitations such as upper and lower bounds on inflow rates of production units, product quality specifications such as permissible concentration ranges of key components of final products, various types of costs, etc. Over the whole scheduling horizon, refinery scheduling determines which operation modes to run and when to switch modes on production units, how many materials are produced or consumed in operations, what are the inventory levels or the backlogging levels of final products. The objective is to minimize the sum of raw material, production, inventory, and backlogging costs. For the detailed problem characteristics, please refer to Shi et al.1

basic formulation is a significantly smaller-sized model and consequently leads to a great reduction in computational effort, and the FCNF-based reformulations together with the valid inequalities contribute to remarkably tighter linear relaxations of MILP models and much-reduced solution times for refinery production scheduling problems. This paper is organized as follows. In section 2, a refinery production scheduling problem similar to that in Shi et al.1 is stated. We present the basic formulation and slightly modify model RPSMiT in Shi et al.1 to describe the problem. In section 3, three FCNF-based reformulations are proposed. In section 4, we derive the lot-sizing relaxations, valid inequalities of which are developed for the FCNF-based reformulations. The computational results are shown in section 5. Some conclusions follow in section 6.

2. PROBLEM STATEMENT AND THE BASIC FORMULATION 2.1. Problem Statement. We first review the flow sheet of the refinery production system as illustrated in Figure 1, which is the same as that in Shi et al.1 The raw material for refinery is assumed to be a kind of mixed crude oil with relatively steady properties and sufficient supply. The production process consists of the following production units: one atmospheric distillation unit (ATM), one vacuum distillation unit (VDU), one fluidized catalytic cracking unit (FCCU), one hydrogen desulfurization unit (HDS), one etherification unit (ETH), two hydro-refining units (HTU1 and HTU2), one catalytic reformer (RF), and one methyl tert-butyl ether unit (MTBE), on which various oil refining tasks run to transform the crude oil into intermediates and finally into five kinds of gasoline blend components (C5, rf, mtbe, HG, and EG) as well as three kinds of diesel blend components (LD, 7873


Article


Figure 2. Illustration of the definitions of binary variables denoting operating states.

physically though. By reclassifying all the materials involved in production process, there are no mixtures of different materials flowing through ports. We only need to define flow rate for each intermediate and blend component, and nonrepeated variable definitions are achieved. (3) In the scheduling, key properties of blend components should be constant so that their proportions in blending operations can be determined. When faced with the five kinds of blend components whose key properties vary with the operating states of production units, Shi et al.1 further classified each of them into certain types according to the operating states under which it is produced. Hence, in model RPSMiT, a set including 47 kinds of blend components was defined, the size of which was much larger than the actual number 8. This can result in more binary variables to formulate the blend operations, so that model RPSMiT is much more difficult to solve. In the basic formulation, with new constraints establishing the correspondence between blend components and current operating states, as long as the operating states are determined, the key properties of blend components are also definite. Therefore, it is considered that there are only eight kinds of blend components, and the redefined set of blend components become much smaller-sized. 2.2.2. Mathematical Formulation. The basic formulation employs the state-task network (STN) representation combined with a common discrete-time grid. It includes five groups of constraints: mode switching constraints express logical relationships in operation mode decisions for production units; capacity constraints specify production, blending and inventory capacity limitations; material balance constraints enforce flow conservation restrictions on all the materials involved; demand constraints describe deliveries of final products; quality constraints ensure quality specifications of final products. The objective function is the total cost to be minimized. An overall description of all the sets, parameters, and variables in the basic formulation is listed in the Notation section. (1). Mode Switching Constraints. We redefine the binary variable ymu,t to denote steady operation modes on production units as follows.

In addition, the refinery production scheduling problem studied in this paper is slightly different from that in Shi et al.1 in the following aspects. (1) With the real ranges of flow rates controlled by valves on pipelines being considered, upper and lower bounds on blending flow rate of each blend component are given as the transfer capacity limitations of pipelines connected to blender units. (2) For more reliable deliveries, we enforce the permissible maximum backorder in schedules. Considering that backorder with no limitation would decrease customer satisfaction to its maximal level, a specified upper bound on the stockout can lessen its serious negative impact on the overall refinery operations. (3) The demand orders are different from those in Shi et al.1 While the orders with delivery time windows for various final products were considered in Shi et al.,1 each demand order here only involves one kind of final product and one single time interval for delivery. This type of demand order has been extensively considered in refinery scheduling problems,16,43 and can be viewed as the decomposition of the multiproduct multiinterval demand order. 2.2. Basic Formulation. 2.2.1. Challenges and Improvements. We first propose the basic formulation with optimized sets and variables, by which the challenges that RPSMiT fails to tackle are overcome. The basic formulation makes the corresponding improvements as follows. (1) Because of the fact that yields of outflows and production costs vary with the operating states of production units, we have to specify the current operating state so as to determine the yield and the production cost. Hence, in model RPSMiT the nonlinear items were derived from the products of different decision variables denoting operating states and inflow rates, respectively. In the basic formulation, the bilinear terms can be avoided by defining new continuous variables to associate inflow rates with operating states and the trilinear terms can be further avoided with the binary variables denoting operation modes redefined. (2) Shi et al.1 considered the one-to-one correspondence between each pipeline and the material transferred by the pipeline. Because of the separators and aggregators on pipelines, the material flow for any unit port may be a mixture of multiple kinds of intermediates and blend components. Hence, in model RPSMiT, flow rate variables were defined for each port as well as each intermediate and blend component. However, this can result in repeated flow rate definitions for the materials transferred by pipelines without connecting to separators or aggregators. In the basic formulation, we view the materials transferred by pipelines with common starting points or terminal points as one kind, which may not be true

yum, t

⎧1, if steady operation mode m runs on unit u ⎪ = ⎨ during time interval t ⎪ ⎩ 0, otherwise

By comparison, it was defined in model RPSMiT as follows. 7874


Article


yum, t

⎧1, if steady operation mode m or operational ⎪ ⎪ transition to mode m runs on unit u =⎨ ⎪ during time interval t ⎪ 0, otherwise ⎩

Necessary Condition. In constraints 4 and 5, if xm,m u,t ′ is equal to 1, the operational transition from mode m to m′ runs on production unit u during time interval t, then the preceding operating state of unit u must be the steady operation mode m, and the left side of constraints 4 and 5 must be at least 1. Similarly, constraints 6 and 7 together make sure that the next operating state of production unit u must be the steady operation mode m′, and they also imply that any two operational transitions cannot run one after another. Whereas the values of time interval t are different in constraints 4−7.

The difference between the definition of ymu,t here and that in Shi et al.1 is illustrated in Figure 2. As shown in Figure 2, production unit u is in steady operation mode C during time interval 1, in steady mode A during time intervals 4, 5, 6, and in steady mode B during time intervals 9, 10. Therefore, the corresponding values of ymu,t in this paper are ⎧1, t = 4, 5, 6 yuA, t = ⎨ ⎩ 0, otherwise

⎧1, t = 9, 10 yuB, t = ⎨ ⎩ 0, otherwise

t−1

∑ t ′= t − TTum , m ′

⎧1, t = 1 yuC, t = ⎨ ⎩ 0, otherwise

m ∈ Mu , m′ ∈ Mu , t > TTum , m ′

yuC, t

∑ yum,t′ ≥ xum,t,m′ ∀ u ∈ U ,

⎧1, t = 7, 8, 9, 10 yuB, t = ⎨ ⎩ 0, otherwise

t ′= 1

m ∈ Mu , m′ ∈ Mu , 2 ≤ t ≤ TTum , m ′

⎧1, t = 1 =⎨ ⎩ 0, otherwise

∑

+

m ∈ Mu

∑

∑ ∑

xum, t, m ′

≤ Tmax − TTum , m ′ Tmax

∑

∑

≥

+

TTum , m ′·(yum, t − 1

+

xum, t, m ′

t + TTum , m ′− 1

∑

− 1) ≤

xum, t,′m ′

t ′= t

∀ u ∈ U , m ∈ M u , m′ ∈ M u , 2 ≤ t < Tmax − TTum , m ′ + 1

−1

(8)

yum, t − 1 + xum, t, m ′ ≤ 1

t ″= t + 1

∀ u ∈ U , m ∈ M u , m′ ∈ M u ,

∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , m ≠ m′

t ≥ Tmax − TTum , m ′ + 1

(2)

(9)

d. State Correlation Constraints. As the upstream unit of HDS and ETH, FCCU changes its outflows with its operating states. Since the states of HDS and ETH are only determined by their inflow properties, the three units FCCU, HDS, and

yum, t + yum, t′+ 1 ≤ 1 ∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , m ≠ m′, t < Tmax

(7)

c. Completeness Constraints. Completeness constraints ensure that the duration of any operational transition must be equal to its transition time. This means that any operational transition cannot be interrupted by a new mode switching. Equations 8 and 9 here impose this restriction as follows. For any u ∈ U, m ∈ Mu, and m′ ∈ Mu, if ymu,t − 1 and ym,m u,t ′ are both equal to 1, the production unit u starts the operational transition from mode m to m′ at time interval t, then the transition must last until time interval t + TTm,m u ′ − 1 to ensure the completeness, and we have constraint 8. However, at the time intervals t ≥ Tmax − TTm,m u ′ + 1 which are so late that operational transitions cannot run completely at the end of the scheduling horizon, the start-up of transition is not allowed, which is guaranteed by constraint 9.

(1)

t ∈ T , t ′ ∈ T , t ′ > t + 1, t ′ < t + 2TTum , m ′ + 2

yum, t ′ ≥ xum, t, m ′

∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , Tmax − TTum , m ′ < t ≤ Tmax − 1

t ′− 1

yum, t′′

(6)

t ′= t + 1

b. Consistency Constraints. Consistency constraints maintain the consistency of an operational transition and its preceding as well as next steady operation modes, which is expressed as a sufficient condition and a necessary condition in eqs 2−7 here. Sufficient Condition. In constraint 2, if ymu,t and ymu,t′′ are both equal to 1, then ∑t″ t=′ −t +1 1 xm,m u,t″ ′ must be greater than 0. This represents that if any two different steady operation modes run on any production unit successively, there must be the corresponding operational transition between these two steady modes. The requirements on time intervals in constraint 2 are t′ > t + 1, t′ < t + 2TTm,m u ′ + 2, which ensure that the steady modes m and m′ run successively. Furthermore, two different steady operation modes cannot run on any production unit during two adjacent time intervals, which is indicated by constraint 3. Hence, at most one of ymu,t and yu,tm′+ 1 can be equal to 1 in constraint 3. yum, t

yum, t ′ ≥ xum, t, m ′ ∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , t

t ′= t + 1

= 1 ∀ u ∈ U, t ∈ T

m ∈ Mu m ′∈ Mu

xum, t,″m ′

(5)

t + TTum , m ′

a. Allocation Constraints. Allocation constraints describe that either a steady operation mode or an operational transition is allocated to any production unit during any time interval. This category of constraints consists of eq 1 here and eqs 5, 6, and 7 from RPSMiT. On the basis of the redefinition of variable ymu,t, for any u ∈ U and t ∈ T, one and only one binary variable ymu,t or xm,m u,t ′ must be equal to 1, which means that a single steady operation mode or operational transition should run on any production unit during any time interval. yum, t

(4)

t−1

Whereas the corresponding values of ymu,t in Shi et al.1 are ⎧1, t = 2, 3, 4, 5, 6 yuA, t = ⎨ ⎩ 0, otherwise

yum, t ′ ≥ xum, t, m ′ ∀ u ∈ U ,

(3) 7875


Article


lower bounds restrict the amount of inflow consumed by unit u under steady mode m during time interval t. The same capacity limitations on production units under operational transitions are enforced by constraint 15.

ETH can always work in the same operating states. Considering that the operation modes of other units are only determined by the outflow reaching its maximal production, with different inflow properties, the operating conditions are adjusted correspondingly so that the units can work in desired operation modes. Specifically, if the properties of the raw material processed by the ATM change, in order to maintain the production performance required by schedule, the operating conditions of the ATM are indeed changed and the ATM would be in a desired mode. All the production units other than HDS and ETH are in the similar situations, and the operating states of them can be changed independently. Constraint 10 ensures that the same kind of steady operation mode runs on HDS and ETH simultaneously, where the mode of HDS denoted m1 corresponds with the mode of ETH denoted m2. Since the known transition times of these two production units are equal, constraint 10 also implies that identical operational transitions run on HDS and ETH simultaneously. Thus, they must be in the same operating states at any time.

QIuminyum, t ≤ QIum, t ≤ QIumaxyum, t ∀ u ∈ U , m ∈ Mu , t ∈ T

QIuminxum, t, m ′ ≤ QIum, t, m ′ ≤ QIumaxxum, t, m ′ ∀ u ∈ U , m ∈ M u , m′ ∈ M u , t ∈ T

(10)

For FCCU and HDS, constraints 11 and 12 enforce that they must start the same kind of operational transition simultaneously, where the modes of FCCU denoted m1 and m1′ correspond with the modes of HDS denoted m2 and m2′, m1,m1 respectively. In constraint 11, if ym1 FCCU,t − 1 and xFCCU,t′ are both equal to 1, FCCU starts the operational transition from mode m1 to m1′ at time interval t, then HDS must also start the identical transition at time interval t, that is, ym2 HDS,t − 1 and xm2,m2 HDS,t ′ are both equal to 1. In a similar way, constraint 12 is developed. Also, since the known transition times of FCCU are longer than those of HDS by one time interval, if steady operation mode m2 runs on HDS during time interval t, then either steady mode m1 or a transition to m1 must run on FCCU simultaneously. This is stated by constraint 13.

M minboc , o , t ≤ Q oc , o , t ≤ M max boc , o , t ∀ oc ∈ OC , o ∈ O , t ∈ T

∀ m1 ∈ MFCCU , m1′ ∈ MFCCU , m2 ∈ MHDS , (11)

∑ ∑ m2 2·(yHDS, t−1

+

m2, m2 ′ x HDS, t

− 1) ≤

m1 yFCCU, t−1

+

m1, m1 ′ x FCCU, t

u

m1 yFCCU, + t

u

∑ ∑ u

∑ ∑ ∑

feed u , oiQIum, t +

m ∈ Mu

yield mu ,,oim ′QIum, t, m ′


∑ ∑ ∑ u

feed u , oiQIum, t, m ′


∀ oi ∈ OI , t ∈ T

(12)

(17)

m1 ′ , m1 m2 ∑ xFCCU, t ≥ yHDS, t

b. Material Balance Constraints for Storage Tanks. In this category of constraints, eqs 18 and 19 are flow conservation restrictions of blend component tanks, and eqs 20 and 21 are those of final product tanks. Constraint 18 represents that the inventory level of blend component oc during time interval t is equal to that during the previous time interval t − 1 adjusted by the production and consumption at production units as well as the consumption in blending operations. Constraint 19 is for the special case in the first time interval, which involves the initial inventory levels of blend components.

m1′

∀ m1 ∈ MFCCU , m1′ ∈ MFCCU , m2 ∈ MHDS , m1 = m2, t ∈ T

yield mu , oiQIum, t +

m ∈ Mu

=

∀ m1 ∈ MFCCU , m1′ ∈ MFCCU , m2 ∈ MHDS , m2′ ∈ MHDS , m1 = m2, m1′ = m2′, t ≥ 2

(16)

c. Inventory Capacity Constraints. This category of constraints specifies the upper and lower bounds on inventory capacities for each storage tank. Equations 21 and 22 from RPSMiT are included in the basic formulation. Especially, for the five kinds of blend components with varied key properties, min the upper bounds INVmax oc and lower bounds INVoc are all set to 0 to represent that no tanks are for their inventories. (3). Material Balance Constraints. a. Material Balance Constraints for Intermediates. The left side of constraint 17 is the production of any intermediate oi, and its consumption is represented by the right side of constraint 17. They must be equal to each other during any time interval t.

m1 m1, m1 ′ m2 m2, m2 ′ 2·(yFCCU, + x FCCU, t − 1) ≤ yHDS, t − 1 + x HDS, t t−1

m2′ ∈ MHDS , m1 = m2, m1′ = m2′, t ≥ 2

(15)

b. Blending Capacity Constraints. Equation 16 imposes the blending capacity restrictions, where binary variables boc,o,t denote whether any blend component is blended to produce any final product during any time interval. If and only if boc,o,t is equal to 1, the amount of component oc blended into final product o during time interval t is limited by the upper and lower bounds on blending flow rates. Otherwise no component oc is consumed in blending operation for final product o during time interval t.

m1 m2 yHDS, = yETH, , ∀ t ∈ T, t t

m1 ∈ MHDS , m2 ∈ METH , m1 = m2

(14)

(13)

(2). Capacity Constraints. a. Production Capacity Constraints. Equations 14 and 15 here enforce the production capacity limitations. For any production unit, new continuous variables QImu,t and QIm,m u,t ′ associate the rates of inflow with its operating states, and they denote the inflow rates under steady operation modes and operational transitions, respectively. In constraint 14, if and only if ymu,t is equal to 1, the upper and 7876


Article

Industrial & Engineering Chemistry Research INVoc , t = INVoc , t − 1 +

u

+∑ u

−∑ u

−

∑ ∑

m ∈ Mu

yield mu ,,ocm ′QIum, t, m ′ −

∑ ∑


∑ ∑

and the flow rates for materials are represented by these m m products involved in eqs 17−19, such as yieldu,oi QIu,t , m,m feedu,ocQIu,t ′, etc. Thus, the reclassification of materials in production process contributes to the nonrepeated variable definitions, and the types of intermediates are reduced as well. (4). Blending Operation Constraints. The upper and lower bounds on blending ratios are given by eq 23 from RPSMiT. The permissible concentration ranges of key components of final products are represented by eq 24 from RPSMiT. In addition, new eqs 22−24 as follows are established for the five kinds of blend components with varied key properties, so that the redefined set of blend components is simplified. Constraints 22 and 23 associate the components in blend operations with the current operating states of production units, hence the operating states under which the blend components are produced can be specified. Among them, voc denotes the five kinds of blend components, subvoc − s denotes the blend component voc produced under the operating state s, and s represents steady operation mode or operational transition in eqs 22 and 23, respectively. Since a single operating state runs on any production unit at any time, constraints 22 and 23 ensure that at most one Qsubvoc−s,o,t can be greater than 0 for any voc ∈ VOC, o ∈ O and t∈T.The relationship between Qvoc,o,t and Qsubvoc−s,o,t are established by constraint 24. On the basis of eqs 22−24, it is considered that there are only eight kinds of blend components with constant key properties, and the simplified blend components set further reduces the number of boc,o,t, which are the binary variables related to blend operations.

yield mu , ocQIum, t

∑ ∑

u

feed u , ocQIum, t

m ∈ Mu

feed u , ocQIum, t, m ′


∑ Q oc ,o,t ∀ oc ∈ OC , t ≥ 2 o

(18)

INVoc , t = INVinioc +

∑ ∑ u

+∑ u

−∑ u

∑ ∑

yield mu , ocQIum, t

m ∈ Mu

yield mu ,,ocm ′QIum, t, m ′

−

∑ ∑


∑ ∑

u

feed u , ocQIum, t, m ′ −


feed u , ocQIum, t

m ∈ Mu

∑ Q oc ,o,t o

∀ oc ∈ OC , t = 1 (19)

For any t ∈ T, the inflow of storage tank for any final product o is the production from blending operations, and the outflow is used for delivery. Constraint 21 is for the special case in the first time interval, and the initial inventory levels of final products are considered. INVo , t = INVo , t − 1 +

∑ Q oc ,o,t − Do,t ∀ o ∈ O , t ≥ 2 oc

(20)

INVo , t = INVinio +

∑ Q oc ,o,t − Do,t ∀ o ∈ O , t = 1

0 ≤ Q subvoc − s, o , t ≤ M max yum, t

oc

(21)

∀ subvoc − s ∈ SUBVOC, o ∈ O , t ∈ T

The material balance constraints benefit from two important improvements made by the basic formulation, which are the avoidance of nonlinear items and the simplification of variable definitions. First, the bilinear and trilinear terms in material balance constraint 12 from RPSMiT are avoided by defining new m m,m variables QImu,t, QIm,m u,t ′ and redefining variables yu,t. QIu,t ′ denotes inflow rate when and only when operational transition from mode m to m′ runs on production unit u at time interval t, and it is equivalent to the product of a binary variable denoting the operational transition and a continuous variable denoting the inflow rate of unit u at time interval t. Thus, the corresponding bilinear term xu,m,m′,tQIu,t in eq 12 from RPSMiT m can be replaced by QIm,m u,t ′. With the redefined yu,t, steady operation mode m runs on production unit u during time interval t if and only if ymu,t is equal to 1. Hence, the trilinear term yu,m′,t · (1 − ∑m xu,m,m′,t) · QIu,t in eq 12 from RPSMiT can be rewritten as the bilinear term yu,m′,t · QIu,t. Similar to the case of m QIm,m u,t ′, yu,m′,t · QIu,t can be further replaced by QIu,t. Thus, all the nonlinear terms in material balance constraint 12 from RPSMiT are avoided. Second, for pipelines with common starting points or terminal points, we reclassify the materials transferred by them into the same kind, so that each inflow and outflow of production units can be a certain kind of intermediate or blend component. Hence, the flow rate variables for ports are no longer needed. We multiply new variables QImu,t or QIm,m u,t ′ by parameters which denote yields of outflows or kinds of inflows,

(22)

0 ≤ Q subvoc − s, o , t ≤ M max xum, t, m ′ ∀ subvoc − s ∈ SUBVOC, o ∈ O , t ∈ T Q voc, o , t =

∑

(23)

Q suboc − s, o , t

suboc − s

∀ voc ∈ VOC, o ∈ O , t ∈ T

(24)

(5). Demand Constraints. Equation 25 represents the deliveries of final products to satisfy the external demands. For any t ∈ T for delivery, two restrictions are enforced by constraint 25: the delivery of any final product is no more than its demand; the backlogging of any final product is limited by its permissible maximum shortage. 0 ≤ Ro , t − Do , t ≤ Shortmax o , t ∀ o ∈ O , t ∈ T

(25)

(6). Objective Function. Same as those in material balance constraints, the bilinear and trilinear terms in objective function 29 of RPSMiT can also be avoided. In the objective function 26 here, the costs summed in sequence are crude oil costs under different operating states, production costs of units under different operating states, inventory costs, and backlogging costs. 7877


Article


Figure 3. “Transferring flow of modes” for a production unit.

min z =

t

+

∑∑ m

m m,m′ QIATM, t Crudeprice)

m′

+∑ ∑ (∑ QIum, t OpCost mu + t

u

(1) With the avoidance of nonlinear items, all the auxiliary variables xyu,m′,t, xQIu,m,m′,t, xQI1u,m,m′,t, xyQIu,m′,t, xyQI1u,m′,t and the auxiliary constraints 30−43 for linearization in m-RPSMiT no longer exist in BF. Specifically, the number of binary variables xyu,m′,t is |U| · |Mu| · |T|, where U denotes the set of production units, Mu denotes the set of operation modes of production unit u, and T denotes the set of time intervals within the scheduling horizon; the total number of continuous variables xQIu,m,m′,t, xQI1u,m,m′,t, xyQIu,m′,t and xyQI1u,m′,t is 2 · |U| · |Mu| · |T| · (1 + |Mu|); the total number of eqs 30−43 in mRPSMiT is |U| · |Mu| · |T| · (9 + 5|Mu|). (2) Owing to the simplification of flow rate variable definitions, we reduce the number of the following continuous variables and constraints. Specifically, flow rate variables QIu,t, QIu,oi,t, QOu,s,t, QOu,oi,t and QIu,oc,t in m-RPSMiT are replaced by QImu,t and QIm,m u,t ′ in BF. The total number of the former is |U| · | T| · (|S| + 2|OI| + |OC′| + 1), where S denotes the set of production unit ports, OI denotes the intermediates set in which the intermediates are almost as many as the pipelines used, and OC′ denotes the nonredefined set of blend components. Whereas the total number of the latter is |U| · | Mu| · |T| · (1 + |Mu|). Since the types of operation modes are much fewer than those of intermediates and blend components, a large number of continuous variables are reduced. Also, we no longer need eq 12 in m-RPSMiT and the constraints which establish the correspondence between flow rates for ports and flow rates for materials, the total number of which is |U| · |T| · (1 + 2|S|). (3) The redefined set of blend components OC reduces its size from |SUBVOC| + |OC| − |VOC| in m-RPSMiT to |OC| in BF, where VOC denotes the set of blend components voc with varied key properties and SUBVOC denotes the set of subcomponents of blend component voc. Hence, the reduction of binary variables boc,o,t is (|SUBVOC| − |VOC|) · |O| · |T|, where O denotes the set of final products. For the refinery production system studied in this paper, the size of set OC is reduced from 47 in m-RPSMiT to 8 in BF, and the corresponding reduction in boc,o,t is 39 · |O| · |T|. So far, the total reduction in binary variables are |T| · (|U| · | Mu| + (|SUBVOC| − |VOC|) · |O|), the total reduction in continuous variables are |U| · |T| · (|Mu| + |Mu|2 + |S| + 2|OI| + | OC′| + 1), and the total reduction in constraints are |U| · |T| · (9|Mu| + 5|Mu|2 + 1 + 2|S|). In conclusion, on the basis of the improvements above, BF is much smaller-sized than mRPSMiT and allows us to solve the refinery production scheduling problems involving transitions with significantly less computational efforts.

m ∑ (∑ QIATM, t Crudeprice

m

∑ ∑ QIum,t,m′ m

m′

tOpCost mu , m ′) +∑ α(∑ INVoc , t + t

oc

∑ INVo,t) + β ∑ ∑ (Ro,t − Do,t ) o

o

t

(26)

Therefore, we obtain the basic formulation denoted BF as follows. BF: min z

s.t. eqs 1−25 in this paper and eqs 5−7, 21−24 from RPSMiT. 2.3. Modified RPSMiT. In order to compare with model BF, the mathematical formulation of RPSMiT is slightly modified as follows to adapt to the scheduling problem here. Equation 16 in this paper is included to impose the blending capacity restrictions. To represent the single-product singleinterval demand orders, eqs 15, 17, and 18 from RPSMiT are replaced by eqs 20 and 21 in this paper, and the objective function 29′ of RPSMiT should be rewritten as the objective function 27 here. Equations 25−28 from RPSMiT are replaced by eq 25 in this paper to involve the given maximum backorder. min z′ =

∑ (QIATM,t OPC t

+

∑ ∑ ∑ xQIu ,m,m′ ,t tOpCost u ,m,m′ u

+

m

∑ ∑ xyQIu ,m′ ,t OpCost u ,m′) u

+

m′

m′

∑ α(∑ INVo,t + ∑ INVoc ,t) t

o

oc

+ β ∑ ∑ (Ro , t − Do , t ) o

t

(27)

Therefore, we obtain the modified RPSMiT denoted mRPSMiT as follows: m‐RPSMiT: min z′

s.t. eqs 16, 20, 21, 25 in this paper and eqs 1−14, 16, 19−24, 30−43 from RPSMiT. 2.4. Summary of Model Improvements. The major improvements in modeling made by BF are summarized as follows. 7878


Article


3. FCNF-BASED REFORMULATIONS From computational tests, we observe that the integrality gap between the optimal solutions of MILP model BF and its linear relaxation is large. Thus, it is a natural choice to strengthen BF for an accelerated solution. Considering the FCNF characteristic of operation mode switching, we would like to reformulate relevant constraints with standard FCNF structure and develop FCNF-based reformulations with tighter linear programming relaxations. The standard FCNF problem determines a minimum cost arc combination to transfer objects from starting nodes to end nodes on a given digraph. The total costs include the fixed charge of arc usage and the transferring charge which varies with the amount transported along arcs. For the refinery production scheduling problem, a definite operating state is determined for each production unit during every time interval, which can be either a steady operation mode or an operational transition between two different steady modes. Thus, as shown in Figure 3, for any production unit, the operating states form a consecutive “transferring flow of modes” over the entire scheduling horizon, which represents the logical relationship between two adjacent states. In Figure 3, over a scheduling horizon with 15 time intervals, there are four realizable operation modes called m1,m2,m3,m4 for a production unit. The bold horizontal arcs represent that a certain steady operation mode runs on the production unit and lasts at least one time interval, whereas the bold oblique arcs represent that an operational transition from one steady mode to another is running within the corresponding time intervals. At first, the unit is in steady mode m1, and its next operating state can only be either the same steady mode or a transition from m1 to other modes. We can see that the unit remains in m1 during time interval 2 and starts the transition from m1 to m2 at time interval 3. Then, it can be deduced that the unit will be in steady mode m2 after the transition. As we can see, the transition lasts three time intervals and the unit indeed enters steady mode m2 at time interval 6. The description and analysis of subsequent operation mode switching are similar. Over the whole scheduling horizon, the mode of this unit starts from m1, then is transferred to m2 and m4 successively, and finally returns to m1 again. At any time in the scheduling horizon, the production cost of any production unit depends on its operating state and the amount of inflow. In order to minimize the total production cost, for each unit, the scheduling chooses an optimal arc combination to transfer the operating state from an initial steady mode to a terminal steady mode and determines the amount of inflow under each operating state. Hence, the FCNF characteristic of operation mode switching is sufficiently indicated by the logical relationship involved in operation mode switching and the production cost related to operating states. The three FCNF-based reformulations are developed as follows. 3.1. First FCNF-Based Reformulation. Since we know the transition times of all production units, if any unit starts an operational transition at a certain time interval, then we can identify the subsequent time intervals when it is still in the transition or enters the destined steady mode. Therefore, only when the production unit is in a certain steady operation mode, the scheduling needs to determine whether it maintains the current steady mode or starts a transition at the next time interval; once a transition is started, no decisions on operating

states of the production unit are needed during the transition. In this way, the decision-making effort can be reduced and the solution process of the scheduling problem can be accelerated. We combine this idea with the FCNF characteristic of operation mode switching to propose the first FCNF-based reformulation. Before reformulating model BF, we redefine the binary variable xm,m u,t ′ to denote the start of operational transitions on production units as follows. xum, t, m ′

⎧1, if operational transition from mode m to m′ ⎪ = ⎨ starts running on unit u at time interval t ⎪ ⎩ 0, otherwise

By comparison, it was defined in RPSMiT and BF as follows. xum, t, m ′

⎧1, if operational transition from mode m to m′ ⎪ = ⎨ runs on unit u during time interval t ⎪ ⎩ 0, otherwise

The difference between the definition of xm,m u,t ′ here and that in Shi et al.1 is also illustrated in Figure 2. We can see that for production unit u, operational transition from mode C to mode A lasts during time intervals 2, 3, and operational transition from mode A to mode B lasts during time intervals 7, 8. Therefore, the corresponding values of xm,m u,t ′ in this paper are ⎧1, t = 2 ⎧1, t = 7 xuA,B xuC,A ,t = ⎨ ,t = ⎨ ⎩ 0, otherwise ⎩ 0, otherwise 1 Whereas the corresponding values of xm,m u,t ′ in Shi et al. are

⎧1, t = 2, 3 A,B ⎧1, t = 7, 8 xuC,A xu , t = ⎨ ,t = ⎨ ⎩ 0, otherwise ⎩ 0, otherwise

Then, we give the following constraints to describe the FCNF characteristic explicitly. (1). Mode Switching Constraints. a. Allocation Constraints. In model BF, this category of constraints consists of eq 1 in this paper and eqs 5, 6, and 7 from RPSMiT. Equation 1 can be rewritten as eq 28, and eqs 5 and 7 from RPSMiT can be rewritten as eqs 29 and 30, respectively.

∑

yum, t +

m ∈ Mu

∑ ∑

xum, t, m ′ ≤ 1


∀ u ∈ U, t ∈ T

(28)

For any u ∈ U and t ∈ T, if any steady operation mode m runs, only the corresponding variable ymu,t is equal to 1; if any operational transition from mode m to m′ starts running, only the corresponding variable xm,m u,t ′ is equal to 1; if any transition has started and is running, all the variables ymu,t and xm,m u,t ′ are equal to 0. Thus, we have constraint 28.

∑

yum, t1 = 1 ∀ u ∈ U

m ∈ Mu

∑ m ∈ Mu

yum, T

max

(29)

=1∀u∈U (30)

Constraints 29 and 30 enforce that the initial and terminal operating states of production units must be steady operation modes. 7879


Article

Industrial & Engineering Chemistry Research t

b. Consistency Constraints and Completeness Constraints. In model BF, consistency constraints and completeness constraints are two categories of constraints involving eqs 2−7 and 8 and 9 in this paper, respectively. They can be rewritten as eqs 31−33 which impose the restrictions of consistency and completeness simultaneously.

∑

m′ ′ xum, t,−mTT = m,m′ + y u,t−1 u

m ∈ Mu

∑

QIumin

t ′= 2

(31)

′ xum, t,−mTT m , m ′ = 0 ∀ u ∈ U , m ∈ M u , m′ ∈ M u , u

yum, t − 1 ≥

∑

FR1: min z

(32)

s.t. eqs 10, 14, 16−25, 28−37 in this paper and eqs 6, 21−24 from RPSMiT. 3.2. Second FCNF-Based Reformulation. On the basis of model FR1, the first time interval and the subsequent time intervals of any steady state are distinguished by defining new variables. Considering the restriction that any production unit can start an operational transition only after its current steady mode has lasted at least one time interval, we define new binary variables yamu,t and ybmu,t to represent the first time interval and the subsequent time intervals of any steady operation mode, respectively, for any u ∈ U, m ∈ Mu, and t ∈ T, yamu,t denotes whether steady mode m starts running on unit u at time interval t, and ybmu,t denotes whether time interval t falls into the subsequent time intervals when unit u maintains steady mode m. In this way, the first time interval of any steady mode is distinguished from the corresponding subsequent time intervals. The mode switching constraints involving ymu,t in model FR1 can be reformulated with yamu,t and ybmu,t. a. Allocation Constraints. In model FR1, allocation constraints consist of eqs 28−30 in this paper and eq 6 from RPSMiT. New constraint 38 is added to this category of constraints, which can build a relationship between yamu,t, ybmu,t, and ymu,t. Equation 29 needs to be rewritten as eq 39.

xum, t, m ′ ∀ u ∈ U , m ∈ Mu , t ≥ 2 (33)

m ′∈ Mu

With the standard FCNF structure, constraint 31 can be viewed as a flow conservation restriction on operation modes at the decision-making nodes. During any time interval t − 1, if any production unit u is in the last time interval of any operational transition to mode m′ (∑m∈muxm,m ′ um,m′ is equal to 1) or in the u,t‑TT m steady operation mode m′ (yu,t′ − 1 is equal to 1), then in the next time interval t, the unit u can start a new transition from mode m′ (∑m∈muxmu,t′,m is equal to 1) or enter/maintain the steady mode m′ (ymu,t′ is equal to 1). However, two operational transitions are not allowed to be one after another in the schedules, and we enforce this restriction by constraint 33 which is complementary to constraint 31. Constraints 31−33 together with all the allocation constraints describe the “transferring flow of modes” over the entire scheduling horizon. c. State Correlation Constraints. In model BF, this category of constraints consists of eqs 10−13 in this paper. Among them, eqs 11 and 12 can be rewritten as eq 34 and eq 13 can be rewritten as eq 35. m1, m1 ′ m2, m2 ′ x FCCU, t = x HDS, t

yaum, t + ybum, t = yum, t ∀ u ∈ U , m ∈ Mu , t ∈ T

∀ t ∈ T , m1 ∈ MFCCU , m1′ ∈ MFCCU , m2 ∈ MHDS ,

∑

m2′ ∈ MHDS , m1 = m2, m1′ = m2′

(35)

Constraint 34 represents that FCCU and HDS must start the same kind of operational transitions simultaneously. Since the transition times of FCCU are longer than those of HDS, if a certain steady operation mode runs on FCCU during time interval t (ym1 FCCU,t is equal to 1), then HDS must also be in the corresponding steady mode simultaneously (ym2 HDS,t is equal to 1) but not vice versa. Thus, we have constraint 35. (2). Capacity Constraints. We focus on the production capacity constraints in model BF for reformulation, and eq 15 in this paper, which involves xm,m u,t ′, can be rewritten as follows. t

∑

∑

t ′= t − TTum , m ′+ 1

∑

m′ m,m′ ′ xum, t,−mTT m , m ′ = yau , t ∀ u ∈ U , m′ ∈ M u , t > TTu u

m ∈ Mu

(40)

yaum, t − 1

+

ybum, t − 1

=

ybum, t

+

∑

xum, t, m ′

∀ u ∈ U , m ∈ Mu , t ≥ 2

m ′∈ Mu

(41)

Constraint 40 states that any production unit must be in the destined steady mode after the corresponding operational transition. Constraint 41 can also be viewed as a flow conservation constraint with standard FCNF structure at the nodes where operating states need to be decided. If and only if the left side of constraint 41 is equal to 1, any steady mode m runs on any unit u during time interval t − 1, then the unit u can either maintain the steady mode m (ybmu,t is equal to 1) or

t

xum, t,′m ′ ≤ QIum, t, m ′ ≤ QIumax

(39)

b. Consistency Constraints and Completeness Constraints. In model FR1, this category of constraints consists of eqs 31−33 in this paper. Equations 31 and 33 can be rewritten as eqs 40 and 41.

m2 m1 yHDS, ≥ yFCCU, ∀ t ∈ T , m1 ∈ MFCCU , m2 ∈ MHDS , m t t

QIumin

(38)

yaum, t1 = 1 ∀ u ∈ U

m ∈ Mu

(34)

1 = m2

(37)

For any u ∈ U, m ∈ Mu, and m′ ∈ Mu, if the sums of xm,m u,t′ ′ in constraints 36 and 37 are equal to 1, the production unit u must be in the operational transition from mode m to m′ during time interval t. Constraints 36 and 37 only differ in the values of time interval t. The objective function of this reformulation is the same as that of model BF. Therefore, we obtain the first FCNF-based reformulation denoted FR1 as follows:

m ∈ Mu

t ≤ TTum , m ′

t ′= 2

∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , 2 ≤ t < TTum , m ′

xum, t′ , m + yum, t′

∀ u ∈ U , m′ ∈ Mu , t > TTum , m ′

t

∑ xum,t,′m ′ ≤ QIum,t,m ′ ≤ QIumax ∑ xum,t,′m ′

xum, t,′m ′

t ′= t − TTum , m ′+ 1

∀ u ∈ U , m ∈ Mu , m′ ∈ Mu , t ≥ TTum , m ′ (36) 7880


Article

Industrial & Engineering Chemistry Research start transiting to another mode (∑m′∈Muxm,m u,t ′ is equal to 1) at time interval t. Furthermore, adding constraints 40 and 41 together can establish constraint 31. c. State Correlation Constraints. In model FR1, this category of constraints consists of eqs 10, 34, 35 in this paper. Equation 35 can be rewritten as follows. m2 yaHDS, t−1

=

them. We use oi6 to denote the intermediate which is produced by FCCU and consumed by HDS. The constraint 17 on oi6 can be rewritten as follows. 1 m1 m2 yield mFCCU, QIAFCCU, t1 = feedHDS, oi6QIAHDS, t1 oi6

∀ m1 ∈ MFCCU , m2 ∈ MHDS , m1 = m2

m1 yaFCCU, t

∀ t ≥ 3, m1 ∈ MFCCU , m2 ∈ MHDS , m1 = m2

1 m2 yield mFCCU, QI m1 = feedHDS, oi6QIBHDS, t oi6 FCCU, t

(42)

∀ t ≥ 2, m1 ∈ MFCCU , m2 ∈ MHDS , m1 = m2

m2 m1 yaHDS, t1 = yaFCCU, t1

∀ m1 ∈ MFCCU , m2 ∈ MHDS , m1 = m2

m1 ∈ MFCCU

∑

∀ t ≥ 2, m1′ ∈ MFCCU , m2′ ∈ MHDS , m1′ = m2′

FR3: min z

s.t. eqs 10, 16−25, 28, 30, 32, 34, 36−49 in this paper and eqs 6, 21−24 from RPSMiT.

4. VALID INEQUALITIES FOR THE LOT-SIZING RELAXATIONS It is an effective solution technique to improve the mathematical formulations of MILP models via valid inequalities. With no integer feasible solutions cut off, valid inequalities tighten the feasible spaces of linear programming relaxations of MILP models, so that the linear relaxations become much closer to the convex hulls of MILP feasible solutions. Therefore, valid inequalities can decrease the nodes explored in branchand-bound trees and speed up solution processes. As mentioned in section 1, high-level relaxations are canonical production models which formulate different characteristics of discrete manufacturing systems such as production batches, sequences, costs, capacities, and resources. A variety of highlevel relaxations exist as submodels in the MILP models of lotsizing problems extensively. A great many valid inequalities for high-level relaxations are reported in the literature and show their powerful effect on enhancing the computational performances of MILP models. In this paper, we derive high-level relaxations from the FCNF-based reformulations of the refinery production scheduling problem and develop valid inequalities based on the relevant results in lot-sizing area. 4.1. Lot-Sizing Relaxations from the FCNF-Based Reformulations. The three FCNF-based reformulations all include eqs 16, 20, 21, and 25 in this paper and eq 22 from RPSMiT, and these constraints represent the restrictions on production and deliveries of final products. We observe that a submodel including all of them is highly similar to the MILP formulation of the single-item uncapacitated lot-sizing (LS-U) problem. With the upper and lower bounds on inventories and blending capacities relaxed, an initial high-level relaxation

QIuminyaum, t ≤ QIA um, t ≤ QIumaxyaum, t ∀ u ∈ U , m ∈ Mu , t ∈ T (44) ≤

≤

∀ u ∈ U , m ∈ Mu , t ∈ T (45)

QIum, t = QIA um, t + QIBum, t ∀ u ∈ U , m ∈ Mu , t ∈ T

(49)

Constraint 47 requires that FCCU and HDS must be in the same kind of initial steady operation mode, whereas constraints 48 and 49 represent that HDS would enter any steady operation mode one time interval earlier than FCCU at the subsequent time intervals in the scheduling horizon. The objective function of this reformulation is still the same as that of model BF. Therefore, we obtain the third FCNFbased reformulation denoted FR3 as follows:

s.t. eqs 10, 14, 16−25, 28, 30, 32, 34, and 36−43 in this paper and eqs 6 and 21−24 from RPSMiT. 3.3. Third FCNF-Based Reformulation. On the basis of model FR2, the inflow rates in the first time interval and the subsequent time intervals of any steady state are represented, respectively, by defining new variables. So far, we have focused on the binary variables denoting operating states and have reformulated relevant constraints to describe the FCNF characteristic of operation mode switching explicitly. Now we further consider the close relationship between the inflow rate of production unit and its current operating state, which is built by the production capacity constraints, and we try to expand the scope of reformulation to the modeling of material flows. m m New continuous variables QIAu,t and QIBu,t are defined m corresponding to binary variables yau,t and ybmu,t, respectively, for any u ∈ U, m ∈ Mu, and t ∈ T, QIAmu,t denotes the inflow rate of unit u when t is the first time interval of steady mode m, and QIBmu,t denotes the inflow rate of unit u when steady mode m continues running in time interval t. Then, we can give the following constraints. (1). Capacity Constraints. We focus on the production capacity constraints in model FR2 for reformulation. Equation 14 which involves the inflow rate under any steady state needs to be rewritten as eqs 44 and 45. New constraint 46 builds a relationship between QIAmu,t, QIBmu,t, and QImu,t.

QIumaxybum, t

m2, m2 ′ m2 ′ feedHDS, oi6QIHDS, t + feedHDS, oi6QIAHDS, t

m2 ∈ MHDS

FR2: min z

QIBum, t

(48)

1, m1 ′ yield mFCCU, QI m1, m1 ′ = oi6 FCCU, t

∑ (43)

Since the transition times of FCCU are one time interval longer than those of HDS, HDS would always enter any steady operation mode earlier than FCCU by one time interval, and we have constraint 42. Constraint 43 indicates that HDS and FCCU must be in the same kind of initial steady operation mode. The objective function of this reformulation is still the same as that of model BF. Therefore, we obtain the second FCNFbased reformulation denoted FR2 as follows:

QIuminybum, t

(47)

(46)

(2). Material Balance Constraints. We focus on the material balance constraints for intermediates in model FR2 for reformulation. Specifically, combining the state correlation restrictions on FCCU and HDS, we reformulate the material balance constraint for the intermediate transferred just between 7881


Article


which are the tailored WW-U valid inequalities 52 and 53 for all the FCNF-based reformulations.

denoted R0 is derived from all the FCNF-based reformulations as follows. R0 is a LS-U relaxation where M is a large positive number approaching infinity. INVo , t = INVo , t − 1 +

t

INVo , k − 1 ≥

∑ Q oc ,o,t − Do,t ∀ o ∈ O , t ≥ 2

l=k

oc

INVo , t = INVinio +

l

∑ [(1 − ∑ bo,p)(Ro,l − Shortmaxo,l)] p=k

∀ o ∈ O , 2 ≤ k ≤ t ≤ Tmax

∑ Q oc ,o,t − Do,t ∀ o ∈ O , t = 1 oc

t

l

∑ [(1 − ∑ bo,p)(Ro,l − Shortmaxo,l)]

0 ≤ Ro , t − Do , t ≤ Shortmax o , t ∀ o ∈ O , t ∈ T

INVinio ≥

0 ≤ Q oc , o , t ≤ Mboc , o , t ∀ oc ∈ OC , o ∈ O , t ∈ T

∀ o ∈ O , t ≤ Tmax

l=1

boc , o , t ∈ {0, 1} ∀ oc ∈ OC , o ∈ O , t ∈ T

Pochet and Wolsey32 mentioned that, when developing valid inequalities for lot-sizing relaxations, it is necessary to consider the compromise between the tightness of valid inequalities and their sizes. As another common high-level relaxation, the singleitem Wagner-Whitin uncapacitated lot-sizing (WW-U) model can be obtained by further relaxing the LS-U model. Hence, the valid inequalities of WW-U model must be also valid for the LSU model. As presented by Pochet and Wolsey,44 there are only O(n 2 ) facet-defining valid inequalities for the WW-U formulation as opposed to O(2n) for the LS-U formulation, and the run-time order of separation algorithm is O(n) under the Wagner-Whitin cost condition, which is better than O(n log n) in the general case. Therefore, with the initial LS-U relaxation R0, we try to take advantage of the valid inequalities for WW-U relaxations to develop effective valid inequalities for the refinery production scheduling problem. First, we reformulate relaxation R0. New binary variable bo,t is defined to represent whether final product o is produced in blending operations during time interval t, and new eqs 50, 51 are established.

NUMOC

≤ bo , t ≤

∀ o ∈ O , 2 ≤ k ≤ Tmax , k ≤ l ≤ Tmax

oc

(54)

For any final product o ∈ O, if it is not produced during the time frame [k,l], ∑lp=k bo,p is equal to 0, then eq 54 represents that the required minimal delivery during time interval l can be satisfied by its inventory. By summing eq 54 over the set l = k,...,t, we further have constraint 55, thus “the single-level inventory restrictions” can be indicated by valid inequality 52. t

INVo , k − 1 ≥

(50)

∑ boc ,o,t ∀ o ∈ O , t ∈ T

∑ bo,p)(Ro,l − Shortmaxo,l) p=k

oc

∑oc boc , o , t

(53)

l

INV lo , k − 1 ≥ (1 −

∑ Q oc ,o,t ≤ M max ·NUMOC·bo,t ∀

o ∈ O, t ∈ T

p=1

Also, valid inequalitites 52 and 53 can be viewed as “the singlelevel inventory constraints” for the refinery production scheduling problem. They describe that if no production runs for any final product during a certain time frame, there must be enough inventory on hand at the beginning of the time frame to ensure the necessary delivery of it within the time frame. A detailed explanation is provided as follows. For any o ∈ O and k ≥ 2, let INVlo,k − 1 represent the inventory of final product o during time interval k − 1 which is to satisfy its necessary delivery during time interval l(l ≥ k). Considering the permissible maximum backorder, we have constraint 54 to enforce the restrictions on “single-level inventories of final products”.

INVo , t ≥ 0 ∀ o ∈ O , t ∈ T

M minbo , t ≤

(52)

t

≥

∑ INV lo,k− 1 l=k

l

∑ [(1 − ∑ bo,p)(Ro,l − Shortmaxo,l)] l=k

p=k

(51)

∀ o ∈ O , 2 ≤ k ≤ t ≤ Tmax

In constraint 50, ∑oc Qoc,o,t is not equal to 0 if and only if bo,t is equal to 1, which means that at least one kind of component is blended into final product o during time interval t. In constraint 51, bo,t is equal to 1 if and only if ∑oc boc,o,t is not equal to 0, which means that final product o is produced by blending at least one kind of component during time interval t. It is obvious that constraints 50 and 51 are exactly equivalent to each other, and either of them is able to introduce bo,t. Then, two reformulated LS-U relaxations are derived from the FCNF-based reformulations, both of which are about the production and deliveries of final products. We denote the two high-level relaxations by R1 and R2, respectively, and they are also equivalent to each other: R1 is composed of R0 and eq 50 in this paper; R2 is composed of R0 and eq 51 in this paper. 4.2. Valid Inequalities and the Tightening Formulations. According to the facet-defining valid inequalities for WW-U model reported in Pochet and Wolsey,44 we can construct the same valid inequalities for relaxations R1 and R2,

(55)

For the first time interval (k is equal to 1), the initial inventory levels of final products INVinio are involved in valid inequality 53 to represent “the single-level inventory restrictions”. Furthermore, we can use valid inequalities 52 and 53 to tighten the FCNF-based reformulations. By including the valid inequalities for the two relaxations R1 and R2 in the three reformulations FR1, FR2, and FR3, we propose nine tightening formulations which are listed in Table 1. For a detailed explanation, we take models FR11 and FR12 as examples. FR11 and FR12 only differ in the high-level relaxations for which the valid inequalities are developed. By adding eq 50 into model FR1, we can derive the LS-U relaxation R1 from FR1, obtain valid inequalities 52 and 53 and include them in FR1, and the tightening formulation FR11 is developed. Whereas for FR12, eq 51 is added to model FR1, so that the LS-U relaxation R2 can be derived to obtain the same valid inequalities 52 and 53, which are also included in FR1 to 7882


Article


Unless stated otherwise, the valid inequalities are added dynamically as cutting planes to tighten the linear relaxation submodels at each node of branch-and-bound trees, which is known as the branch-and-cut algorithm. The parameter values involved in the instances are given in the Supporting Information (see Tables S1−S6). 5.1. Comparison between Models BF and m-RPSMiT. The model and solution statistics of BF and m-RPSMiT for instances I4, I6, and I8 are shown in Table 3. The column “0-1 Vars” indicates the number of binary variables, the column “Cont. Vars” indicates the number of continuous variables, and the total number of constraints is given in the column “Cons.”. The number of binary variables in BF is less than a third of that in m-RPSMiT, and the numbers of continuous variables and constraints reduce by about half in BF. It is obvious that the size of model BF is much smaller than that of model m-RPSMiT. Columns “Obj. Value (¥)” and “Opt. Gap (%)” indicate the objective value and the optimality gap of final solution, respectively. The number of nodes and the computational time are given in columns “nodes” and “CPU time (s)”, respectively. As a consequence of reduction in model size, BF decreases the number of nodes and solution times greatly. As the scales of instances get larger, the reductions in computational efforts for BF are more remarkable. For instance I6, the nodes explored by m-RPSMiT are more than 5555 times as many as those by BF, and the CPU time taken by m-RPSMiT is over 706 times as long as that by BF. Since m-RPSMiT cannot obtain the optimal solution of instance I8 within the time limit, the corresponding reductions in nodes and solution time are only lower bounds on the true values. 5.2. Comparison between Model BF and the Nine Tightening Formulations. We solve model BF and the 9 tightening formulations on the 30 instances of I14, I20, and I28. The computational results are comprehensively analyzed as follows. 5.2.1. Average CPU Times and Numbers of Nodes over I14. The 10 instances of I14 are solved to optimality with model BF and the nine tightening formulations. The average computational times and the average numbers of nodes over the 10 instances are compared in Figure 4 and Figure 5, respectively. From Figure 4, we can see that all the tightening formulations perform much better than BF. On average, the CPU time consumed by the nine tightening formulations is less than two-fifths as long as that by model BF. The average CPU times for FR1 and FR3 are both shorter than that for FR2, which implies that the reformulating technique employed in FR2 makes the least effect on speedup. With the inclusion of valid inequalities, the reductions in CPU times consumed by FR2 and FR3 are obvious; FR22 reduces by up to 20% the CPU time for FR. Whereas the CPU times consumed by FR11

Table 1. All the Tightening Formulations denotation

explanation the first FCNF-based reformulation FR1 tightened with VIs eqs 52 and 53 FR1 tightened with VIs eqs 52 and 53 the second FCNF-based reformulation FR2 tightened with VIs eqs 52 and 53 FR2 tightened with VIs eqs 52 and 53 the third FCNF-based reformulation FR3 tightened with VIs eqs 52 and 53 FR3 tightened with VIs eqs 52 and 53

FR1 FR11 FR12 FR2 FR21 FR22 FR3 FR31 FR32

for relaxation R1 for relaxation R2 for relaxation R1 for relaxation R2 for relaxation R1 for relaxation R2

develop the tightening formulation FR12. The tightening models FR21, FR22, FR31, and FR32 are developed in the similar manner.

5. COMPUTATIONAL RESULTS We randomly generate multiple instances and select the feasible ones to validate the effectiveness of all the alternative formulations and valid inequalities. The scheduling horizons are uniformly discretized with 8 h time intervals, and the instances vary in demand orders and time horizons. All the instances involved in the computational experiments are one instance I4, one instance I6, one instance I8, 10 instances of I14, 10 instances of I20, and 10 instances of I28, and we list their numbers of time intervals in Table 2. Also, we specifically Table 2. Summary of Instances in Tests denotation

I4

I6

I8

I14

I20

I28

no. of time intervals

4

6

8

14

20

28

name the 10 instances with 28 time intervals as I28-1, I28-2, I28-3, and so on. In the same manner, each of the instances with 14 time intervals and 20 time intervals are named. First, with the solution results of the three instances I4, I6, and I8, we compare the formulations BF and m-RPSMiT on model size as well as computational effort. Second, we focus on testing the effectiveness of the nine tightening formulations, that is, the three FCNF-based reformulations with/without the inclusion of the valid inequalities are compared with model BF. The test contents consist of computational time, number of nodes explored in branch-and-bound tree, optimality gap which indicates the percentage deviation of final solution from global optimum within the time limit, and integrality gap. We solve the nine tightening formulations and model BF on all the instances I14, I20, and I28. Every 10 instances with the same scheduling horizons only differ in demand orders. All the computational tests are carried out using GAMS 24.0.2/CPLEX 12.5 on a computer with 2 GB of RAM and a 2.50 GHz Intel Core (i5-2340M) processor running on Windows 7. A time limit of 5000 s is set for each instance. Table 3. Model and Solution Statistics of BF and m-RPSMiT Ins.

model

0-1 Vars

Cont. Vars

Cons.

Obj. value (¥)

nodes

CPU time (s)

Opt. Gap (%)

I4

BF m-RPSMiT BF m-RPSMiT BF m-RPSMiT

304 1092 500 1682 696 2272

1565 2976 2301 4418 3037 5860

2705 5904 4317 9032 5957 12160

69858910.4 69858910.4 122015573.6 122015573.6 250170994.6 250170994.6

5 2735 29 161105 515 135316

0.44 9.20 1.48 1045.74 11.90 5000

0 0 0 0 0 31.4

I6 I8

7883


Article


Figure 4. Average CPU times over I14.

Figure 5. Average numbers of nodes over I14.

Figure 6. Average CPU times over I20.

problems with relatively short time horizons (14 time intervals). As shown in Figure 5, all the tightening formulations except FR2 can reduce by at least 20% and at most 65% the average number of nodes explored by BF. On average, more nodes are explored by FR2 whereas less CPU time is consumed when compared with BF. This means that FR2 contributes to more

and FR12 are slightly longer than that by FR1. FR31 and FR32 are the most effective models; FR31 can reduce by up to 65% the average CPU time consumed by BF, and FR32 can reduce by about 64%. This indicates that the combination of the valid inequalities and all the reformulating techniques based on FCNF can accelerate the solution processes of the scheduling 7884


Article


Figure 7. Average numbers of nodes over I20.

efficient LP submodels so that much less computational time is consumed at each node of branch-and-bound trees. Also, this verifies that the reformulating techniques based on FCNF can significantly accelerate the solution of the linear relaxation of MILP submodels. Same as the trend in Figure 4, FR31 and FR32 still perform best; FR1 and FR3 perform much better in reducing nodes than FR2. When the valid inequalities are added to them, the numbers of nodes are further decreased and FR22 reduces by up to 63% the number of nodes explored by FR2. It is obvious that the valid inequalities can help obtain optimal solutions with remarkably less computational efforts. 5.2.2. Average CPU Times and Numbers of Nodes over I20. The 10 instances of I20 are also solved to optimality by the 9 tightening formulations and model BF. We show the average computational times and the average numbers of nodes over the 10 instances in Figure 6 and Figure 7, respectively. In Figure 6, as the scheduling horizons become longer, the acceleration of solutions achieved by the tightening formulations is more remarkable. The CPU time consumed by model BF is more than 6.8 times as long as the average CPU time over the 9 tightening formulations. FR1 reduces by nearly 90% the average computational time for BF, of which the computational performance is better than FR2 and FR3. The average CPU time for FR1 is further reduced with the inclusion of valid inequalities for relaxation R1, whereas the valid inequalities make little effect on accelerating solutions when added to FR2 and FR3. Thus, FR11 and FR1 are the most effective models. In Figure 7, all the tightening formulations except FR21 can reduce by at least 13% and at most 48% the average number of nodes explored by BF. The performance of FR21 is similar to that of FR2 in Figure 5, and the explanation is the same as that over I14. When the scheduling horizons become relatively medium (20 time intervals), FR11 and FR32 perform best in decreasing nodes. FR3 provides less effect on shortening the search routes for global optimums than FR1 and FR2. The valid inequalities work well in cutting down the numbers of nodes when included in FR1 and FR3, except that FR21 explores more nodes on average than FR2. 5.2.3. Average CPU Times and Numbers of Nodes over I28. For each of the 10 instances of I28, the formulations with which the optimal solutions can be obtained within the time limit are listed in Table 4. Only the five instances I28-4, I28-5, I28-6, I28-7, and I28-8 can be solved to optimality with all the 9 tightening formulations and BF, while the other five instances

Table 4. Summary of the Obtainment of Global Optimums for I28 Ins.

models that can reach optimality

I28-1 I28-2

the nine tightening models FR1, FR11, FR12, FR22, FR3, FR31, FR32 FR1, FR11, FR12, FR2, FR3, FR31, FR32 BF, the nine tightening models BF, the nine tightening models BF, the nine tightening models BF, the nine tightening models BF, the nine tightening models FR1, FR11, FR12, FR22, FR3, FR31, FR32 FR1, FR11, FR12, FR21, FR3, FR31, FR32

I28-3 I28-4 I28-5 I28-6 I28-7 I28-8 I28-9 I28-10

models that cannot reach optimality BF BF, FR2, FR21 BF, FR21, FR22

BF, FR2, FR21 BF, FR2, FR22

I28-1, I28-2, I28-3, I28-9, and I28-10 are solved to optimality with some tightening models. We show the average computational time over the five instances I28-4, I28-5, I28-6, I28-7, and I28-8 for each of the 10 formulations in Figure 8. The CPU time consumed by model BF is up to more than 15.3 times the average CPU time over the nine tightening formulations, the speedup effects of which become much more obvious. FR3 performs much better than FR1 and FR2 w.r.t. the average CPU time. The valid inequalities work well in reducing computational times when included in FR1 and FR2, and FR12 further reduces by more than 40% the CPU time consumed by FR1. Among all the 10 formulations, FR3 and FR31 perform best with reductions by more than 95% in the average computational time for BF. We compare the average numbers of nodes over the five instances I28-4 to I28-8 for the 10 formulations in Figure 9. Similar to the trend in Figure 7, all the tightening formulations other than FR21 can reduce by at least 10% and at most 61% the average number of nodes explored by BF, and the explanation about FR21 is the same as that about FR2 over I14. When the scheduling horizons are relatively long (28 time intervals), FR3 and FR31 become the most effective models on reducing nodes. However, the valid inequalities do not always assist in searching for global optimums with less computational efforts. Only FR12 and FR22 contribute to reductions by 40% and 33% in the number of nodes explored by FR1 and FR2, respectively. 7885


Article


Figure 8. Average CPU times over I28-4 to I28-8.

Figure 9. Average numbers of nodes over I28-4 to I28-8.

For the other five instances with 28 time intervals I28-1, I282, I28-3, I28-9, and I28-10, we summarize the computational results about CPU times and optimality gaps of final solutions in Table 5. BF cannot obtain the optimal solutions of all the five

solutions much closer to global optimums within the time limit. The optimality gaps reached by BF are reduced by at least more than 55% for instance I28-9 and at most nearly 90% for instance I28-2. For each of the five instances in Table 5, at least seven tightening formulations can obtain the optimal solution within CPU times much less than the time limit. This demonstrates that the computational performances of most tightening models are steady and excellent. 5.2.4. Closed Integrality Gaps of the Tightening Formulations. In Figure 10, we compare the effects on tightening linear relaxations of MILP models made by the nine tightening formulations. For each of them, we use eq 56 to calculate the percentage of reduction in integrality gap (PRIG) on each instance, and the average percentage of reduction in integrality gap (APRIG) over all the 30 instances of I14, I20, and I28 is illustrated by the corresponding bar in Figure 10. In eq 56, RMILP (RMILPtighten) denotes the best objective value of the linear relaxation model BF (the linear relaxation tightening formulation), and MILP denotes the objective value of the optimal integer solution of the instance.

Table 5. Summary of Solution Results for I28-1 to I28-3 and I28-9 and I28-10

Ins. I28-1 I28-2 I28-3 I28-9 I2810

optimality gap of final solution for BF within 5000s (%)

avg optimality gap of final solutions for tightening models that cannot reach optimality (%)

avg CPU time for tightening models that reach optimality (s)

0.60 1.74 0.61 0.56 0.49

0 0.08 0.08 0.24 0.09

1001.2 2385.1 3356.4 2786.9 2740.5

instances within the time limit, and the second column lists the final gap reached by BF for each of them. The third column lists the average final gaps reached by the corresponding tightening formulations shown in the last column of Table 4. The fourth column lists the average CPU times consumed by the corresponding tightening formulations shown in the second column of Table 4. The instance I28-1 is solved to optimality with all the tightening formulations. For the other four instances, the tightening models can either reach optimality or obtain final

PRIG = (RMILPtighten − RMILP)/(MILP − RMILP) × 100%

(56)

The solution results in Figure 10 show that the linear relaxations of all the tightening formulations are tighter than that of model BF. The FCNF-based reformulating technique employed in FR3 contributes most to tightening and more than 7886


Article


Figure 10. Average percentages of reduction in integrality gaps of the nine tightening formulations.

45% the integrality gap of BF is closed, whereas the effect of valid inequalities is not so obvious. 5.3. Summary. First, compared with m-RPSMiT, model BF has a much smaller size and quite better computational performance. The scheduling problems that m-RPSMiT fails to address in reasonable times can be quickly solved to optimality by BF with rather less efforts. Second, compared with BF, the tightening formulations can reach optimality in much less computational times or obtain smaller optimality gaps within the time limit, and the numbers of nodes can also be cut down in most cases. Integrality gaps of the FCNF-based reformulations are smaller than that of model BF. The valid inequalities can further reduce the numbers of nodes explored by the FCNF-based reformulations. Both the FCNF-based reformulating techniques and the valid inequalities can enhance the computational performances of MILP models remarkably. As the scheduling horizons become longer, the computational efficiency of the tightening formulations are more remarkable. Therefore, with model BF, the FCNF-based reformulations, and the valid inequalities, medium-scale or even large-scale refinery production scheduling problems of industrial relevance can be addressed with significantly increased possibilities.

definitely reduced. Also, developing new valid inequalities involving refinery characteristics is a desirable direction of the future work.

■

ASSOCIATED CONTENT

S Supporting Information *

Complete details of all instances such as operation modes, operational transitions, production costs, yields for outflows and transition times of production units, various upper and lower bounds, inventory and backlogging costs, initial holdup values, and crude oil price. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b00600.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the National Natural Science Foundation of China (Grant No. 61273039 and Grant No. 21276137).

6. CONCLUSIONS A refinery production scheduling problem involving operational transitions is studied in this paper. We propose several alternative MILP formulations. One is the basic formulation with optimized sets and variables, and the others are three FCNF-based reformulations with constraints describing the FCNF characteristic of operation mode switching. In addition, valid inequalities for two lot-sizing relaxations are developed to tighten the FCNF-based reformulations. Numerical experiments show that all the alternative formulations and the valid inequalities can enhance the computational performances of MILP models remarkably, and the effectiveness is especially apparent for large-scale problems with long scheduling horizons. We can conclude that the standard FCNF structure contributes to improving MILP models by closing integrality gaps. It is a valuable choice to formulate models with the standard FCNF structure for scheduling problems in chemical production. Valid inequalities for lot-sizing relaxations have an effect on reducing the number of nodes explored in the branchand-cut algorithm, but the computational time cannot be

■

NOTATION

Indices/Sets

t,t′ ∈ T = set of time intervals within the scheduling horizon u ∈ U = set of production units m,m′ ∈ Mu = set of operation modes of production unit u oi ∈ OI = set of intermediates oc ∈ OC = set of blend components voc ∈ VOC = set of blend components with varied key properties subvoc − s ∈ SUBVOC = set of subcomponents of blend component voc o ∈ O = set of final products p ∈ P = set of key product properties

Parameters

TTm,m u ′ = transition time from operation mode m to m′ of production unit u yieldmu,oi = yield for intermediate oi of production unit u in its steady operation mode m

7887


Article


xm,m u,t ′ (in all the tightening models) = One if operational transition from mode m to m′ starts running on production unit u at time interval t, 0 otherwise boc,o,t = One if blend component oc is processed in blending operation for final product o during time interval t, 0 otherwise bo,t = One if final product o is produced in blending operations during time interval t, 0 otherwise

feedu,oi = proportion of intermediate oi consumed by production unit u yieldmu,oc = yield for blend component oc of production unit u in its steady operation mode m feedu,oc = proportion of blend component oc consumed by production unit u yieldm,m u,oi ′ = yield for intermediate oi of production unit u in its operational transition from operation mode m to m′ yieldm,m u,oc ′ = yield for blend component oc of production unit u in its operational transition from operation mode m to m′ QImin = lower bound on inflow rate for production unit u u QImax = upper bound on inflow rate for production unit u u INVmin oc = lower bound on inventory level in storage tank for blend component oc INVmax oc = upper bound on inventory level in storage tank for blend component oc INVmin o = lower bound on inventory level in storage tank for final product o INVmax = upper bound on inventory level in storage tank for o final product o Tmax = terminal time interval of the scheduling horizon Ro,t = demand for final product o during time interval t Shortmaxo,t = permissible maximum shortage for final product o during time interval t INVinioc = initial holdup value in storage tank for blend component oc INVinio = initial holdup value in storage tank for final product o OpCostmu = production cost of production unit u in its steady operation mode m for per ton of inflow tOpCostm,m u ′ = production cost of production unit u in its operational transition from mode m to m′ for per ton of inflow Mmin = lower bound on amount of blend component processed in blending operation Mmax = upper bound on amount of blend component processed in blending operation NUMOC = upper bound on the number of blend components used to blend gasolines or diesels Crudeprice = price of crude oil per ton α = inventory cost for per ton of blend component and final product β = backlogging cost per ton PROmin o,p = lower bound on value of key property p for final product o PROmax o,p = upper bound on value of key property p for final product o PROoc,p = value of key property p for blend component oc rmin oc,o = lower bound on ratio of blend component oc to blend final product o rmax oc,o = upper bound on ratio of blend component oc to blend final product o

Continuous Nonnegative Variables

■

QImu,t = inflow rate of production unit u when steady operation mode m runs on it during time interval t QIAmu,t = inflow rate of production unit u when steady operation mode m starts running on it at time interval t QIBmu,t = inflow rate of production unit u when it remains steady operation mode m during time interval t after it is in steady mode m during time interval t − 1 QIm,m u,t ′ = inflow rate of production unit u when operational transition from mode m to m′ runs on it during time interval t INVoc,t = inventory level of blend component oc in storage tank during time interval t INVo,t = inventory level of final product o in storage tank during time interval t Qoc,o,t = amount of blend component oc processed in blending operation for final product o during time interval t Do,t = delivery of final product o during time interval t

REFERENCES

(1) Shi, L.; Jiang, Y.; Wang, L.; Huang, D. Refinery Production Scheduling Involving Operational Transitions of Mode Switching under Predictive Control System. Ind. Eng. Chem. Res. 2014, 53 (19), 8155−8170. (2) Lu, W.; Zhu, Y.; Huang, D.; Jiang, Y.; Jin, Y. A New Strategy of Integrated Control and On-Line Optimization on High-Purity Distillation Process. Chin. J. Chem. Eng. 2010, 18 (1), 66−79. (3) Chen, X.; Grossmann, I.; Zheng, L. A Comparative Study of Continuous-Time Models for Scheduling of Crude Oil Operations in Inland Refineries. Comput. Chem. Eng. 2012, 44, 141−167. (4) Wu, N.; Chu, F.; Chu, C.; Zhou, M. Hybrid Petri Net Modeling and Schedulability Analysis of High Fusion Point Oil Transportation Under Tank Grouping Strategy for Crude Oil Operations in Refinery. IEEE Trans. Syst., Man, Cybern. C 2010, 40 (2), 159−175. (5) Wu, N.; Chu, F.; Chu, C.; Zhou, M. Tank Cycling and Scheduling Analysis of High Fusion Point Oil Transportation for Crude Oil Operations in Refinery. Comput. Chem. Eng. 2010, 34 (4), 529−543. (6) Menezes, B. C.; Kelly, J. D.; Grossmann, I. E. Improved SwingCut Modeling for Planning and Scheduling of Oil-Refinery Distillation Units. Ind. Eng. Chem. Res. 2013, 52 (51), 18324−18333. (7) Wu, N.; Chu, F.; Chu, C.; Zhou, M. Short-Term Schedulability Analysis of Crude Oil Operations in Refinery with Oil Residency Time Constraint Using Petri Nets. IEEE Trans. Syst., Man, Cybern. C 2008, 38 (6), 765−778. (8) Jia, Z. Y.; Ierapetritou, M. Efficient Short-Term Scheduling of Refinery Operations Based on a Continuous Time Formulation. Comput. Chem. Eng. 2004, 28 (6−7), 1001−1019. (9) Zhang, N.; Zhu, X. A Novel Modelling and Decomposition Strategy for Overall Refinery Optimization. Comput. Chem. Eng. 2000, 24 (2−7), 1543−1548. (10) Mouret, S.; Grossmann, I. E.; Pestiaux, P. Time Representations and Mathematical Models for Process Scheduling Problems. Comput. Chem. Eng. 2011, 35 (6), 1038−1063. (11) Cao, C.; Gu, X.; Xin, Z. A Data-Driven Rolling-Horizon Online Scheduling Model for Diesel Production of a Real-World Refinery. AIChE J. 2013, 59 (4), 1160−1174.

Binary Variables

ymu,t = One if steady operation mode m runs on production unit u during time interval t, 0 otherwise yamu,t = One if steady operation mode m starts running on production unit u at time interval t, 0 otherwise ybmu,t = One if production unit u remains steady operation mode m during time interval t after it is in steady mode m during time interval t − 1, 0 otherwise xm,m u,t ′ (in BF) = One if operational transition from mode m to m′ runs on production unit u during time interval t, 0 otherwise 7888


Article

Industrial & Engineering Chemistry Research (12) Jiao, Y.; Su, H.; Hou, W.; Liao, Z. A Multiperiod Optimization Model for Hydrogen System Scheduling in Refinery. Ind. Eng. Chem. Res. 2012, 51 (17), 6085−6098. (13) Mouret, S.; Grossmann, I. E.; Pestiaux, P. A Novel Priority-Slot Based Continuous-Time Formulation for Crude-Oil Scheduling Problems. Ind. Eng. Chem. Res. 2009, 48 (18), 8515−8528. (14) Wu, N.; Chu, F.; Chu, C.; Zhou, M. Short-Term Schedulability Analysis of Multiple Distiller Crude Oil Operations in Refinery with Oil Residency Time Constraint. IEEE Trans. Syst., Man, Cybern. C 2009, 39 (1), 1−16. (15) Wu, N.; Chu, C.; Chu, F.; Zhou, M. Schedulability Analysis of Short-Term Scheduling for Crude Oil Operations in Refinery with Oil Residency Time and Charging-Tank-Switch-Overlap Constraints. IEEE Trans. Automat. Sci. Eng. 2011, 8 (1), 190−204. (16) 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 (4), 614−634. (17) Wu, N.; Zhou, M.; Chu, F. A Petri Net-Based Heuristic Algorithm for Realizability of Target Refining Schedule for Oil Refinery. IEEE Trans. Automat. Sci. Eng. 2008, 5 (4), 661−676. (18) Wu, N.; Bai, L.; Zhou, M.; Chu, F.; Mammar, S. A Novel Approach to Optimization of Refining Schedules for Crude Oil Operations in Refinery. IEEE Trans. Syst., Man, Cybern. C 2012, 42 (6), 1042−1053. (19) Castillo, P. A. C.; Mahalec, V.; Kelly, J. D. Inventory Pinch Algorithm for Gasoline Blend Planning. AIChE J. 2013, 59 (10), 3748−3766. (20) Velez, S.; Maravelias, C. T. Advances in Mixed-Integer Programming Methods for Chemical Production Scheduling. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 97−121. (21) Joly, M.; Moro, L. F. L.; Pinto, J. M. Planning and Scheduling for Petroleum Refineries Using Mathematical Programming. Braz. J. Chem. Eng. 2002, 19 (2), 207−228. (22) Hamisu, A. A.; Kabantiok, S.; Wang, M. Refinery Scheduling of Crude Oil Unloading with Tank Inventory Management. Comput. Chem. Eng. 2013, 55, 134−147. (23) Fetterolf, P. C.; Anandalingam, G. A Lagrangian Relaxation Technique for Optimizing Interconnection of Local Area Networks. Oper. Res. 1992, 40 (4), 678−688. (24) Van Vyve, M. Fixed-Charge Transportation on a Path: Linear Programming Formulations. In Proceedings of the 15th Conference on Integer Programming and Combinatorial Optimization (IPCO XV), Armonk, NY, June 15−17, 2011; pp417−29. (25) Eksioglu, S. D.; Romeijn, H. E.; Pardalos, P. M. Cross-Facility Management of Production and Transportation Planning Problem. Comput. Oper. Res. 2006, 33 (11), 3231−3251. (26) Pochet, Y.; Wolsey, L. A. Lot-Size Models with Backlogging: Strong Reformulations and Cutting Planes. Math. Program. 1988, 40 (3), 317−335. (27) Rardin, R. L.; Wolsey, L. A. Valid Inequalities and Projecting the Multi-Commodity Extended Formulation for Uncapacitated Fixed Charge Network Flow Problems. Eur. J. Oper. Res. 1993, 71 (1), 95− 109. (28) Ortega, F.; Wolsey, L. A. A branch-and-cut algorithm uncapacitated, fixed-charge for the single-commodity, network flow problem. Networks 2003, 41 (3), 143−158. (29) Miller, A. J.; Nemhauser, G. L.; Savelsbergh, M. W. P. On the Polyhedral Structure of a Multi-Item Production Planning Model with Setup Times. Math. Program. 2003, 94 (2−3), 375−405. (30) Miller, A. J.; Nemhauser, G. L.; Savelsbergh, M. W. P. A MultiItem Production Planning Model with Setup Times: Algorithms, Reformulations, and Polyhedral Characterizations for a Special Case. Math. Program. 2003, 95 (1), 71−90. (31) Agra, A.; Andersson, H.; Christiansen, M.; Wolsey, L. A Maritime Inventory Routing Problem: Discrete Time Formulations and Valid Inequalities. Networks 2013, 62 (4), 297−314. (32) Pochet, Y.; Wolsey, L. A. Production Planning by Mixed Integer Programming; Springer: New York, 2006.

(33) Sanjeevi, S.; Kianfar, K. Mixed N-Step MIR Inequalities: Facets for the N-Mixing Set. Discrete Optimization 2012, 9 (4), 216−235. (34) Di Summa, M.; Wolsey, L. A. Mixing Sets Linked by Bidirected Paths. Siam Journal on Optimization 2011, 21 (4), 1594−1613. (35) Zhao, Y.; Klabjan, D. A Polyhedral Study of Lot-Sizing with Supplier Selection. Discrete Optimization 2012, 9 (2), 65−76. (36) Constantino, M. A Polyhedral Approach to a Production Planning Problem. Annals of Operations Research 2000, 96, 75−95. (37) Zhang, M.; Kuecuekyavuz, S.; Yaman, H. A Polyhedral Study of Multi-Echelon Lot Sizing with Intermediate Demands. Oper. Res. 2012, 60 (4), 918−935. (38) Gaglioppa, F.; Miller, L. A.; Benjaafar, S. Multitask and Multistage Production Planning and Scheduling for Process Industries. Oper. Res. 2008, 56 (4), 1010−1025. (39) Merchan, A. F.; Velez, S.; Maravelias, C. T. Tightening Methods for Continuous-Time Mixed-Integer Programming Models for Chemical Production Scheduling. AIChE J. 2013, 59 (12), 4461− 4467. (40) Velez, S.; Sundaramoorthy, A.; Maravelias, C. T. Valid Inequalities Based on Demand Propagation for Chemical Production Scheduling MIP Models. AIChE J. 2013, 59 (3), 872−887. (41) Saharidis, G. K. D.; Minoux, M.; Dallery, Y. Scheduling of Loading and Unloading of Crude Oil in a Refinery Using Event-Based Discrete Time Formulation. Comput. Chem. Eng. 2009, 33 (8), 1413− 1426. (42) Shah, N. K.; Ierapetritou, M. G. Short-Term Scheduling of a Large-Scale Oil-Refinery Operations: Incorporating Logistics Details. AIChE J. 2011, 57 (6), 1570−1584. (43) Gothe-Lundgren, M.; Lundgren, J. T.; Persson, J. A. An Optimization Model for Refinery Production Scheduling. Int. J. Prod. Econ. 2002, 78 (3), 255−270. (44) Pochet, Y.; Wolsey, L. A. Polyhedra for Lot-Sizing with WagnerWhitin Costs. Math. Program. 1994, 67 (3), 297−323.

7889


Alternative Formulations and Improvements with Valid Inequalities for

Recommend Documents