A Novel Two-Level Optimization Framework Based on Constrained

Article pubs.acs.org/IECR

A Novel Two-Level Optimization Framework Based on Constrained Ordinal Optimization and Evolutionary Algorithms for Scheduling of Multipipeline Crude Oil Blending Bai Liang,†,‡ Jiang Yongheng,†,‡ and Huang Dexian*,†,‡ †

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

‡

ABSTRACT: This paper introduces a practical scheduling of multipipeline crude oil blending (SMCOB) problem. It is formulated as a complex mixed integer nonlinear programming (MINLP) model, taking the charging sequence and flow rates of oil tanks as decision variables, which cannot be efficiently solved by traditional deterministic methods and solvers. Then, a novel two-level optimization framework based on constrained ordinal optimization (COO) and evolutionary algorithms (EA) is proposed. The solution methodology has two stages based on the main procedures of COO. At the crude evaluation stage, discrete EA are used to search for sequence solutions in the outer level. It evolves the sequence solutions on the basis of their rough evaluation of the feasibility and objective value obtained from the inner level and keeps certain number of probably best sequence solutions. At the accurate evaluation stage, the probably best sequence solutions kept by the crude evaluation stage are accurately evaluated by inner-level continuous EA. The COO approach ensures that some true good enough sequence and flow rate solutions can be obtained from the accurate evaluation stage with high probability. COO-based EA are compared with mixed-coding EA to verify the framework’s efficiency and robustness.

1. INTRODUCTION Nowadays, as a result of the globalization of international crude oil market, lots of refineries, located near seas or rivers, often process many different kinds of crude oil imported from around the world. The significant difference of crude oil property causes severe fluctuation to crude distillation units and secondary processing units and makes advanced control unachievable or ineffective. Moreover, in order to produce more valuable products according to market demand, plant-level optimization has moved to refinery-wide optimization along with the rapid development of modeling and optimization techniques. However, successful implementation of refinery-wide scheduling optimization needs a suitable feed crude oil with desired property as the prerequisite. Scheduling of crude oil blending is a necessary solution to obtain a quality feedstock.1,2 Since the 1970s, relevant literature has been published.3−8 However, most of the literature mainly focused on minimizing scheduling costs rather than property-optimized crude oil blending. Bai et al.9 established a scheduling model for a practical two-pipeline crude oil blending process, where a certain type of crude (named “major crude”) is always transferred in one of the pipelines and other crudes are transferred in another pipeline, in turn. However, in practical production, there is also the problem of multipipeline blending without major crude, which will be further studied in this paper. Scheduling of blending problem is often formulated as mixed integer nonlinear programming (MINLP) models.9−12 With respect to the optimization of MINLP models, some deterministic methods have been available since years ago.13 These methods require the prior step of identification and elimination of nonconvexities and decompose the MINLP models into different kinds of nonlinear programming (NLP) and mixed integer linear programming (MILP) subproblems that have to be iteratively solved. The most common ones are branch and bound, outer-approximation, © 2012 American Chemical Society

generalized benders decomposition, etc. Based on these methods, some commercial MINLP solvers have been developed,14 wherein the MINLPs with special properties can be handled. Evolutionary and stochastic methods, because of their simple concept, easy scheme, and the global search capability independent of gradient information, have been developed rapidly since the last half-century.15 Evolutionary and stochastic optimization methods for MINLPs can be categorized into two classes in terms of treating discrete and continuous variables together or separately. As to the “together” class, some mixedcoding methods were proposed, which include mixed-coding genetic algorithm (GA),16 information-guided genetic algorithm (IGA),17 mixed-integer hybrid differential evolution (MIHDE),18 modified differential evolution (MDE),19 and improved particle swarm optimization (PSO),20 etc. Compared with their original versions, mixed-coding methods often adopt hybrid encoding and operation schemes or some transformation mechanism mapping between the continuous space and discrete space.16−21 As to the “separately” class, some methods based on a two-level optimization framework were presented. Ponce-Ortega et al.22 proposed a two-level approach based on GA to optimize the heat exchanger networks (HENs). The outer level is used to perform the structural optimization, for which a binary GA is used. To evaluate the minimum cost of each structure, each candidate or individual in the binary GA is sent to inner-level NLPs optimization for heat load, where a continuous GA is applied. Based on the cost obtained for each structure, binary GA ranks the HENs structure and produces Received: Revised: Accepted: Published: 9078

September 28, 2011 April 1, 2012 April 17, 2012 April 17, 2012 dx.doi.org/10.1021/ie202224w | Ind. Eng. Chem. Res. 2012, 51, 9078−9093

Industrial & Engineering Chemistry Research

Article

Figure 1. Flow-sheet of the practical P-pipeline crude oil blending process.

Figure 2. Gantt chart of SMCOB problem with P pipelines and I tanks.

new structures until algorithm converges to optimal solution. Thereafter, Chen et al.,23 Khorasany and Fesanghary,24 and Gorji-Bandpy et al.25 successively presented similar two-level solving strategies for synthesis of cost-optimal HENs, the difference of which is the subalgorithms they used in the outer and inner level optimization within the two-level framework. Specifically, Chen et al.23 combined tabu search (TS) and sequential quadratic programming (SQP), and Khorasany and Fesanghary24 combined harmony search (HS) and SQP, while Gorji-Bandpy et al.25 combined GA and SQP. In 2011, GarciaHerreros and Gomez26 solved a MINLP model maximizing an economic criterion for the industrial production of bioethanol with a two-level approach, where SA is used to optimize design variables in the outer level and interior point algorithm is applied to optimize operating variables in the inner level. Being “time-consuming” is the main challenge for the above two-level methodologies when the problem size is increasing. Most of the two-level approaches adopt stochastic−deterministic strategy,23−26 and only a few use stochastic−stochastic strategy.22

The reason is that, for the optimization of inner-level submodels, if the gradient exists and a satisfied solution can be obtained, deterministic methods converge much faster than stochastic methods. For discrete event dynamic system (DEDS) optimization problems, ordinal optimization (OO), put forward by Professor Yu-Chi Ho,27 suggests using a crude model to estimate the performance of solution candidates roughly but fast and then seek the good enough ones with high probability. The OObased optimization techniques will cut down computation without significant loss of precision. They are effective tools for DEDS optimization problems. Derived from OO, COO28 further suggests using a feasibility model with some mistake to fast screen out the feasible solutions before crude-model evaluation, which rightly reduces redundant computation on infeasible solutions and enables COO to efficiently handle constrained optimization problems. The COO approach has been successfully applied in complex and time-consuming planning optimization problems of a remanufacturing system.28,29 9079

dx.doi.org/10.1021/ie202224w | Ind. Eng. Chem. Res. 2012, 51, 9078−9093


Article

stored in each tank; (4) amount of crude stored in each tank; (5) trace elements concentration of each blending crude and the concentration limitations in feedstock; (6) process load of CDU or flow rate in pipelinemix ; (7) flow rate limitations in pipelines; (8) yield of light component limitations in feedstock; (9) due emptying time for each tank; (10) time horizon limit. 2.2. Mathematical Formulation. The current research on crude oil scheduling mainly focuses on reducing operating costs rather than feedstock property optimization. Although many studies refer to the issue of crude oil blending according to concentration of some trace elements, such as sulfur and heavy metals, these trace elements do not have direct reflection on product yield and operational stability. In this paper, TBP data (yield of narrow fractions for full range distillation) are used as the main quality indexes for crude oil evaluation. Based on these data, a scheduling model is established. Before the formulation, the following assumptions are made regarding the refinery operations: (1) The charging tanks, which are feeding through pipelines, are not emptied exactly at the same time. Then, when a real oil tank is emptied and finishes feeding, one scheduling time slot is defined; (2) The changeover time for tank switch is neglected. (3) As to the evaluating index of crude oil property, the TBP data are taken as the main index, which is going to be kept stable, while the trace element content is taken as the assistant index, which is going to be kept below the certain limits. (4) Perfect mixing and linear blending relation are assumed. Due to the assumptions, tanks switch one by one. Moreover, as expressed by operation rules and shown in Figure 2, the SMCOB process is different from the scheduling of batch process, as in gasoline blending, in which exist production breaks. In SMCOB, the blending status, namely, the types of crude charging through pipelines, will not change until the next tank switch happens; thus, only one scheduling period (time slot) is needed between two switches. Therefore, when a charging tank is emptied and switched with a minimal blending period, a scheduling period or a time slot is defined. Then, a total of T = Ib + Is = Iv + Is time slots are considered. In the following, we take two groups of decision variables to model and optimize the SMCOB problem as a two-level problem, where the outer level searches the sequence field of the charging tanks and the inner level optimizes the blending flow rates under the certain sequence:

To the best of our knowledge, there is not much work dealing with the property-optimized SMCOB problem,30,31 which is a complex optimization problem rising in the practice. In this paper, we model the SMCOB problem as a twolevel problem, where the outer level searches the sequence field of the charging tanks and the inner level optimizes the blending flow rates under the certain sequence. The COO idea is introduced to solve the two-level problem by taking blending flow rates optimizing as the evaluation of the sequence. The rest of this paper is organized in the following way. In section 2, the SMCOB problem is described and the mathematical formulation is established. In section 3, the basics of EA and COO are briefly reviewed. In section 4, a novel solution methodology is proposed for the SMCOB model. In section 5, without loss of generality, as the outstanding representatives of EA, DE, and GA are in turn embedded into the COO-based framework for solving a real-case SMCOB model and the results are compared with their mixed-coding solving strategies from the literature. In section 6, we make a brief conclusion and point out the future directions.

