Emission Constrained Dynamic Scheduling for Ethylene Cracking

Article pubs.acs.org/IECR

Emission Constrained Dynamic Scheduling for Ethylene Cracking Furnace System Shujing Zhang,† Sujing Wang,‡ and Qiang Xu*,† †

Dan F. Smith Department of Chemical Engineering and ‡Department of Computer Science Lamar University, Beaumont, Texas 77710, United States ABSTRACT: The scheduling of an ethylene cracking furnace system should not only focus on the plant profitability but also address the concern of environmental sustainability. On the basis of the previous study (Zhao et al.1), this Article has developed a new scheduling model for the optimal operation of the cracking furnace system in ethylene plants by considering more practical operating features, such as nonhomogeneous last-batch ending times and furnace load makeups, as well as the background air-quality conscious decoking. With the newly developed model, the simultaneous furnace shutdown can be successfully avoided at the end of the scheduling time horizon; meanwhile, certain manufacturing capacity losses and upsets can be alleviated during furnace decoking operations. Furthermore, the obtained scheduling solution will be responsible for background air-quality concerns by wisely allocating furnace decoking emission peaks into time windows that would cause less adverse environmental impacts. The efficacy of the developed model has been demonstrated by a comprehensive case study.

1. INTRODUCTION Ethylene is the most important organic compound of petrochemical/chemical industry with the largest bulk production that significantly outweighs other petrochemical products. A broad spectrum of ethylene derivatives, such as ethylene oxide, vinyl acetate, linear alcohols, ethylene dichloride, and high/low density polyethylene, play extremely important roles in people’s daily life. According to the IHS report in 2014,2 during the year of 2009 through 2014, the consumption rate of ethylene climbed at a yearly rate of nearly 4.5%, with a capacity increase at an annual rate of 3.5%. Conceivably, with both increments of the global population and economy, ethylene demand and manufacturing will be continuously growing worldwide. The steam thermal cracking is still the dominating way for ethylene manufacturing in current ethylene plants. In an ethylene plant, the cracking operation serves as the first operating sector, which largely determines downstream product yields and energy consumption of the entire plant. In such an operating section, multiple cracking furnaces are working in parallel to crack different types of hydrocarbon feedstock into the olefin enriched cracked gas, mostly ethylene and propylene. After that, the cracked gas will go through operating sections of quenching, compression, chilling, and separation to recover products, such as ethylene and propylene, from the cracked gas. As well established, thermal reactions inside ethylene cracking furnaces are highly endothermic and accompanied by many side reactions. Thus, both high reaction temperature and short resident time are generally required. During such thermal reactions, the carbon coke will be continuously formed and accumulated on furnace radiant tubes. The accumulated coke adds to the resistance of heat transfer and finally leads to © XXXX American Chemical Society

prohibitively high TMT (tube metal temperature) and pressure drop, which degrades both the reaction selectivity and poses danger to furnaces. Therefore, periodic decoking operations with steam/air mixtures have to be conducted for furnace regenerations. From the plant operating point of view, decoking operations will incur upsets to the plant operation and, thus, potentially affect the plant operability and productivity. Therefore, ethylene plants will normally follow the rule of single furnace decoking at a time. Nowadays, with the increasingly tighter margins and volatile feedstock market, ethylene plants are encouraged to diversify their feedstock portfolio and increase their operation resilience in face of raw material uncertainties. Moreover, the shale gas booming has brought renewed competitive edge to petrochemical/chemical industry and spurred the recent wave of cracker expansions worldwide. There has never been a better time for ethylene producers to review the importance of furnace system scheduling, which could provide cost-effective optimal decisions, such as how to optimally allocate multiple feeds to different types of furnaces, and determine the optimum processing time of each batch operation, as well as the decoking schedule of the entire furnace system. In light of the underlying opportunities for ethylene cracking furnace system, there have been some research endeavors devoted to optimizing its operations for the past decade. Jain and Grossmann3 developed an MINLP (mixed integer nonlinear programming) for cyclic scheduling of multiple Received: Revised: Accepted: Published: A

July 24, 2016 October 16, 2016 January 17, 2017 January 17, 2017 DOI: 10.1021/acs.iecr.6b02822 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Graphical representation of the scheduling problem for a cracking furnace system.

Facing the challenge of volatile feed supply, Zhao et al.1 developed a dynamic scheduling strategy to dynamically and seamlessly generate reschedules to accommodate new feed deliveries. Recently, Wang et al.9 developed a synchronized model to integrate upstream naphtha inventory planning and downstream furnace scheduling, which involves a subcooperative model to achieve alignment between these two subsections. As noted, all of the previous works have been centered on maximizing the average net profits for the scheduling operations of ethylene cracking furnaces. However, several practical concerned aspects have not drawn sufficient attentions yet. To temporally mitigate the capacity loss and process upset during a decoking operation, plant schedulers usually tune up the feed load of other working furnaces to make up for the capacity loss, which also helps reduce the decoking upsets to the downstream operations. Moreover, nonsimultaneous decoking operations should not only be maintained for intermediate batch operations but also for the last batch operation, which necessitates the realization of nonhomogeneous ending times between different furnaces. Another concern related to the furnace scheduling but still lacking systematic studies is the emission issue during furnace decoking operations. Large amounts of emissions such as nitrogen oxides (NOx), carbon monoxide (CO), and particulate matter (PM) in the flue gas are generated from ethylene cracking furnaces during decoking operations.10 These emissions certainly would cause negative impacts to local air quality. For example, it has been reported that CO and NOx can contribute to the ground-level ozone production through photolysis reactions,11 the high concentration of which can seriously hurt human lungs, airways and other tissues under hot and little-wind weather.12 In an attempt to mitigate the formation of ground level-ozone, decoking operations are better scheduled when there is less solar radiation which serves as the incentive for photolysis reactions. The major source of

feeds on parallel ethylene furnaces, where the objective was to maximize the average net profits through defining an exponential decaying yield model. Also, a solving algorithm was proposed to reach the global optimality. However, their model was oversimplified in the sense that it enforced the same processing duration under the same type of feed cracked in the same furnace. Most importantly, nonsimultaneous decoking constrains were not incorporated into the model. On the basis of the previous work, Schulz et al.4 proposed an extended model with recycled ethane as the internal feed. The recycling mean value was estimated through the addition of a simplified plant model. The developed MINLP model employed a timedependent empirical variable indicating coil internal roughness, which determined whether a furnace needed to be shut down for decoking operation. However, this model still ignored the nonsimultaneous decoking constraints and only considered a single type of feedstock. Lim et al.5 put forth a furnace scheduling model which introduced a neural network-based cracker model deduced from dynamic data of product yields, S/H (steam to hydrocarbon) ratio, COT (coil outlet temperature), and coke thickness to predict the time-varying process variables, such as ethylene/propylene fraction, tube skin temperature, pressure drop, and coke thickness. In their work, three alternative solution strategies were investigated and compared for handling this MINLP problem. As a further study, they proposed a proactive scheduling model6 which triggered re-estimation of model parameters whenever the gap between model prediction and actual measurement exceeded a threshold value. The merit of the proactive method has been manifested with comparison to the heuristic and reactive methods. However, their model did not consider multiple feedstock. In 2009, Liu et al.7 developed a cyclic scheduling model with consideration of multiple feedstock, multiple furnaces, and nonsimultaneous decoking. Later, Zhao et al.8 introduced new feature of secondary ethane cracking to enhance the applicability of scheduling results. B

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

Article

