Efficient Two-Level Hybrid Algorithm for the Refinery Production

Article pubs.acs.org/IECR

Efficient Two-Level Hybrid Algorithm for the Refinery Production Scheduling Problem Involving Operational Transitions Lu Zhang,† Yongheng Jiang,*,† Xiaoyong Gao,† Dexian Huang,†,‡ and Ling Wang†,‡ †

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: An overall refinery production scheduling problem involving operational transitions was studied by Shi et al. (Ind. Eng. Chem. Res. 2014, 53 (19), 8155−8170), which is an intractable large-scale mixed-integer linear programming (MILP) problem. To deal with this challenge, an efficient two-level hybrid algorithm based on the hierarchy of decisions is proposed. In the outer level, the key techniques of a discrete particle swarm optimization algorithm are designed to search for the operating states assignment of production units. A queue-based solution representation is offered by an elaborate encoding scheme for the outer-level problem so that much fewer discrete decision variables are involved. All the constraints on operating states assignment are represented with the encoding scheme and can be satisfied by all the particles throughout the evolutionary process. In the nested inner level, under the assignment specified in the outer level, the detailed production process is optimized by dual simplex method, which includes the production, consumption, inventories of materials, and deliveries of products. Considering that Zhang et al. (Ind. Eng. Chem. Res. 2015, 54 (32), 7871−7889) studied a similar problem and provided an improved MILP model with much better performance, computational comparisons are made between the two-level hybrid algorithm and a full-space MILP model which is mainly adapted from the work of Zhang et al. The results show that the proposed method takes much less time to obtain near-optimal or even optimal solutions robustly. The advantage of the hybrid algorithm is especially obvious to largescale problems, and better solutions can be obtained in minutes when GAMS cannot reach the optimality within reasonable times.

1. INTRODUCTION Production scheduling plays a key role in the refinery industry for great challenges from the competitive market. Refinery scheduling allocates materials and operations for production units to produce final products, which has been viewed as an important approach to improving economic and social benefits for refineries. The research work on short-term scheduling of refineries has received increasing attention for decades. Various optimization models have been formulated for the crude oil scheduling,1 essential production units,2 or overall refinery optimization.3 Based on the state-task network4 (STN) or the resource-task network5 (RTN) problem representations, the mixedinteger linear programming6,7 (MILP) models or the mixed-integer nonlinear programming8,9 (MINLP) models employed discretetime,10 continuous-time,11 single time grid,12 or unit-specific time grids13 representations. Different solution methods have also been proposed to solve the formulations efficiently. Jia and Ierapetritou14 developed a comprehensive mathematical programming formulation for the scheduling of overall refinery operations. A spatial decomposition strategy was applied to generate three subproblems which were solved separately. Mendez et al.15 presented an iterative MILP-based method to deal with the simultaneous optimization of short-term scheduling and off-line blending problem. A sequential MILP approximation was implemented to replace the solution of a complex MINLP model. Luo and Rong16 proposed a © XXXX American Chemical Society

hierarchical approach for short-term refinery scheduling. The optimization model involved in the upper level, the simulation system involved in the lower level, and the iteration procedure between them were presented. To address a nonconvex MINLP model for a crude oil operations scheduling problem, Mouret et al.17 explored a method involving constraint programming (CP) to generate a tighter linear relaxation for the MINLP model. An efficient Lagrangian decomposition approach was introduced by Mouret et al.18 to solve a large-scale MINLP model. The refinery planning and crude oil scheduling problems were integrated in this model and solved separately by the decomposition approach. For a gasoline blending scheduling problem, Castillo-Castillo et al.19 proposed an efficient multiperiod inventory pinch-based algorithm, the three levels of which corresponded to inventory pinchbased optimization, discrete-time approximate scheduling, and continuous-time detailed scheduling, respectively. With constrained ordinal optimization combined, a bilevel hybrid algorithm was developed by Bai et al.20 for a multipipeline crude oil blending problem, and different evolutionary algorithms were applied in the two levels. To optimize an integrated refinery scheduling problem, Received: February 16, 2016 Revised: May 21, 2016 Accepted: June 23, 2016

A

DOI: 10.1021/acs.iecr.6b00631 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Shah et al.21 first decomposed it into a production unit scheduling subproblem and a final product blending and delivery scheduling subproblem, then applied an iterative heuristic algorithm to obtain feasible solutions. Cerda et al.22 proposed a sequential solution strategy for the crude oil blending and scheduling problem, which consisted of solving a very tight MILP model and a nonlinear programming (NLP) model that employed the MILP-solution as its starting point. In addition, some extensive reviews23−26 on the models and solution methods for refinery production scheduling have been published, in which the advances, key issues, and opportunities presented in the refinery industry were discussed. Undoubtedly, all the research work promotes the considerable development of refinery scheduling. With the unit-level advanced control technology implemented by our research group,27 a finite number of realizable operation modes are defined for production units. Due to the remarkable dynamic characteristic of oil refining process, mode switching would lead to long operational transitions. With both the operation modes and operational transitions considered, Shi et al.28 developed a large-scale discrete-time MILP formulation for an overall refinery production scheduling problem, which can reflect the process dynamics and offer implementable schedules to improve economic benefits. The decision variables contain two parts: the operating states assignment of production units; the detailed production process under a certain assignment. Thus, the introduction of operation mode brings about the hierarchy of decisions. Since a great many binary variables are needed to represent the operating state for each unit during each time interval, the model becomes very computationally expensive. Then, Zhang et al.29 studied a similar problem with some subtle differences only in the lower and upper bounds on blending flow rates, the forms of demand orders, and the permissible maximum backorder. By redefining some variables and parameters, they proposed an improved MILP model which was much smallersized and computationally more efficient. However, as the scheduling horizons increase, the optimal solutions still cannot be obtained within reasonable times. Hierarchical strategy is a valuable choice to address large-scale chemical production scheduling problems, especially when the variables can be divided and optimized in different levels, respectively. The decisions in the outer and inner levels can be disaggregating products into groups and scheduling of each group,30 configuration of the system and complete operating conditions under the configuration,31 allocation of utilities and detailed process of operation,3 planning and scheduling,32 etc. Based on the hierarchical framework, the hybrid algorithm with various subalgorithms used for the two-level optimization has been a research hotspot for decades. Different hybrid algorithms33−36 with meta-heuristics and exact methods employed were developed for large-scale chemical production scheduling problems. As a category of stochastic optimization algorithms, meta-heuristics are effective for problems with large search space and can obtain satisfactory solutions within reasonable times. The popular metaheuristics include genetic algorithm,37 particle swarm optimization (PSO),38 ant colony optimization,39 cuckoo search algorithm,40 etc. However, meta-heuristics cannot guarantee the quality of final solutions obtained and have difficulty in representing complex constraints to find feasible solutions. Exact algorithms are the most common solution methods for mathematical programming models, such as dual simplex method for linear programming (LP) model,30 interior point method for NLP model,2 branch-and-bound method for MILP model,7 outer approximation method for MINLP model,1 etc. All these exact algorithms