2. SCHEDULING OF MULTIPIPELINE CRUDE OIL BLENDING PROBLEM 2.1. Multipipeline Crude Oil Blending Process. In this paper, the SMCOB problem is to optimize the tanks charging schedule for crude distillation unit (CDU).4 As shown in Figure 1, a multipipeline crude oil blending system is composed of several charging tanks, P transfer pipelines, a mixing pipeline, a flow control system, and a CDU. Several types of crude oil, each of which has a certain amount of reserves, are stored in oil tanks. The tanks charge through P pipelines in a certain sequence, and then, oils are mixed in pipelinemix to charge CDU. The flow control system will keep the bending flow rates (ratios) stable. The Gantt chart of SMCOB problem with total P pipelines and I tanks is illustrated Figure 2. The variables above the lines denote the tanks' allocation to pipelines in different time slots. The full lines represent the charging process of real tanks, whereas the long dashes represent the charging process of virtual tanks which store the crude oil to arrive in the future of the scheduling period and are added at the end of scheduling. The scheduling task is to study how to allocate tanks to charge through pipelines properly, so as to achieve the objective of providing CDU with a mixed feedstock of qualified and stable property. During the scheduling process, some operating rules have to be followed: (1) Each charging tank can charge, at most, one transfer pipeline at each time slot. (2) Each transfer pipeline must be connected to one and only one oil tank at each time slot. (3) The oil tanks should feed pipelines one by one in a certain sequence, and the next charging tank cannot deliver oil before the preceding one is emptied. (4) Pipelines should be charged all the time during the scheduling horizon. In other words, the state of “empty pipeline” is not allowed. (5) Each charging tank should be emptied before a specified due date. (6) The crude oils to arrive in the future scheduling periods are considered. Moreover, data available from the refinery will be the following: (1) yield of narrow fractions for blending crudes and desired mixed crude; (2) optimization weights on different narrow fractions in the objective function; (3) types of crude

(1) s = (s1, ..., st, ... sIs; vIs+1, ..., vt, ..., vT), the blending sequence of charging tanks (2) f = f p,t, p = 1, ..., P; t = 1, ..., T, the corresponding blending flow rate in each pipeline and time slot. where st ∈ i = Ib + 1, ..., Ib + Is represents the scheduling tank beginning to participate in blending at the start of time slot t, t = 1, ..., Is; vt ∈ i = Ib + Is + 1, ..., Ib + Is + Iv represents the virtual tank beginning to charge at the start of slot t, t = Is + 1, ..., T; f p,t represents the blending flow rate in pipeline p during time slot t, p = 1, ..., P, t = 1, ..., T. The model with sequence and flow rate variables describes the blending process more intuitionisticly and characteristically and also avoid many operation constraints which are necessary in the model with binary variables representation. (a). Objective Function. In this model, TBP data are taken as the main quality index for evaluating crude oil property. des ηmix specify the yield of the lth narrow fraction of mixed l,t , ηl feedstock in time slot t and that of desired feedstock, respectively. On the right side of eq 1, the first term in 9080



Article

brackets is used to limit the yield change in lth narrow fraction of mixed feedstock between the scheduling time slot and its succeeding one, while the second term is used to keep the yield in lth narrow fraction of mixed feedstock around the objective ηdes l . The sum of the two terms, namely, the feedstock property deviations in different narrow fractions, are weighted by wl. Then, the objective of the function formulated as eq 1 is to reduce the total property deviation Δη. Because it is more important to keep the property of light components stable, we set larger wl for smaller l, namely, more weight on lighter fractions. Equation 2 is a blending equation by linear relation, where ηip,t,o,l denotes the yield of the lth narrow fraction of crude o stored in tank ip,t and ip,t denotes the certain tank that charges through pipeline p in time slot t. The relation between the composition of feedstock and the product yield is approximated by linear relation, considering that the CDUs are not so sensitive. L

∑

P

∑ p = 1 ηi

(i , o) ∈ IO;

t = Is + 1, ..., T (7)

τp , t = Vip,t , t − 1/f p , t

ip , t ∈ i = 1, ..., I

p = 1, ..., P

t = 1, ..., T

(8)

p,t , o , l

·f p, t

(1)

S1 = 0

l = 1, ..., L ;

t = 1, ..., T

p=i

i = 1, ..., Ib

ip , t = st

st ∈ i = Ib + 1, ..., Ib + Is

(2)

t=1

Fip,t , t = f p , t (Et − St )

t=1

Vi , t = Vi , t − 1

p = 1, ..., P (13)

i = {1, ..., I } − {ip , t , p = 1, ..., P} (14)

(f). Flow Rate Constraints. Constraint 15 indicates the blending flow rate in each pipeline is subject to the upper and lower charging limits. Constraint 16 indicates the flow rate of mixed feedstock is specified by throughput of CDU. f pmin ≤ f p , t ≤ f pmax

p = 1, ..., P

t = 1, ..., T

(15)

P

f mix,min ≤

(5)

t = 2, ..., Is

ip , t ∈ 1, ..., I

t = 1, ..., T

st ∈ i = Ib + 1, ...Ib + Is

p = 1, ..., P

p = 1, ..., P

t = 1, ..., T

∑ fp,t p=1

if Vip,t−1, t − 1 = 0

ip , t ∈ i = 1, ..., I

(12)

Vip,t , t = Vip,t , t − 1 − Fip,t , t

ip , t − 1 ∈ i = 1, ..., I

t = 2, ..., T

(11)

t = 1, ..., T

Constraint 5 expresses that if tank ip,t−1 is not emptied at the end of time slot t − 1, t = 2, ..., T, it will go on charging through pipeline p in time slot t. Constraint 6 indicates that if tank ip,t−1 is emptied at the end of time slot t − 1, t = 2, ..., Is, the scheduling tank st ∈ i = Ib + 1, ..., Ib + Is, in accordance with blending sequence (s1, ..., st, ..., sIs), will begin to charge through pipeline p in slot t. Similarly, constraint 7 indicates that if tank ip,t−1 is emptied at the end of time slot t − 1, t = Is + 1, ..., T, virtual tanks vt ∈ i = Ib + Is + 1, ..., I are added to blend one by one according to sequence (vIs+1, ..., vt, ..., vT) until all the real tanks finish blending. if Vip,t−1, t − 1 ≠ 0

t = 1, ..., T − 1

(e). Mass Balance Constraints. Equation 12 indicates that Fip,t,t, the volume of crude tank ip,t chargeing through pipeline p during slot t, can be calculated by flow rates f p,t and time slot duration (Et − St). Then, Vip,t,t, the volume of crude stored in tank ip,t at the end of time slot t, can be further calculated by eq 13; the volume of other tanks remain unchanged, which is expressed by eq 14.

(3)

p=P

(10)

St + 1 = Et

(4)

ip , t − 1 ∈ i = 1, ...I

(9)

(d). Time Slot Sequencing Constraints. Based on the assumptions and operating rules, once a tank is emptied and switched, a time slot is defined. Then, a new tank takes over the emptied tank to charge and a new time slot starts. Constraint 11 establishes a sequence between consecutive time slot t and t + 1.

ip , t ∈ i = 1, ..., I ;

ip , t = i

p = 1, ..., P

t = 1, 2, ..., T

p

(b). Tanks Allocation Constraints. The time point when a charging tank is just emptied can be defined as the beginning of the whole scheduling. Without loss of generality, as shown in Figure 2, it is assumed that the tank connected to pipeline P just finished charging and initial tanks i = 1, ..., Ib (Ib = P − 1) are, in turn, charging through pipelines p = 1, ...,P − 1 with remaining storage at the beginning of scheduling, which is expressed by eqs 3 and 4.

ip , t = st

p = 1, ..., P

(c). Period Duration Constraints. Equation 8 indicates that τp,t, the tank charging duration in pipeline p and slot t, can be calculated with Vip,t,t − 1, the volumetric storage of tank ip,t remaining from time slot t − 1 and f p,t, the charging flow rates in pipeline p and time slot t. Moreover, based on the assumption (1), when a tank finishes blending, one time slot is defined. Then, it can be inferred that the duration of slot t, namely (Et − St), is the minimum τp,t with respect to p, which is expressed in eq 9.

Et − St = min τp , t

− ηl mix )2 ] ,t

P ∑p=1 fp,t

ip , t = ip , t − 1

vt ∈ i = Ib + Is + 1, ..., I

T