Industrial & Engineering Chemistry Research PM emissions from an ethylene plant stems from decoking operations of ethylene furnaces,10 for which the uncontrolled intermittent emissions are 41.6 lb/h given an ethylene plant with a capacity of 1.1 × 109 lb/yr, and the controlled emissions are 0.9 lb/h. The fine PM usually contains microscopic solids or liquids that can penetrate into human lungs and leads to serious health consequences. On the basis of studied PM sources, the background PM concentration will reach the peak value during day time with intensive solar radiation, industrial and human activities.13 To help mitigate PM level, it is advised to schedule plant decoking operations during the time periods with lower level of solar radiations, so that PM peak emissions from furnace decoking operations can avoid the background PM emission peaks. In this Article, a new dynamic scheduling model for cracking furnace system operations has been developed by considering the aforementioned operational and environmental concerns. Dynamic scheduling is preferred over cyclic scheduling as it is more closely aligned with practical plant operations, and plant personnel can readily specify the scheduling horizon of their interest. With the new model, simultaneous furnace shutdowns can be avoided at the end of the scheduling time horizon, and manufacturing capacity loss and process upsets because of some furnace decoking can be alleviated by adjusting the feed load of active furnaces. Meanwhile, the obtained scheduling solution will be responsible for background air-quality concerns by wisely allocating furnace decoking emission peaks into time windows that would cause less adverse environmental impacts. The efficacy of the new scheduling model has been demonstrated by a comprehensive case study.

Assumptions: (1) Each batch operation in any furnaces can process only one type of feed; (2) when a furnace is under decoking operation, every leftover furnaces will additionally process a certain amount of makeup feed with their fixed production yield rates; and (3) the emission profile during decoking operation roughly follows a normal distribution with the emission peak falling at the middle point of the decoking duration.

3. DYNAMIC SCHEDULING MODEL 3.1. Objective Function. The objective function of the scheduling model is to maximize the average net profit per day over the customized scheduling horizon, which is shown in eq 1. Altogether there are three terms involved, the first of which represents the sales revenue from various products, and the second term denotes the feedstock purchase cost, while the last term indicates the decoking cost. Note that SPl,i,j,k is the total amount of product l produced during batch k of furnace j cracking feedstock i and SCi,j,k refers to the total amount of feedstock i consumed by furnace j during batch k, which will be calculated in the following constraints: max profits = ( ∑ ∑ ∑

∑ PrSP l l ,i,j,k

l ∈ L i ∈ Ij j ∈ J k ∈ K

−

∑ ∑ ∑ CrSC i i , j , k − ∑ ∑ ∑ Csi , jyi , j , k )/ H i ∈ Ij j ∈ J k ∈ K

i ∈ Ij j ∈ J k ∈ K

(1)

3.2. Mass Balance Constraints. Because the product yield is decaying with respect to cracking time and usually follows a nonlinear trend, an exponential formula is adopted to describe the dynamic product yield during cracking process. As shown from eq 2, the exponential yield model is characterized by three coefficients, where ci,j,l is the constant coefficient of product l’s yield model for feed i cracked in furnace j, ai,j,l is the preexponential factor, and bi,j,l is the power factor. Using this exponential yield model, the total amount of products during any batch operation can be calculated, as shown in the first term of eq 2. As observed, the second part in eq 2 represents another source of product yield from load increase, where pj is the percentage of load increase, and ydi,j,l is the constant yield parameter during the load increase period which takes the baseline product yield in the yield function with the time variable in the exponential term to be zero since the decaying rate is relatively small and the time duration of load increase (Tsmi,j,k) only takes a small portion of the total batch processing time. To further corroborate this assumption, the recalculated profits after considering the time effect in the yield function for the load increase period are 0.03% more than the profits generated by neglecting the time effect. Likewise, eq 3 calculates the total amount of feedstock i consumed by furnace j during batch operation k, where the first part corresponds to the feedstock consumed by a normal batch operation and the second part comes from the load increase.

2. PROBLEM STATEMENT To facilitate the understanding of the scheduling problem for a cracking furnace system, Figure 1 gives a detailed graphical illustration of batch processing time, decoking operation, nonsimultaneous decoking, scheduling horizon, load increase, and dynamic scheduling. Various types of feedstock are assigned to different batch operations as labeled by alphabetic letters (A−D). Decoking operations among different furnaces are forbidden to be overlapping each other. As seen from Figure 1, the ending time of the last batch operation of each furnace is allowed to cross the customized scheduling time horizon. Additionally, whenever a furnace is under a decoking operation, other furnaces will pick up additional feed loads to make up for the capacity loss and mitigate the downstream upsets. Note that the developed scheduling model should be implemented in a dynamic approach, which can initiate a rescheduling whenever it is necessary. Given information: (1) Average feedstock supply flow rate, (2) feasible types of feed that can be processed in each furnace, (3) cracking yield information on different types of feedstock in different furnaces, (4) decoking time expense for each furnace, (5) lower/upper bounds on batch processing time, (6) decoking cost, (7) feedstock purchase cost and product sale price, and (8) customized scheduling time horizon. Information to be determined: (1) The number of batches assigned for each furnace, (2) feed type processed in each batch operation, (3) the starting and ending time points for each batch operation, (4) makeup feed amount of every leftover furnaces when a furnace is under decoking operation, and (5) the sequence of decoking operations for the entire furnace system.

SPl , i , j , k =

∫0

ti , j , k

FLiin, j(ci , j , l + ai , j , lebi ,j ,lt ) dt

+ pj FLiin, jydi , j , lTsmi , j , k , ∀ i ∈ Ij ; j ∈ J ; l ∈ L ; k ∈ K (2)

SCi , j , k =

FLiin, jti , j , k

+

pj FLiin, jTsmi , j , k ,

∀ i ∈ Ij ; j ∈ J ; k ∈ K (3)

C


Article


designated batch slot. Equation 15 denotes that yej,k is equal to 1 only if the specific batch k is the last batch operation of furnace j.

According to the principle of mass conservation, the total amount of products will be equal to the total amount of feedstock consumed as shown in eq 4. Also, if a certain batch operation is not activated, then the amounts of related products and feedstock during that batch will be enforced to zero, as formulated by eqs 5 and 6.

ylj , k =

(5)

up SCi , j , k ≤ tiup , j Fi , j yi , j , k , ∀ i ∈ Ij ; j ∈ J ; k ∈ K

(6)

ylj , k ≥

(8)

3.3. Integrality Constraints for Batch Selection. The total number of scheduled batches for each furnace is unknown beforehand, thus, some logical constraints are configured to define the batch selection sequence. Equation 9 indicates that a single batch operation can process at most one type of feed, and eq 10 enforces that the first batch of any furnaces will be selected for sure; while eq 11 defines the batch sequence that the earlier batch slots have to be activated before the later ones can be scheduled.

Sj , k ≤ H + UB(1 − ylj , k), ∀ j ∈ J ; k ∈ K

j ∈ J; k ∈ K

(18)

∑ yi ,j ,k = 1, ∀ j ∈ J ; k ∈ K , k = 1

H − Ej , k ≤

(10)

≥

∑ τi ,jyi ,j ,k

+ UByej , k + UB(1 − ylj , k), ∀

i ∈ Ij

∑ yi ,j ,k′ , ∀ j ∈ J ; k , k′ ∈ K , k < k′

j ∈ J; k ∈ K

i ∈ Ij

(11)

(19)