have been implemented in various commercial solvers to address large-scale scheduling problems. Exact algorithms can provide the optimality gap indicating the difference of current solution from the optimum, and obtain the optimum if run for a long time. But low-quality final solutions may be inevitable owing to reasonable time limitations. As mentioned in several reviews,25,26 optimization methods with complementary strengths, such as metaheuristics and exact methods, are preferred to be combined in the hybrid algorithm to solve scheduling problems efficiently. In this paper, with the hierarchy of decisions exploited, an efficient two-level hybrid algorithm is developed for the refinery production scheduling problem proposed by Shi et al.28 In the outer level, with the key techniques designed, a discrete particle swarm optimization (DPSO) algorithm is applied to effectively search the discrete solution space of operating states assignment. We present an elaborate queue-based encoding scheme for solution representation, so that much fewer discrete variables are involved in DPSO. All the constraints on the operating states assignment are also represented with the encoding scheme. In addition, we design the population initialization and update operations adaptively to guarantee the feasibility of all the candidates throughout the evolutionary process. In the nested inner level, with the commercial solver CPLEX used, the continuous variables are determined by dual simplex method to optimize the detailed production process under a certain assignment. Specifically, the production, consumption, inventories of materials, and deliveries of products are optimized. The objective value of inner level is taken as the evaluation of the assignment provided by outer level. This paper is organized as follows. In section 2, the refinery production scheduling problem from Shi et al.28 is stated. To adapt to the problem here, several modifications are made on the improved MILP model in Zhang et al.29 and all of them are presented. In section 3, we analyze the hierarchical structure of original problem and propose the two-level hybrid algorithm. The outer- and inner-level optimizations are introduced respectively and the designed DPSO is discussed in detail. In section 4, the experimental results are shown and analyzed. Finally, some conclusions follow in section 5.

2. PROBLEM STATEMENT AND MATHEMATICAL FORMULATION 2.1. Problem Statement. First, we review the flow sheet of a typical refinery production system as depicted in Figure 1, which is the same as that in the work of Shi et al.28 The raw material for the refinery is a kind of blended crude oil with relatively stable properties and sufficient supply. The production units include an atmospheric distillation unit (ATM), a vacuum distillation unit (VDU), a fluidized catalytic cracking unit (FCCU), a hydrogen desulfurization unit (HDS), a etherification unit (ETH), two hydro-refining units (HTU1 and HTU2), a catalytic reformer (RF), and a methyl tert-butyl ether unit (MTBE). Under different operating states, the production units transform the crude oil into various intermediates and finally into several blend components. Then, blender units blend these components into final products, including five types of gasoline (JIV93, JIV97, GIII90, GIII93, and GIII97) and three types of diesel (GIII0, GIIIM10, and GIV0). A storage tank is dedicated for each of the final products and three types of gasoline blend components (C5, rf, and mtbe), while the other two types of gasoline blend components (HG and EG) and three types of diesel blend components (LD, RD1, and RD2) are directly transferred for blending operation. B


Article

Industrial & Engineering Chemistry Research

Figure 1. Flow sheet of the refinery production system.

formulation and notation of P are both given in the Supporting Information (See sections A and B). (1) Since no lower and upper bounds are imposed on blending flow rates, the binary variable boc,o,t and eq 16 from BF are excluded. Also, eqs 22 and 23 from BF are modified as eqs 1 and 2 here, where M denotes a large positive value approaching infinity.

The available information about the refinery production scheduling problem are as follows: operation modes and operational transitions with fixed transition times, lower and upper bounds on inflow rates, yields of outflows for production units under different operating states, maximum and minimum ratios of components in blending, product quality specifications about key component concentration ranges, demands for final products and delivery times of orders, various types of costs, etc. During the scheduling horizon, the overall refinery scheduling determines the sequencing and timing of operation modes on each production unit, production and consumption of each material in operations, inventory levels of storage tanks, and deliveries of final products. The objective is to minimize the total costs of raw material, production, inventory, and backorder penalties. For the detailed problem characteristics, please refer to the work of Shi et al.28 2.2. Mathematical Formulation. A discrete-time MILP model (denoted RPSMiT) was first developed by Shi et al.28 for the refinery production scheduling problem involving operational transitions. By redefining some variables and parameters, Zhang et al.29 then developed an alternative discrete-time MILP model (denoted BF) for a similar problem, in which the nonlinearity was avoided and the variables and elements in sets were simplified. Compared with RPSMiT, BF is much smaller and can be solved with significantly less computational effort. Thus, by modifying BF slightly, an improved discrete-time MILP model (denoted P) can be developed for the problem here. There are two groups of constraints in P: mode switching constraints describe logical relationships in operating state assignments of production units; production process constraints specify limitations of capacity, flow conservation restrictions on materials, deliveries, and quality specifications of final products to describe the detailed production process. Only all the modifications we make on BF are listed as follows. The complete mathematical

0 ≤ Q subvoc‐s, o , t ≤ Myum, t ∀ subvoc‐s ∈ SUBVOC, o ∈ O , t ∈ T , u ∈ U , m ∈ Mu (1) 0 ≤ Q subvoc‐s, o , t ≤

Mxum, t, m ′

∀ subvoc‐s ∈ SUBVOC, o ∈ O , t ∈ T , u ∈ U , m ∈ Mu , m′ ∈ Mu

(2)

(2) Since the demand orders are given in the form of delivery time windows for all kinds of final products and the permissible maximum backorder is no longer considered, eq 25 from BF is replaced by eqs 3−5 here. Also, eqs 20−21 from BF are modified as eqs 6 and 7 here, which are the other constraints involving product delivery. Tl1− 1

∑ Dl ,o,t = 0

∀ l ∈ L, o ∈ O (3)

t=1

Tmax

∑

Dl , o , t = 0

∀ l ∈ L, o ∈ O (4)

t = Tl2 + 1

∑ Dl ,o,t ≤ Rl ,o

∀ l ∈ L, o ∈ O (5)

t

INVo , t = INVo , t − 1 +

∑ Q oc ,o,t − ∑ Dl ,o,t oc

∀ o ∈ O, t ≥ 2 C

l

(6) DOI: 10.1021/acs.iecr.6b00631 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research INVo , t = INVinio +

only continuous variables left to be optimized. By calling innerlevel LP algorithm to solve this submodel to optimality, we obtain the optimal production process c opt|d̅ under the operating states assignment d̅ and the optimal objective value of this LP submodel z(d̅, c opt|d̅ ). Obviously, the latter is the minimum total costs under the corresponding operating states assignment and can be viewed as the evaluation value of d̅. On the other hand, after the inner-level LP submodels of multiple outer-level candidates are optimized one by one, each feasible operating states assignment d̅ gets its evaluation value. By updating the values of decision variables, the discrete optimization algorithm tries to improve the quality of outer-level candidates and reduce the minimum total costs achieved further. Finally, the best operating states assignment dbest with minimum evaluation value z(dbest, copt|dbest) and the corresponding optimal production process copt|dbest are selected as the final solution. For the discrete optimization in the outer level, the combinatorial complexity increases with the number of scheduling time intervals significantly; we select PSO, which is an outstanding meta-heuristic algorithm, as the optimization method of outer level. Considering the continuous nature of solution representation and updating mechanism in standard PSO, we apply DPSO which arises from expanding the standard PSO to discrete optimization, in which the key techniques including encoding scheme, population initialization, and update operations are adaptively designed to optimize the outer-level decision variables effectively. For the continuous optimization in the inner level, since exact LP methods can address LP problems efficiently, we select dual simplex method as the optimization algorithm of inner level, which has been extensively used and successfully implemented in commercial solvers. Thus, based on the twolevel optimization structure, a hybrid algorithm combing DPSO and LP solver is developed. All the mode switching constraints in model P (see eqs S1−S17 in the Supporting Information) are the constraints on operating states assignment, which only involve discrete variables and should be handled in the outer-level optimization. While all the other constraints in model P (see eqs S18−S34 in the Supporting Information) enforce restrictions on the flows of materials, inventories and deliveries, which are included in the inner-level LP submodel. Specifically, DPSO searches the feasible solution space of operating states assignments and generates multiple candidates; to evaluate the minimum total costs under any feasible assignment, each candidate in DPSO is sent to the inner level and the corresponding detailed production process is optimized by calling LP solver to solve an inner-level LP submodel; based on the evaluation values obtained for all feasible operating states assignments, DPSO continues its optimization and performs update operations to generate the next generation of outer-level candidates, which also should be evaluated in the same way as above; before the stopping criterion is met, the update and evaluation of feasible operating states assignments would be repeated so that DPSO adequately evolves the outerlevel candidates and searches the outer-level solution space for optimality. The algorithm would be terminated when the predetermined maximum number of iterations is reached. The above two-level hybrid algorithm is illustrated in Figure 2. 3.2. Outer-Level Optimization. 3.2.1. Overview of PSO. Particle swarm optimization41,42 (PSO) is a population-based evolutionary algorithm derived from the social behavior of bird flocking. Each particle represents a candidate and multiple particles coexist as a swarm. Benefiting from the cooperation

