Development of Optimal Decoking Scheduling Strategies for an

Jul 4, 2006 - Development of Optimal Decoking Scheduling Strategies for an Industrial Naphtha Cracking Furnace System ... Proactive Scheduling Strateg...
48 downloads 13 Views 412KB Size
5738

Ind. Eng. Chem. Res. 2006, 45, 5738-5747

Development of Optimal Decoking Scheduling Strategies for an Industrial Naphtha Cracking Furnace System Heejin Lim,§ Jaein Choi,†,| Matthew Realff,† Jay H. Lee,*,† and Sunwon Park‡ KT&G Central Research Institute, 302 Shinseong-dong, Yuseong-gu, Daejeon, 305-805 Korea, Department of Chemical and Biomolecular Engineering, KAIST, Deajeon, Korea, and School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332

A naphtha cracking furnace produces various petrochemical products ranging from ethylene to pitch. Continuous operation of the furnace leads to coke formation, in increasing amounts, on the inner surface of the cracking coils. The coke deposits inside the coils decrease the productivity of the furnace, and therefore, its operation is periodically stopped for decoking to restore the productivity. An industrial furnace system has many furnaces in parallel, which share a common naphtha feed flow. This leads to an optimization problem in which the decoking schedule for a set of parallel furnaces must be decided simultaneously along with the distribution of the inlet naphtha feed flow among them. The optimization problem can be formulated as an MINLP, which may be computationally intractable for an industrially sized problem. Three alternative solution strategies have been developed to circumvent the inevitable nonlinear terms in the objective function. The performance of each strategy is presented in terms of the solution quality and computational time. 1. Introduction Naphtha cracking furnaces are used to convert the naphtha of a refinery into smaller hydrocarbon molecules, mostly ethylene and propylene. Cracking is an endothermic reaction and is carried out in “cracking coils”, which typically have short residence times. Continuous operation of the furnace leads to coke formation on the inner surface of the cracking coils. The amount of coke inside the coils increases with operation time. The coke hinders heat transfer through the coil wall and, thus, decreases the productivity of the furnace. To maintain the productivity, the input energy must be continually increased, and this raises the tube skin temperature (the coil surface temperature). Furthermore, coke deposition in excessive amounts jeopardizes the plant safety by increasing the pressure drop inside the coils. In practice, excessive coke deposits can plug the coil, which can result in a catastrophic failure of the furnace system. For the sake of production efficiency and plant safety, the furnace is periodically removed from production to be stripped of the deposited coke inside its coils using steam. This coke removal process is called “decoking”. Typically, this is done if the tube skin temperature or the pressure drop of the coil reaches its maximum limit.1 The operation time between two successive decoking operations is termed a “decoking interval”. If the decoking interval is too long, deposited coke in excessive amount lowers the furnace production efficiency. Conversely, if it is too short, the actual production time would decrease.2 Therefore, optimization of the decoking schedule is highly desirable to maximize the overall profit of the furnace system. To optimize the decoking schedule, a furnace model approximating the furnace operation is required. There are several commercial software packages for naphtha cracking furnace simulation, such as SPYRO3 and CRACKER,4 but these * To whom all correspondence should be addressed. Tel: (404) 3852148. Fax: (404) 894-2866. E-mail: [email protected]. † Georgia Institute of Technology. ‡ KAIST. § KT&G Central Research Institute. | Present address: Process Systems Team, Process Technology R&D, LG Chem. Ltd. Research Park, Daejeon, 305-380, Korea.

Figure 1. One furnace system of 10 furnaces in parallel.

packages are not suitable for direct use in solving the scheduling problem. When naphtha is thermally cracked in the furnace, thousands of reactions occur.5 Due to the huge number of reactions, furnace simulation, which has to be embedded in the decoking scheduling problem, has a high computational load. As a replacement for the rigorous simulation model for the decoking scheduling problem, we have developed a regression model using data obtained from the rigorous simulation. A furnace system has many furnaces in parallel that share a common naphtha feed to the system, as shown in Figure 1. The naphtha feed flow rate of a furnace varies within a given operating range. Over a given operating horizon, the inlet naphtha feed flow rate to each furnace is controlled to accommodate a given decoking schedule; however, the total operating capacity of the furnace system has to be kept constant to stabilize the downstream processing units, such as coolers and separators. Thus, the total flow rate in Figure 1, Ftotal, should be kept constant despite the intermittent decoking operations of the individual furnaces. This means that, if one furnace is not in operation due to decoking, the extra load must be shared by the other furnaces in the same furnace system. Determining the decoking schedule is very important to ethylene producers, but modeling and optimization present a difficult challenge due to the complexity of thermal cracking reactions. Only a few researchers have tried to solve the decoking scheduling problem. Edwin and Balchen6 solved a dynamic optimization problem of the decoking interval to obtain time-dependent operation trajectories of a single thermal cracking furnace with respect to the coke formation. They formulated a rigorous furnace model, including the mass, energy, and

10.1021/ie050129n CCC: $33.50 © 2006 American Chemical Society Published on Web 07/04/2006

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5739

Figure 2. Neural network model of the cracking furnace for the decoking scheduling problem.