3.6. Feed Allocation Constraints. Since only one furnace has the capability to process ethane feed and the other furnaces will crack other types of feedstock, the following feed allocation constraints are established:

3.4. Integrality Constraints for the Last Batch Operation. To allow the last batch operation to pass across the customized scheduling horizon, two new binary variables (ylj,k and yej,k) are defined. ylj,k is used to identify the last active batch operation of a furnace, which is equal to 1 if the k-th batch operation of furnace j counts as the last one; and yej,k indicates whether or not the last batch operation of furnace j exceeds the customized scheduling horizon. Thus, eq 12 implies that if ∑i ∈ I yi , j , k is equal to 1 and ∑i ∈ I yi , j , k + 1 is equal to j

(17)

Ej , k − H ≤ UByej , k + UB(1 − ylj , k), ∀ j ∈ J ; k ∈ K

(9)

i ∈ Ij

(16)

Ej , k − H ≥ −UB(2 − yej , k − ylj , k) + 1E − 5, ∀

≤ 1, ∀ j ∈ J ; k ∈ K

i ∈ Ij

(15)

3.5. Timing Constraints for the Last Batch Operation. after defining the binary variables, it is necessary to utilize them to restrict the timing variables. Equation 16 enforces that the starting time variable of the last active batch operation should be less than the customized scheduling horizon. Equations 17 and 18 define that only when both ylj,k and yej,k are equal to 1, the ending time of the last batch operation will pass the boundary of customized scheduling horizon, otherwise, the last batch operation will be bounded within the customized scheduling horizon. Equation 19 implies that for the scenario that the last batch operation is kept within the customized scheduling horizon, the allowance between the last batch ending time and the scheduling horizon should not be bigger than the corresponding decoking time.

(7)

k∈K

(14)

yej , k ≤ ylj , k , ∀ j ∈ J ; k ∈ K

∑ (FLiin,j ∑ ti ,j ,k), ∀ i ∈ I

Gi ≤ (Fiup − Filo)H , ∀ i ∈ I

i ∈ Ij

∑ yi ,j ,k , ∀ j ∈ J ; k = |K | i ∈ Ij

Equations 7 and 8 are established with respect to the fact that the total amount of each feedstock consumed by all cracking furnaces should be kept within a supply capability. As shown, eq 7 calculates the consumption amount of feed i during the customized scheduling horizon H, where Gi represents the excess amount of feedstock i which goes above the lower bound. And eq 8 imposes an upper bound on Gi.

∑ yi ,j ,k

(13)

k∈K

l∈L

∑ yi ,j ,k

i ∈ Ij

∑ ylj ,k = 1, ∀ j ∈ J

(4)

∑ SPl ,i ,j , k ≤ tiup, j Fiup,j yi ,j ,k , ∀ i ∈ Ij ; j ∈ J ; k ∈ K

j∈J

∑ yi ,j ,k+ 1 , ∀ j ∈ J ; k ∈ K , k ≠ |K | (12)

i ∈ Ij l ∈ L

FiloH + Gi =

−

i ∈ Ij

∑ SCi ,j ,k = ∑ ∑ SPl ,i ,j ,k i ∈ Ij

∑ yi ,j ,k

yi ,NC, k = 0, ∀ i = 1, ..., NF − 1; k = 1, ..., NB

(20)

yNF, j , k = 0, ∀ j = 1, ..., NC − 1; k = 1, ..., NB

(21)

3.7. General Timing Constraints. There are some straightforward relations between the timing variables, which are reflected in eqs 22−24. Equation 22 indicates that the starting time variable of current batch is equal to the ending time variable of the previous batch plus its decoking duration. Equation 23 constrains that the ending time variable of a batch operation is equal to its starting time variable plus the batch duration. Equation 24 relates the batch duration with the binary variable for batch selection, which means that the batch duration of a not-selected batch is forced to zero.

j

0, then the k-th batch operation must be the last active batch of furnace j (ylj,k = 1). Otherwise, if both ∑i ∈ I yi , j , k and j

∑i ∈ I yi , j , k + 1 are equal to 1, the k-th batch operation will not j

be the last one of furnace j. Equation 13 enforces the logic that each furnace will have only one batch as the last active batch. Equation 14 indicates that if a furnace uses all the assigned batch slots, the last active batch operation will be the last D


Article


i ∈ Ij

Ej , k = Sj , k +

3.10. Integrality Constraints for Load Increase. Before calculating the time duration for load increase, a new binary variable (yvj,k,j′,k′) is defined that denotes whether the k-th decoking operation of furnace j falls within the k′-th batch operation of furnace j′. Equation 30 indicates that if the k-th decoking operation of furnace j is nonoverlap behind the k′-th decoking operation of furnace j′ (xj,k,j′,k′ = 1), the k-th decoking operation of furnace j will not fall within the k″-th batch operation of furnace j′ given that k″ ≤ k′, which suggests yvj,k,j′,k″ = 0. As shown in eq 31, if the k-th decoking operation of furnace j is nonoverlap ahead of the k′-th decoking operation of furnace j′ (xj,k,j′,k′ = 0), the k-th decoking operation of furnace j will not fall within the k″-th batch operation (k′ < k″) of furnace j′, which means yvj,k,j′,k″ = 0. Figure 2 depicts the physical meaning of these two equations, where Figure 2a corresponds to eq 30, and Figure 2b reflects eq 31.

∑ τi ,jyi ,j ,k− 1 , ∀ j ∈ J ; k ∈ K , k ≠ 1

Sj , k = Ej , k − 1 +

(22)

∑ ti ,j ,k , ∀ j ∈ J ; k ∈ K i ∈ Ij

(23)

tilo, j yi , j , k ≤ ti , j , k ≤ tiup , j yi , j , k , ∀ i ∈ Ij ; j ∈ J ; k ∈ K

(24)

3.8. Nonsimultaneous Decoking Constraints. As stated earlier, the decoking operations of different furnaces are not allowed to be overlapping with each other. In order to achieve this idea, a binary variable xj,k,j′,k′ is defined. It equals 1, if the kth decoking operation of furnace j is nonoverlap behind the k′th decoking operation of furnace j′; otherwise, it is equal to 0 if the k-th decoking operation of furnace j is nonoverlap ahead of the k′-th decoking operation of furnace j′. Equation 25 corresponds to the scenario that xj,k,j′,k′ = 0, and eq 26 indicates the other scenario that xj,k,j′,k′ = 1. Ej , k +

∑ τi , jyi ,j , k

− UBxj , k , j ′ , k ′ − UB(1 −

i ∈ Ij

∑ yi ,j , k ) i ∈ Ij

≤ Ej ′ , k ′ + UB(1 −

∑

yi ′ , j ′ , k ′), ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K

i ′∈ Ij ′

(25)

Ej , k + UB(1 −

∑ yi ,j ,k ) ≥ Ej ′ ,k ′ + ∑ i ∈ Ij

− UB(1 −

∑

τi ′ , j ′yi ′ , j ′ , k ′

i ′∈ Ij ′

yi ′ , j ′ , k ′) − UB(1 − xj , k , j ′ , k ′), ∀ j , j′

i ′∈ Ij ′

∈ J , j ≠ j′ ; k , k ′ ∈ K

(26)