∑ Q oc ,o,t − ∑ Dl ,o,t oc

l

∀ o ∈ O, t = 1

(7)

(3) Since no upper bound is imposed on stockout, the calculation of backorder penalties should be changed and the objective function from BF is modified as follows. min z =

m m,m′ ∑ (∑ QIATM, t Crudeprice + ∑ ∑ QIATM, t Crudeprice) t

+

m

∑ ∑ (∑ t

+

m

u

QIum, t OpCost mu

m

+

m′

∑ ∑ QIum,t,m′tOpCostmu ,m′) m

m′

∑ α(∑ INVoc ,t + ∑ INVo,t) + ∑ ∑ βl(Rl ,o − ∑ Dl ,o,t ) t

oc

o

l

o

t

(8)

(4) Given the total length of scheduling horizon and the fixed transition time of any operational transition, we can calculate the maximum number of time intervals when the operating states of any unit are operational transitions. Since any unit must be in any steady operation mode during the first and the last time intervals, which are enforced by eqs S2 and S4 in the Supporting Information, we introduce the following constraint to impose this upper bound on the number of time intervals for operational transitions, where ⌊x⌋ refers to the maximum integer which is no more than x.

∑ ∑ ∑ t ∈ T m ∈ Mu m ′∈ Mu

⎢T − 1⎥ xum, t, m ′ ≤ TTu·⎢ max ⎥ ⎣ TTu + 1 ⎦

∀u∈U (9)

3. TWO-LEVEL HYBRID ALGORITHM 3.1. Two-Level Optimization Structure Involving DPSO and LP Solver. For this scheduling problem, the decision variables can be divided into two parts. One part represents the operating states assignment of production units, which refers to the sequencing and timing of operation modes on each production unit. The other part represents the detailed production process involving all devices under a certain operating states assignment, which contains the production, consumption, inventories of materials and deliveries of products. Based on the hierarchy of decisions, a two-level optimization structure where the inner level is nested in the outer level can be designed. For any feasible operating states assignment provided by the outer-level algorithm, the detailed production process is decided by calling the inner-level algorithm to minimize the total costs under the corresponding assignment. With the optimal objective value provided by the inner level, the operating states assignment is improved gradually in the outer level to reduce the total costs further. Thus, the optimization on the two levels are interdependent and work together to obtain the optimal schedule for refinery production. In addition, the operating states assignment is represented by all the discrete decision variables while the detailed production process is represented by all the continuous decision variables. The two-level optimization structure can be expressed as follows, where c and d denote all the continuous and discrete variables in model P, respectively; z() denotes the objective function of model P. min z(d , c) ⇔ min(min z(d̅ , c)|d = d̅ )

c ∈ Cd ∈ D

d∈D c∈C

(10)

On the one hand, with any feasible operating states assignment d̅ generated in the outer level, the original MILP model P becomes a computationally efficient LP submodel min z(d̅ , c)|d = d̅ with c∈C

D


Article


The basic procedures for standard PSO are as follows. Step 1: Initialize the velocity and position of each particle in a population randomly. Step 2: Evaluate the objective value of each particle. Set the current position and objective value of each particle as its pbest and the objective value of pbest. For the particle with the best objective value in the initial swarm, set its current position and objective value as gbest and the objective value of gbest. Step 3: Use eqs 11 and 12 to update the velocities and positions of particles. Step 4: Evaluate the objective values of all particles. Step 5: Compare the current objective value of each particle with the objective value of its pbest, if the current objective value is better, update its pbest and the objective value of pbest with its current position and objective value. Step 6: For the pbest with the best objective value in the current swarm, if its objective value is better than the objective value of gbest, update gbest, and the objective value of gbest with its position and objective value. Step 7: If a stopping criterion is reached, output gbest and the objective value of gbest, and the PSO algorithm ends; otherwise, go back to step 3. Due to the continuous nature of solution representation and updating mechanism, the standard PSO is generally used for continuous optimization problems. In order to address combinational optimization problems, based on the framework of PSO, the positions and velocities are discretized to develop discrete particle swarm optimization (DPSO) algorithms.45,46 The encoding scheme, population initialization and update equations are the key techniques of DPSO. The encoding schemes determine the composition of particles for different problems, and the update equations are adaptively formulated to update the integer values of particles. 3.2.2. Design of Key Techniques in DPSO. (1) Encoding Scheme. (a) Representation of Operating States Assignment. For any production unit with more than one operation mode, mode switching is inevitable so that various desired steady modes can run on units to satisfy the final demands. The operating state in any time interval is either a steady operation mode or an operational transition arising from mode switching. An example of the operating states assignment for such a unit is illustrated in Figure 3. This unit has three operation modes numbered 1, 2, and 3, and its fixed transition time is two time intervals. Figure 3 shows that the steady modes 2, 1, and 3 run on the unit sequentially, and the operational transitions from mode 2 to mode 1 and then from mode 1 to mode 3 start at time intervals 2 and 7, respectively. We define “mode queue” and “switch moment queue” to describe the operating states assignment for any unit with several operation modes. Throughout the scheduling horizon, the former is a list of different modes which run on the unit in sequence, while the latter is a list of time intervals when all the transitions start in sequence. Thus, the mode queue is 2−1−3 and the switch moment queue is 2−7 in Figure 3. For any unit, if mode queue and switch moment queue are specified, its operating states assignment are definite, and vice versa. In this way, the operating states assignment of any unit can be represented by its mode queue and switch moment queue together. It is noted that in all the figures of this paper, the serial numbers of operation modes are denoted by integers in bold while those of time intervals are denoted by integers in italic. (b) Constraints Representation. In order to represent a feasible operating states assignment, the following properties should be

Figure 2. Scheme of the two-level hybrid algorithm.

and competition mechanism in population, particles fly in the solution space of optimization problem and look for a desired position with the optimal objective value. PSO has the memory of best positions encountered which contributes to information sharing among particles. Specifically, each particle remembers the best position (denoted pbest) it encountered during the past; each generation of swarm remember the best position (denoted gbest) they encountered during the past; both of them guide the search directions of particles. Thus, the experiences of individual and swarm are combined to balance exploration and exploitation. In PSO, a particle has two attributes: its velocity and its position. Vki = [νki,1, νki,2, ..., νki,D] and Xki = [xki,1, xki,2, ..., xki,D] denote the velocity and the position of the ith particle in the D-dimensional search space at iteration k, respectively. pbestki = [pbestki,1, pbestki,2, ..., pbestki,D] and gbestk = [gbestk1, gbestk2, ..., gbestkD] denote the best position corresponding to the best objective value reached by the ith particle and the whole swarm until iteration k, respectively. Assuming that the size of swarm is NP and the maximum number of iterations is Imax, the velocity of any particle is updated as follows:43 νik, d+ 1 = ωνik, d + c1r1(pbest ik, d − xik, d) + c 2r2(gbestkd − xik, d) ∀ i = 1, 2, ..., NP, d = 1, 2, ..., D , k = 0, 1, ..., Imax − 1 (11)