l=1 t=1 T−1 + (ηl mix ,t+1 t=1

=

if Vip,t−1, t − 1 = 0

ip , t − 1 ∈ i = 1, ..., I

∑ wl[∑ (ηlmix − ηldes)2 ,t

min Δη =

ηl mix ,t

ip , t = vt

≤ f mix,max

t = 1, ..., T (16)

(g). Trace Elements Concentration Constraints. High concentration of some trace elements, such as sulfur, nitrogen, nickel, etc., will cause equipment corrosion and catalyst

(6) 9081



Article

and operator representation. Generally, the basic algorithm procedures for standard EA are as follows: Step 1: Set parameters adopted by EA, such as population size, probabilities of crossover and mutation, stopping criterion, etc. Step 2: Represent the decision variables by individuals of EA using effective coding scheme. For the sake of facilitating the expression, floating point coding is often used for continuous variables and binary coding or sequence coding is usually adopted for discrete variables. Step 3: Randomly initialize the population within the high and low limits of decision variables. Step 4: Calculate the fitness of all the individuals in the population. A constraint handling mechanism, such as penalty function method, is often used to change the constrained optimization problems into unconstrained ones before calculation. Step 5: Update the best so far individual and corresponding best so far objective value. Step 6: A selection operator is used to replace some individuals in the current population with the children by fitness comparison. Step 7: Create the children from parents through evolutionary operators, such as crossover, mutation, etc. The parents are certain individuals chosen to reproduce. Step 8: If a stopping criterion is met, output the best so far individual and corresponding best so far objective value; otherwise, go back to Step 4. 3.2. Core Ideas of COO. Before introducing COO, we review OO first.29 OO is originally applied to stochastic simulation optimization in discrete event dynamic systems (DEDS), which can be described as following model:

poisoning. Constraint 17 indicates the concentration limitations of some trace elements in mixed feedstock should satisfy the requirements. Specifically, Bip,t,o,k represents the content of trace element k in crude o, which is stored in tank ip,t, Bmix,max k represents the maximum content limit of trace element k in mixed feedstock. P

P

∑ Bi

f ≤ Bkmix,max ∑ f p , t p,t , o , k p , t

p=1

ip , t ∈ i = 1, ..., I

p=1

(i , o) ∈ IO

k = 1, ..., K

t = 1, ..., T

(17)

(h). Yield of Light Components Constraints. As is mentioned in ref 7, if the mixed feedstock is too heavy, there will not be enough products in the distillation tower top section to produce a proper internal reflux flow rate, or if the crude is too light, it will cause difficulties in pressure control. Therefore, constraint 18 is set to limit the yield of light distillate of mixed feedstock. ηlgt,max and ηlgt,min denote the upper and lower yield limits of light components, Llgt is maximum number of light narrow fractions, and Llgt < L. Llgt

η

lgt,min

≤

≤ η lgt,max ∑ ηlmix ,t

t = 1, ..., T (18)

l=1

(i). Tank Emptying Due Date Constraints. It can be noted that, due to assumption (1), Et, the ending time of time slot t, also denotes the emptying time of the tank which finishes charging and is switched at the end of time slot t. Constraint 20 expresses that, according to production plan, each tank should be emptied out before a due date Hdue ut so as to receive incoming new crudes. Indices variable ut ∈ i = 1,..., Ib + Is denotes the tank finishing delivering in period t, and it is determined by eq 19. Moreover, constraint 21 indicates the whole scheduling should end before the scheduling horizon H. ut = i p , t

if Vip,t , t = 0

ip , t ∈ i = 1, ..., I

ut ∈ i = 1, ..., Ib + Is

ET ≤ H

θi ∈Θ

p = 1, ..., P

t = 1, ..., T Et ≤ Hudue t

min J(θi) ≡ min E[L(θi , ξi)] ≈ min

(20) (21)

If we nominate the search space of the sequence variable s and flow rate variable f as S and F, the SMCOB optimization problem can be shorted as following formula:

min J(f , s)

f ∈ F, s ∈ S

θi ∈Θ

Ni

∑ L(θi , ξij) j=1

(23)

where θi ∈ Θ represents the various system design parameters; Θ is the design space for θi; ξi represents all the randomness that may occur in the system when evaluating θi; L is a function of θi that defines the performance metric of the system in which we are interested. The expectation is taken with respect to the distribution of all the randomness, ξi. Computationally, it is calculated by sample mean. It is well-known that the accurate estimate of J(θi), namely the accurate evaluation of θi, needs running the simulation model sufficient times (even Ni→∞). Since each sample run of the simulation model of a complex system may consume considerable time, accurate evaluation may cause a heavy computation burden. Therefore, sample mean with finite Ni is used as a crude model to give θi a crude evaluation, and the less Ni is, the cruder the evaluation value is. Based on this concept, the OO approach on the optimization of DEDS has two stages: the crude evaluation stage and the accurate evaluation stage. At the crude evaluation stage, the major task applying OO is using the fast-computing but crude model to construct the crudely good set U, which includes several probably good θi selected based on the crude evaluation values. The size of U is far less than that of Θ, namely |U| ≪ |Θ|. Then, at the accurate evaluation stage, sample mean with sufficient simulation runs are used as accurate model to evaluate each θi ∈ U. OO theory will ensure that there are some true good enough solutions of Θ exist in U with high probability as

(19)

t = 1, ..., T

θi ∈Θ

1 Ni

(22)

that is, the following variables are going to be optimized: (1) s = (s1, ..., st, ... sIs; vIs+1,...,vt,..., vT), the blending sequence of charging tanks (2) f = f p,t, p = 1,...,P; t = 1,...,T, the corresponding blending flow rates in each pipeline and time slot.

3. PRINCIPLE OF EA AND COO 3.1. Basics of EA. EA take advantage of population evolution mechanism to solve the optimization problems through collaboration and competition with individuals in the population. This evolution process continues until a termination criterion is met. All the EA basically share a common set of underlying assumptions and the evolution of individual structures, but they differ in the offspring generation strategy 9082



Article

model with only flow rate variable f ∈F left to optimize. This can be formulated as follows:

long as the crude model has proper accuracy. Expressed in mathematical form, it is defined as

⎛ ⎞ min J(f , s) ⇔ min⎜min J(f , s ̅ )| s ̅ = s⎟ ⎠ s∈S ⎝ f ∈F

Construct U based on the crude evaluation value s.t. AP = Pro[|U ∩ G| ≥ 1] ≥ PA

(24)

Equation 26 shows if the submodel minf ∈ F J( f, s)| is ̅ s=s ̅ solved optimally by continuous-optimization solvers in the inner level, the optimal flow rate solution with respect to s, namely f opt|s, is obtained and the corresponding optimal value can be considered as the accurate evaluation value of s ∈ S, namely eopt|s. Conversely, when each s ∈ S gets an accurate evaluation value, the global best sequence solution sopt ∈ S with minimum value can be picked up by discrete optimization solvers in the outer level and the global best flow rate solution f opt|sopt is also obtained correspondingly. 4.2. Embedding of EA and COO. As reviewed in the introduction, when the two-level optimization structure is applied for MINLPs, stochastic−deterministic or stochastic− stochastic strategy is usually adopted. In our SMCOB model, the inner-level sub model minf ∈ F J(f,s)| is still complex to ̅ s=s ̅ solve with deterministic methods and solvers. Moreover, when a large number of tanks are prepared for scheduling, there is also high combination level in the outer-level discrete optimization. Therefore, in this paper, evolutionary algorithms, because of their simple concept and easy scheme, are selected as both of the discrete-optimization solver in outer level (denoted as EAs) and the continuous-optimization solver in inner level (denoted as EAf); that is, stochastic−stochastic strategy is used. However, according to stochastic optimization theory, minf∈F J( f, s)|̅ s=s can be solved optimally with probability 1; that ̅ is, the optimal flow rate variables f opt|s is surely obtained or s is accurately evaluated, if the algorithm parameters are properly set and iterations are infinite. Therefore, similar to the case in DEDS, there is also randomness in EA optimization and accurate evaluation of s is also computation-costing. To improve the optimization efficiency and robustness, the idea of COO28 is introduced here to seek a true good enough sequence solution sg ∈ S with small computation and high probability. Then, two evaluation stages, namely, the crude evaluation stage and the accurate evaluation stage, will be successively adopted in our methodology for SMCOB model. At crude evaluation stage, the main task is

where, G is the true good enough set, usually the top-|G| good solutions in the search space Θ. AP is called alignment probability, which is the probability representing that there are truly at least one good enough solution in U. PA is a constant high probability value (i.e. 0.95). Furthermore, if constraints are considered, they are defined as the following eq 25, 1 s.t. gm(θi) ≡ E[Lm(θi , ξi)] ≈ Ni , m m = 1, ..., M

(26)

f ∈ F, s ∈ S

Ni , m

∑ Lm(θi , ξij) ≤ 0 j=1

(25)

where M is the number of constraints. Similar to the objective function, it is also time-consuming to precisely determine the feasibility of θi beforehand, and too much needless computation will be used on infeasible designs. Therefore, Guan28 proposed a feasibility model-based approach, called constrained ordinal optimization. The core idea of COO is to first use a feasibility model, which is fast but probably with some mistakes, to roughly estimate the feasibility of θi. For instance, the feasibility model in DEDS can be sample mean with much less Ni,m. Then, the infeasible solutions can be quickly screened out before the crude evaluation process so as to further save computation. The illustration of COO on the DEDS optimization is showed in Figure 3.

finding estimated top‐N g ̂ sequence solutions s g (̂ j), j = 1, 2, ..., N g ̂ to keep in S g ̂

Figure 3. Simple illustration of COO on the optimization of DEDS.

s.t. e fealID |s g (̂ j) = feasible ̂

In summary, COO has two evaluation stages, crude evaluation stage and accurate evaluation stage, and the spirit of COO can be summarized as two models, the feasibility model and crude model.

j = 1, 2, ..., N g ̂

and

̂ e opt < ... < e opt ̂ |s g (l) ̂ |s g (̂ j) < ... < e opt ̂ |s g (̂ N g )̂ < e opt ̂ |s ∉ s g ̂

(27) fealID

opt

where, ê expresses the estimated feasibility; ê expresses the estimated optimal value (or crude evaluation value); Sĝ represents crudely optimal set; s ∉ Sĝ represents all the sequence solutions that do not belong to Sĝ. Specifically, in inner level, EAf is adopted to crudely optimize

4. TWO-LEVEL OPTIMIZATION FRAMEWORK BASED ON COO AND EA FOR SMCOB MODEL From the point view of optimization efficiency and robustness, a novel two-level optimization framework based on COO and EA is proposed for SMCOB model in this section. 4.1. Two-Level Optimization Structure. As the foundation of the framework, a two-level optimization structure is first introduced. On the basis of eq 26, once sequence variable s ∈ S is determined, minf∈F, s∈S J(f,s) becomes a relatively simpler

min J( f , s ̅ )| s ̅ = s f ∈F

(28)

for each s from outer level. It performs the traditional operations of continuous EA and iterates only a few steps. During the iterative process, the best so far evaluation data 9083



Article

(1) A feasibility model that is used to quickly, but probably with some mistakes, screen out the feasible solutions before crude-model evaluation, so as to avoid needless computation on infeasible solutions. (2) A crude model that is used to quickly but roughly evaluate the estimated feasible solutions. Then, true good enough solutions can be picked up with certain probability on the basis of crude evaluation value. In this section, the two models of our methodology for the SMCOB problem will be introduced. If penalty function method is adopted, when EAf iterates to kf th step, the best so far expanded objective value e(kf)|s, namely, the kf th-step evaluation value of s ∈ S, can be formulated as the sum of objective value o(kf)|s and penalty value p(kf)|s:

(objective function value plus penalty value) and penalty data in each step are recorded. Then the êfealID|s and êopt|s can be obtained by feasibility model and crude model with the data. The design of the feasibility model and crude model in our methodology will be described in detail in section 4.3(a). The other hand, in outer level, EAs is used to optimize min(e opt ̂ |s ) s∈S

fealID

s.t. e ̂

|s = feasible

(29)

It iterates based on êfealID|s and êopt|s generated from the inner level and update Sĝ, which includes the estimated top-Nĝ sequence solutions sĝ (j) ∈ Sĝ, j = 1, 2, ..., Nĝ. At accurate evaluation stage, the accurate model, namely EAf with sufficient iterations, is used to optimize minf∈F J( f, s)|̅ s=s̅ ĝ(j),sĝ(j) ∈ Sĝ, j = 1 ,2, ..., Nĝ. Then, f opt|sĝ(j) and eopt|sĝ(j), the optimal flow rate and the accurate evaluation value with respect to sĝ(j), can be obtained. Finally, a globally true good enough solution for the SMCOB model will be available with high probability by eqs 30 and 31: s g = arg min( e opt|s ĝ(j)) ĝ s (j)

s g ̂ ( j) ∈ S g ̂

e(k f )|s = o(k f )|s + p(k f )|s

Based on the stochastic optimization theory, if the parameters of EAf are set properly, it will probably converge to the optimum with infinite iterations. Then, e(∞)|s will be eopt|s, namely the optimal value of minf∈F J( f, s)| or the ̅ s=s ̅ accurate evaluation value of s ∈ S. Similarly, p′(∞)|s denotes the minimum penalty value of s ∈ S, where p′(kf)|s, kf = 1, 2, ..., Kf are defined as follows:

j = 1, 2, ..., N ĝ (30)

f g = f opt |s g

(32)

⎧ p(k f )|s if p(k f )|s < p′(k f − 1)|s ⎪ p′(k f )|s = ⎨ ⎪ otherwise ⎩ p′(k f − 1)|s

(31)

The above whole methodology is simply illustrated in Figure 4. The detailed steps will be introduced in section 4.4.

(33)

Then, the following judgment holds: ⎧ infeasible, if p′(∞)|s > 0 feasibility of s = ⎨ otherwise ⎩ feasible, ⎪

⎪

(34)

Therefore, for time-limited optimization, the task of feasibility model and crude model in the two-level framework is to fast predict p′(∞)|s and e(∞)|s with a small amount of data p′(kf)|s and e(kf)|s, kf = 1, 2, ..., Kf (Kf ≪ ∞). As shown in Figures 5 and 6, with respect to an infeasible sequence solution

Figure 4. Simple description of the two-level optimization structure based on COO and EA. Figure 5. Evaluation and penalty curves of an infeasible sequence solution s ∈ Sinfea fitted by function 35 with small amount of data.

4.3. Specific Instructions Related. (a). Feasibility Model and Crude Model in COO-Based Framework. As briefly introduced in section 3.2, the spirit of COO can be summarized as two models:

s ∈ Sinfea or a feasible one s ∈ Sfea, both e(kf)|s and p′(kf)|s, kf =1, 2, ..., ∞ are plotted as decline curves. Moreover, the penalty 9084



Article

accurate evaluation process of s illustrated in Figure 5, exhausted computation can be avoided. ⎧ ∂2 θkf ≈ θkf − 1 + ⎨ 2 Jk (θ ) ⎩ ∂θ f

⎫−1 ⎬ ⎭

θk f

⎡ × ⎢ −2ϖkf (z(k f )|s − z(̂ θ , k f )|s ) ⎢⎣ ×

∂(z(k f )|s − z(̂ θ , k f )|s ) ∂θ

θk f − 1

θk f − 1

⎤ ⎥ ⎥⎦

k f = 1, 2, ..., K f (36)

⎧∂ ⎨ 2 Jk (θ ) ⎩ ∂θ f

Figure 6. Evaluation and penalty curves of a feasible sequence solution s ∈ Sfea fitted by function 35 with small amount of data.

×

⎫ ⎬ ⎭

−1

2

θk f

⎧ ⎧ ⎪ ⎡ ∂2 ⎪ ⎪ = ⎨I − ⎨2ϖkf ⎢ 2 Jk − 1(θ ) ⎣ ∂θ f ⎪ ⎪ ⎩ ⎪ ⎩

θk f − 1

⎤−1 ⎥ ⎦

∂(z(k f )|s − z(̂ θ , k f )|s ) ∂θ θk f − 1

⎡ ∂(z(k f )|s − z(̂ θ , k f )|s ) ×⎢ ⎢⎣ ∂θ

⎫ ⎤τ ⎪ ⎥⎬ ⎥⎦ ⎪ θk f − 1 ⎭

⎧ ⎡ ∂(z(k f )|s − z(̂ θ , k f )|s ) ⎪ /⎨I + 2ϖkf ⎢ ⎢⎣ ∂θ ⎪ ⎩ ⎡ ∂2 × ⎢ 2 Jk − 1(θ ) ⎣ ∂θ f

⎡ ∂2 × ⎢ 2 Jk − 1(θ ) ⎣ ∂θ f

curve Cp′ is stabilized faster than the evaluation curve Ce. Therefore, when EAf iterates on, p′(∞)|s can be predicted with fewer steps than e(∞)s; that is, feasibility can be determined faster than crude evaluation. This feature fits well with the idea of COO. Through the simulation analysis, the entire trajectory of both Ce and Cp′ can be approximately described by the following parametrized power function, which also has a good prediction performance. k f = 1, ..., ∞

⎤ ⎥ ⎦

−1

θk f − 1

⎫ ⎫⎪ ∂(z(k f )|s − z(̂ θ , k f )|s ) ⎪⎪ ⎬⎬ × ∂θ ⎪⎪ ⎭ θk f − 1 ⎪ ⎭

Figure 7. Convergence of parameters in function 35 by iterative eqs 36 and 37 corresponding to the first picture in Figure 5.

z(̂ θ , k f )|s = Ak f −B + C

θk f − 1

⎤τ ⎥ ⎥⎦

θk f − 1

⎤−1 ⎥ ⎦

k f = 1, 2, ..., K f (37)

where Jkf(θ) is the mean square error (MSE) of fitting results of function 35 when EAf iterates to the kf th step. It is defined as eq 38, where ϖj is the weighting coefficient: kf

Jk (θ ) = f

(35)

∑ ϖj(z(j)|s − z(̂ θ , j)|s )2

k f = 1, 2, ..., K f

j=1

where ẑ(θ,kf)|s denotes the fitted value of kf th step with respect to θ, which can be replaced by p̂′(θ,kf)|s or ê(θ,kf)|s; θ = (A, B, C) represents the model parameters. Figure 7 shows that they can be fitted by the following iterative eqs 36 and 37 (the derivation is formulated in Appendix A) with a small amount of data p′(kf)|s or e(kf)|s, kf = 1, 2, ..., Kf (Kf ≪ ∞). Then p′(∞)|s and e(∞)|s can be further predicted using the fitted models; that is, p̂′(∞)|s and ê(∞)|s can be roughly determined. Although there is probably inaccuracy due to the fitting results, compared with the

(38)

For the SMCOB optimization, the initial values of the iterative eqs 36 and 37 can be set as eqs 39 and 40, where ε is a small value:

9085

θ0 = (A 0 , B0 , C0) = (0, 0.5, z(1)|s )

(39)

⎧ ∂2 ⎨ 2 J0 (θ ) ⎩∂ θ

(40)

θ0

⎫−1 ⎬ = εI ⎭



Article

Figure 8. Flowchart for the two-level optimization framework based on COO and DE.

(1) The feasibility model: When minf∈F J(f,s)| is solved by ̅ s=s ̅ EAf, the algorithm first iterate to Kf FMth step (Kf FM ≪ ∞). During the optimization process, θ of p̂′(θ,kf)|s as function 35

Based on the points introduced above, at the crude evolution stage of our methodology, the two models for SMCOB problem can be summarized as follows. 9086



Article

Table 1. Information about Yield of Crude Oil Narrow Fractions and Their Corresponding Weights yield of narrow fractions for different crude oil (%) boiling range (°C) L = 18 Llgt = 14 (500

wt coefficients wl

OA

OB

OC

OD

OE

OF

Odes

50 50 50 50 50 50 30 30 30 30 30 10 10 10 2 2 2 0

10.34 5.39 4.42 1.93 1.32 3.8 2.98 2.75 3.23 4.81 5.19 5.54 5.9 4.85 6.11 11.02 6.94 12.68

7.05 4.1 2.04 2.03 3.84 5.3 3.53 3.4 3.71 5.71 3.89 5.07 5.02 3.8 4.77 10.61 6.64 18.72

4.06 3.78 3.87 1.86 1.92 3.08 3.19 3.14 3.46 4.06 4.24 3.35 5.32 5.59 8.08 9.93 7.78 23.36

1.08 1.1 1.61 0.89 0.8 1.99 1.98 2.34 2.2 3.09 3.93 4.99 4.22 4.03 4.78 10.35 10.13 41.13

0.89 0.31 1.05 0.92 0.56 1.93 1.88 1.77 2.46 3.24 4.73 4.25 5.22 4.52 4.53 17 12.88 31.44

0.8 0.43 0.05 0.03 0.07 0.33 0.56 0.75 0.95 1.23 2.12 2.91 3.71 4.46 3.88 9.63 11.65 54.22

2.35 1.83 2.5 1.01 1.36 2.72 2.53 3.19 3.13 4.06 4.93 4.98 5.15 5.22 4.76 15.75 9.87 23.77

are iteratively fitted with data p′(kf)|s, kf = 1, 2, ..., Kf FM by eqs 36 and 37 and p′(∞)|s is further predicted. Then, êfealID|s, the estimated feasibility of s can be determined by eq 34. If êfealID|s = infeasible, p̂′(∞)|s is recorded as the estimated minimum penalty value of s. (2) The crude model: If êfealID|s = feasible, EAf go on evolving to Kf CDth step (Kf FM < Kf CD ≪ ∞). Similarly, θ of ê′(θ,kf)|s as function 35 are iteratively fitted with data e(kf)|s , kf =1, 2, ..., Kf CD by eqs 36 and 37 and e(∞)|s is predicted. Then, ê (∞)|s (or êopt|s) is recorded as the crude evaluation value of s. (b). Comparison of Candidate Solutions in EAs at the Crude Evaluation Stage. When sequence individuals and Sĝ are updated in EAs at the crude evaluation stage of our methodology, comparison between individuals should be made. Deb proposed an individual comparison technique called feasibility rules32 for EA: (1) Any feasible individual is preferred to any infeasible individual. (2) Between two feasible individuals, the one having better objective function value is preferred. (3) Between two infeasible individuals, the one having smaller constraint violation is preferred.

M

p(nf , k f ) = ω·v sum(nf , k f ) = ω· ∑ v norm(nf , k f , m) m=1 M

v (n f , k f , m )

= ω· ∑

N

m=1

∑n f= 1 v(nf , k f , m) f

(41)

where v(nf,kf,m) denotes the mth constraint violation of nf th candidate individual when EAf iterates to kf th step, which is further normalized as vnorm(nf,kf,m); vsum(nf,kf) represents the sum of total normalized constraint violation of nf th candidate individual in kf th step; ω is penalty factor. Ideally, the penalty factor setting should follow the minimum penalty rule.34,35 However, it is not necessarily easy to implement. The reason is that the exact location of the boundary between the feasible and infeasible regions is unknown. Therefore, the penalty factor is empirically set as moderate to balance the constraints and the objective.33 4.4. Detailed Steps of Two-level Optimization Framework Based on COO and EA. In this section, we take differential evolution (DE) as the representative of EA. The detailed algorithm steps of COO-based DE for SMCOB model are introduced in Figure 8. The main features of the solution methodology can be summarized as follows: (1) To solve SMCOB model more efficiently and robustly, COO-based DE is proposed. It is implemented in a twolevel optimization framework. Discrete DE, namely DEs, is used in outer level to optimize s, the blending sequence of charging tanks. Continuous DE, namely DEf (bold in Figure 8), is used in inner level to optimize f, the corresponding blending flow rates in each pipeline and time slot. (2) The solution methodology has two stages based on the idea of COO. At the crude evaluation stage, DEs in the outer level evolves and keeps several probably best so far sequence solutions in the crudely best set based on their roughly estimated feasibility and evaluation value obtained from the inner level. At the accurate evaluation stage, each sequence solution within the crudely best set is further accurately evaluated, that is, the corresponding

Motivated by this technique, modified feasibility rules (MFR) are proposed in this paper based on êfealID|s, ê(∞)|s (or êopt|s), and p̂(∞)|s. Supposing that sa ∈ S and sb ∈ S are two sequence individuals for SMCOB model, sa is considered better than sb at any of the following scenarios: (1) êfealID|sb = infeasible, but êfealID|sa = feasible (2) Both êfealID|sa and êfealID|sb= feasible, but ê(∞)|sa < ê(∞)|sb (3) Both êfealID|sa and êfealID|sb = infeasible, but p̂′(∞)|sa < p̂′(∞)|sb (c). Definition of Penalty in EAf. The penalty function method is one of the most popular techniques to handle constraints in EA.33 In this paper, when minf∈F J(f,s)| is ̅ s=s ̅ solved by EAf, penalty function is defined as follows: 9087



Article

Table 2. Information about Sulfur Concentration and Flow Rate Constraints crude oil

sulfur concentration (ppm)

flow rate (min, max) (m3/h)

crude oil

sulfur concentration (ppm)


OA OB OC

1087 848 2833

(100, 400)

OD OE OF

552 794 2040

(100, 400)

Table 3. Information about Tank Storage, Due Date for Emptying tank

oil stored

amount (m3)

due emptying time (day)

tank

oil stored

amount (m3)

due emptying time (day)

1 2 3 4 5 6 7 8

OA OA OB OB OB OC OC OC

15 000 (remaining) 50 000 20 000 20 000 10 000 10 000 20 000 10 000

10 10 30 30 30 30 30 30

9 10 11 12 13 14 15 16

OC OD OD OE OE OF OF OA

20 000 20 000 50 000 20 000 10 000 10 000 50 000 MV

30 30 30 30 30 20 20 MH

flow rate is optimally (or fully) solved by inner-level DEf to find true good enough solutions. (3) At the crude evaluation stage, a feasibility model is proposed to quickly estimate the feasibility of each sequence solution so as to avoid further computation on infeasible sequence solutions, and a crude model is proposed to provide each estimated feasible sequence solution a rough but computationally fast evaluation value instead of an accurate evaluation value. (4) Best so far crudely optimal set with a proper size instead of the single best so far optimum will ensure that at least a pair of true good enough sequence and corresponding flow rate solutions can be obtained for the SMCOB model from the accurate evaluation stage with high probability.

Table 5. Key Parameters for Mixed-Coding EA and COOBased EA mixedcoding DE key parameters population size crossover probabilitiy scaling factor mutation probabilitiy size of Sĝ (computing time, size)

5. EXPERIMENTAL STUDY As the outstanding representatives of EA, DE, and GA are successively embedded into the COO-based framework (hereinafter COO-based DE, GA) to solve a real-case SMCOB model. It has been mentioned at the beginning of section 4.2, with respect to two-level methodology for MINLPs, stochastic−deterministic strategy23−26 is hard to work in our problem, and traditional stochastic−stochastic strategy22 consumes computing time severely. Therefore, in this section, our approach is compared with mixed-coding EA.16,18 The encoding and operators of discrete and continuous optimization used in two methodologies are set the same, so as to verify the effectiveness of the COO-based hierarchical framework. Other modified schemes, such as migration operation,18 etc., are not considered here, because when they are used in mixed-coding methodology to further enhance exploitation or exploration, they can also be applied to COO-based methodology for improvement. Specifically, in the DE case, floating point encoding is used for both the flow rate variable and the sequence variable, and values corresponding to the sequence variable are limited within range (0, 1). Then the DE/rand/l/bin version36 is adopted for evolutionary operations. When the fitness of individuals is

sulfur (max) (ppm)

yield of light components (min) (%)

yield of light components (max) (%)


Omix

2000

35

50

(500, 500)

DEs

DEf

128 0.9 0.24

128 0.82 0.36

350 s, 3 500 s, 4 650 s, 5 800 s, 6 950 s, 7 1200 s, 9 1500s, 11 1800s, 13 2400s, 16

Mixedcoding GA

COO-based DE GAs

GAf

256 0.78

128 0.82

128 0.72

0.22

0.64 0.02 350 s, 3 500 s, 4 650 s, 5 800 s, 6 950 s, 7 1200 s, 9 1500s, 11 1800s, 13 2400s, 16

evaluated, random keys scheme21 is used to map the floating values to permutation for sequence variable temporarily. In the GA case, permutation encoding, partially mapped crossover, inversion mutation, tournament selection are used for sequence variable, and floating point encoding, arithmetic crossover, uniform mutation, tournament selection are used for flow rate variable.37−39 In the real-case SMCOB problem, 16 oil tanks (14 scheduling tanks, 1 initial tank, and 1 virtual tank) stored with 6 kinds of crude (TBP curves are shown in Figure 12) are charging through 2 pipelines (p1 and p2). All the source information is real data from a certain refinery and tabulated in Tables 1−4. With respect to the information, some points have to be specified: (1) wl are normalized before using and set from high to low, indicating that more attention is paid to light components of crude oil; (2) as for trace elements, only sulfur is considered. Moreover, key parameters for both methodologies are adjusted properly in advance and tabulated in Table 5. Under the Matlab2007a environment and with Intel Core 2 Quad CPU Q6600 @2.40 GHz and 1.96 GB RAM, the performance comparison with different computing time, such as 350 s, 500 s, ..., 2400 s, is made, and 35 independent runs are executed in each case. Two indexes, namely, the average objective values and variance values, are used to statistically analyze the

Table 4. Information about Specified Property and Flow Rate Requirements for Mixed Feed Oil mixed crude oil

256 0.86 0.38

COO-based DE

9088



Article

Figure 9. Optimization results of two methodologies using DE.

Figure 10. Optimization results of two methodologies using GA.

accuracy is gradually improved when the learning data increase along with the inner-level optimization iterating on. Besides, the best so far crudely optimal set Sĝ with a proper size instead of the single optimum will help reduce the risk of trapping in a local minimal sequence solution. In a word, the COO-based hierarchical

optimization results. The performance comparison between two methodologies using DE (GA) is illustrated in Figure 9 (Figure 10), which shows that the COO-based methodology will statistically outperform the mixed-coding counterpart within 1200 s in terms of both indexes. When the computing time lasts longer, the advantage of efficiency and robustness over mixedcoding EA would decrease; nevertheless, the computing time is impractical, and the solution obtained within 1200 is good enough. By tracking the evolution process of these two kinds of methodology, we find that the mixed-coding EA are easy to stick in a local minimal sequence solution. This situation only can be improved through increasing the mutation scaling factor. However, this may result in a hard convergence, unless sufficient iterations are implemented. As for the COO-based methodology, the optimization process of sequence variable and flow rate variable are separated. The performance of the whole methodology mainly depends on the feasibility estimation and crude evaluation from the inner-level optimization for each sequence solution, and this process is based on the use of a predictive model and an online parameters learning process, which is equivalent to introducing new structure information into the optimization process, and its

Figure 11. Gantt chart for an optimized schedule of SMCOB problem. 9089



Article

Table 6. Optimized Schedule of SMCOB Problem scheduling time slots start time(day) end time(day) added tanks tanks for p1 tanks for p2 emptied tanks crude oil in p1 crude oil in p1 flow in p1(m3/h) flow in p2 (m3/h) sulfur (ppm)

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

0 2.02 15 1 15 1 OA OF 309 191 1450

2.02 8.75 2 2 15 2 OA OF 309 191 1450

8.75 10.4 4 4 15 15 OB OF 240 260 1467

10.4 12.0 14 4 14 14 OB OF 240 260 1467

12.0 12.2 12 4 12 4 OB OE 231 269 819

12.2 14.0 5 5 12 5 OB OE 231 269 819

14.0 15.1 3 3 12 12 OB OE 231 269 819

15.1 17.9 10 3 10 3 OB OD 213 287 678

17.9 18.0 8 8 10 10 OC OD 251 249 1698

18.0 19.5 11 8 11 8 OC OD 251 249 1698

19.5 22.9 9 9 11 9 OC OD 251 249 1698

22.9 24.5 6 6 11 6 OC OD 251 249 1698

24.5 26.4 7 7 11 11 OC OD 251 249 1698

26.4 27.7 13 7 13 7 OC OE 264 236 1871

27.7 28.1 16 16 13 13 OA OE 215 285 920

Figure 12. Improvement of TBP curves of crude oils under the optimized schedule.

(1) In the case study, the outer level and inner level of the COO-based framework are embedded with the same EA for the sake of making comparison in the same conditions with mixed-coding methodology. However, Because COO is a soft stochastic optimization technique,29 other stochastic optimization methods with combinatorial optimization-suited nature, such as tabu search, simulated annealing, harmony search, etc., can also be adopted in the high level of the COObased framework and the performance can be further studied. (2) According to no free lunch theorem,40 uniform sampling performs as good as any search algorithm on average without some knowledge of the structural information of a problem. In this paper, the proposed feasibility model and crude model based on iteratively fitted predictive curve model can be considered as new self-learning structure information introduced into the optimization process. However, the accuracy and efficiency with respect to this self-learning mechanism are contradictory, because more accuracy needs more fitting data, which causes less efficiency, and vice versa. Therefore, how to make computation balance between accuracy and efficiency should be further studied. (3) The generality of the methodology on other complex MINLPs, such as the optimization of chemical design and synthesis problems,22−26 etc., will also be our future work.

framework is an effective and robust tool to solve the SMCOB problem, when efficiency and reliability are emphasized, which are also the common requirements for practical engineering application. Moreover, one of the detailed optimized schedules by COObased DE with 350 s are shown by Gantt chart in Figure 11 and Table 6 (the decision variable solutions are in bold type), which represents that heavy and light, sour and sweet crudes are reasonably arranged to blend, and all the specified production requirements are satisfied. Moreover, the first picture of Figure 12 shows the TBP curves of crude oils before scheduling, where the discrepancy of property between different crudes is obvious. Then, the second image in Figure 12 illustrates that, after scheduling, the mixed crudes in different scheduling time slots all have a similar property and close to the desired feedstock. Therefore, all the results indicate the effectiveness of the proposed model and solving strategy in this paper.

6. CONCLUSION In this paper, a novel solution methodology is proposed to solve a practical SMCOB problem. First, for the purpose of improving the quality and stability of feedstock property, the general mathematical formulation of SMCOB problem is established, which is a complex MINLP model. Then, from the point view of optimization efficiency and robustness, a novel two-level optimization framework based on COO and EA is proposed. Finally, the COO-based framework is successively embedded with DE and GA, and applied to a real-case SMCOB problem. The optimization results indicate both the effectiveness of the model and the efficiency as well as the robustness of the solution methodology. As for the solution methodology, there are some other directions that deserve further study.

■

APPENDIX

A. Derivation of Iterative Fitting Algorithm

Supposing that the downward trend of evaluation curve can be approximately fitted by a power function model: z(̂ θ , n) = A ·n−B + C 9090

n = 1, 2, ...

(A1)



Article

⎧ 2 ⎪ ∂ θn ≈ θn − 1 + ⎨ 2 Jn − 1(θ) ⎪ ⎩ ∂θ

where θ = (A, B, C) denotes the model parameters; ẑ(θ,n) denotes the fitted value of nth step with respect to θ. The mean square error (MSE) can be defined as n

(A2)

j=1

where ϖj is the weighting factor and z(j), j = 1, 2, ..., n are real data. Model parameters θ can be fitted by solving the following unconstrained optimization problem:

θn − 1

⎡ × ⎢ − 2ωn(z(n) − z(̂ θ , n)) ⎣ θ

×

n−1

(A3)

θ

−1 τ⎫

⎤⎪ ⎬ ⎥⎪ ⎦ θn − 1 ⎭

∂(z(n) − z(̂ θ , n)) ∂θ

θn − 1

⎤ ⎥ ⎦

Based on eqs A9 and A10, further we get

It can be noted that ẑ(θ,n) and Jn(θ) are once and twice continuously differentiable for θ, respectively. If ∃θn ∈ R3 satisfies formula A3, the following equality holds. ∂ J (θ) = 0 ∂θ n θn

⎧ 2 ⎪ ∂ θn ≈ θn − 1 + ⎨ 2 Jn − 1(θ ) ⎪ ⎩ ∂θ

2

Jn (θ) = Jn − 1(θ) + ϖn(z(n) − z(̂ θ , n))

(A5)

n−1

(A7)

× (θn − θn − 1) θn − 1

⎤τ ⎥ ⎦

≈

∂2 J (θ) ∂θ 2 n − 1

θn

⎡ ∂(z(n) − z(̂ θ , n)) ×⎢ ⎣ ∂θ

+ 2ϖn

θn − 1

⎤ ⎥ ⎦

θn − 1 −1 τ⎫

⎤⎪ ⎬ ⎥⎪ ⎦ θn − 1 ⎭

θn − 1

⎡ ∂2 × ⎢ 2 Jn − 1(θ) ⎣ ∂θ θ

⎤−1 ⎥ ⎦

θn − 1

θn − 1

⎤τ ⎥ ⎦

⎫⎫ ⎪⎪ ⎬⎬ ⎪⎪ ⎭⎭

(A12)

When initial value θ0 = (A0, B0, C0) and {∂2/∂θ2J0(θ)|θ0}−1 = εI are given properly, θ can be obtained iteratively through formulas A11 and A12 without calculating the second derivative of Jn(θ) with respect to θ and its inverse matrix.

■

∂(z(n) − z(̂ θ , n)) ∂θ

θn − 1 τ

θn − 1

⎤−1 ∂(z(n) − z(̂ θ , n)) ⎥ × ∂θ ⎦

n−1

(A8)

θn − 1

∂(z(n) − z(̂ θ , n)) ∂θ

⎡ ∂2 × ⎢ 2 Jn − 1(θ) ⎣ ∂θ θ

n−1

⎫ ⎪ ∂(z(n) − z(̂ θ , n)) × (θn − θn − 1)⎬ × ⎪ ∂θ ⎭

+ 2ϖn

⎧ ⎡ ∂(z(n) − z(̂ θ , n)) ⎤⎪ ⎪ ⎡ ∂(z(n) − z(̂ θ , n)) ×⎢ ⎥ ⎬/⎨I + 2ϖn⎢ ⎣ ⎦ ⎣ ⎪ ⎪ ∂θ ∂θ θn − 1 ⎭ ⎩

2

θn − 1

(A11)

τ⎫

We derive once and twice on both sides of eq A7 with respect to θ and then we make θ = θn. It is easy to see eq A8 and eq A9, respectively:

θn − 1

θn − 1

⎧ ⎧ ⎡ 2 ⎤−1 ∂(z(n) − z(̂ θ , n)) ⎪ ⎪ ∂ = ⎨I − ⎨2ϖn⎢ 2 Jn − 1(θ) ⎥ × ⎪ ⎣ ∂θ ∂θ ⎦ ⎪ ⎩ θn − 1 ⎩

⎧ ⎫2 ⎡ ∂(z(n) − z(̂ θ , n)) ⎤τ + ϖn⎨(z(n) − z(̂ θ , n))|θn − 1 + ⎢ |θn − 1⎥ × (θ − θn − 1)⎬ ⎣ ⎦ ∂θ ⎩ ⎭

⎡ ∂(z(n) − z(̂ θ , n)) +⎢ ⎣ ∂θ

⎤ ⎥ ⎦

n−1

⎡ ∂(z(n) − z(̂ θ , n)) ×⎢ ⎣ ∂θ

∂2 1 (θ − θn − 1)τ × 2 Jn − 1(θ)|θn − 1 × (θ − θn − 1) 2 ∂θ

⎧ ⎪ + 2ϖn⎨(z(n) − z(̂ θ , n)) ⎪ ⎩

θn − 1

⎤ ⎥ ⎦

⎧ 2 ⎧ ∂2 ⎫−1 ⎪ ∂ ⎨ 2 Jn (θ) ⎬ = ⎨ 2 Jn − 1(θ) ⎪ ⎩ ∂θ ⎭ ⎩ ∂θ θn

Given eq A4, the above equation changes into

θn

∂(z(n) − z(̂ θ , n)) ∂θ

Moreover, considering formula A9, when the Reversion Formula is applied, such that

(A6)

∂ J (θ) ∂θ 2 n − 1

×

⎤⎪ ⎥⎬ ⎦⎪ θn − 1 ⎭

n

∂(z(n) − z(̂ θ , n)) × ∂θ

∂2 × 2 Jn − 1(θ)|θn − 1 × (θ − θn − 1) ∂θ ⎫2 ⎧ ⎤τ ⎡ ∂(z(n) − z(̂ θ , n)) |θn − 1⎥ × (θ − θn − 1)⎬ + ϖn⎨(z(n) − z(̂ θ , n))|θn − 1 + ⎢ ⎦ ⎣ ∂θ ⎭ ⎩

≈

θn − 1 −1 τ⎫

⎧ ∂2 ⎫−1 ⎡ ≈ θn − 1 + ⎨ 2 Jn (θ ) ⎬ × ⎢− 2ϖn(z(n) − z(̂ θ , n)) ⎣ ⎩ ∂θ ⎭ θ θ

1 ∂ J (θ)|θn − 1 × (θ − θn − 1) + (θ − θn − 1)τ 2 ∂θ n − 1

Jn (θ) ≈ Jn − 1(θn − 1) +

∂(z(n) − z(̂ θ , n)) ∂θ

θn − 1

⎡ × ⎢− 2ϖn(z(n) − z(̂ θ , n)) ⎣ θ

We further apply second-order and first-order Taylor expansion to Jn−1(θ) and (z(n)) − ẑ(θ,n)) on θn−1 ∈ R3, respectively. Then, we have the following: Jn (θ) ≈ Jn − 1(θn − 1) +

+ 2ϖn

⎡ ∂(z(n) − z(̂ θ , n)) ×⎢ ⎣ ∂θ

(A4)

Rewrite the formula A2 as a recursive form:

∂2 J (θ) ∂θ 2 n

θn − 1

(A10)

min Jn (θ)

∂ J (θ) ∂θ n

∂(z(n) − z(̂ θ , n)) ∂θ

⎡ ∂(z(n) − z(̂ θ , n)) ×⎢ ⎣ ∂θ

∑ ϖj(z(j) − z(̂ θ , j))2

Jn (θ) =

+ 2ϖn

AUTHOR INFORMATION

Corresponding Author

θn − 1

*Phone: +86-10-62784964. Fax:+86-10-62786911. E-mail: [email protected].

(A9)

Notes

The authors declare no competing financial interest.

Given eq A4 and eq A8, it follows that 9091



■

Article

Superscripts Model

ACKNOWLEDGMENTS We would like to thank National Basic Research Program of China (No. 2012CB720500), National Natural Science Foundation of China (No. 60974008).

■

des = desired mixed feedstock due = due time lgt = light components mix = mixed feedstock min = minimum limit max = maximum limit norm = normalized sum = summarized

NOTATION

Acronym

AP = alignment probability Cp′ = penalty curve Ce = evaluation curve COO = constrained ordinal optimization DE = differential evolution DEDS = discrete event dynamic system EA = evolutionary algorithms EAf = evolutionary algorithm selected for optimizing the flow rate variable in the inner level EAs = the evolutionary algorithm selected for optimizing the sequence variable in the outer level GA = genetic algorithm MFR = modified feasibility rules OO = ordinal optimization Pro = the probability SMCOB = scheduling of multipipeline crude oil blending TBP = true boiling point

Methodology

feas = feasible feaID = the index of feasibility g = truly good ĝ = estimated good infea = infeasible opt = optimal CD = crude model FM = feasibility model Parameters Model

Bi,o,k = content of trace element k in curde o stored in tank i Bmix,max = maximum content limit of trace element k in mixed k feedstock Llgt = maximum number of narrow fractions of light components in mixed feed oil f max = maximum flow rate limit in pipeline p p f min = minimum flow rate limit in pipeline p p f mix,max = maximum flow rate limit of mixed feedstock f mix,min = minimum flow rate limit of mixed feedstock Hdue = specified emptying time of real tank i i H = scheduling horizon MV = most storage volume in charging tanks wl = weight coefficient for deviation in lth narrow fraction (normalized first) ηi,o,l = yield of lth narrow fraction of crude o stored in tank i ηdes = yield of lth narrow fraction of desired mixed feedstock l ηlgt,min,ηlgt,max = minimum and maximum yield limit of light fractions in mixed feedstock

Indices and Set Model

i i i i k l o p t IO F S

Sfea Sinfea Sĝ

1, ..., Ib (Ib = P − 1) = charging tanks feeding at the beginning time Ib + 1,...Ib + Is = charging tanks for scheduling Ib+ Is +1, ..., I (I = Ib + Is + Iv = Is + 2(P − 1)) = virtual charging tanks 1, ..., Ib, Ib + 1, ...Ib + Is, Ib + Is + 1, ..., I (I = Ib + Is + Iv = Is + 2(P − 1)) = all the charging tanks (i.e. real tanks plus virtual tanks) 1, ..., K = trace element 1, ..., L = narrow fractions oA, oB...= crude oils 1, ..., P = pipelines 1,...,T (T = Ib + Is = Iv + Is = P − 1 + Is) = scheduling period or time slot {(i,o)| set of pairs such that charging tank i is stored with crude oil o} max {f = f p,t, p = 1,...,P; t = 1,...,T | f min p ≤ f p,t ≤ f p , p = 1, ..., P; t = 1, ..., T } = search space of blending flow rate {s = (s1, ..., st, ... sIs; vIs+1,...,vt,..., vT)| s1∪...∪st∪...∪sIs = Ib + 1, ..., Ib+ Is; s1 ∩...∩ st ∩...∩ sIs = ○̷ ; vIs+1∪...∪ vt ∪...∪ vIs+Iv = Ib + Is + 1, ..., Ib + Is + Iv; vIs+1 ∩...∩ vt ∩...∩ vIs+Iv = ○̷ } = charging sequences of tanks feasible charging sequences, Sinfea ∩ Sfea = ○̷ , Sinfea ∪ Sfea= S infeasible charging sequences, Sinfea ∩ Sfea = ○̷ , Sinfea ∪ Sfea= S crudely optimal sequences, Sĝ ⊂ S

Methodology

Nĝ = number of estimated top sequences kept in Sĝ Ns, Nf = population size in DEs and DEf (or GAs and GAf) Pcs, Pcf = crossover probability in DEs and DEf (or GAs and GAf) Pms, Pmf = mutation probability in DEs and DEf (or GAs and GAf) COO

Ni = number of simulation runs in DEDS for design θi PA = a constant high probability value (i.e. 0.95)

Variables Model

ip,t = the tank charges through pipeline p in time slot t. ut = the tank finishes charging at the end of slot t ηl,tmix = yield of lth narrow fraction of mixed feedstock during slot t Fi,t = volumetric amount of crude tank i has charged during slot t Vi,t = volumetric amount of crude stored in tank i at the end of time slot t St = starting time of time slot t Et = ending time of time slot t

COO

U = crudely good set G = true good set θi = θ1, ..., θ|Θ| = Θ =designs for DEDS, where | · | represents the size of set 9092



Article

τp,t = the tank charging duration in pipeline p and slot t s = (s1, ..., st, ... sIs; vIs+1, ..., vt, ..., vT) = the charging sequence of tanks, where st and vt represent the scheduling tank and virtual tank participating in charging at the beginning of slot t, respectively f = f p,t, p = 1, ..., P; t = 1, ..., T = blending flow rates in pipeline p during the time slot t

(18) Lin, Y.; Hwang, K.; Wang, F. A mixed-coding scheme of evolutionary algorithms to solve mixed-integer nonlinear programming problems. Comput. Math. Appl. 2004, 47, 1295. (19) Angira, R.; Babu, B. V. Optimization of process synthesis and design problems: A modified differential evolution approach. Chem. Eng. Sci. 2006, 61, 4707. (20) Luo, Y.; Yuan, X.; Liu, Y. An improved PSO algorithm for solving non-convex NLP/MINLP problems with equality constraints. Comput. Chem. Eng. 2007, 31, 153. (21) Bean, J. C. Genetic algorithms and random keys for sequencing and optimization. ORSA J. Comput. 1994, 6, 154. (22) Ponce-Ortega, J. M.; Serna-Gonzalez, M.; Jimenez-Gutierrez, A. Heat exchanger network synthesis including detailed heat exchanger design using genetic algorithms. Ind. Eng. Chem. Res. 2007, 46, 8767. (23) 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, 2320. (24) Khorasany, R. M.; Fesanghary, M. A novel approach for synthesis of cost-optimal heat exchanger networks. Comput. Chem. Eng. 2009, 33, 1363. (25) Gorji-Bandpy, M.; Yahyazadeh-Jelodar, H.; Khalili, M. Optimization of heat exchanger network. Appl. Therm. Eng. 2011, 31, 779. (26) Garcia-Herreros, P.; Gomez, J. M. 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, 3977. (27) Ho, Y. C.; Sreenivas, R.; Vakili, P. Ordinal optimization of discrete event dynamic systems. J. Discrete Event Dyn. Syst. 1992, 2, 61. (28) Guan, X.; Song, C.; Ho, Y.; Zhao, Q. Constrained ordinal optimization−A feasibility model based approach. Discrete Event Dyn. Syst. 2006, 16, 279. (29) Ho, Y. C.; Zhao, Q. C.; Jia, Q. S. Ordinal Optimization: Soft Optimization for Hard Problems; Springer: New York, 2007. (30) Jiang, Y.; Cai, Y.; Huang, D. Crude oil blending scheduling problem based on order. Huagong Xuebao (Chin. Ed.) 2010, 61, 365. (31) Jiang, Y.; Cai, Y.; Huang, D. Crude oil blending schedule problem without primary oil. Huagong Xuebao (Chin. Ed.) 2010, 61, 2015. (32) Deb, K. An efficient constraint handling method for genetic algorithms. Comput. Methods Appl. Mech. Eng. 2000, 186, 311. (33) Coello Coello, C. A. Theoretical and numerical constrainthandling techniques used with evolutionary algorithms: a survey of the state of the art. Comput. Methods Appl. Mech. Eng. 2002, 191, 1245. (34) Richardson, J. T.; Palmer, M. R.; Liepins, G. E.; Hilliard, M. Some guidelines for genetic algorithms with penalty functions. Proceedings of the Third International Conference on Genetic Algorithms; Morgan Kaufmann: San Mateo, 1989; p 191. (35) Davis, L. Genetic Algorithms and Simulated Annealing; Pitman: London, 1987. (36) Storn, R.; Price, K. Differential evolutionA simple and efficient heuristic for global optimization over continuous spaces. J. Global Optim. 1997, 11, 341. (37) Goldberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison-Wesley: MA, 1989. (38) Michalewicz, Z.; Janikow, C. Z.; Krawczyk, J. B. A modified genetic algorithm for optimal-control problems. Comput. Math. Appl. 1992, 23, 83. (39) Srinivas, M.; Patnaik, L. M. Genetic algorithms: A survey. Comput. J. 1994, 27, 17. (40) Wolpert, D. H.; Macready, W. G. No free lunch theorems for optimization. IEEE Trans. Evol. Comput. 1997, 1, 67.

Methodology

êfealID|s = estimated feasibility of s eopt |s (or e(∞)|s) = optimal evaluation value (or accurate evaluation value) of s êopt|s (or ê(∞)|s) = estimated optimal evaluation value (or crude evaluation value) of s p′(∞)|s = minimum penalty value of s p̂′(∞)|s = estimated minimum penalty value of s f opt|s = optimal flow rate solution with respect to sequence s sopt = global best sequence solution sĝ (j) = estimated top-jth best so far sequence solution sg = global good enough sequence solution f g = global good enough flow rate solution

■

REFERENCES

(1) Kelly, J. D.; Mann, J. L. Crude oil blend scheduling optimization: an application with multimillion dollar benefits. Part 1. The ability to schedule the crude oil blendshop more effectively provides substantial downstream benefits. Hydrocarb. Process. 2003, 82, 47. (2) Kelly, J. D.; Mann, J. L. Crude oil blend scheduling optimization: an application with multimillion dollar benefits. Part 2. The ability to schedule the crude oil blendshop more effectively provides substantial downstream benefits. Hydrocarb. Process. 2003, 82, 72. (3) Shah, N. Mathematical programming techniques for crude oil scheduling. Comput. Chem. Eng. 1996, 20, S1227. (4) Lee, H. M.; Pinto, J. M.; Grossmann, I. E.; Park, S. Mixed-integer linear programming model for refinery short-term scheduling of crude oil unloading with inventory management. Ind. Eng. Chem. Res. 1996, 35, 1630. (5) 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, 3085. (6) Reddy, P.; Karimi, I. A.; Srinivasan, R. Novel solution approach for optimizing crude oil operations. AIChE J. 2004, 50, 1177. (7) Moro, L. F. L.; Pinto, J. M. Mixed-integer programming approach for short-term crude oil scheduling. Ind. Eng. Chem. Res. 2004, 43, 85. (8) Wang, J.; Rong, G. Robust optimization model for crude oil scheduling under uncertainty. Ind. Eng. Chem. Res. 2010, 49, 1737. (9) Bai, L.; Jiang, Y. H.; Huang, D. X.; Liu, X. G. A Novel scheduling strategy for crude oil blending. Chinese J. Chem. Eng. 2010, 18, 777. (10) Mendez, C. A.; Grossmann, I. E.; Harjunkoski, I.; Kabor, P. A Simultaneous optimization approach for off-line blending and scheduling of oil-refinery operations. Comput. Chem. Eng. 2006, 30, 614. (11) Pan, M.; Li, X. X.; Qian, Y. New approach for scheduling crude oil operations. Chem. Eng. Sci. 2009, 64, 965. (12) Li, J.; Karimi, I. A.; Srinivasan, R. Recipe determination and scheduling of gasoline blending operations. AIChE J. 2010, 56, 441. (13) Grossmann, I. E. Review of nonlinear mixed-integer and disjunctive programming techniques. Optim. Eng. 2002, 3, 227. (14) Ambrosio, C. D.; Lodi, A. Mixed integer nonlinear programming tools: A practical overview. 4OR: Q. J. Oper. Res 2011, 1. (15) Bäck, T. Evolutionary Algorithms in Theory and Practice; Oxford University Press: New York, 1996. (16) Costa, L.; Oliveira, P. Evolutionary algorithms approach to the solution of mixed integer nonlinear programming problems. Comput. Chem. Eng. 2001, 25, 257. (17) Young, C. T.; Zheng, Y.; Yeh, C. W.; Jang, S. S. Informationguided genetic algorithm approach to the solution of MINLP problems. Ind. Eng. Chem. Res. 2007, 46, 1527. 9093


A Novel Two-Level Optimization Framework Based on Constrained

Recommend Documents