3.9. Additional Logic Constraints. There are some logic constraints constructed to reduce the solution space as shown in eqs 27−29. Equation 27 indicates that if the k″-th decoking operation of furnace j is nonoverlap ahead of the k′-th decoking operation of furnace j′ (xj,k″,j′,k′ = 0), then the earlier decoking operation k (k < k″) of furnace j will definitely be nonoverlap ahead of the k′-th decoking operation of furnace j′ (xj,k,j′,k′ = 0). Similarly, eq 28 implies that if the k-th decoking operation of furnace j is nonoverlap behind the k″-th decoking operation of furnace j′ (xj,k,j′,k″ = 1), it must be nonoverlap behind the earlier decoking operations (k″ > k′, xj,k,j′,k′ = 1). Equation 29 suggests that if the k-th decoking operation of furnace j is sandwiched by two adjacent decoking operations (k′-th and k′-th +1) of another furnace j′, the k′-th +1 batch operation of furnace j′ must have been selected (∑i ∈ I yi , j ′ , k ′+ 1 = 1).

Figure 2. Illustration for relation between x and yv as shown in eqs 30 and 31: (a) yvj,k,j′,k″ ≤ 1 − xj,k,j′,k′ and (b) yvj,k,j′,k″ ≤ xj,k,j′,k′.

Equation 32 implies that if the k-th decoking operation of furnace j lies between two adjacent decoking operations (k′, k′ + 1) of furnace j′, the k-th decoking operation of furnace j must fall within the k′-th + 1 batch operation of j′ as shown in Figure 3. Equations 33 and 34 restrict that if a batch is not selected (∑i ∈ I yi , j , k = 0), its accompanied decoking operation will not j

fall within any batch operation of any other furnaces (yvj,k,j′,k′ = 0) and vice versa, if a batch is not selected (∑i ∈ I yi , j ′ , k ′ = 0),

j′

xj , k , j ′ , k ′ ≤ xj , k ″ , j ′ , k ′ , ∀ j , j′ ∈ J , j ≠ j′; k , k′, k″ ∈ K , k < k″

j′

(27)

there will not be any decoking operations residing within its batch operation (yvj,k,j′,k′ = 0).

xj , k , j ′ , k ′ ≥ xj , k , j ′ , k ″ , ∀ j , j′ ∈ J , j ≠ j′; k , k′, k″ ∈ K , k″ > k′ (28)

xj , k , j ′ , k ′ − xj , k , j ′ , k ′+ 1 = 1 ⇒

yvj , k , j ′ , k ″ ≤ 1 − xj , k , j ′ , k ′ , ∀ j , j′ ∈ J , j ≠ j′; k , k′, k″

∑ yi ,j ′ ,k′+ 1 = 1, ∀ j , j′

∈ K , k″ ≤ k′

(30)

i ∈ Ij ′

yvj , k , j ′ , k ″ ≤ xj , k , j ′ , k ′ , ∀ j , j′ ∈ J , j ≠ j′; k , k′, k″ ∈ K , k′

∈ J , j ≠ j′ ; k , k ′ ∈ K , k ′ < | K |

< k″

∑ yi ,j ′ ,k′+ 1 ≥ xj ,k ,j ′ ,k′ − xj ,k ,j ′ ,k′+ 1, ∀ j , j′ ∈ J , j ≠ j′; k , k′

yvj , k , j ′ , k ′+ 1 ≥ xj , k , j ′ , k ′ − xj , k , j ′ , k ′+ 1 , ∀ j , j′ ∈ J , j ≠ j′; k , k′

i ∈ Ij ′

∈ K , k ′ < |K |

(31)

∈ K , k ′ ≠ |K |

(29) E

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

Article

Industrial & Engineering Chemistry Research Tcl j , k , j ′ , k ′ ≤

∑ τi ,jyi ,j ,k , ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K i ∈ Ij

(38)

Tcl j , k , j ′ , k ′ ≥

∑ τi ,jyi ,j ,k

− UB(1 − yvj , k , j ′ , k ′), ∀ j , j′ ∈ J , j

i ∈ Ij

≠ j′ ; k , k ′ ∈ K

(39)

Tcl j , k , j ′ , k ′ ≤ UByvj , k , j ′ , k ′, ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K (40)

Tcl j , k , j ′ , k ′ ≤ UB ∑ yi , j ′ , k ′ , ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K i ∈ Ij ′

(41)

Equations 42 and 43 are employed to calculate the starting time variable of load increase period within the k′-th batch operation of furnace j′ caused by the k-th decoking operation of furnace j. Similarly, eqs 44 and 45 are built for the ending time variable. In order to calculate the gain in feedstock and products associated with the load increase period, a new time variable Tsmi,j,k is established, which represents the aggregated time that the k-th batch operation of furnace j has a raised feed load, as shown in eqs 46 and 47. Equation 48 restricts that if a batch operation is not scheduled, its total time for load increase will be enforced to 0.

Figure 3. Illustration for the relation between x and yv as shown in eq 32: (a) yvj,k,j′,k′+1 ≥ xj,k,j′,k′ − xj,k,j′,k′+1 (if xj,k,j′,k′ − xj,k,j′,k′+1 = 1, then yvj,k,j′,k′+1 = 1) and (b) yvj′,k′,j,k+1 ≥ xj,k+1,j′,k′ − xj,k,j′,k′ (if xj,k+1,j′,k′ − xj,k,j′,k′ = 1, then yvj′,k′,j, k+1 = 1).

yvj , k , j ′ , k ′ ≤

∑ yi ,j ,k , ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K i ∈ Ij

yvj , k , j ′ , k ′ ≤

(33)

Tic sj , k , j ′ , k ′ ≤ Ej , k + UB(1 − yvj , k , j ′ , k ′), ∀ j , j′ ∈ J , j

∑ yi ,j ′ ,k′ , ∀ j , j′ ∈ J , j ≠ j′; k , k′ ∈ K i ∈ Ij ′

(34)

≠ j′ ; k , k ′ ∈ K

With aforementioned constraints, the value of yvj,k,j′,k′ cannot be determined when a certain decoking operation of a furnace falls within the first batch operation of another furnace. Therefore, a new binary variable zj,k,j′,1 is proposed to help identify this situation, which assumes the value of one when the ending time of the k-th batch operation of furnace j is not less than the starting time variable of the first batch operation of another furnace j′; otherwise, it is equal to zero. Equations 35 and 36 are constructed to fulfill this idea. Equation 37 enforces that if zj,k,j′,1 = 1 and xj,k,j′,1 = 0 simultaneously, the k-th decoking operation of furnace j must fall within the first batch operation of furnace j′ (yvj,k,j′,1 = 1). Ej , k − Sj ′ ,1 ≤ UBzj , k , j ′ ,1 , ∀ j , j′ ∈ J , j ≠ j′; k ∈ K

(42)

Tic sj , k , j ′ , k ′ ≥ Ej , k − UB(1 − yvj , k , j ′ , k ′), ∀ j , j′ ∈ J , j ≠ j′ ; k , k ′ ∈ K Tic jf, k , j ′ , k ′ ≤ Ej , k +

∑ τi ,jyi ,j ,k

(43)

+ UB(1 − yvj , k , j ′ , k ′), ∀ j , j′

i ∈ Ij

∈ J, j ≠ j′ ; k , k ′ ∈ K Tic jf, k , j ′ , k ′ ≥ Ej , k +

∑ τi ,jyi ,j ,k

(44)

− UB(1 − yvj , k , j ′ , k ′), ∀ j , j′

i ∈ Ij

(35)

∈ J, j ≠ j′ ; k , k ′ ∈ K

Ej , k − Sj ′ ,1 ≥ UB(zj , k , j ′ ,1 − 1), ∀ j , j′ ∈ J , j ≠ j′; k ∈ K (36)