where ω is inertia weight to control the effect of the previous velocity on the current velocity; c1 and c2 are acceleration coefficients to represent the effects of self-experience and information sharing in population, respectively; r1 and r2 are two independent random values uniformly distributed in [0, 1]. Then, the position of any particle is updated as follows:44 xik, d+ 1 = xik, d + vik, d+ 1 ∀ i = 1, 2, ..., NP, d = 1, 2, ..., D , k = 0, 1, ..., Imax − 1 (12) E


Article


Figure 3. Operating states assignment for an illustrative example.

Figure 4. Unique switch moment queue with the maximum number of mode switching.

If Tmax − 1 is divisible by TTu + 1, when MNMSu is reached, the feasible switch moment queue of unit u is unique; that is, each transition can only start in a unique time interval which can be determined. An example of this case is shown in Figure 4, where the scheduling horizon lasts for ten time intervals and the fixed transition time of unit u is two time intervals, so that Tmax − 1/TTu + 1 is equal to 3. As we can see, the first transition starts in time interval 2 as early as possible; the last transition ends in time interval 9 as late as possible; the start time of any two adjacent transitions is TTu + 1 time intervals apart. Thus, Figure 4 shows the unique feasible operating states assignment and the corresponding switch moment queue of unit u. However, If Tmax − 1 is not divisible by TTu + 1 and the remainder is r, when MNMSu is reached, the feasible switch moment queue of unit u is nonunique; that is, each transition can start in any of several time intervals which can be determined. Based on the analysis of the former case, it is reasonable to infer that the nonuniqueness comes from the various allocations of the extra r time intervals. An example of this case is shown in Figure 5, where the scheduling horizon extends to 11 time intervals and the fixed transition time of unit u is still two time intervals, so that MNMSu is still equal to 3 but the remainder of Tmax − 1/TTu + 1 becomes 1. Compared with the operating states assignment shown in Figure 4, the extra one time interval can be allocated before the start time or after the finish time of any transition. Thus, Figure 5 shows as many as four feasible operating states assignments and the corresponding switch moment queues of unit u. Property 4: There are minimum and maximum values for any moment in switch moment queue. The upper and lower bounds on the start time of any transition can be deduced from mode switching constraints S1−S2, S4, and S9−S12 in the Supporting Information. If each transition starts as early as possible, the value of each moment in switch moment queue is minimum. That is, the first moment is 2; the second moment is 2 + (TTu + 1); the third moment is 2 + 2(TTu + 1) and so on. Similarly, if each transition starts as late as possible, the value of each moment in switch moment queue is maximum. That is, the last moment is Tmax − TTu, the penultimate moment is Tmax − TTu − (TTu + 1), the antepenultimate moment is Tmax − TTu − 2(TTu + 1), and so on. (c) Position Representation of Particles. The position of each particle represents a feasible operating states assignment, which consists of a mode queue subposition and a switch moment queue subposition. In this problem, the seven units ATM, VDU,

satisfied by the mode queue and switch moment queue of any unit: Property 1: Any two adjacent modes in mode queue are different, and the length of mode queue is one more than that of switch moment queue. We can infer from the definition of mode queue that any two adjacent modes must be different. According to mode switching constraints S1, S3, and S5−S10 in the Supporting Information, over the scheduling horizon, if m steady modes in which any two adjacent ones are different run on the unit sequentially, m − 1 transitions must be between these steady modes. In this case, the m operation modes in sequence constitute the mode queue, and the m − 1 time intervals when mode switching start in sequence constitute the switch moment queue. Obviously, there is a correspondence between the mode queue and the switch moment queue: the transition from the first mode to the second mode in mode queue starts at the first moment in switch moment queue; the transition from the second mode to the third mode in mode queue starts at the second moment in switch moment queue, and so on. This can be verified by the example of operating states assignment in Figure 3. Property 2: The difference between any two adjacent moments in switch moment queue is no less than TTu + 1. According to mode switching constraints S9−S12 in the Supporting Information, on the one hand, the duration of any operational transition for unit u is a fixed transition time TTu, on the other hand, only after the new steady mode lasts for at least one time interval can another operational transition starts. Thus, the minimum difference between any two adjacent moments is TTu + 1. Property 3: The maximum length of switch moment queue can be determined with the length of scheduling horizon and the fixed transition time. The maximum length of switch moment queue, that is, the maximum number of mode switching can be inferred from eq 9 and can also be deduced from mode switching constraints S1, S2, S4, and S7−S10 in the Supporting Information. The formula for calculating the maximum number of mode switching for any unit u (denoted MNMSu) is as follows, where ⌊x⌋ refers to the maximum integer which is no more than x. ⎢T − 1⎥ MNMSu = ⎢ max ⎥ ⎣ TTu + 1 ⎦

∀u∈U (13) F


Article


Figure 5. Nonunique switch moment queues with the maximum number of mode switching.

Figure 6. Position representation for an illustrative example.

units FCCU, HDS, and ETH are identical, and so are their switch moment queues. To avoid needless repetition, the mode queue and switch moment queue of FCCU are considered in position

FCCU, HDS, ETH, HTU1, and HTU2 have several operation modes, respectively. According to mode switching constraints S14−S17 in the Supporting Information, the mode queues of the G


Article


randomly until a different value is obtained, so that the difference in any two adjacent modes specified by property 1 can be satisfied. Second, we initialize the switch moment queue subposition of each particle, in which the coding segment of unit u is denoted by smqu,i (i = 1 ∼ MNMSu). For any unit u, only when there are mode switching over the scheduling horizon should the switch moment queue be initialized; the length of switch moment queue, that is, the actual number of mode switching is determined by the corresponding preamble mqu,0 in mode queue subposition and is one less than the length of mode queue, so that property 1 is entirely satisfied; property 3 is also satisfied naturally because the initialized length of mode queue is no more than MNMSu + 1 and that of switch moment queue must be no more than MNMSu; with the initialized actual number of mode switching mqu,0 − 1, the optional time intervals when each transition starts can be calculated according to properties 2 and 4; then from the far front of switch moment queue (denoted smqu,i (i = 1 ∼ (mqu,0 − 1))), each element is initialized as a random integer uniformly distributed among its optional values; while the redundant elements (denoted smqu,i (i = mqu,0 ∼ MNMSu)) are left to be blank. In this way, the complete position of any particle is initialized. With the preamble set for any unit, the initialized actual number of mode switching can be easily obtained as key information during the above initialization process. Also, the population initialization guarantees that the four properties described above are satisfied by all the mode queues and switch moment queues, so that all the particles in the initialized swarm can represent various feasible operating states assignments. The pseudocode of population initialization is as follows.