momentum balances. The trajectories of feed flow rate, steamto-hydrocarbon ratio, and energy consumption were optimized. The effects of operating condition changes on the determined decoking interval were shown for a single furnace operation. In real situations, however, multiple cracking furnaces are operated in parallel. If the decoking interval or operating conditions (especially the flow rate) of a furnace are changed, it affects the operating conditions of the other furnaces in the same furnace system. For multiple furnaces, Jain and Grossmann2 formulated and solved the decoking scheduling problems as an MINLP to determine the optimal decoking interval of a multiple furnace system. They assumed that conversion to ethylene decreases exponentially with on-stream time and then determined optimal decoking intervals for different types of feeds when the inlet furnace feed flow rate and tube skin temperature are kept constant. As previously explained, to keep constant product yields, the total naphtha feed flow rate into the furnace system has to be kept constant, but the feed flow rate to each furnace may vary within a given range. In addition, the tube skin temperature needed for a constant yield increases as the coke deposits inside the coils. In this study, we solve the decoking scheduling problem for a multiple tube system with the constraint of a constant total throughput. To obtain a full-scale solution to the problem, the inlet naphtha feed flow rate to each furnace has to be treated as a continuous decision variable and the decoking decision, as a binary variable. This problem is naturally formulated as an MINLP, which is difficult to solve for practically sized furnace systems. Instead of solving the MINLP, we propose three alternative solution strategies that circumvent the nonlinear terms in the objective function. The performance of each strategy will be presented in terms of the solution quality and computational time. The best strategy among the three alternatives is then selected to solve the decoking scheduling problem for a 120-day time horizon. The decoking schedule formed by the selected strategy is also compared with the one generated by simple heuristics. 2. Furnace Modeling Using Simulation Data The furnace simulators such as SPYRO3 and CRACKER4 simulate ethylene furnaces by solving the mass, energy, and momentum balance equations with their own thermal cracking reaction sets, including the coke deposition reactions. These furnace simulators are useful tools for developing process operating rules and optimizing operating conditions, but they are not convenient for use in scheduling due to their complexity and their implicit computational structure. As an alternative, for the purpose of solving the decoking scheduling problem, we derived a regression model from simulation data obtained with CRACKER.4 The model calculates the coke thickness, yields of the products, and furnace conditions for a given operating condition. The thickness of the coke inside the coil is the most important variable for the decoking scheduling. The

coke deposition dynamics are of integrating type, and the current coke layer thickness and operating conditions combine to determine the further growth of the coke layer.7-9 To represent the integrating dynamics, the coke thickness is fed back from the output node to the input node in the proposed model, as shown in Figure 2. In this study, the sampling interval of the model data is single, chosen as 1 day, which is the same as the time required for a decoking operation. Hence, the model gives a new prediction of the output every 1 day. When the furnace is on line, it is very important to keep the yields of the products constant. Fluctuations in the yields act as disturbances to downstream processing units of the furnace system, such as the coolers and separators. To minimize the fluctuation, four main operating variables are controlled: naphtha feed composition (naphtha feed type), coil outlet temperature (COT), dilution steam to inlet naphtha feed flow rate ratio (dilution steam ratio), and inlet naphtha feed flow rate. In the proposed databased model, COT, feed type, the flow rates of steam and the feed, and the previous coke thickness is applied as inputs. The naphtha feed consists of hundreds of components, and a general model would require the composition of each component as input. In our study, we used as input the type of blended naphtha feeds, the classification used by the company that provided the data. A furnace produces various products ranging from ethylene to pitch. To see the effect of operation optimization, two major products of the furnace, ethylene and propylene, are selected. Outputs of the model are the yields of ethylene and propylene, tube skin temperature, pressure drop, and current coke thickness, as shown in Figure 2. The rate of thermal cracking reactions is dependent on the inputs of the proposed model.3,5,7,8 If the mass flow rate of each product is used as the model output instead of the product yield, which is related to the reaction rate directly, the output becomes too sensitive to feed flow rate changes. Thus, the fractions of ethylene and propylene in the outlet flow are selected as outputs of the model. Coke affects the tube skin temperature and the pressure drop. Coke insulates the inside of the coils and increases the heat transfer resistance of the coils. To keep the yields of products constant, the input energy to the coils would need to be increased, and this raises the tube skin temperature. Reduction of the cross-sectional area inside the coils by coke deposition causes the residence time to decrease and the pressure drop to increase. Due to these effects, the production rate is changed with changing coke thickness. After trying several different types of activation function, the linear type was found to be the best. Applying a linear activation function, the network model in Figure 2 is mathematically represented as follows: J

yˆ k,i,t )

∑j wk,jxˆ j,i,t + bk

(1)

5740

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Figure 3. Comparison of CRACKER data with neural network model data (+: training data, O: validation data).