Tsmi , j ′ , k ′ ≤

(45)

∑ ∑ Tclj ,k ,j ′ ,k′ + UB(1 − yi ,j ′ ,k′), ∀ j∈J ,j≠j′ k∈K

yvj , k , j ′ ,1 ≥ zj , k , j ′ ,1 − xj , k , j ′ ,1 , ∀ j , j′ ∈ J , j ≠ j′; k ∈ K

i ∈ I j ′ ; j′ ∈ J ; k ′ ∈ K

(37)

Tsmi , j ′ , k ′ ≥

3.11. Timing Constraints for Load Increase. On the basis of the definition of yvj,k,j′,k′, the time duration of load increase (Tclj,k,j′,k′) can be determined. Accordingly, eqs 38 and 39 are built to calculate Tclj,k,j′,k′ of each load increase when there is a decoking operation. Note that Tclj,k,j′,k′ represents the time duration of load increase within the k′-th batch operation of furnace j′ due to the k-th decoking operation of furnace j. And eqs 40 and 41 enforce that Tclj,k,j′,k′ is equal to 0 if the k-th decoking operation of furnace j does not fall within the k′-th batch operation of furnace j′.

(46)

∑ ∑ Tclj ,k ,j ′ ,k′ − UB(1 − yi ,j ′ ,k′), ∀ j∈J ,j≠j′ k∈K

i ∈ I j ′ ; j′ ∈ J ; k ′ ∈ K Tsmi , j , k ≤ ti , j , k , ∀ i ∈ Ij ; j ∈ J ; k ∈ K

(47) (48)

3.12. Constraints for Decoking Emission Time Window. To avoid the daytime period, which has more sunlight to promote the formation of ground-level ozone and possesses higher background PM level, the emission peak of F


Article


batch ending time (Ej,Kj) plus its corresponding decoking time (∑i ∈ I τi , jyi , j , K ). This leftover time (Tlj) actually marks the

decoking operations should be scheduled inside some designated time windows (roughly from the late night of a day to the early morning of the next day). A new binary variable xej,k,d is defined which implies whether or not the k-th decoking operation of furnace j lies within day d. Equation 49 is used to express the logic that if a batch operation k is selected, its decoking operation must fall within a specific day. Subsequently, eqs 50 and 51 identify the particular day that a decoking operation falls into. After capturing which day the decoking operation resides, eqs 52 and 53 are used to schedule the decoking emission peak within a designated time widow. As seen, the left-hand side terms of eqs 52 and 53 calculate the time point of the decoking emission peak, and Atw/Btw correspond to the boundary points of a designated emission

j

Kj = {k|min(t R − Sj , k) ≥ 0}, ∀ j ∈ J Tl j = Ej , K j − t R +

(49)

≥ Atw + (ord(d) − 1)

i ∈ Ij

Ej , k +

∑ 0.5τi ,jyi ,j ,k

(52)

≤ Btw + ord(d) + UB(1 − xej , k , d ), ∀

i ∈ Ij

j ∈ J; k ∈ K; d ∈ D

furnace j

leftover feed type

Tlj (day)

1 2 3 4 5 6

Fb Fb Fb Fb − Fd

16 2 27 21 0 24

horizon is 110 days (H), and all related parameters are summarized in Table 2 which covers the furnace decoking time, designated feed flow rate, lower and upper bounds on batch processing time, as well as the yield parameters for furnace model. Generally speaking, ethylene yield declines with time, while propylene and ethane yields possess an opposite trend, thus an exponential yield model is deployed to capture the decaying behavior, for which the yield parameters are obtained through regression of the industrial historical data. The associated decoking and material costs as well as product prices have been given in Tables 3 and 4. As a comparison, plant heuristics have been implemented to schedule the ethylene cracking operations. The rationale behind the heuristics is that LPG has the highest ethylene yield in furnaces 2 and 5, and furnace 3 generates the highest propylene yield for LD, while naphtha has the largest supply amount. Thus, furnaces 2 and 5 are assigned to crack LPG, furnace 3 is used to process LD, and the remaining furnaces except ethane furnace 6 are allocated to crack naphtha. Because decoking operations are considered as nonprofit, they should be avoided unless it is mandatory to conduct them, which makes the scheduler to set the batch processing times as close as to the imposed upper bounds. Figure 5 shows the heuristic scheduling

(50)

(51)

− UB(1 − xej , k , d ), ∀ j ∈ J ; k ∈ K ; d ∈ D

(56)

Table 1. Initial Condition of the Case Study

Ej , k ≤ ord(d) + UB(1 − xej , k , d ), ∀ j ∈ J ; k ∈ K ; d ∈ D

∑ 0.5τi ,jyi ,j ,k

(55)

4. CASE STUDY To demonstrate the efficacy of the developed dynamic scheduling model, a comprehensive case study has been conducted with industry data. In the investigated case study, four different types of feedstock including liquefied petroleum gas (LPG, Fa), naphtha (Fb), light diesel (LD, Fc), and ethane (Fd) are processed by the cracking furnace system, which generates four kinds of products such as ethylene, propylene, ethane and other lumped products designated respectively as Pa, Pb, Pd, and Pc. The studied system comprises six cracking furnaces working in parallel, among which five are conventional furnaces processing LPG, naphtha, and LD; only one furnace is ethane furnace which cracks the freshly fed ethane. To carry out the dynamic scheduling model, the leftover time carried over from the previous schedule has to be provided as initial values for current schedule, as listed in Table 1. The given scheduling

Ej , k ≥ (ord(d) − 1) − UB(1 − xej , k , d ), ∀

Ej , k +

j

Sj ,1 = Tl j, ∀ j ∈ J

∑ xej ,k ,d ≥ ∑ yi ,j ,k , ∀ j ∈ J ; k ∈ K

j ∈ J; k ∈ K; d ∈ D

∑ τi ,jyi ,j ,K , ∀ j ∈ J i ∈ Ij

time window. Figure 4 has presented the physical meaning clearly. i ∈ Ij

(54)

k

Figure 4. Illustrative scheduling of time windows for the decoking emission peak.

d∈D

j

starting time to carry out rescheduling operations as shown in eq 56.

(53)

3.13. Initial Condition Identification for Rescheduling. Since the model is implemented dynamically, initial conditions need to be identified first. When the delivery of new feed has been confirmed, the delivery time tR along with feed type and quantity are all disclosed, which helps establish the initial conditions for rescheduling. Here, eq 54 is used to identify the batch index of the ongoing cracking operation when the new feed supply is coming in. Next, eq 55 calculates the leftover time which refers to the time duration starting from tR to the G


Article

Industrial & Engineering Chemistry Research Table 2. Parameter Values for the Cracking Furnace System furnace j

feed parameter τi,j (day)

Di,j (ton/day)

tloi,j (day)

tupi,j (day)

ai,j,Pa

ai,j,Pb

ai,j,Pd

bi,j,Pa (1/day)

bi,j,Pb (1/day)

bi,j,Pd (1/day)

ci,j,Pa

ci,j,Pb

ci,j,Pd

ydi,j,Pa

ydi,j,Pb

i

1

2

3

4

5

6

Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc Fd Fa Fb Fc

1 1 1 0 204.4 202.2 201.0 0 20 19 18 0 70 53 41 0 −0.0539 −0.0539 −0.0129 0 −0.0267 −0.0267 −0.0149 0 −0.2445 −0.0783 −0.0436 0 0.0050 0.0060 0.0070 0.0010 −0.0041 −0.0041 −0.0105 −0.0040 −0.0007 −0.0012 −0.0014 −0.0010 0.5900 0.3600 0.3000 0 0.0650 0.1600 0.1600 0 0.5701 0.1136 0.0686 0 0.5361 0.3061 0.2871 0 0.0383 0.1333 0.1451