while those of HDS and ETH are not, and the operating states assignments of HDS and ETH can be specified according to that of FCCU during decoding. Thus, the mode queue subposition and switch moment queue subposition consist of the mode queues and switch moment queues of the remaining five units, respectively. However, the variable actual numbers of mode switching can result in the variable numbers of elements in mode queue and switch moment queue, which would complicate the update operations on particles. To overcome this challenge, fixed-length encodings are applied to both mode queue subposition and switch moment queue subposition. For any unit u, the fixed length of switch moment queue is MNMSu while that of mode queue is MNMSu + 1. If the number of mode switching reaches its upper bound, both the mode queue and the switch moment queue of unit u can be filled; otherwise, a list of values would be arranged from the far left of the two queues while the redundant elements are blank or viewed to be meaningless. In addition, in the mode queue subposition, we set a preamble for each unit and place it before the mode queue. The preamble recodes the current actual length of corresponding mode queue, which contributes to the implementation of both population initialization and update operations. On the one hand, with the preambles initialized and updated randomly, the solution space can be fully searched by particles. On the other hand, owing to the necessary information provided by preambles, only the meaningful elements in mode queues and switch moment queues are operated to update particles, while the ineffective operations on the redundant elements are avoided, so that the efficiency of DPSO is optimized. We use bold squares to denote the preambles in all the figures of this paper. Therefore, the total length of a switch moment queue subposition (denoted TLSMQS) is the sum of MNMSu for the five units, and the total length of a mode queue subposition (denoted TLMQS) is the sum of MNMSu + 1 for the five units and the five preambles, which are listed as follows. TLSMQS = MNMSATM + MNMS VDU + MNMSFCCU + MNMSHTU1 + MNMSHTU2 TLMQS = TLSMQS + 10

(14) (15)

An illustrative example of the position of a particle and the operating states assignment represented by it are shown in Figure 6. (2) Population Initialization. With the encoding scheme for positions described above, only when the variable actual lengths of mode queue and switch moment queue for each unit are known can the two corresponding queues be generated. First, we initialize the mode queue subposition of each particle, in which the coding segment of unit u is denoted by mqu,i (i = 0 ∼ (MNMSu + 1)). For any unit u, the preamble (denoted mqu,0) is initialized as a random integer uniformly distributed in [1, MNMSu + 1], which indicates the initialized actual number of elements in mode queue; with the length given by its preamble mqu,0, the initialization of mode queue (denoted mqu,i (i = 1 ∼ mqu,0)) starts from the front end and each element is a random integer uniformly distributed between 1 and the total number of operation modes (denoted Mu) for unit u, which denotes the serial number of any mode; while the redundant elements (denoted mqu,i (i = (mqu,0 + 1) ∼ (MNMSu + 1))) are left to be blank. Also, if the values of any two adjacent elements in mode queue are same, the latter should be reinitialized

(3) Update Operations. (a) Crossover Operation. The crossover operation involves two particles which are named parent 1 and parent 2, respectively. First, among the five units ATM, VDU, FCCU, HTU1 and HTU2, we randomly select one H


Article


Figure 7. Crossover operation for an illustrative example.