Hence, the model fitted has a simple affine structure. All the input and output variables of this model (xˆ j,i,t and yˆ k,i,t) are normalized. xˆ coke,i,t is the recurrent variable of yˆ coke,i,t-1. Initial thickness of furnace i, xˆ coke,i,1, is represented by Ci,0. Data for estimation and validation of model parameters wk,j and bk are obtained from simulation data obtained with CRACKER over a wide operating range. All the data are used after normalization. 960 data sets are used for the parameter estimation, and 120 data sets are used for the validation of the model. The validation sets are selected randomly from the whole data sets. Using the back-propagation with the gradient descent method, model parameters are estimated with 0.1% relative error for the estimation data sets. This model gave 0.075% relative error for the validation sets. The variance of normalized simulation data is 0.0004. The simulation data from CRACKER and the predictions from the neural network model are plotted in Figure 3. The figure shows that the simple affine model predicts the furnace conditions very well. The proposed neural network model was developed for general purposes, but it is simplified on the basis of the plant operating rules we are familiar with. The following three assumptions are introduced to reduce the number of variables for the furnace decoking scheduling problem: (1) COT is kept constant over the scheduling horizon. In normal operation of a furnace, COT is kept constant because COT affects the yields of the products significantly. (2) The naphtha feed composition is also kept constant over the scheduling horizon. The feed composition affects the yields of products. In practice, to minimize the variation, the composition is controlled by blending various naphtha inventories. (3) The dilution steam ratio is kept constant. It is possible to change this ratio, but it is not preferred due to the expensive cost of steam and its effect on the downstream processes. With these assumptions, the feed flow rate remains as a key operating variable of the model. 3. Problem Formulation The decoking scheduling problem of a system composed of multiple parallel furnaces is formulated as a MINLP. The nonlinearity arises from the inlet naphtha feed flow rate multiplying the product weight fraction to give the product amount. However, the MINLP strategy is not a practical choice

in solving industrial decoking scheduling problems due to its heavy computational load and potential problems with local optimality. Computational speed and robustness of the scheduling solution method is crucial in industrial applications where the solution may need to be updated frequently on the basis of updated plant information. The MINLP model can be reformulated into an MILP or NLP by making reasonable assumptions. In this section, we develop the MINLP model and propose three alternative solution strategies that give significant computational advantages over the original MINLP. 3.1. MINLP Formulation. The objective of this problem is to maximize the profit, which is the weighted sum of product amounts. For the sake of simplicity, the products we consider are ethylene and propylene only, of which the produced amounts are weighted by cost parameters CC2 and CC3. In commercial naphtha crackers, the production rates of other product types can also be significant, and they can easily be added to the formulation if this is the case. The objective function is represented as follows: T

obj ) max

I

∑ ∑ (CC2pC2,i,t + CC3pC3,i,t) t)1 i)1

(2)

pC2,i,t ) Fi,tyC2,i,t

(3)

pC3,i,t ) Fi,tyC3,i,t

(4)

T is the total number of sample intervals within the scheduling horizon, and I is the number of furnaces. In the objective function (2), production amounts pC2,i,t and pC3,i,t are represented as flow rate Fi,t () xfeed,i,t) multiplied by weight fractions of ethylene (yC2,i,t) and propylene (yC3,i,t), as in eqs 3 and 4. The flow rate, Fi,t, of furnace i at time t is a continuous variable in the range of Fmin e Fi,t e Fmax. As explained previously, the sum of the feed flow rate for each furnace at time t is kept constant at the system capacity. Thus, the total flow rate is constrained as follows: I

Fi,t ) Ftotal ∑ i)1

(5)

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5741

nonlinear terms. The quadratic terms can be linearized by a piecewise linear approximation. The default form of the objective function is the same as eq 2.

Table 1. MINLP Equation Summary objective function equalities

objective function, 2 furnace model, 1 production amount, 3, 4, 6, 7 flowrate constraint, 5 coke thickness, 9 decoking constraint, 8 tube skin temperature, 10 pressure drop, 11

inequalities

T

objS ) max

yC2,i,t ) (yˆ C2,i,t(yC2max - yC2min) + yC2min)(1 - di,t)

(6)

yC3,i,t ) (yˆ C3,i,t(yC3max - yC3min) + yC3min)(1 - di,t)

(7)

In eqs 6-7, to set the production amounts to zero when the decoking is performed for furnace i at time t, (1 - di,t) is multiplied. di,t is a binary variable that represents whether the decoking process is performed or not for furnace i at time t. If a furnace i is decoked at time t, di,t is 1. Otherwise, di,t is 0. During operation, only one furnace can be decoked at any given time due to the flow rate constraint of eq 5 and furnace capacity. Therefore, the binary variable di,t is constrained as follows: I

(8)

In the same manner as in eqs 6 and 7, coke thickness is calculated from the furnace model of eq 1. xˆ coke,i,t, a normalized value of current coke thickness, is a recurrent variable from the previously calculated coke thickness, yˆ coke,i,t-1, as below in eq 9.

xˆ coke,i,t ) yˆ coke,i,t-1(1 - di,t)

S S (CC2pˆ C2,i,t + CC3pˆ C3,i,t ) ∑ ∑ t)1 i)1

S S ) pC2,i,t (1 - di,t) pˆ C2,i,t