1 1 1 0 214.9 213.7 211.1 0 23 21 19 0 72 54 43 0 −0.0500 −0.0560 −0.0140 0 −0.0243 −0.0287 −0.0149 0 −0.1675 −0.0855 −0.0485 0 0.0054 0.0062 0.0070 0.0010 −0.0043 −0.0042 −0.0100 −0.0040 −0.0010 −0.0011 −0.0013 −0.0010 0.6378 0.3698 0.3109 0 0.0610 0.1581 0.1523 0 0.4430 0.1208 0.0735 0 0.5878 0.3138 0.2969 0 0.0367 0.1294 0.1374

1 1 1 0 188.2 187.2 185.8 0 18 17 16 0 65 49 39 0 −0.0574 −0.0534 −0.0139 0 −0.0288 −0.0269 −0.0143 0 −0.3466 −0.0940 −0.0525 0 0.0053 0.0054 0.0069 0.0010 −0.0038 −0.0038 −0.0115 −0.0040 −0.0006 −0.0010 −0.0012 −0.0010 0.5637 0.3548 0.2979 0 0.0747 0.1663 0.1631 0 0.6944 0.1293 0.0775 0 0.5063 0.3014 0.284 0 0.0459 0.1394 0.1488

1 1 1 0 202.3 199.6 197.6 0 20 18 17 0 69 52 43 0 −0.0510 −0.0528 −0.0125 0 −0.0242 −0.0244 −0.0146 0 −0.2318 −0.0723 −0.0420 0 0.0049 0.0058 0.0073 0.0010 −0.0037 −0.0038 −0.0112 −0.0040 −0.0007 −0.0013 −0.0015 −0.0010 0.5987 0.3626 0.3011 0 0.0724 0.1620 0.1572 0 0.5359 0.1076 0.0670 0 0.5477 0.3098 0.2886 0 0.0482 0.1376 0.1426

1 1 1 0 212.7 210.6 209.8 0 23 22 20 0 73 56 45 0 −0.0585 −0.0584 −0.0134 0 −0.0278 −0.0257 −0.0148 0 −0.2901 −0.0671 −0.0394 0 0.0050 0.0057 0.0068 0.0010 −0.0038 −0.0043 −0.0100 −0.0040 −0.0006 −0.0014 −0.0016 −0.0010 0.6432 0.3628 0.2985 0 0.0599 0.1524 0.1513 0 0.5733 0.1024 0.0644 0 0.5847 0.3044 0.2851 0 0.0321 0.1267 0.1365

0 0 0 1 0 0 0 210.1 0 0 0 59 0 0 0 101 0 0 0 −0.0542 0 0 0 −0.0264 0 0 0 −0.1659 0.0010 0.0010 0.0010 0.0051 −0.0040 −0.0040 −0.0040 −0.0039 −0.0010 −0.0010 −0.0010 −0.0010 0 0 0 0.6067 0 0 0 0.0666 0 0 0 0.4732 0 0 0 0.5525 0 0 0

H


Article

Industrial & Engineering Chemistry Research Table 2. continued furnace j

feed i

parameter

1

Fd Fa Fb Fc Fd

ydi,j,Pd

2

0 0.3256 0.0353 0.025 0

3

0 0.2755 0.0353 0.025 0

0 0.3478 0.0353 0.025 0

4 0 0.3041 0.0353 0.025 0

5 0 0.2832 0.0353 0.025 0

6 0.0402 0 0 0 0.3073

Table 3. Furnace Decoking Cost Values furnace j

feed parameter

i

1

2

3

4

5

6

Csi,j ($)

Fa Fb Fc Fd

86700.0 87361.0 89154.0 0

92758.3 90964.4 87879.6 0

81960.6 79695.6 84239.9 0

94973.8 80087.9 81375.3 0

86733.3 82357.1 92398.7 0

0 0 0 88625.2

proposed in this paper. Figure 6 presents the optimization scheduling results for both previous and new models. It can be seen from Figure 6a that simultaneous decoking operations are not eliminated at the end of the scheduling horizon. However, the new model as shown in Figure 6b has achieved nonuniform ultimate ending times for different furnaces by allowing the ending time of the last batch to go over the ending time point of customized scheduling horizon. Scrutinizing the scheduling results reveals that feedstock selection differs a lot for furnace 2, where the previous model picks Fa for all the last three batches, but the new model chooses the series: Fb-Fb-Fa. Also, because of the incorporation of automatic load increase at the time of decoking operation, there are many red lines lying above the baseline of feed flow rate, which represents the portion of boosted yield due to automatic load increase. As a result, the overall feed load of new model is greater than that of the previous model. The normal tune-up ratio of feedstock flow

Table 4. Additional Parameter Values of Feeds and Products feed i Fa (gas) Fb (naphtha) Fc (light diesel) Fd (ethane)

Cri ($/ton)

Floi (ton/ day)

Fupi (ton/ day)

444.13 415.49

150 380

450 1140

411.01

70

444.13

110

product l

Pl ($/ton) 866.67 718.53

210

Pa (ethylene) Pb (propylene) Pc (others)

410

Pd (ethane)

444.13

356.67

results, which completely complies to the aforementioned rules. In this figure, the leftover batches carried over from the previous schedule have been labeled as green, while the batch operations in current schedule have been colored in blue. On the basis of heuristics, the resulting profit per day is $194 150, which is equivalent to 70.9 million dollars annually. The same problem has also been solved successfully with the previous model developed by Zhao et al.1 and the new model

Figure 5. Heuristic scheduling results for the ethylene production. I


Article


Figure 6. Optimal scheduling results of C2H4 productivity based on (a) the previous scheduling model and (b) the new scheduling model.

model (487.07 ton/day), as labeled by the red horizontal lines in Figure 7. Likewise, the variations of C3H6 flow rate for the previous and new models have been compared in Figure 8. As observed, there are much less flow rate upsets for the new model comparing to the previous model. In summary, the standard deviation of the previous model is 12.6 comparing to the value of 3.87 for the new model. Regarding the C3H6 average flow rate, the previous model reaches at 102.48 ton/day and the new model has a value of 124.9 ton/day. In addition, the detailed scheduling results for the previous and new model have been summarized in Tables 5 and 6 respectively in the form of cash flow. Since during the period after sunset and before sunrise, the formation of ground-level ozone is suppressed and the background particulate level is comparatively low, the characteristic emission window is specified as from 8:00 pm today to 5:00 am of the next day.

rate ranges from 10% to 20%, which can be adjusted according to plant experience. Furthermore, flow rate variations of C2H4 and C3H6 have been explicitly shown in Figures 7 and 8. From Figure 7, it can be seen that the ethylene flow rate fluctuation along the scheduling horizon based on the previous model is more drastic than that from the new model. More quantitatively, the standard deviation of flow rate based on the previous model is 59.04 and for the new model, it is 15.02. The underlying reasons should be 2-fold: (1) the introduction of automatic load increase makes up for the capacity loss during decoking operations and (2) the new model stretches the schedule beyond the customized scheduling horizon, which leads to longer batch duration within the customized scheduling horizon. However, the previous model achieves a larger average flow rate of C2H4 (522.68 ton/day) comparing to the new J


Article