of them as the crossover unit u. Second, we copy the preamble (denoted mq2u,0), mode queue (denoted mq2u,i (i = 1 ∼ mq2u,0)) and switch moment queue (denoted smq2u,i (i = 1 ∼ (mq2u,0 − 1))) of the crossover unit in parent 2 to replace the corresponding 1 1 2 1 , mqu,i (i = 1 ∼ mqu,0 ), and smqu,i elements (denoted mqu,0 2 (i = 1 ∼ (mqu,0 − 1)), respectively) of the crossover unit in parent 1, while the coding segments of the other units in parent 1 remain unchanged, in this way the offspring is generated. Since the preamble, mode queue and switch moment queue of the crossover unit are entirely copied from parent 2 to parent 1 and crossover operation does not destroy the feasibility, the four properties are satisfied by all the mode queues and switch moment queues in offspring, which still represents a feasible operating states assignment. For the crossover unit, owing to the actual lengths of mode queue and switch moment queue obtained from its preamble mq2u,0, only the meaningful elements in parent 2 are copied, so that the ineffective replacement of redundant elements (denoted mq2u,i (i = (mq2u,0 + 1) ∼ (MNMSu + 1)) and smq2u,i (i = (mq2u,0 ∼ MNMTu), respectively) is avoided. The pseudocode of crossover operation is as follows, and this operation is illustrated in Figure 7 where FCCU is the crossover unit.

(b) Mutation Operation. The mutation operation only involves one particle which is regarded as the parent. First, among the five units ATM, VDU, FCCU, HTU1, and HTU2, we randomly select one of them as the mutation unit u. Second, the preamble (denoted mqu,0), mode queue (denoted mqu,i (i = 1 ∼ mqu,0)) and switch moment queue (denoted smqu,i (i = 1 ∼ (mqu,0 − 1))) of the mutation unit in this parent are all reinitialized in the same way as population initialization, while the coding segments of the other units remain unchanged, in this way the offspring is generated. Since the preamble, mode queue and switch moment queue of the mutation unit are entirely reinitialized which guarantees the satisfaction of all the properties, the feasibility of offspring is I


Article


Figure 8. Mutation operation for an illustrative example.

maintained by the mutation operation. For the mutation unit, owing to the actual number of mode switching obtained from its updated preamble, the corresponding mode queue and switch moment queue can be reinitialized correctly in length. The pseudocode of mutation operation is as follows, and this operation is illustrated in Figure 8 where HTU2 is the mutation unit.

Xik + 1 = c1 ⊕ C(Λik + 1, pbestik) ∀ i = 1, 2, ··· , NP, k = 0, 1, ···, Imax − 1

where ⊕ denotes an operation is performed with a certain probability. First, the velocity of particle Vk+1 is calculated by i eq 16. M() denotes the mutation operation which is performed with the probability ω to get Vk+1 i . Specifically, a random number r uniformly distributed in [0, 1] is generated. If r is no more than ω, the mutation operation is performed to obtain a perturbed operating states assignment by Vk+1 = M(Xki ); otherwise the i current operating states assignment remains unchanged by Vk+1 i = Xki . Second, the intermediate item Λk+1 is calculated by eq 17. The i crossover operation denoted C() is performed on the first parent Vk+1 and the second parent gbestk with the probability c2, so that i Λk+1 is either C(Vki + 1, gbestk) or Vk+1 according to a random i i number uniformly distributed in [0, 1]. This reflects the information sharing among particles. Third, with the probability of c1, the crossover operation is applied again while the first and the second parent becomes pbestki in parent becomes Λk+1 i eq 18. Based on a random number uniformly distributed in [0, 1], the updated position of the ith particle Xk+1 is either C(Λk+1 i i , k k+1 pbesti ) or Λi . This reflects the individual experience of any particle. (4) Summary of DPSO Improvements. The major improvements in DPSO made by the key techniques designed are summarized as follows: (a) The operating states assignments of all the units are optimized by DPSO, and the corresponding decision variables are ymu,t and xm,m u,t ′ in model P. If DPSO searches for the optimal values of these variables directly, that is, the position of each particle consists of ymu,t and xm,m u,t ′, the total length of position would be |U||Mu||T| + |U||Mu|2|T|, where U is the set of units, Mu is the set of operation modes of unit u, and T is the set of time intervals over the scheduling horizon. While in our proposed DPSO, with mode queue and switch moment queue encoded, a queue-based encoding scheme is designed and the total length of position is

(c) Update Equations. The position and velocity of the ith particle at iteration k can be updated as follows: V ik + 1 = ω ⊕ M(Xik )

∀ i = 1, 2, ..., NP, k = 0, 1, ..., Imax − 1 (16)

Λik + 1 = c 2 ⊕ C(V ik + 1 , gbestk) ∀ i = 1, 2, ..., NP, k = 0, 1, ..., Imax − 1

(18)

(17) J


Article


Step 2: Initialization Step 2.1: Set the iteration k = 0. Initialize the position of the ith particle Xki , ∀ i = 1, 2, ..., NP Step 2.2: For each particle i, evaluate its objective value zki by calling the dual simplex optimizer of CPLEX to solve the innerlevel submodel, ∀ i = 1, 2, ..., NP Step 2.3: Initialize gbestk as the particle with the best evaluation value in the population Step 2.4: Initialize pbestki as a copy of Xki , ∀ i = 1, 2, ..., NP Step 3: Repeat until the iteration k = Imax Step 3.1: Set the iteration k = k + 1. For each particle i, update its velocity Vki and position Xki by eqs 16−18, ∀ i = 1, 2, ..., NP Step 3.2: For each particle i, evaluate its objective value zki by calling the dual simplex optimizer of CPLEX to solve the innerlevel submodel, ∀ i = 1, 2, ..., NP k Step 3.3: For each particle i, if zki < z(pbestk−1 i ), set pbesti = k Xi , ∀ i = 1, 2, ..., NP k k Step 3.4: Set minpbest = pbest i where pbest i satisfies

2(MNMSATM + MNMSVDU + MNMSFCCU + MNMSHTU1 + MNMSHTU2) + 10, which represents the maximum number of decision variables in DPSO. Since the maximum number of mode switching for any unit is much less than the total number of scheduling time intervals, and only the five units ATM, VDU, FCCU, HTU1, and HTU2 are considered in position, the total length of position becomes much shorter with significantly fewer decision variables involved, and considerable computational effort can be reduced. (b) All the constraints on operating states assignment are considered to design the encoding scheme. Only when the four properties described above are satisfied by the mode queues and switch moment queues of all the units can any particle represent a feasible assignment. Population initialization guarantees the feasibility of each initial candidate. With the update operations including both crossover operation and mutation operation designed, the feasibility of all the candidates can be maintained throughout the evolutionary process. Population initialization and update mechanism contribute to the desired searching behavior of DPSO together, and the search efforts are focused on the feasible solution space of operating states assignments. Therefore, an effective and efficient optimization of DPSO can be achieved. 3.3. Inner-Level Optimization. With the integers provided by the outer-level optimization and the constraints S18−S34 (see the Supporting Information) included, the LP submodel min(min z(d̅ , c)|d = d̅ ) is solved to optimize all the continuous

k

z(pbest i ) = min z(pbestik) ∀ i = 1, 2, ..., NP. If z(gbestk−1) > i

z(minpbest), set gbestk = minpbest Step 4: Output the best-known solution given by gbestk and the optimal solution of its corresponding inner-level submodel

4. COMPUTATIONAL RESULTS With the demand orders generated randomly, several cases in different scales are used to test the proposed two-level hybrid algorithm. The scheduling horizons are discretized uniformly with 8 h time intervals, and the cases vary in time horizons and demand orders. For each case, the numbers of time intervals and orders are listed in Table 1 while the specific demands and delivery times are given in the Supporting Information (See Tables S1−S12).On a CPU Intel Xeon E5−2609 at 2.5 GHz with RAM 32.0 GB, the two-level hybrid algorithm is implemented in Visual C++ 2012. With the CPLEX Concert Technology, the inner-level optimization modeled in C++ calls CPLEX 12.6.1, in which the dual simplex optimizer is selected with appropriate parameter settings. Except for the inner-level solver, we code all the other programs including the DPSO algorithm. Each case is independently solved 20 times by the two-level hybrid algorithm from different random initialization. For comparison, the cases are also solved with the MILP model P by using GAMS 24.2.2/ CPLEX 12.6.0 on the same computer. The time limit of GAMS is

d∈D c∈C

variables. We call the commercial solver CPLEX in which dual simplex method is applied to this submodel. Under a certain operating states assignment provided by outer level, the optimal solution of inner level represents the optimal detailed production process, which contains the production, consumption, inventories of materials and deliveries of products. The optimal objective value is the minimum total costs under the corresponding assignment. 3.4. Complete Two-level Hybrid Algorithm. In summary, the detailed steps of the two-level hybrid algorithm are presented as follows. Step 1: Setting parameters Set parameters of DPSO, such as population size NP, the maximum number of iterations Imax, mutation probability ω, crossover probabilities c1 and c2, etc.

Table 1. Sizes of Model P and Maximum Number of Discrete Variables in DPSO on All Cases model P case 1 2 3 4 5 6 7 8 9 10 11 12

DPSO

no. time intervals

no. orders

no. discrete variables

no. continuous variables

no. constraints

maximum no. discrete variables

reduction of discrete variables (%)

12

1 2 3 2 3 4 3 4 5 3 4 5

680

4424 4520 4616 6728 6872 7016 9128 9320 9512 11384 11624 11864

9460 9484 9492 16660 16684 16708 25444 25468 25492 35788 35812 35836

44

93.53

68

93.61

88

93.92

112

93.89

18

24

30

1064

1448

1832

K


Article


best, worst values, and the standard deviation of objective values obtained by the hybrid algorithm, respectively; for GAMS, columns obj. value (¥) and opt. gap (%) denote the objective value and optimality gap of final solution within the time limit, respectively. The optimality gap is provided by GAMS and refers to the percentage difference of a solution from the best known lower bound. Compared with the final objective values obtained by GAMS, the relative percentage differences of those obtained by the hybrid algorithm are listed in Table 4, where columns avg, best, worst, and sd denote the average, best, worst values, and the standard deviation of relative percentage differences, respectively. The relative percentage difference is calculated by (zhybrid − zGAMS)/zGAMS × 100%, where zhybrid and zGAMS denote the objective values obtained by the hybrid algorithm and GAMS, respectively. The computational times consumed by the two approaches are also listed in Table 4, and the time for the hybrid algorithm refers to the average value over 20 runs. As shown in Table 4, the qualities of solutions obtained by the hybrid algorithm are very close to or even better than those obtained by GAMS. For all the cases, the average relative percentage differences are less than 0.4%, and even the worst relative percentage differences remain below 1%, both of which indicate the robustness of the hybrid algorithm. For the three small-sized cases with only 12 time intervals, GAMS can quickly find the global optima within about 1 min. During the 20 independent runs, the average computational times of the hybrid algorithm are less than 15 s, and the global optima can also be obtained. When a longer scheduling horizon over 18 time intervals is involved, GAMS consumes much more computational times which are up to 4818 s to find the global optima. In addition, there is a wide variation in the time requirements of GAMS ranging from 353 to 4818 s, which reflects the lack of steadiness and is an undesirable character for practical applications. While for the hybrid algorithm, all the global optima can still be obtained in at least one run, and the average computational times are less than 85 s with little variation among the three cases. For the three relatively medium-sized cases with 24 time intervals, the computational times of GAMS are further increased remarkably, and only case 9 can be solved to optimality within 3 h. However, with less than 5 min averagely consumed by the hybrid algorithm, the best objective value obtained for case 8 is as good as the final solution of GAMS, and the best objective value obtained for case 7 is even better with a reduction by about 0.13%. When the three large-sized cases with 30 time intervals are tackled,

set to 10800 s. All the parameter values involved in cases are given in the Supporting Information (See Tables S13−S18). For each case, the size of model P and the maximum number of discrete variables in DPSO are listed in Table 1. It is obvious that the sizes of P enlarge as the scheduling horizons increase. For the cases with the same number of time intervals, the number of discrete variables in P and the maximum number of discrete variables in DPSO remain unchanged, while the numbers of continuous variables and constraints in P increase slightly when more orders are involved. The maximum number of discrete variables in DPSO equals to the fixed total length of position, and the last column of Table 1 lists the relative percentage differences between the number of discrete variables in P and the maximum number of discrete variables in DPSO. As discussed in Summary of DPSO Improvements, compared with encoding the discrete variables of P in particles directly, the proposed encoding scheme has a great effect on shortening the lengths of particles, so that more than 93% decision variables are no longer needed in DPSO for all the cases in different scales. Since the major computational effort of MILP optimization is for optimizing discrete variables, an efficient outer-level optimization with significantly less computational effort can be achieved. All the DPSO parameters are listed in Table 2, which include population size NP, the maximum number of iterations Imax, Table 2. DPSO Parameters case

NP

Imax

ω

c1

c2

1−3 4−6 7−9 10−12

20 50 80 100

50 80 120 130

0.9

0.4

0.4

mutation probability ω, crossover probabilities c1 and c2. The values of ω, c1, and c2 are same in all the cases. For NP and Imax, considering the number of discrete variables increases with the total length of scheduling horizon, we set larger values in cases with longer scheduling horizons so that the expanded outer-level solution spaces can be fully searched by DPSO, while same values are set in cases with identical scheduling horizons. The values of all the parameters are selected by trial-and-error based on their typical values in standard PSO and the guidelines of parameter settings. For all the cases, the solution statistics of the two-level hybrid algorithm and GAMS are summarized in Tables 3 and 4. In Table 3, columns avg, best, worst, and sd denote the average,

Table 3. Final Objective Values for the Two Approaches and Final Optimality Gaps for GAMS hybrid algorithm obj value (¥)

GAMS

case

avg

best

worst

sd

obj. value (¥)

opt. gap (%)

1 2 3 4 5 6 7 8 9 10 11 12

146539624.5 89939917.9 170990765.2 174594210.6 151833794.0 314029059.6 231141648.2 178878502.0 249008340.8 485979631.9 226775925.6 301820817.9

146434913.0 89913283.9 170676548.4 173966601.0 151615859.5 313559098.8 230793992.3 178438999.6 248878756.3 485656559.5 226547898.0 301730303.9

147753893.5 90413963.5 171605320.7 175390135.6 152076699.7 315617024.4 232517329.3 179789455.3 249138331.8 488219487.2 227259407.8 302130361.5

296754.4 108772.8 348992.4 426737.5 160446.5 553592.2 434658.0 293152.0 79425.9 686840.4 177009.3 111865.9

146434913.0 89913283.9 170676548.4 173966601.0 151615859.5 313559098.8 231100965.8 178438999.6 248878356.3 487132467.4 227015719.5 302514400.3

0.000 0.000 0.000 0.000 0.000 0.000 1.098 0.400 0.000 2.029 1.341 1.871

L


Article

Industrial & Engineering Chemistry Research Table 4. Relative Percentage Differences for Hybrid Algorithm and Computational Times for the Two Approaches relative percentage difference for hybrid algorithm (%) case

avg

best

worst

sd

average computational time for hybrid algorithm (s)

computational time for GAMS (s)

1 2 3 4 5 6 7 8 9 10 11 12

0.072 0.030 0.184 0.361 0.144 0.150 0.018 0.246 0.052 −0.237 −0.106 −0.229

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 −0.1328 0.0000 0.0002 −0.3030 −0.2061 −0.2592

0.9007 0.5568 0.5442 0.8183 0.3040 0.6563 0.6129 0.7568 0.1045 0.2231 0.1073 −0.1269

0.203 0.121 0.204 0.245 0.106 0.177 0.188 0.164 0.032 0.141 0.078 0.037

9.58 14.55 9.27 84.35 70.77 67.54 253.24 292.21 279.38 381.30 553.17 471.12

22.50 63.68 47.53 4818.51 412.72 353.70 10800.00 10800.00 6339.32 10800.00 10800.00 10800.00

the variables in different levels. In this way, the computational performance of MILP model can be remarkably enhanced. Second, the following measures can improve the optimization performance of meta-heuristics significantly: (1) the encoding scheme should be designed adaptively and effectively to make the size of individual much smaller, so that much fewer variables need to be optimized; (2) with the constraints considered by the solution representation and evolutionary operations, all the individuals are guaranteed to be feasible, so that considerable computational effort for the detection of infeasibility can be reduced and the feasible solution space can be fully searched.

none of them can be solved to optimality by GAMS within 3 h, which indicates that using GAMS to pursue optimality in reasonable times is unpractical for cases in such a scale. By contrast, for the hybrid algorithm, with less than 10 min averagely consumed, both the best and the average objective values obtained are better than the final solutions of GAMS, and even the worst objective value obtained is better for case 12. Compared with GAMS, the hybrid algorithm can obtain the same or better solutions of all cases except case 9, where the best objective value obtained is extremely close to the global optimum with a relatively percentage difference of 0.0002%. All the improvements in qualities of solutions made by the hybrid algorithm are up to more than 0.3%. In summary, when GAMS can obtain the global optima of small-sized cases in reasonable times, the hybrid algorithm takes less computational times to find near-optimal solutions robustly and is able to obtain the global optima during several independent runs; when GAMS cannot obtain the global optima of large-sized cases within 3 h, the hybrid algorithm takes only minutes to find satisfactory solutions which are better than the final solutions of GAMS. As the scheduling horizons become longer, the computational times of GAMS increase significantly, while those of the hybrid algorithm increase much more slowly and are less than 10 min even for the large-sized cases.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b00631. Complete mathematical formulation and notation of MILP model P and complete details of all instances such as demands and delivery times for orders, operation modes, operational transitions, production costs, yields for outflows and transition times of production units, various upper and lower bounds, inventory and backorder penalty costs, initial holdup values, crude oil prices, etc. (PDF)

■

5. CONCLUSIONS In this paper, by optimizing the operating states assignment and the detailed production process in a nested optimization structure, a two-level hybrid algorithm combining DPSO and LP solver is proposed for a refinery production scheduling problem. In the outer level, with the encoding scheme, population initialization and update operations adaptively designed for DPSO, much fewer discrete variables are involved and all the outerlevel constraints are considered in the solution representation of operating states assignment. In the inner level, under the specified outer-level assignment, the detailed production process represented by all the continuous variables can be efficiently optimized by dual simplex optimizer to evaluate the candidate from outer level. Numerical experiments show the two-level hybrid algorithm is effective and efficient in terms of both solution quality and computational time, and the outstanding performance is especially obvious for large-scale problems. Two constructive conclusions can be drawn from this paper. First, instead of using a full-space method to address the largescale MILP model in long computational time, we can develop a hierarchical framework according to the hierarchy of decisions, and combine appropriate algorithms to respectively optimize

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 (No. 61273039, No.21276137) and the National Science Fund for Distinguished Young Scholars of China (No. 61525304).

■

REFERENCES

(1) 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. (2) Menezes, B. C.; Kelly, J. D.; Grossmann, I. E. Improved Swing-Cut Modeling for Planning and Scheduling of Oil-Refinery Distillation Units. Ind. Eng. Chem. Res. 2013, 52 (51), 18324−18333. (3) Zhang, N.; Zhu, X. X. A novel modelling and decomposition strategy for overall refinery optimization. Comput. Chem. Eng. 2000, 24 (2−7), 1543−1548. M


Article

Industrial & Engineering Chemistry Research (4) Jia, Z. Y.; Ierapetritou, M.; Kelly, J. D. Refinery short-term scheduling using continuous time formulation: Crude-oil operations. Ind. Eng. Chem. Res. 2003, 42 (13), 3085−3097. (5) Castro, P. M. Optimal Scheduling of Pipeline Systems with a Resource-Task Network Continuous-Time Formulation. Ind. Eng. Chem. Res. 2010, 49 (22), 11491−11505. (6) Mouret, S.; Grossmann, I. E.; Pestiaux, P. Time representations and mathematical models for process scheduling problems. Comput. Chem. Eng. 2011, 35 (6), 1038−1063. (7) Hamisu, A. A.; Kabantiok, S.; Wang, M. Refinery scheduling of crude oil unloading with tank inventory management. Comput. Chem. Eng. 2013, 55, 134−147. (8) 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. (9) Moro, L. F. L.; Pinto, J. M. Mixed-integer programming approach for short-term crude oil scheduling. Ind. Eng. Chem. Res. 2004, 43 (1), 85−94. (10) 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. (11) 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. (12) Castro, P. M.; Grossmann, I. E. Global Optimal Scheduling of Crude Oil Blending Operations with RTN Continuous-time and Multiparametric Disaggregation. Ind. Eng. Chem. Res. 2014, 53 (39), 15127−15145. (13) Yadav, S.; Shaik, M. A. Short-Term Scheduling of Refinery Crude Oil Operations. Ind. Eng. Chem. Res. 2012, 51 (27), 9287−9299. (14) 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. (15) 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. (16) Luo, C.; Rong, G. Hierarchical approach for short-term scheduling in refineries. Ind. Eng. Chem. Res. 2007, 46 (11), 3656−3668. (17) Mouret, S.; Grossmann, I. E.; Pestiaux, P. Tightening the Linear Relaxation of a Mixed Integer Nonlinear Program Using Constraint Programming. Integration of Ai and or Techniques in Constraint Programming for Combinatorial Optimization Problems, Proceedings 2009, 5547, 208−222. (18) Mouret, S.; Grossmann, I. E.; Pestiaux, P. A new Lagrangian decomposition approach applied to the integration of refinery planning and crude-oil scheduling. Comput. Chem. Eng. 2011, 35 (12), 2750− 2766. (19) Castillo-Castillo, P. A.; Mahalec, V. Inventory pinch gasoline blend scheduling algorithm combining discrete- and continuous-time models. Comput. Chem. Eng. 2016, 84, 611−626. (20) Bai, L.; Jiang, Y.; Huang, D. A Novel Two-Level Optimization Framework Based on Constrained Ordinal Optimization and Evolutionary Algorithms for Scheduling of Multipipeline Crude Oil Blending. Ind. Eng. Chem. Res. 2012, 51 (26), 9078−9093. (21) Shah, N. K.; Sahay, N.; Ierapetritou, M. G. Efficient Decomposition Approach for Large-Scale Refinery Scheduling. Ind. Eng. Chem. Res. 2015, 54 (41), 9964−9991. (22) Cerda, J.; Pautasso, P. C.; Cafaro, D. C. Efficient Approach for Scheduling Crude Oil Operations in Marine-Access Refineries. Ind. Eng. Chem. Res. 2015, 54 (33), 8219−8238. (23) Shah, N. K.; Li, Z.; Ierapetritou, M. G. Petroleum Refining Operations: Key Issues, Advances, and Opportunities. Ind. Eng. Chem. Res. 2011, 50 (3), 1161−1170. (24) Grossmann, I. E. Advances in mathematical programming models for enterprise-wide optimization. Comput. Chem. Eng. 2012, 47, 2−18. (25) Velez, S.; Maravelias, C. T. Advances in Mixed-Integer Programming Methods for Chemical Production Scheduling. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 97−121.

(26) Harjunkoski, I.; Maravelias, C. T.; Bongers, P.; Castro, P. M.; Engell, S.; Grossmann, I. E.; Hooker, J.; Mendez, C.; Sand, G.; Wassick, J. Scope for industrial applications of production scheduling models and solution methods. Comput. Chem. Eng. 2014, 62, 161−193. (27) Lue, 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. (28) 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. (29) Zhang, L.; Jiang, Y.; Huang, D. Alternative Formulations and Improvements with Valid Inequalities for the Refinery Production Scheduling Problem Involving Operational Transitions. Ind. Eng. Chem. Res. 2015, 54 (32), 7871−7889. (30) Harjunkoski, I.; Grossmann, I. E. A decomposition approach for the scheduling of a steel plant production. Comput. Chem. Eng. 2001, 25 (11−12), 1647−1660. (31) Garcia-Herreros, P.; Gomez, J. M.; et al. Optimization of the Design and Operation of an Extractive Distillation System for the Production of Fuel Grade Ethanol Using Glycerol as Entrainer. Ind. Eng. Chem. Res. 2011, 50 (7), 3977−3985. (32) Erdirik-Dogan, M.; Grossmann, I. E. A decomposition method for the simultaneous planning and scheduling of single-stage continuous multiproduct plants. Ind. Eng. Chem. Res. 2006, 45 (1), 299−315. (33) Tian, W.; Sun, S. Application of Hybrid Simulated Annealing Algorithm of Crude Inventory Management in Refinery. Comput. Eng. Appl. 2005, 41 (26), 220−223. (34) Sand, G.; Till, J.; Tometzki, T.; Urselmann, M.; Engell, S.; Emmerich, M. Engineered versus standard evolutionary algorithms: A case study in batch scheduling with recourse. Comput. Chem. Eng. 2008, 32 (11), 2706−2722. (35) Chen, X.; Li, Z.; Yang, J.; Shao, Z.; Zhu, L. Nested tabu search (TS) and sequential quadratic programming (SQP) method, combined with adaptive model reformulation for heat exchanger network synthesis (HENS). Ind. Eng. Chem. Res. 2008, 47 (7), 2320−2330. (36) You, X.; Rodriguez-Donis, I.; Gerbaud, V. Investigation of Separation Efficiency Indicator for the Optimization of the AcetoneMethanol Extractive Distillation with Water. Ind. Eng. Chem. Res. 2015, 54 (43), 10863−10875. (37) Ramteke, M.; Srinivasan, R. Large-Scale Refinery Crude Oil Scheduling by Integrating Graph Representation and Genetic Algorithm. Ind. Eng. Chem. Res. 2012, 51 (14), 5256−5272. (38) Pan, H.; Wang, L. Blending scheduling under uncertainty based on particle swarm optimization with hypothesis test. Computational Intelligence and Bioinformatics, Pt 3, Proceedings 2006, 4115, 109−120. (39) Chan, F. T. S.; Shekhar, P.; Tiwari, M. K. Dynamic scheduling of oil tankers with splitting of cargo at pickup and delivery locations: a Multi-objective Ant Colony-based approach. Int. J. Prod. Res. 2014, 52 (24), 7436−7453. (40) Yang, H.; Ma, W.; Zhang, X.; Li, H.; Tian, S. A Method for Crude Oil Selection and Blending Optimization Based on Improved Cuckoo Search Algorithm. China Pet. Process. & Pet. Technol. 2014, 16 (4), 70−78. (41) Kennedy, J.; Eberhart, R. Particle swarm optimization. IEEE International Conference on Neural Networks Proceedings 1995, 4, 1942−1948. (42) Eberhart, R. C.; Shi, Y.; Kennedy, J. Swarm Intelligence; Elsevier Science, 2001. (43) Shi, Y.; Eberhart, R. A modified particle swarm optimizer. IEEE International Conference on Evolutionary Computation Proceedings 1998, 69−73. (44) Shi, Y.; Eberhart, R. C. Fuzzy adaptive particle swarm optimization. Proceedings of the IEEE Conference on Evolutionary Computation 2001, 1, 101−106. (45) Lian, Z.; Jiao, B.; Gu, X. A similar particle swarm optimization algorithm for job-shop scheduling to minimize makespan. Applied Mathematics and Computation 2006, 183 (2), 1008−1017. (46) Lian, Z.; Gu, X.; Jiao, B. A similar particle swarm optimization algorithm for permutation flowshop scheduling to minimize makespan. Applied Mathematics and Computation 2006, 175 (1), 773−785. N


Efficient Two-Level Hybrid Algorithm for the Refinery Production

Recommend Documents