In eqs 3 and 4, the weight fractions yC2,i,t and yC3,i,t are values of the model outputs from eq 1 but scaled back to their original units.The model outputs yˆ C2,i,t and yˆ C3,i,t can be scaled back to the original units as follows:

di,t e 1 ∑ i)1

I

(9)

To express the coke removal after each decoking, xˆ coke,i,t is set to zero if di,t is 1. The tube skin temperature, Tskin,i,t, and the pressure drop, ∆Pi,t, are calculated from the model of eq 1 according to

- yTmin ) + yTmin ) e TUskin (yˆ Tskin,i,t(yTmax skin skin skin

(10)

max min min - y∆P ) + y∆P ) e ∆PU (yˆ ∆P,i,t(y∆P

(11)

Whether a decoking is performed or not for furnace i at time t, the tube skin temperature and the pressure drop should not exceed their upper limits. Thus, the denormalized values of the tube skin temperature and pressure are not multiplied by (1 - di,t). The MINLP is summarized in Table 1. During normal operation of the given furnace system, the decoking interval of each furnace is ∼50-60 days, depending on the operating conditions. To see at least one whole decoking interval of every furnace, a minimum 120-day time horizon is required. In this scheduling horizon, the MINLP problem has 1200 binary and 16 920 continuous variables with 12 120 equality and 2520 inequality constraints, including the furnace model. It is well-known that any large scale MINLP has computational limitations in practice. When solving this MINLP for a 120-day time horizon by CONOPT of GAMS, the relative gap of its integer solution does not decrease under 30% in 24 h on a machine with 1.7 GHz of P4 CPU and 1 GB of RAM. To reduce the computational load of this problem, we next propose several alternative strategies. 3.2. MILP Formulation with Piecewise Linear Functions. The objective function of eq 2 in the MINLP has quadratic

S pˆ C3,i,t

)

S pC3,i,t (1

(13)

- di,t)

(14)

S S S ) Fi,t yC2,i,t pC2,i,t

pC3,i,tS S yC2,i,t

) )

(12)

(15)

S Fi,t

S yC3,i,t (yˆ C2,i,t(ymax C2

min - ymin C2 ) + yC2 )

(16) (17)

S min min yC3,i,t ) (yˆ C3,i,t(ymax C3 - yC3 ) + yC3 )

(18)

To set the production amounts to zero when the decoking is S performed in furnace i at time t, (1 - di,t) is multiplied by pC2,i,t S and pC3,i,t in eqs 13 and 14. These equations are linearized using the big-M method as follows: S S e pC2,i,t + MC2di,t pˆ C2,i,t

(19)

S S e -pC2,i,t + MC2di,t -pˆ C2,i,t

(20)

S pˆ C2,i,t

e MC2(1 - di,t)

(21)

S -pˆ C2,i,t

e MC2(1 - di,t)

(22)

For the nonlinear terms in eqs 15 and 16, a separable programming technique, a piecewise linearization-based method, is applied. To apply this method, all nonlinear terms should be in separable forms, which are functions of a single variable; S S for example, x2 or 1/x.10 The nonlinear terms of pC2,i,t and pC3,i,t in eqs 15 and 16 are nonseparable since they are functions of two variables, F and y. The nonseparable term of eq 15 can be modified to separable functions as follows: S S S1 S2 Fi,t yC2,i,t ) (uC2,i,t )2 - (uC2,i,t )2

1 S S1 S ) (Fi,t + yC2,i,t ) uC2,i,t 2 1 S S2 S ) (Fi,t - yC2,i,t ) uC2,i,t 2

(23) (24) (25)

S In eqs 23-25, the nonlinear term Fi,tyC2,i,t is replaced with S1 S2 S1 single variable functions (uC2,i,t)2 and (uC2,i,t)2. Let pC2,i,t ) ( S1 uC2,i,t)2. Then, its piecewise linear function can be applied to S1 S1 and uC2,i,t as follows: pC2,i,t

m5

S1 ) pC2,i,t

∑ m)m

S1 AS1 C2,m rC2,m,i,t

(26)

S1 BS1 C2,m rC2,m,i,t

(27)

1

m5

S1 uC2,i,t )

∑ m)m

1

m5

∑ m)m

S1 rC2,m,i,t )1

(28)

1

S1 AS1 C2,m and BC2,m of eqs 26 and 27 are the vertices of the S1 S1 ) (uC2,i,t )2. In this particular piecewise linear function of pC2,i,t S1 case, a five-vertex linearization has been used. To make pC2,i,t

5742

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Table 2. Equation Summary of MILP with Piecewise Linear Functions objective function equalities

objective function, 12 furnace model, 1 production amount, 19-22, 26-28, 29-34 flowrate constraint, 5 coke thickness, 35-38 decoking constraint, 8 tube skin temperature, 10 pressure drop, 11

inequalities

F1 or F2. When the furnace system has Nf furnaces and one furnace is being decoked, the remaining (Nf - 1) furnaces have the same flow rate of F1 ) Ftotal/(Nf - 1). If all furnaces are operated and none is being decoked, they all have the flow rate of F2 ) Ftotal/Nf. Because the flow rates are switched if one of the furnaces is decoked at time t, eqs 15 and 17 with appropriately replaced variables, representing the ethylene production, are changed as follows using the big-M method: I