Figure 7. Variation of C2H4 production flow rate based on (a) the previous scheduling model and (b) the new scheduling model. Figure 8. Variation of C3H6 production flow rate based on (a) the previous scheduling model and (b) the new scheduling model.

This implies that on a daily basis, the emission peaks of decoking operations should be scheduled later than 8:00 pm (0.83 day) of current day and meanwhile, earlier than 5:00 am (0.21 day) of the next day. Assuming that the emission profile takes a normal distribution with the peak value at the middle point of decoking duration, with the previous model, only the first batch of furnace 1, the first batch of furnace 3, and the first batch of furnace 5 have decoking emission peaks residing inside the designated emission window. However, for the new model, all decoking emission peaks have been scheduled inside the assigned emission window. With respect to profits, the new model has produced a significantly larger amount of total product revenues comparing to the previous model due to longer time horizon introduced. As seen from Tables 5 and 6, the total net profit for the previous model is 23.88 M$ and the new model achieves a total profit of 38.01 M$. Note that the first batches of all furnaces except furnace 5 come from the leftover batches. If setting the comparison base to be H (110 days), the profit harvested by the new model is 23.14 M$ which is equivalent to $210 364 per day, comparing to the profit of $217 091 per day from the previous model. This has shown that the average profit of new model within H is slightly less than that of the previous model, which should be a trade-off brought in by the emission window scheduling constraints. Recalling the average profit from

heuristic scheduling ($194 150 per day), both optimization models have accomplished higher profit gains. Certainly, the proposed scheduling model has improved downstream upsets and provided better environmental benefits. Alternative values for the length of scheduling horizon are tested to investigate the impact of scheduling horizon on model complexity and solution times. The results are summarized in Table 7. As observed, when the scheduling horizon is increased to 240 days, the required number of batches should be increased correspondingly, which leads to larger model size in terms of variables and constraints and subsequently longer solution time needed. Finally, it should be pointed out that although the development of environmental emission constraints is a good addition to furnace scheduling problems, the potential real applications based on this study should be very cautious. The constraints of emission time window employed in this study are built upon most common background and meteorological conditions. Thus, the solution results from this study are limited to less adverse environmental impacts, that said, uncertainties which stem from more complicated scenarios of K


Article

Industrial & Engineering Chemistry Research Table 5. Detailed Optimization Results Based on the Previous Scheduling Model furnace 1

2

3

4

5

6

batch number

feed type

batch starting time (day)

batch ending time (day)

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

Fb Fb Fb Fb Fb Fa Fa Fa Fb Fb Fb Fb Fb Fb Fb Fb Fa Fa Fa Fd Fd

0 16 48.63 81.27 0 2 38.09 74.17 0 27 53.63 80.27 0 21 50.13 79.27 0 36.59 73.17 0 24

15 47.63 80.27 109 1 37.09 73.17 109 26 52.63 79.27 109 20 49.13 78.27 109 35.59 72.17 109 23 109

total

ethylene revenue (M $)

propylene revenue (M $)

ethane revenue (M $)

other product revenue (M$)

feed cost (M$)

decoking cost (M$)

batch net profit (M$)

0.80 1.67 1.67 1.46 0.06 3.80 3.80 3.78 1.25 1.24 1.24 1.38 1.06 1.49 1.49 1.57 3.80 3.80 3.82 2.30 8.34 49.82

0.29 0.62 0.62 0.54 0.02 0.21 0.21 0.21 0.49 0.49 0.49 0.54 0.40 0.56 0.56 0.59 0.18 0.18 0.19 0.14 0.57 8.10

0.05 0.10 0.10 0.09 0.003 0.93 0.93 0.93 0.08 0.08 0.08 0.09 0.06 0.09 0.09 0.10 0.96 0.96 0.97 0.66 2.49 9.84

0.57 1.20 1.20 1.05 0.04 0.27 0.27 0.27 0.91 0.90 0.90 1.01 0.74 1.04 1.04 1.10 0.27 0.27 0.27 0.17 0.66 14.15

−1.26 −2.66 −2.66 −2.33 −0.09 −3.35 −3.35 −3.32 −2.02 −1.99 −1.99 −2.23 −1.66 −2.33 −2.33 −2.47 −3.36 −3.36 −3.38 −2.15 −7.93 −56.22

−0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.09 −0.09 −0.09 −0.09 −0.09 −1.81

0.36 0.84 0.84 0.72 −0.06 1.77 1.77 1.78 0.63 0.64 0.64 0.71 0.52 0.77 0.77 0.81 1.76 1.76 1.78 1.03 4.04 23.88

Table 6. Detailed Optimization Results of the New Scheduling Model furnace 1

2

3

4

5

6

batch number

feed type

batch starting time (day)

batch ending time (day)

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

Fb Fb Fb Fb Fb Fb Fb Fa Fb Fb Fb Fb Fb Fb Fb Fb Fa Fa Fa Fd Fd Fd

0 16 63.33 107.71 0 2 56.33 108.71 0 27 68.33 105.71 0 21 62.33 106.71 0 53.33 104.71 0 24 109.71

15 62.33 106.71 160.71 1 55.33 107.71 179.71 26 67.33 104.71 154.71 20 61.33 105.71 158.71 52.33 103.71 159.71 23 108.71 210.71

total

ethylene revenue (M $)

propylene revenue (M $)

ethane revenue (M $)

other product revenue (M$)

feed cost (M$)

decoking cost (M$)

batch net profit (M$)

0.81 2.49 2.32 2.82 0.06 3.07 2.99 7.75 1.3 1.98 1.76 2.38 1.09 2.16 2.32 2.76 5.63 5.45 6.00 2.37 8.56 9.96 76.03

0.30 0.94 0.87 1.07 0.02 1.10 1.08 0.45 0.51 0.78 0.70 0.95 0.41 0.82 0.89 1.06 0.30 0.27 0.31 0.15 0.59 0.69 14.26

0.05 0.17 0.15 0.19 0.003 0.20 0.20 1.95 0.09 0.13 0.11 0.16 0.07 0.14 0.15 0.18 1.43 1.41 1.57 0.70 2.58 3.01 14.64

0.57 1.78 1.68 2.05 0.04 2.15 2.10 0.54 0.93 1.44 1.28 1.74 0.75 1.50 1.63 1.95 0.40 0.37 0.39 0.17 0.63 0.78 24.87

−1.28 −3.94 −3.71 −4.54 −0.09 −4.75 −4.69 −6.87 −2.08 −3.20 −2.85 −3.87 −1.69 −3.38 −3.66 −4.38 −4.94 −4.83 −5.33 −2.20 −8.09 −9.52 −89.89

−0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.08 −0.09 −0.09 −0.09 −0.09 −0.09 −0.09 −1.90

0.36 1.35 1.22 1.50 −0.06 1.68 1.59 3.73 0.67 1.05 0.92 1.28 0.55 1.16 1.25 1.49 2.73 2.58 2.85 1.10 4.18 4.83 38.01

multiple feed, multiple furnaces, and nonsimultaneous decoking, this Article extends the previous dynamic scheduling model for ethylene furnaces by considering more practical operating features such as nonhomogeneous last-batch ending and furnace load makeup, as well as the new environmental concern of background air-quality conscious decoking. With the newly developed model, the simultaneous furnace shutdowns can be successfully avoided at the end of the customized

meteorological and background air-quality should be further investigated to strengthen its applicability in future studies.

5. CONCLUDING REMARKS This Article has developed a new scheduling model for the optimal operation of the cracking furnace system in ethylene plants. In addition to the major operational aspects, such as L