S1 and uC2,i,t lie on one of the straight-line segments of the S1 in the piecewise functions, the continuous variables rC2,m,i,t range of [0, 1] are introduced as weights to the vertices of the S1 piecewise linear functions. Two adjacent rC2,m,i,t are allowed be nonzero for each furnace i at time t; the others must be forced to be zero.10 A new set of binary variables, aC2,m,i,t, is introduced to represent this stipulation. The size of the set and number of constraints are dependent on the number of vertices used in the linearization, in the following manner:

F F min min pC2,i,t e F1(yˆ C2,i,t (ymax C2 - yC2 ) + yC2 ) + MC2(1 -

(39)

I

F F min min e -F1(yˆ C2,i,t (ymax -pC2,i,t C2 - yC2 ) + yC2 ) + MC2(1 -

(40)

di,t ∑ i)1

F F min min pC2,i,t e F2(yˆ C2,i,t (ymax C2 - yC2 ) + yC2 ) + MC2

(29)

S1 rC2,m - aC2,m1,i,t e 0 1,i,t

(30)

S1 - aC2,m1,i,t - aC2,m2,i,t e 0 rC2,m 2,i,t

(31)

S1 - aC2,m2,i,t - aC2,m3,i,t e 0 rC2,m 3,i,t

(32)

S1 - aC2,m3,i,t - aC2,m4,i,t e 0 rC2,m 4,i,t

(33)

S1 - aC2,m4,i,t e 0 rC2,m 5,i,t

(34)

F F min min -pC2,i,t e -F2(yˆ C2,i,t (ymax C2 - yC2 ) + yC2 ) + MC2

1

S2 S2 Likewise, uC2,i,t and pC2,i,t in eq 25 are represented in the S of eq 14 can be replaced same manner as in eqs 26-34. pC3,i,t with piecewise linear functions in the same manner as in eqs 23-34. The coke thickness recursion eq 9 can be linearized by the big-M method as follows:

S S e yˆ coke,i,t + Mcokedi,t xˆ coke,i,t

(35)

S -yˆ coke,i,t

(36)

+ Mcokedi,t

S xˆ coke,i,t e Mcoke(1 - di,t)

(37)

e Mcoke(1 - di,t)

(38)

S -xˆ coke,i,t

(41)

I

∑ aC2,m,i,t ) 1 m)m

e

di,t) ∑ i)1

I

m4

S -xˆ coke,i,t

di,t) ∑ i)1

Other constraints on the total flow rate, the tube skin temperature, and the pressure drop are the same as in eqs 5, 10, and S 11 by replacing Fi,t and yˆ k,i,t with Fi,t and yˆ Sk,i,t. This MILP with piecewise linear functions is summarized in Table 2. 3.3. MILP Formulation with Fixed Flow Rates. The separable programming technique enables the nonlinear relationships to be made linear, but it increases the problem complexity by increasing the number of binary variables. As an alternative, one can remove the nonlinear terms of the objective function by fixing the values of the naphtha feed flow rate, Fi,t , a priori. For this alternative, we assume that all furnaces except the one under the decoking equally share the total flow rate, Ftotal, to satisfy the constraint of eq 5. The formulation is the same as the piecewise linear one, except that the pˆ S, pS variables are replaced by pˆ F, pF to represent fixed product quantities. Similarly, yˆ S, yS variables are replaced by yˆ F, yF. By the assumption of the equal share of total feed, Ftotal, the flow rate of each furnace FF,i,t is fixed as one of two flow rates

di,t ∑ i)1

(42)

The analogous constraints should be written for the C3 variables to replace eqs 16 and 18. Other constraints for the total flow rate, coke thickness, tube skin temperature, and pressure drop are the same as in eqs 5, 35-38, 10, and 11 by replacing Fi,t and yˆ k,i,t with FF,i,t and yˆ Fk,i,t. 3.4. Sequential Solution Strategy. As an extension of the previous strategy, this strategy solves MILP as a master problem and then solves NLP as a slave problem. In general, a similar sequential solution strategy has been applied to solve the MINLP problems for the heat exchanger network synthesis (HENS) problems.11 In these HENS problems, a master problem and a slave problem are iteratively solved until their solutions converge. But the sequential strategy for the decoking scheduling problem is different from the HENS problems due to certain operational characteristic of the furnace system. At the master problem level, the decoking schedule di,t is determined by solving the MILP with fixed flow rates. At the subproblem level, NLP is solved to optimize the flow rate, Fi,t, under the decoking schedule determined from the master problem. If the decoking schedule is fixed by the master problem, no naphtha feed is fed to the furnace during the decoking process. Because there is no information to improve or change the zero flow rate of the previously determined decoking days, the proposed strategy sequentially solves the master problem and the subproblem only once. The master problem is exactly the fixed flow rate formulation discussed in Section 3.3, and the subproblem is congruent to the original MINLP problem of eqs 1-11, except for the binary variables di,t. At the subproblem level, di,t is replaced with a parameter dˆ i,t, which has the same value as di,t obtained from the master problem. As in the previous alternatives, the objective function of the subproblem can be written as follows: T

objN ) max

I

N N + CC3pC3,i,t ) ∑ ∑ (CC2pC2,i,t t)1 i)1

(43)

N N N min min ) Fi,t (yˆ C2,i,t(ymax ˆ i,t) pC2,i,t C2 - yC2 ) + yC2 )(1 - d

(44)

N N N min min ) Fi,t (yˆ C3,i,t(ymax ˆ i,t) pC3,i,t C3 - yC3 ) + yC3 )(1 - d

(45)

N In this NLP, the flow rate, Fi,t , is a continuous variable having the range of Fmax e Fi,t e Fmin. Equation 5 is applied to this

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5743 Table 3. Equation Summary of Sequential Strategy eq type

eq name

master problem

subproblem

objective function equalities

objective function furnace model production amount flowrate constraint coke thickness decoking constraint tube skin temperature pressure drop

2 1 39-42 5 35-38 8 10 11

43 43 44-45 46 47 none 48 49

inequalities

problem as follows:

Figure 5. Profit of each solution strategy. I

N Fi,t ) Ftotal ∑ i)1

(46)

Other constraints of eqs 9-11 are replaced with the following equations. N N ) yˆ coke,i,t-1 (1 - dˆ i,t) xˆ coke,i,t

(47)

- yTmin ) + yTmin ) e TUskin (yˆ TNskin,i,t(yTmax skin skin skin

(48)

N max min min (y∆P - y∆P ) + y∆P ) e ∆PU (yˆ ∆P,i,t

(49)

The proposed sequential strategy formulation is summarized in Table 3. 4. Computational Results 4.1. Comparison of the Three Solution Strategies. By solving the furnace optimization problems with a 120-day horizon, we cannot obtain a reasonable solution with the piecewise linear strategy due to degeneracy in reasonable computational time. Thus, to compare the computational load, the horizon is reduced to 20 days. The three proposed solution strategies are compared for 20 sets of randomly generated initial coke thicknesses using a 20-day scheduling horizon. GAMS was used to formulate the problem, and CPLEX and CONOPT solvers are used for the MILP and NLP, respectively, on a machine with 1.7 GHz of P4 CPU and 1 GB of RAM. The performance of each strategy is presented in terms of the solution quality and computational time. Figure 4 shows the computational time of the three solution strategies for the 20 sets of initial coke thickness to achieve 5% optimality gap. The MILP formulation with the piecewise linear functions shows significantly poorer computational performance as compared to the other two strategies due to the large number of binary variables introduced with the use of piecewise linear approximations. The sequential strategy takes slightly more computational time than the fixed flow rate

Figure 4. Computational time comparison of each solution strategy.

strategy because it further solves the flow rate optimization problem (NLP). Figure 5 shows the solution quality of each strategy. The profit in Figure 5 is the objective function value (sum of the weighted products amounts) minus assumed fixed operating cost. The specific profits are dependent on the initial coke thickness; however, the profits of the fixed flow rate strategy solutions are smaller than those of the other strategies. To compare the solution qualities, the result of each strategy for 20 cases is reported in Figure 6. The solutions of the three strategies for one of the 20 cases are classified as “worst”, “mid”, and “best” according to the profits. This figure represents the solutions from the fixed flow rate strategy that are, on average, better than those from the other strategies. The fixed flow rate strategy and the sequential one have the same decoking schedule because the MILP model of the fixed flow rate strategy is the first part of the sequential strategy. The improved quality of the sequential strategy demonstrates the importance of the flow rate optimization. Despite the much heavier computational load, the solution strategy with piecewise linear functions shows little improvement as compared to the sequential strategy, as shown in Figure 4, because its solution quality is limited by the piecewise linearization. In summary, the sequential strategy appears to be the best strategy in terms of solution quality and computational efficiency. 4.2. Comparison of the Sequential Strategy to a Simple Heuristic Strategy. We compared the three computationally feasible solution strategies of the decoking scheduling problem in the previous section and concluded that the sequential strategy is the best. In this section, the sequential strategy is tested for a more realistically sized problem by extending the optimization horizon from 20 to 120 days. Decoking scheduling with a 20day time horizon is not meaningful when the schedule is applied to real operation because the actual decoking interval of our furnace system is over 50 days. The scheduling solution obtained with the limited time horizon tends to leave the furnaces heavily coked at the end of the time horizon to maximize the total profit only in the given time horizon. Furthermore, during real operation, the furnaces may have some uncertain or unexpected events, such as extremely fast coke deposition by hot spots on the furnace coil wall. With this significant uncertainty in real operations, applying the deterministic solution may cause serious operational problems, such as violating the furnace operation limits. Thus, we suggest applying only the first half of the deterministic solution. To make this strategy practicable, the deterministic solution strategy has to be computationally efficient to obtain a solution for a sufficient time horizon, 120 days, in which at least two decoking decisions are made for each furnace regardless of the initial coke thickness of each furnace. In the previous section, we demonstrated the benefit of feed flow rate optimization for each furnace along with the decoking decision; however, in the current operating policies, the inlet

5744

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Figure 6. Comparison of solution qualities of the three strategies.

Figure 7. Profit of 20 random initial cases for 120-day time horizon.

Figure 8. Profit vs decoking day for the heuristic solution of one case (case 5 in Figure 7).

naphtha feed flow rate for each furnace is often kept fixed under the constant total feed flow rate requirement. Furthermore, the decoking interval is also fixed on the basis of prior operation experiences so as not to exceed the maximum values of the tube skin temperature and pressure drop. In our furnace system, under a deterministic assumption of the coke growth rate, a simple heuristic scheduling method can be proposed to compare the performance of our solution method with the current practice. Because there is no information to determine the feed flow rate in this heuristic method, we optimized the constant operating period under a fixed flow rate assumption. By trying different constant periods to the fixed flow rate strategy, we determine the 56-day constant operating period to be optimal. In this study,

we consider the decoking point as well as the decoking interval and flow rates as decisions. The simple heuristic strategy is applied to the decoking schedule problem with 20 sets of random initial coke thicknesses. The decoking scheduling solutions by the sequential strategy are compared with solutions by the simple heuristic decoking strategy. GAMS is used for the sequential strategy and MATLAB for the simple heuristic strategy. To see more than one decoking interval for all furnaces, a 120-day time horizon is applied. Figure 7 compares the profits achieved by the sequential strategy vs the simple heuristic strategy. The profits achieved by the sequential strategy are consistently higher than those by the heuristic strategy. Figure

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5745

Figure 9. Gantt chart of decoking schedule of case 5 by the simple heuristic and the sequential strategies.

Figure 10. Accumulated profit difference of the sequential strategy and the heuristic strategy (case 5).

8 shows the profit vs decoking interval for case 5 in Figure 7. The simple heuristic method has various profits, depending on the decoking intervals. From Figure 8, the 56-day decoking interval case gives a maximum profit of 1537.6 in the given 120-day horizon. With the same initial coke thickness distribution, the sequential method gives a profit of 1537.9, which is slightly better. To compare the solutions in detail, the Gantt chart for the decoking schedules of case 5 by the heuristic strategy and the sequential strategy is shown in Figure 9. In this Gantt chart, the decoking schedule of the sequential strategy has a trend that is very similar to that of the heuristic solution, but the decoking interval of the sequential strategy is slightly shorter than that of the heuristic solution. This difference arises from optimizing the decoking schedule and feed flow rate of the furnaces. As previously explained, the furnace recovers its production efficiency after the coke removal by the periodic decoking process. Even though the decoking process sacrifices valuable production days, the decoking increases the total production amount of the furnace by allowing the productivity to recover. Additionally, the operation time loss is compensated by the flow rate increase in the sequential optimization strategy. The average profit of a furnace during normal operation in the sequential strategy is 0.09% higher than that in the heuristic method. Due to the characteristic of our problem, the effect of optimization by reducing one decoking increases 1/120 of profit in one furnace. Even though it is a very low increase in profit, it can have an effect because the ethylene plant has a huge production size.

To compare the production difference between the sequential and heuristic schedules, the difference in their accumulated profits is plotted in Figure 10. The accumulated profit difference fluctuates depending on the decoking schedule, but it has an ascending trend. The profit by the sequential optimization strategy becomes higher and higher than the profit by the heuristic strategy as the operation goes on. There is no significant difference in the product yields of the sequential and heuristic methods; however, their solutions have different profits due to the flow rate differences, which are shown in Figure 11. To maximize the profit in the time horizon, the sequential strategy increases the flow rate as high as possible during the early phase of a decoking interval. The increased flow rate encourages the growth of coke, and the fast coke growth decreases the decoking interval. Even though the possible production days decrease due to the more frequent decoking, the increased productivity by decoking more than makes up for the loss. 5. Conclusion Optimal scheduling of decoking of individual furances in an industrial cracking furnace system is highly desirable because the decoking operation sacrifices valuable production time. The decoking scheduling problem is difficult to solve due to the complexity of the reactions occurring in industrial cracking furnaces. Thus, a much simpler input-output model of a naphtha cracking furnace was developed by data-fitting the results from a comprehensive simulator called CRACKER. Then, an MINLP

5746

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

suffers from inevitable uncertainties, such as those in the coke growth rates and measurement errors. The required extension for application to real operations with safety constraints is to develop a reactive or proactive rescheduling strategy based on the proposed sequential solution algorithm. Acknowledgment This work was supported by the BK21 Project. J.H.L. gratefully acknowledges financial support from the U.S. National Science Foundation (CTS-#0301993). Nomenclature General Notations

Figure 11. Operating conditions and production yields of furnace 1 in case 5 by sequential and heuristic strategies (solid line, sequential solution; dotted line, heuristic solution).

was formulated on the basis of the simplified model, and three solution strategies for the MINLP were proposed. The three solution strategies were compared with each other in 20 cases with randomly chosen initial coke thickness distributions. In terms of solution quality and computational time, the sequential solution strategy of optimizing the decoking schedule through MILP and then the flow rates through NLP were significantly better than the others. This sequential solution strategy was applied to the decoking scheduling problem for a 120-day time horizon. The 120-day decoking scheduling solutions by the sequential strategy were compared with heuristic solutions. The heuristic decoking schedules were similar to the sequential strategy solutions, but the flow rates were different. Despite the decrease in the production days for the individual furnaces in the sequential strategy solution, the increased productivity afforded by the decoking process more than made up for the loss. In terms of the plant safety, the heuristic strategy is worse than the sequential strategy because the tube skin temperature and pressure drop of the heuristic strategy were seen to be closer to their limit values. The sequential strategy shows performance that is superior to the other proposed solution strategies in terms of solution quality and computational load. But to apply it to the real furnace operation, the proposed sequential strategy needs to be extended to consider uncertainties in the actual operation. Because the coke thickness inside coils cannot be measured directly during the operation, applying the scheduling solution to the real plant

a ) binary variable constraining separable programming variable, r d ) binary variable for decoking decision p ) [ton/h] production amount pˆ ) production amount multiplied by (1 - di,t) r ) intermediate discrete variable for separable programming u ) intermediate continuous variable for separable programming x ) denormalized neural network model input y ) denormalized neural network model output xˆ ) normalized neural network model input yˆ ) normalized neural network model output A, B ) vertex for piecewise linear function of p ) u2 F ) [ton/h] naphtha feed flow rate F1 ) [ton/h] naphtha feed flow rate (F1 ) Ftotal/(Nf - 1)) F2 ) [ton/h] naphtha feed flow rate (F2 ) Ftotal/Nf) Ftotal ) [ton/h] total naphtha feed flow rate to a furnace system M ) large parameter for big-M method Nf ) total number of furnaces in a furnace system obj ) objective function value Subscripts i ) furnace number j ) jth node of normalized neural network model input k ) kth node of normalized neural network model output m ) index of vertices in piecewise linear function t ) [day] discrete time coke ) coke thickness COT ) [K] coil outlet temperature C2 ) ethylene fraction C3 ) propylene fraction ∆P ) [kPa] pressure drop through cracking coil feed ) [ton/h] feed flow rate type ) feed type number Tskin ) [K] tube skin temperature steam ) [ton/h] steam flow rate Superscripts max ) maximum bound min ) minimum bound F ) superscript for MILP with fixed flow rate variables N ) superscript for sequential strategy variables S ) superscript for MILP with piecewise linear functions using separable programming variables S1 ) variable for separable programming in MILP with piecewise linear functions S2 ) variable for separable programming in MILP with piecewise linear functions

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5747