Article

Industrial & Engineering Chemistry Research tR

Table 7. Impact of Scheduling Horizon Length on Model Statistics

tloi,j tup i,j

H (days) no. of batches no. of continuous variables no. of binary variables no. of constraints solution time (s)

120

170

240

3 7451 6114 27 729 31.64

3 7451 6114 27 729 40.98

5 13 911 10 790 54 931 288.07

ydi,j,l

Variables

Ej,k Gi

scheduling time horizon; meanwhile, certain manufacturing capacity losses and upsets can be alleviated during furnace decoking operations. Additionally, the obtained scheduling solution is background air-quality conscious, which optimally allocates the decoking emission peak within appropriate time windows that will cause less adverse environmental impacts on ground-level ozone and PM.

■

feed delivery time for the next schedule lower bound of batch processing time for feed i cracked in furnace j upper bound of batch processing time for feed i cracked in furnace j the yield parameter of product l for feedstock i cracked by furnace j during load increase period

Kj pj Sj,k SCi,j,k SPl,i,j,k Tclj,k,j′,k′

AUTHOR INFORMATION

Corresponding Author

*Phone: 409-880-7818. Fax: 409-880-2197. E-mail: Qiang.xu@ lamar.edu.

Ticsj,k,j′,k′

Notes

The authors declare no competing financial interest.

Ticfj,k,j′,k′

■

ACKNOWLEDGMENTS This work was partially supported by Texas Air Research Centre (TARC) headquartered at Lamar University, the Centre for Advances in Water and Air Quality at Lamar University, and Graduate Student Scholarship from Lamar University.

■

TlRj Tsmi,j,k ti,j,k

NOMENCLATURE

xj,k,j′,k′

Sets

K = {k|k = 1, ..., NB} I = {i|i = 1, ..., NF} J = {j|j = 1, ..., NC} L = {l|l = 1, ..., NP} Ij

set of operating batches set of feeds set of cracking furnaces set of considered products feedstock i that can be cracked in furnace j

xej,k,d

Parameters

yi,j,k

τi,j ai,j,l

decoking time used after feed i cracked in furnace j pre-exponential factor of product l’s yield formula for feed i cracked in furnace j Atw the lower boundary point of the emission time window bi,j,l power factor of product l’s yield formula for feed i cracked in furnace j Btw the upper boundary point of the emission time window ci,j,l constant coefficient of product l’s yield formula for feed i cracked in furnace j Cri raw material cost for feed i Csi,j unit decoking cost for feed i cracked in furnace j Floi lower bound on flow rate of feedstock i Fup lower bound on flow rate of feedstock i i Fup upper bound on flow rate of feedstock i cracked by i,j furnace j FLini,j flow rate of feed i cracked in furnace j H total time of current schedule UB upper bound of the scheduling time period Prl market price for product l Tlj leftover batch processing time for the schedule in furnace j

yej,k ylj,k yvj, k, j′, k′ zj,k,j′,1

■

ending time point of the k-th batch in furnace j the amount of feed flow rate above the lower bound Floi index of current processing batches for furnace j the percentage of load increase for furnace j starting time point of the k-th batch in furnace j total amount of feedstock i that is consumed during the k-th batch of furnace j total amount of product l that is produced during the k-th batch of furnace j cracking feedstock i the time duration of load increase during the k′-th batch operation of furnace j′ due to the k-th decoking operation of furnace j the starting time point of the load increase period within the k′-th batch of furnace j′ caused by the k-th decoking operation of furnace j the ending time point of the load increase period within the k′-th batch of furnace j′ caused by the k-th decoking operation of furnace j leftover batch processing time for the next schedule in furnace j The total time for load increase during the k-th batch operation of furnace j cracking feed i processing time for feed i cracked in the k-th batch of furnace j binary variable which is 1 if the k-th decoking in furnace j is no-overlap behind the k′-th decoking in furnace j′; otherwise, if the k-th decoking in furnace j is no-overlap ahead of the k′-th decoking in furnace j′, it is 0 binary variable which is 1 if the k-th decoking operation of furnace j falls within day d, otherwise, it is equal to 0 binary variable which is 1 if feed i is processed in the k-th batch of furnace j; otherwise, it is 0 binary variable which is 1 if the last batch operation k of furnace j exceeds the customized scheduling horizon (H); otherwise, it is equal to 0 binary variable which is 1 if the k-th batch operation of furnace j is the last active batch; otherwise, it is 0 binary variable which is 1 if the k-th decoking operation of furnace j falls within the k′-th batch operation of furnace j′, otherwise, it is equal to 0 binary variable which is 1 if the k-th batch ending time of furnace j is not less than the starting time of the first batch of furnace j′, otherwise, it is equal to 0

REFERENCES

(1) Zhao, C.; Liu, C.; Xu, Q. Dynamic scheduling for ethylene cracking furnace system. Ind. Eng. Chem. Res. 2011, 50 (21), 12026− 12040. (2) Chemical Economics Handbook Ethylene. https://www.ihs.com/ products/ethylene-chemical-economics-handbook.html (accessed 2014 October). M


Article

Industrial & Engineering Chemistry Research (3) Jain, V.; Grossmann, I. E. Cyclic scheduling of continuous parallel-process units with decaying performance. AIChE J. 1998, 44 (7), 1623−1636. (4) Schulz, E.; Diaz, S.; Bandoni, A. Interaction between process plant operation and cracking furnaces maintenance policy in an ethylene plant. Comput.-Aided Chem. Eng. 2000, 8, 487−492. (5) Lim, H.; Choi, J.; Realff, M.; Lee, J. H.; Park, S. Development of optimal decoking scheduling strategies for an industrial naphtha cracking furnace system. Ind. Eng. Chem. Res. 2006, 45 (16), 5738− 5747. (6) Lim, H.; Choi, J.; Realff, M.; Lee, J. H.; Park, S. Proactive scheduling strategy applied to decoking operations of an industrial naphtha cracking furnace system. Ind. Eng. Chem. Res. 2009, 48 (6), 3024−3032. (7) Liu, C.; Zhang, J.; Xu, Q.; Li, K. Cyclic scheduling for best profitability of industrial cracking furnace system. Comput. Chem. Eng. 2010, 34 (4), 544−554. (8) Zhao, C.; Liu, C.; Xu, Q. Cyclic scheduling for ethylene cracking furnace system with consideration of secondary ethane cracking. Ind. Eng. Chem. Res. 2010, 49 (12), 5765−5774. (9) Wang, Z.; Feng, Y.; Rong, G. Synchronized Scheduling Approach of Ethylene Plant Production and Naphtha Oil Inventory Management. Ind. Eng. Chem. Res. 2014, 53 (15), 6477−6499. (10) Little, A. D., Particulates Summary Report. In Environmental Considerations of Selected Energy-Conserving Manufacturing Process Options; EPA, Aug 1979; Vol. XVIII. (11) Stockwell, W. R.; Lawson, C. V.; Saunders, E.; Goliff, W. S. A review of tropospheric atmospheric chemistry and gas-phase chemical mechanisms for air quality modeling. Atmosphere 2012, 3 (1), 1−32. (12) EPA https://www.epa.gov/ozone-pollution. (13) Ho, T. C. Source Identification of Houston and Beaumont Aerosols using Positive Matrix Factorization. In Advances in Texas Air Quality, April 19 2011.

N


Emission Constrained Dynamic Scheduling for Ethylene Cracking

Recommend Documents