U ) upper limit Literature Cited (1) Wauters, S.; Marin, G. B. Kinetic modeling of coke formation during steam cracking. Ind. Eng. Chem. Res. 2002, 41, 2379-2391. (2) Jain, V.; Grossmann, I. E. Cyclic scheduling of continuous parallelprocess units with decaying performance. AIChE J. 1998, 44, 1623-1636. (3) Dente, M.; Ranzi, E.; Goossens, A. G. Detailed prediction of olefin yields from hydrocarbon pyrolysis through a fundamental simulation model (spyro). Comput. Chem. Eng. 1979, 3, 61-75. (4) Joo, E.; Lee, K.; Lee, M.; Park, S. Cracker - a pc based simulator for industrial cracking furnaces. Comput. Chem. Eng. 2000, 24, 15231528. (5) Joo, E.; Park, S.; Lee, M. Pyrolysis reaction mechanism for industrial naphtha cracking furnaces. Ind. Eng. Chem. Res. 2001, 40, 2409-2415. (6) Edwin, E. H.; Balchen, J. G. Dynamic optimization and production planning of thermal cracking operation. Chem. Eng. Sci. 2001, 56, 989997. (7) Kopinke, F. D.; Zimmermann, G.; Reyniers, G. C.; Froment, G. F. Relative rates of coke formation from hydrocarbons in steam cracking of

naphtha. 2. paraffins; naphthenes; mono-, di-, and cycloolefins; and acetylenes. Ind. Eng. Chem. Res. 1993, 32, 56-61. (8) Kopinke, F. D.; Zimmermann, G.; Reyniers, G. C.; Froment, G. F. Relative rates of coke formation from hydrocarbons in steam cracking of naphtha. 3. aromatic hydrocarbons. Ind. Eng. Chem. Res. 1993, 32, 26202625. (9) Reyniers, G. C.; Froment, G. F.; Kopinke, F. D.; Zimmermann, G. Coke formation in the thermal cracking of hydrocarbons. 4. modeling of coke formation in naphtha cracking. Ind. Eng. Chem. Res. 1994, 33, 25842590. (10) Williams, H. P. Model Building in Mathematical Programming, 3rd ed.; John Wiley and Sons Ltd.: Chichester, England, 1991. (11) Furman, K. C.; Sahinidis, N. V. A critical review and annotated bibliography for heat exchanger network synthesis in the 20th century. Ind. Eng. Chem. Res. 2002, 41, 2335-2370.

ReceiVed for reView February 2, 2005 ReVised manuscript receiVed May 19, 2006 Accepted May 25, 2006 IE050129N