Article pubs.acs.org/IECR
Synchronized Scheduling Approach of Ethylene Plant Production and Naphtha Oil Inventory Management Zihao Wang, Yiping Feng, and Gang Rong* State Key Laboratory of Industrial Control Technology, Institute of Cyber-Systems and Control, Zhejiang University, Hangzhou 310027, China S Supporting Information *
ABSTRACT: Ethylene production process and the feedstock management are the critical activities in an olefin plant. Better analysis of these two problems gives rise to better use of the plant feeding material and enhanced production capabilities. It is important to ensure the feeding material loaded and unloaded continuously and consistently from upstream to downstream in case of any insufficient supply or increased production burden. Besides, considering the changing yield characteristic of the furnace cracking operation and the naphtha blending constraints, how to incorporate the real world operational characteristics into the integrated full space problem from upstream to downstream and improve its economic performance in the strategic level is quite significant. Facing this challenge, this paper develops a novel synchronized decision-supporting framework for the longterm planning and scheduling problem of an ethylene plant by simultaneously considering the upstream naphtha inventory management and the downstream ethylene furnace operation with a subcooperative model, which also incorporate more logistic and operational characteristics of the subproblems in practice. Moreover, the solution performance of both two subproblems is enhanced by the new solving approach. The resulting framework is shown to be robust and seamlessly integrated. The efficacy and the economic potential of this proposed integration framework are illustrated through a comprehensive case study from the real world ethylene plant.
1. INTRODUCTION Ethylene is an important chemical compound in everyone’s daily life. In recent years, the petrochemical plants have faced a rapid growth of ethylene demand due to the rapid development of the world economy and increasing population. To meet this increasing demand and fierce competition in the global market, significant efforts have been made in improving the operation performance of ethylene plants. Among all of these achievements, scheduling is one of the most important operational activities to manage the daily procurement, production and product shipment of a plant over a long or short-term time horizon. The goal of scheduling in process industrial plant is to achieve a long-term or short-term operation decisions and to increase the plant margin while minimize the operational cost at the same time. In chemical industries, the former researches of long-term or short scheduling have been focused on batch process. The comprehensive explanation and work on scheduling and planning can be found in the early works of Reklaitis and Rippin.1,2 They summarized and reviewed the development of batch process systems engineering with the area of planning and scheduling. Since scheduling formulations are often expressed as mixed integer programming models, Pinto and Grossmann3 presented an overview of assignment and sequencing models used in process scheduling with mathematical programming techniques. In general, scheduling formulations usually can be classified into two formulations: batch based and network-based. The networked-based approaches can be further classified into the continuous and discrete time formulation. The early work discretize the time horizon into a number of time intervals, which © 2014 American Chemical Society
may lead to the model inaccuracy, suboptimal solution, and unnecessary increase of the overall size of the mathematical model.4,5 The continuous-time representation can provide great potential for the development of more accurate and efficient modeling and solution approaches. Additionally, the continuoustime formulations are further subdivided into approaches that employ slot-based,6−8 global event-based,9−11 and unit-specific event-based formulations.12−15 Though many studies have been conducted for scheduling and operation optimization in chemical and petroleum plants, limited research can be found for ethylene production. The structure of an ethylene plant includes the upstream feedstock area, ethylene production unit, and downstream secondary production. Among all the processing operations in the ethylene unit, thermal cracking operation in the cracking furnace system is of great significance since the major product yields of an ethylene plant are mainly determined in this operation.16 Thermal cracking operation is a semicontinuous dynamic process, due to fact that the accumulating coke on the reaction coil will affect the heat transfer performance and result in react pressure drop. 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.17 As the limited resource of decoking facility and the safety concern of the disturbance from decoking operation, the simultaneous decoking process is not allowed in the system. Jain and Grossman first Received: Revised: Accepted: Published: 6477
January 7, 2014 March 9, 2014 March 18, 2014 March 18, 2014 dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
the final result is still infeasible due to the shortage of qualified raw material during the a subcycle time or in a month period. According to the problem analysis, we proposed a novel synchronized scheduling framework to solve this integrated scheduling problem of upstream inventory management and the related downstream furnace cracking operation. In this paper, we only considered naphtha oil feeding flow in the tank area since it takes up more than 80% of the total feedstock and comes with different grades of various qualities from the upper level supplier. The naphtha feeding process for the ethylene plants are composed of ship unloading, storage tanks assignment, naphtha transfer between tanks and transferring the charging tanks to the furnaces. All these decisions will affect naphtha composition on charging tanks, therefore the product composition.25 In the proposed synchronized integration framework, the furnace cracking schedule can be modified and accommodate with the detailed feedstock inventory and supply information. On the other hand, the upstream monthly inventory management scheduling problem is also solved month after month in a rolling horizon method with the naphtha demand information of the cracking furnace cyclic scheduling result. The content of this paper is organized as follows. A brief introduction about the problem we study and the basic structure of the integration model is presented in section 2. In section 3, we proposed a new MINLP model with real world operational characteristics for cyclic scheduling in furnace under the operation style with all furnaces running online to minimize the product flow fluctuation. Given the complexity of this MINLP model and increased problem size, a heuristic binary searching algorithm with integer cut is also introduced. The improved upstream scheduling model of naphtha inventory management is presented in section 4 with a continuous time formulation under the specific event demand orders from downstream scheduling result of furnace cracking operation. In section 5, In order to keep the successive and consistent naphtha oil feeding flow in this integration model, a subcooperative model is formulated in the rolling iterative step. The proposed framework is tested through a case study from the real world ethylene plant, and the synchronized approach is also compared with the classical managerial heuristic method in section 6 to prove its efficiency. At last, section 7 draws conclusions of the paper and indicates prospective for future research.
proposed an MINLP model describe the decaying performance of the cracking furnace for a cycling scheduling.18 Schulz et al. considered the recycled ethane and formulated mathematical model of an entire plant at each time interval to carry out production planning for meeting varying demands. Besides, an MINLP model is also proposed by them based on the their former research to increase the net profit.19,20 In the work of Lim et al., they solved the decoking scheduling problem of multiple tube system with the constraint of a constant total throughput and also proposed three alternative solution strategies to solve the MINLP problem.17 However, the situation of multiple feeds is still not considered in Lim’s model. Facing this problem, Zhao et al.16 developed a new MINLP model with multiple feeds and considering secondary ethane cracking. In their model, nonsimultaneous cleanup constraints were incorporated, and an industrial heuristic method was compared with their model to demonstrate the better economic performance. After that, Zhao et al. also extended the model in a rescheduling frame to deal with dynamic information.21 Nevertheless, in their model, the decoking operation could result in a serious impact on the furnace product flow which will lead to a sharp decline while the total operating capacity of the furnace system has to be kept constant to stabilize the downstream processing units, such as coolers and separators. Nowadays, the fierce and volatile raw-material markets have forced ethylene plants to use various feedstocks. In the cracking furnaces system, different grades of naphtha (NP) and other kinds of raw material like liquefied petroleum gas (LPG), light diesel (LD), propylene, or ethane are consumed from upstream tank area or pipelines of refinery. The purchase of the feedstock and its delivery are critical issues in the scheduling problem of ethylene plant since the product composition of pyrolysis depends on the feedstock composition and the operating conditions in cracking furnaces. Moreover, the related reports also indicate that the scheduling operation cost comprises more than 50% of the total production cost. Based on the above analysis, research on the scheduling of the ethylene furnace cracking operation and considering the economic trade off with its upstream inventory management is of highly significant importance and complex depending on the raw material provision, production demands, inventory levels and the operational cost. However, the study of this integration problem is still lacking in the literature. Li and Ierapetritou developed an augmented lagrangian method to solve a large scale integrated problem.22 To reduce the problem size and improve the solution quality of rolling horizon method, Li and Ierapetritou considered the production capacity information into the planning model based on parametric programming technique.23 Sylvain et al. also studied the integrated problem of refinery planning and crude oil scheduling with a new Lagrangian decomposition approach to solve this problem separately.24 Even though some works has noticed the importance of integration planning and scheduling of ethylene furnace operation and the feedstock inventory management, how to manage the upstream logistic and feedstock delivery to cooperate with the furnace production is still unsolved. In practice, the long-term scheduling problem of feedstock inventory management is divided into month period while the furnace cracking operation is a cyclic scheduling problem. Thus, the time scale turns out to be very different which will lead to information inconsistency for the raw material inventory management and furnace scheduling. It is highly likely that without the detailed feedstock information, even though we optimize the furnace cracking during a long-term time horizon,
2. PROBLEM STATEMENT As can be seen in Figure 1, the ethylene production system consists of two parts, the cyclic scheduling problem in the
Figure 1. Flowchart of the integrated ethylene production topology. 6478
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
cracking furnace operation and the naphtha oil feedstock inventory management. Except for the naphtha oil, other feeds such as heavy naphtha, light diesel, liquefied petroleum gas, and ethane (recycled ethane) can be processed in separate furnaces or cocracking with naphtha. Since these supplementary cracking feeds are mainly sent by the pipelines from other petroleum plants or other separate storage tanks, we will not incorporate these delivery logistic operations into the integration model. The goal of the integrated scheduling problem is to gain the total ethylene production profit as much as possible while minimizing the operation cost of the two subproblems simultaneously under the monthly feedstock information and the outside market demands. In the real world ethylene plant, these two problems are solved separately with estimated information from historical reports. However, the asymmetric information could generate infeasible scheduling orders of both sides that may result in safety problems of the ethylene production like increasing the inventory level fluctuation and over running burdens of furnaces. To overcome this, a synchronized scheduling approach is proposed in this integration framework. The basic flow logic of the integrated synchronized scheduling problem is illustrated in Figure 2.
(3) increased running capacity time, (4) increased level for each varying running capacity. Then, the cyclic scheduling result is changed into event-based orders in each month and is used as the demand information for the upstream scheduling problem of each month. In the second stage, based on the demand event orders, we are able to formulate an event-based continuous time formulation model for the naphtha inventory management problem. The objective of naphtha scheduling problem is to minimize the operational cost while satisfy the demand orders at the same time. The given information is as follows: (1) naphtha raw material supply schedule, (2) constraints for component concentration in each tank, (3) quality constraints for cracking naphtha after blending, (4) capacities for tanks, (5) demand orders, (6) operational cost. The variables to be deterimined are as follows: (1) shipment schedule in the port, (2) sequence and starting/ending time for each type of naphtha oil in the tank area, (3) blending transfer in the charging tank, (4) amount of total flow volume for each material transfer, (5) the new demand order after possible adjustment. Since it is desirable to ensure the feasibility of the full-space scheduling problem and keep the consistency of the material flow from upstream to downstream, a subcooperative model is proposed in the third stage. The objective of this model is to adjust the actual furnace cracking operation with the real status of feedstock supply capability and the shipment arriving schedule by changing the time duration of cracking operation or substituting with other feeds. The given information is as follows: (1) original schedule in this month including: raw material assignment in different furnace, starting/ending time of cracking operation in current month, (2) capabilities for additional supply, (3) maximum/minimum time duration for cracking substitute material, (4) feeding flow rate of added substitute cracking material. The variables to be determined are as follows: (1) the sequence of substitute material in each batch of different furnace, (2) time duration for every added serial cracking, (3) new starting or ending time for cracking operation in the current month, (4) average flow rate and time duration of each demand order. Consequently, after solving the subcooperative model if the infeasible situation happens in the current month period, the reoptimized result will be returned to the upstream scheduling model and the three stages problem will be solved iteratively until we get the feasible solution within certain variation tolerance gap. At last, the cyclic scheduling model of furnace cracking operation will also be reoptimized again with updated conditions month by month to form a rolling horizon framework.
Figure 2. Framework of synchronized scheduling approach.
As shown in Figure 2, the scheduling approach can be categorized into three stages. At the first stage, we build a MINLP cyclic scheduling model of furnace groups and optimize the best average profit per day under the approximated supply information. The objective of cyclic scheduling in furnace cracking operation is to maximize the average profit per day. The given information is as follows: (1) average supply feed flow rate, (2) operational cost for each furnace (cleanup cost), (3) time duration for cleanup of different furnace with different feeds, (4) purchase cost and demand price, (5) cracking yield information of different products with different feeds, (6) feeding flow rate of different raw material into different furnaces. The variables to be determined are as follows: (1) feed assignment for each furnace group, (2) starting time and ending time of each batch,
3. LONG-TERM CYCLIC SCHEDULING PROBLEM OF FURNACE CRACKING OPERATION As the first stage of the synchronized scheduling framework, we will introduce the cyclic scheduling problem in this section and 6479
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
Figure 3. Basic operation structure of furnace cracking system.
propose a new MINLP model considering more operational characteristic and detailed structures in the cracking furnace system. The furnace operation system has multiple parallel furnaces to process the multiple feeds from the upstream supplier. These feeds are cracked in the pipe groups to produce the final products for sale or to produce the secondary processing material. The proposed model is a cyclic scheduling problem due to the periodical decoking operation, which keeps the conversion rate above the minimum level of the final products and protects the furnace coil pipes. Normally, there are several furnaces for one ethylene unit, among these furnaces there is one furnace specially designed for cracking ethane. The cracking furnace system usually runs in two operation modes: furnace system with one furnace used as backup furnace (backup furnace running mode) and furnace system with all furnaces online (no backup furnace running mode). In the first mode, all of the furnaces will be running at their maximum designed capacity, which will lead to a sharp decline of the ethylene production yield and increase the coking velocity. If one of the furnaces is shutdown, the backup furnace will be setup for cracking operation to make up for the total product amount decrease of the furnace system. In practice, the no-backup furnace style is widely used in the real plant operations nowadays. In this mode, all of the furnaces are simultaneously running online at 70−85% of the maximum designed capacity. When one furnace is ready for cleanup, the rest of the furnace groups will increase its running capacity temporally to make up for the loss in the ethylene production. The basic structure of the furnace operation system is shown in Figure 3. The minimum producing unit of a cracking furnace is the pipe group. Generally, there are 4−6 pipe groups in the SRT type furnace. The different pipe groups could crack different feeds with similar pyrolysis characteristics of the same furnace, but they should be set up for decoking together at the same time. This kind of operation mode is cocracking or concurrent cracking. On the other hand, if there is a feed flow fluctuation or some uncertain upstream supply condition happens, the furnace will crack different feedstock during successively in the same batch. This operation style is serial cracking, which is not very common in the real world operation. In this full-space problem, the serial cracking is not considered in the first stage cyclic scheduling problem since it is a prescheduling process. When adjusting the original plan in the subcooperative model, serial cracking is introduced as it is always used to cooperate with
emergency and situation for small adjustment in the furnace groups. However, it remains a much more difficult scheduling problem of how to arrange the furnace feeding material in each pipe group and set the best starting and ending time to avoid simultaneous cracking and change the running capacity when one furnace is under decoking. In this model, we improved the furnace cyclic scheduling model from Zhao et al.16 and introduce the aforementioned real world cracking operation constraints. The concepts of the model are illustrated by the Figure 4, which shows the concept of batch
Figure 4. Concepts of cyclic scheduling.
processing time, total cycle time, nonsimultaneous clean up, cocracking. and serial cracking operation, varying running capacity under the no backup furnace operation style. Before we introduce the model, some assumptions should be emphasized: (a) There is enough feedstock supply under specific average flow rate range. (b) The generated yield model used for each pipe groups in the furnace are accurate and has the same characteristics in the same furnace. (c) The recycled ethane is fully reused for the secondary cracking. (d) The running capacity is only variable when other furnace is under decoking. (e) The increased running capacity during the decoking operation will not affect the original decaying trend of the furnace production yield. 3.1. Mathematical Model of Furnace Cracking Operation. The input and output flow of the material balance in the cracking furnace is presented in eqs 1.1 and 1.2. In these 6480
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
Equation 1.9 means the average feeding flow rate of material l should be within the estimated supply capacity.
equations, the operation mode i, i ∈ Isj represents that all of the pipe groups in the same furnace is cracking the same material, which means only one type of feeds can be cracked in furnace during one batch. Thus, mode i stands for material l in this operation mode Isj , and the feeding material l can be connected with mode i in the set l ∈ LiF, which is predetermined set in the mass balance equation. The term FLinl′,i,j represents the feeding flow rate of material l′ of the operation mode i in furnace j. As pointed out in the work of Jain and Grossmann,18 the term cl,l′,j + al,l′,j ebl,l′,jt describes the changing yield of of the product l when cracking the material l′ of furnace j. Hence, the total amount of product l produced by furnace j when running in mode i during batch n is represented by the integral function 1.1. The total amount of feeding material is in calculated by multiplying the flow rate FLl′,i,j and the time duration ti,j,n.
∫0
SPl , i , j , n = j ∈ J,
ti , j , n
FLlin, i , j(cl , l , j ′
n ∈ N,
SCl , i , j , n = FLlin, i , jti , j , n ,
′
+
al , l , j ebl ,l′,jt )dt ,
∀i∈
′
l′ ∈ L Fi ,
l ∈ L Pi
∀ i ∈ I js ,
j ∈ J,
∑ SPl ,i ,j ,n ≤ H × CAPupj × wvi ,j ,n,
∀ i ∈ Ij ,
j ∈ J,
l ∈ Lp
(1.5)
n∈N
∑ SCl ,i ,j ,n ≤ H × CAPupj × wvi ,j ,n,
∀ i ∈ Ij ,
l ∈ L Fi
j ∈ J, (1.6)
n∈N
∑ ∑ SPl ,i ,j ,n = ∑ ∑ SCl ,i ,j ,n ,
∀ j ∈ J,
n∈N
l ∈ L Fi i ∈ Ij
l ∈ L p i ∈ Ij
(1.7)
∀ i ∈ I js ,
∑ ∑ ∑ SPl ,i′,j′,n ≤ FLlin′,i ,j ∑ ti ,j ,n ,
I js ,
i ′∈ I j ′∈ J n ∈ N
j ∈ J,
n∈N
l ∈ L Pi ,
ET
j=J ,
(1.1)
l′ ∈ L Fi
l = LET ,
(1.8)
H × FLa ll ≤
n ∈ N,
∑ ∑ ∑ SCl ,i ,j ,n ≤ H × FLa lu,
l ∈ L Fi
i ∈ Ij j ∈ J n ∈ N
l ∈ L Fi
(1.2)
(1.9)
Co-cracking operation means two or more feeds are concurrently cracked in one furnace. In this model, different feeds are cracked separately in different pipe groups such as naphtha with light diesel, naphtha with ethane or propylene, and naphtha cocracked with liquefied petroleum gas. Equations 1.3 and 1.4 are calculated similarly under the cocracking mode i, i ∈ Icj , and this mode i will represents a set of material l, l ∈ LiF can be cocracked together. We assume that the cocracking mode m represents the proportion of the mixed feeding materials decided by the number of the total the pipe groups in SRT furnace. For example, if there are 4 pipe groups in the each furnace and the cocracking materials are A and B, it is reasonable to assume that we have three modes: the feed A comes into one, two, or three pipe groups while the feed B goes to the left. In addition, the actual destination-pipes of the feed flow is not decided beforehand, since it will not have much influence to the final scheduling plan. Thus, the parameter FLGinl,i,j.m gives the average flow rate of feeding material l under cocracking mode i with combination mode m. Then, the actual feeds flow into the furnace will be decided also along with the set l ∈ LiF
The total number of the batches is a heuristic integer during one cycle for different furnaces. Equations 1.10 through 1.12 give the logical constraints of the batch utilization. Equation 1.10 explains that one batch slot could only be assigned to one cracking mode at most, and eq 1.11 suggests that the first batch should always be used for cracking. According to the industrial experience, every purchased material l will be used according to supply capability arrangement in the long term enterprise planning. Given the stability and higher efficiency of the operation mode when the pipe groups of the furnace are cracking the same material, eq 1.12 is used to ensure this operation mode i, i ∈ Isj will be used at least once during the total cycle time to ensure the usage of material supply.
SPl , i , j , n =
∑ ∑∫ 0
l ′∈ L Fi m ∈ M
∀ i ∈ I jc ,
SCl , i , j , n =
∑
l ∈ L Pi ,
tgi , j , m , n
FLGlin, i , j , m(cl , l , j ′
′
+ al , l
j∈J
FLGlin, i , j , mtgi , j , m , n ,
′,j
e
bl , l , jt
′
∑ wvi ,j ,n ≤ 1,
n>1
∀ j ∈ J,
n ∈ N,
n=1
i ∈ Ij
(1.11)
∑ ∑ wvi ,j ,n ≥ 1,
)dt ,
∀ i ∈ I js (1.12)
j∈J n∈N
mode i, i ∈ Icj , the binary
If one furnace is under the cocracking variable wvgi,j,m,n is used to present that only one mixed combination type m is chosen to be aligned with the furnace operation i, this constraint is explained in eq 1.13. Equation 1.14 and eq 1.15 give the time bounds for the cocracking operation mode. Equation 1.16 is an additional constraint indicates that when one batch slot is not used, the following batches of that furnace will not be utilized either.
l ∈ L Fi ,
m∈M
j∈J
n ∈ N,
(1.10)
∑ wvi ,j ,n = 1,
(1.3)
∀ i ∈ I jc ,
∀ j ∈ J,
i ∈ Ij
(1.4)
The binary variable wvi,j,n is used to assign different operation mode i to the furnaces j during time slot n. This constraint can be modeled as shown in eqs 1.5 and 1.6. Equation 1.7 is the mass balance of furnace operation. Besides, in order to prevent the recycled ethane accumulating within the system, the total amount of ethane produced by the furnace groups during the total cycle time should be less than the processing capacity of the ethane furnace. This relationship is expressed by eq 1.8.
∑
i ∈ I jc ,
wvgi , j , m , n = wvi , j , n,
j ∈ J,
n∈N
m∈M
(1.13)
THli , jwvgi , j , m , n j ∈ J, 6481
≤ tgi , j , m , n ≤
m ∈ M,
THiu, jwvgi , j , m , n,
n∈N
∀i∈
I jc , (1.14)
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research i ∈ I jc ,
∑ tgi ,j ,m,n = ti ,j ,n ,
j ∈ J,
Article
(1.15)
m
∑ wvi ,j ,n ≥ ∑ wvi ,j ,n′, i ∈ Ij
information such as how many cleanups of other furnaces will occur during the batch processing time n of furnace j is necessarily to be collected. Here, a binary variable yvj,n,j′,n′ is introduced to describe this constraint. If yvj,n,j′,n′ is designated as 1, which means the cleanup time window of n′th cleanup of furnace j′ is within the processing time duration window of furnace j in batch n. Otherwise, if yvj,n,j′,n′ = 0, then the referred condition is not satisfied. The value of the binary variable yvj,n,j′,n′ can be deduced from xj,n,j′,n′, for example, if the nth cleanup of furnace j is within the processing time duration window of furnace j′ in batch slot n′, j < j′ then the n′−1th cleanup of furnace j′ is ahead of the nth cleanup of furnace j, meanwhile n′+1th cleanup of furnace j′ is behind the nth cleanup of furnace j. This concept can be written in mathematical forms in eqs 1.24−1.31.
n∈N
n ∈ N,
n < n′ < N ,
i ∈ Ij
∀j∈J
(1.16)
The nonsimultaneous cleanup constraints can be referred to the work of Zhao et al.16 A binary variable xj,n,j′,n′ is used to control the cleanup sequence of different batch slot n, n′ of two different furnaces j, j′, if it has the value of 1, which means the nth cleanup in furnace j is no-overlap behind the n′th cleanup in furnace j′, if it is designated as 0, the nth cleanup in furnace j is no-overlap ahead the n′th cleanup in furnace j′. In the cyclic scheduling model of the furnace system, the cracking process is a round table problem and the first batch may start from the previous cycle or from the current cycle. A binary variable θj is introduced to determine the difference. This constraint is realized by big-M method in eq 1.17. In eqs 1.18 and 1.19, we are able to formulate the constraints which explain the relationship of the cleanup time βi,j between the two adjacent batches. The time duration ti,j,n for each batch is calculated in eqs 1.20 and 1.21. −θjM ≤ Tfi , j ,1 − Tsi , j ,1 ≤ (1 − θj)M ,
∀ i ∈ Ij ,
yvj
n , n′ ∈ N ,
∑ βi ,j wvi ,j ,n′ − (1 − θj)H ,
Tsi , j , n = Tfi , j , n − 1 +
j∈J
n′ = N
j , j ′ ∈ J , j < j′ ,
,
n , n′ ∈ N (1.25)
j , j ′ ∈ J , j < j′ , (1.26)
j , j′ ∈ J ,
j < j′ , (1.27)
∑ βi ,j wvi ,j ,n− 1,
∀ j ∈ J,
j , j′ ∈ J ,
≥ xj , n , j , n − 1 − xj , n , j , n , ′,n′,j,n ′ ′ ′ ′ n , n′ ∈ N
j , j′ ∈ J ,
∀ j ∈ J,
n ∈ N,
(1.21)
zvj , n , j
′,n′
≤ yvj , n , j
zvj , n , j
′,n′
≥ yvj , n , j
(1.22)
∀ i ∈ Ij ,
(1.30)
However, in some furnaces the last used batch slot might not be the maximum batch number N. Thus, the binary variable yvj,n,j′,n′ deduced from the above equations cannot give the right information we need. The following constraints use the information of the batch slot utilization variable wvi,j,n to get the true status of the aforementioned relationship with a positive continuous variable zvj,n,j′,n′ which is inherently 0−1 integer and has the same meaning with variable yvj,n,j′,n′ . This relationship is written in eqs 1.32−1.35.
∀j∈J
TLHi , jwvi , j , n ≤ ti , j , n ≤ TUHi , jwvi , j , n,
j < j′ ,
= 1 − xj , n , j , n + xj , n , j , n , j , j′ ∈ J , j < j′, ′,n′,j,n ′ ′ ′ ″ n , n′, n″ ∈ N , n′ = 1, n″ = N (1.31)
Timing constraints restrict the time variables of all the furnaces. The duration of an exact batch processing time usually takes 20−90 days depending on the furnace type, feed characteristic and cracking severity. Thus, eq 1.22 gives the constraint that the total cycle time of each furnace is the summation of the batch duration time and the cleanup time. Equation 1.23 provides the range of processing time for each batch. i ∈ Ij n ∈ N
(1.29)
yvj
n>1
i ∈ Ij
∑ ∑ (ti ,j ,n + βi ,j wvi ,j ,n) = H ,
j < j′ ,
= 1 − xj , n , j , n + xj , n , j , n , j , j′ ∈ J , ′,n′ ″ ′ ′ ′ ′ n , n′, n″ ∈ N , n = 1, n″ = N
(1.20)
∀ j ∈ J,
(1.28)
yvj , n , j
i ∈ Ij
n=1
j < j′ ,
n ∈ N, yvj
∑ ti ,j ,n − H × θj ,
∑ ti ,j ,n ,
≥ xj , n , j , n − xj , n − 1, j , n , ′,n′ ′ ′ ′ ′ n , n′ ∈ N
yvj , n , j
(1.18)
n>1
n∈N
′,n′−1
≤ 1 − xj , n , j , n , ′,n′,j,n ′ ′ n , n′ ∈ N
∀ j ∈ J,
n , n′ ∈ N (1.24)
yvj
(1.19)
Tfi , j , n = Tsi , j , n +
≤ xj , n , j
j , j′ ∈ J , j < j′ ,
,
≤ 1 − xj , n − 1, j , n , ′,n′ ′ ′ n , n′ ∈ N
i ∈ Ij
n ∈ N,
′,n′,j,n
′,n′
i ∈ Ij
n = 1,
Tfi , j , n = Tsi , j , n +
≤ xj , n , j
yvj , n , j
(1.17)
Tsi , j , n = Tfi , j , n + ′
′,n′
yvj , n , j
j ∈ J, (1.23)
∀ j , j′ ∈ J ,
, ′,n′
′,n′
+
n , n′ ∈ N
∑ wvi ,j ′ ,n′ − 1,
∀ j , j′ ∈ J ,
i ∈ Ii
n , n′ ∈ N
As we have mentioned before, when one furnace shuts down for decoking in the no backup furnace operation mode, the current running capacity of the rest furnaces would rise by 10− 15% to make up for the loss of the product amount. Thus, the
zvj , n , j
′,n′
≤
∑ wvi ,j ′ ,n′,
(1.32)
(1.33)
∀ j , j′ ∈ J ,
n , n′ ∈ N
i ∈ Ij ′
(1.34) 6482
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research zvj , n , j
′,n′
∑ wvi ,j ,n ,
≤
∀ j , j′ ∈ J ,
Article
Tadi,j,n or Tadmi,j,m,n only takes a small proportion of the total processing time in the current batch n. The rectified equation of the consuming and production function is represented in eqs 1.45−1.48.
n , n′ ∈ N
i ∈ Ij
(1.35)
With the variable zvj,n,j′,n′, the total time of increased running capacity during the current batch is given in eqs 1.36−1.44. Equations 1.36 through 1.38 give the exact time of increased running capacity in furnace j of each batch n when the batch n′ of furnace j′ is under decoking operation. Given the difference between cocracking operation mode and processing mode cracking single material, the information of total time length for the increased running capacity in each furnace and each batch is presented separately by Tadi,j,n and Tadmi,j,m,n, which is calculated by eqs 1.39−1.44. Tcl j , n , j
′,n′
≤
∑ βi ,j wvi ,j ′ ,n′,
∀ j , j′ ∈ J ,
′
i ∈ Ij ′
SCl , i , j , n = FLlin, i , jti , j , n + pj FLinl , i , jTad i , j , n, j ∈ J,
Tcl j , n , j
,n
′ ′
≤ Thczvj , n , j
∀ j , j′ ∈ J ,
, ′ ′ ,n
n , n′ ∈ N
SCl , i , j , n =
SPl , i , j , n =
j′
Tad i , j , n ≤
∑ ∑
Tcl j , n , j
j ′∈ J n ′∈ N
∀ i ∈ I js , Tad i , j , n ≥
j ∈ J,
∑ ∑
∀ i ∈ I js ,
j ∈ J,
l∈
− H(1 − wvi , j , n), (1.40)
n∈N
(1.41)
∀i∈
∑ Tcl j ,n,j ′ ,n′ − H(1 − wvgi ,j ,m,n),
Tadmi , j , m , n ≤
j ∈ J,
m ∈ M,
n∈N
∑ ∑∫ 0
L Pi ,
m∈M
j ∈ J,
Tadmi , j , m , n ≤ tgi , j , m , n,
m ∈ M,
∀ i ∈ I jc ,
n∈N
j ∈ J,
(1.47)
FLGlin, i , j , m(cl , l , j + al , l , j ebl ,l′ ,jt )dt ′ ′ ′ ∀ i ∈ I jc , (1.48)
∑ ∑ ∑ ∑ axlSPl ,i ,j ,n l ∈ L p i ∈ Ij j ∈ J n ∈ N
−
∑ ∑ ∑ ∑ bxlSCl ,i ,j ,n
l ∈ L Fi i ∈ Ij j ∈ J n ∈ N
(1.42)
−
∑ Tcl j ,n,j ′ ,n′ + H(1 − wvgi ,j ,m,n), j ∈ J,
n∈N
n∈N
∑ ∑ ∑ cxi ,jwvi ,j , n i ∈ Ij j ∈ J n ∈ N
j′,n′
∀ i ∈ I jc ,
tgi , j , m , n
j ∈ J,
pj FLGlin, i , j , mystl , l , jTadmi , j , m , n ′ ′
max Profit × H =
j′,n′
I jc ,
(1.46)
Objective Function. As shown in eq 1.49, the objective of scheduling problem in furnace cracking system is to maximize the average daily profit over a long cycle period H. The first term in the equation is the total profit earned from the cracking products, the second term is the cost of raw materials, and the third term represents the operational cost. The last two terms are the extra operational cost. It is used to minimize the running burden of different furnaces when they increased their running capacity.
(1.39)
j ∈ J,
j ∈ J,
FLGlin, i , j , mtgi , j , m , n + pj FLlin, i , j , mTadmi , j , m , n ,
l ′∈ L Fi m ∈ M
l ′∈ L Fi
+ H(1 − wvi , j , n),
n∈N
∀ i ∈ I js ,
Tad i , j , n ≤ ti , j , n ,
Tadmi , j , m , n ≥
′,n′
FLlin, i , j(cl , l , j + al , l , j ebl ,l′,jt )dt ′ ′ ′
l ∈ L Fi ,
∑ ∑
+
(1.38)
n∈N
Tcl j , n , j
j ′∈ J n ′∈ N
′,n′
∑
∀ i ∈ I jc ,
(1.37)
n , n′ ∈ N
ti , j , n
m∈M
n , n′ ∈ N
≥ βi , j ( ∑ wvi , j , n + yvj , n , j , n − 1), ′,n′ ′ ′ ′ ′ ′ i∈I
∀ j , j′ ∈ J ,
∫0
(1.45)
+ pj FLlin, i , jystl , l , jTad i , j , n, ∀ i ∈ I js , ′ ′ n ∈ N , l′ ∈ L Fi , l ∈ L Pi
(1.36)
Tcl j , n , j
l ∈ L Fi
n ∈ N,
SPl , i , j , n =
∀ i ∈ I js ,
− (1.43)
∑ ∑ ∑ di ,jTadi ,j , n i ∈ I js j ∈ J n ∈ N
−
m ∈ M,
∑ ∑ ∑ di ,j ∑ i ∈ I jc j ∈ J n ∈ N
Tadmi , j , m , n
m∈M
(1.44)
(1.49)
The mass balance equation for material consuming and production should be modified either. In eq 1.45, the term in pjFLl,i,j Tadi,j,n represents the increased part of the feeding material when furnace j is running under the single-material cracking mode i of batch slot n. In eq 1.46, the term pjFLinl′,i,jystl,l′,jTadi,j,n represents the increased amount of product l when cracking feed l′. The variable pj is a heuristic parameter set for the increased percentage of running capacity, we discretize the variation level into set points such as 10% or 15%, and then this parameter will be iteratively chosen for the optimization result. Considering the complexity caused by the production function with the nonlinear integral term, it is reasonable to assume that the yield of the production material is a constant value ystl,l′,j as the time duration
3.2. Solution Method. The nonbackup cyclic scheduling formulation of ethylene furnace cracking operation is a challenging MINLP model. The nonlinear yield constraints and some multidimensional integer variables make this MINLP problem even harder to solve, let alone getting the optimal solution. To overcome this challenge, we proposed an iterative heuristic solution strategy with integer cut to improve the solution performance of DICOPT solver in the GAMS. The scheduling model is divided into two stage subproblems. The first stage is a cyclic scheduling problem without varying running capacity which is solved by the DICOPT solver in GAMS with CONOPT to handle the NLP subproblem and CPLEX to handle the MIP subproblem. The integer variables such as wvi,j,n in the
n∈N
6483
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
first stage problem is used as the starting point of the second stage problem. Moreover, the heuristic searching rules is employed as the additional constraints to cut significant branches in the searching path to save the calculation time. For example, (i) the maximum batch slot number |FN|j will be added in the constraints and affects the batch duration time for each furnace processing operation. Thus, the solver would be able to jump out of the local optimal solution and choose the best subcycle time. (ii) The cleanup operation with longer time duration should be better within the time window when the system that has a higher ethylene production yield. After each heuristic searching process, if the result of second stage is not better than the previous solution of the scheduling model, integer cut is set to exclude the former tested solutions. This iterative solution repeats until the possible solution is not getting better in certain iterative steps. The mathematical expression is in eq 1.50, in which PI = {(i,j,n)|wv′i.j.n = 1}, NI = {(i,j,n)|wv′i.j.n = 0}, PG = {(i,j,m,n)|wvg′i,j,m,n = 1}, NG = {(i,j,m,n)|wvg′i,j,m,n = 0}, and wv′i.j.n and wvg′i,j,m,n are integer solutions from the previous iterative step.
∑ −
∑
wvi , j , n −
(i , j , n) ∈ PI
(i , j , n) ∈ NI
∑
wvi , j , n +
∑
4. NAPHTHA OIL INVENTORY MANAGEMENT MODEL The operation process for naphtha oil inventory management consists of unloading naphtha oil in berth, transferring naphtha oil from vessels to storage tanks, blending raw material in naphtha charging tanks, and sending the qualified naphtha for pyrolysis in the cracking furnaces. In the tank area site, the raw material usually comes with different amount and various qualities in prescheduled date. Since it is imperative to control the furnaces’ pyrolysis yield and the coke formation that can be affected by the sulfur content and specific gravity in the raw material, the feeds should be blended first in the charging tanks to meet the specific component concentration for the cracking operation. In our proposed scheduling model, two key component concentrations of C4 hydrocarbon fraction and C6−C8 hydrocarbon fraction (paraffin, olefins, and aromatics) are considered in the blending process. The computational complexity comes from the operation constraints in the process of naphtha delivery, storage, blending, and transfer. A large amount of the literature has been reported in the area of crude oil scheduling in refinery. The early work of Lee et al. first address the crude oil inventory management problem using a discrete-time formulation and linearize the bilinear term in the MILP model.26 In the work of Reddy et al.,27 they propose a hybrid approach that is simultaneously discrete and continuous and use an iterative procedure based on the resolution of MILP. Nevertheless, the proposed method in their research can only get a suboptimal solution, while the former work of Li et al.28 avoid the additional resolution of a NLP model. The recent work based on continuous-time formulation is reported by Furman et al.;29 they proposed a novel method for representing the flow in and out of a tank with generalized event-based model and reduced the number of necessary time events and decrease the size of the model. In this paper, we use the event-based unit-specific continuous time formulation for the naphtha oil inventory management problem with the monthly demand orders from the first stage scheduling results. The novel idea in this model is that we use the specific demand orders rather than the normal accumulated total amount orders for the scheduling horizon, and we also use an iterative procedure to deal with the composition discrepancy in bilinear terms while the work of Furman et al.29 only use the linear approximation method. Besides, as we have discussed in section 2, infeasible solution sometimes could happen due to the incomplete supply information in solving the cyclic scheduling problem. Thus, some slack variables should be introduced into this model to provide the feasible schedules with certain variation level for the subcooperative model. Considering the length of the research paper, the mathematical model below is a brief introduction and simplification for the naphtha inventory scheduling problem; more detailed description of the mathematical formulation and the related nomenclature can be referred to the Supporting Information. 4.1. Mathematical Model of Naphtha Inventory Management. Mass Balance Equations. In the scheduling model of naphtha inventory management, the mass balance equation includes the input and output relations in vessels, storage tanks, and charging tanks, which can be explained in eq 2.1.
wvgi , j , m , n
(i , j , m , n) ∈ PG
wvgi , j , m , n ≤ |PI | + |NI | + |PG| + |NG|
(i , j , m , n) ∈ NG
(1.50)
3.3. Calculate the Naphtha Feeding Event Series. From the above cyclic scheduling of furnace cracking operation, the scheduling results includes the best cycle horizon time H, feed flow assignment of every pipe group in each furnace, the starting and ending time of each batch period and the sequence of decoking operation. As a matter of fact, the batch window of the cyclic scheduling result for each furnace is always not in accordance with the scheduling period we executed the upstream inventory management or downstream production. In the manufacturing execution system (MES), the short-term planning model is usually divided into monthly periods. The upper level’s schedule like the provision of naphtha oil is usually given by the supply chain manager in each month for 6−12 months’ plan. In this respect, the optimized result has to be transformed into the monthly operation orders in event series. The illustration of the transformation process is presented in Figure 5, in which the
Vb , t − 1 + Figure 5. Illustration of the transformed naphtha scheduling results.
∑ FLa ,b ,t − ∑ a ∈ As
event slot is mainly separated by the decoking operation in each furnace. According to the flow rate variation, the scheduling result in the 30 days can be transformed into nine events.
a′ ∈ A s ,
b = a′ ,
1 < t < Tn 6484
b ′∈ Bd
FLa
′,b′,t
= Vb , t ,
(ab), (a′b′) ∈ CN,
∀ b ∈ Bd , t ∈ T, (2.1)
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
Component concentration is illustrated in eq 2.2, which is a bilinear term and can be linearized by the formula proposed by McCormick30 and keep this model a MILP problem. The mass balance of the component composition is presented in eq 2.3. The final blending products of each charging tank should meet the quality constraints of the cracking operation; thus, the blending products can be seen as the same type of material. For more details of the blending constraints, the reader can find them in the Supporting Information. FLa , b , t Pa , p = FLPa , b , p , t ,
∀ a ∈ Aa ,
(ab), (a′b′) ∈ CN,
CVb , p , t = CVb , p , t − 1 +
p ∈ P,
t∈T
(2.2)
∑ FLPa ,b ,p,t − ∑ a ∈ As
∀ b ∈ Bd ,
b ∈ Bd ,
a′ ∈ A s ,
(ab), (a′b′) ∈ CN ,
The rest equations are used to control the time sequence of the material transportation. Since the continuous time-formulation is event based, the time variables should also be aligned with the events of material flow with different sources and destinations. The other time sequence and operation logistic constraints include inventory safety level constraints with simultaneous input and output flow, nonsimultaneous input and output constraints in the charging tanks, demand order sequence constraints etc. All of these equations can be found in the Supporting Information eqs 48−66. The objective function minimizes the total cost of the operation in tank area. The first term is the cost of unloading operation, and the second term represents the total the sea waiting cost of the arriving vessel. The following three terms are used to punish variation of the demand charging flow rate. The last two terms consider the inventory cost at two points in per time event of the storage tanks and charging tanks.
b ′∈ Bd
FLPa
′,b′,p,t
,
b = a′ , p ∈ P,
t ∈ T,
min ∑ CLi(Tvfi − Tvsi) +
t>1
The summation of the material flows from the charging tank to the cracking furnace system during time slot t for the demand event f is calculated in eq 2.4. DMf represents the average feeding flow rate to the furnace system in the demand event f and the term Tfc,f,t − Tsc,f,t represents time duration of the charging flow from charging tank c for demand event f.
∑ FLc ,f ,t = DMf ∑ (Tfc ,f ,t c∈C
− Tsc , f , t ),
+ +
∑ FLi ,s ,t ) + ∑
t ∈ T , t < Tn
H (2NT + 1)
Vcc , t − 1 +
2Vccinit
∑ Cinvcc(∑ Vcc ,t
+
c
t
∑ FLs ,c ,t ) t ,s
(2.10)
4.2. Solution Framework. In order to give a precise result of the blending process and minimize the composition discrepancy in the charging tanks, we also propose an iterative framework with an NLP model which generated from the original relaxed MILP model by fixing the binary variables and reformulate the bilinear terms. The objective of the NLP model is minimizing the summation of these slack variables with fixed scheduling orders from the original MILP problem’s result, and the variation variable Var, which stands for the total variation level should be within a certain range. min Var =
∑
(FLPsp, c , p , t + FLPsn, c , p , t )
s ,c ,p,t
− Tsc , f , t ),
+
∑
(FLPcp, f , p , t + FLPcn, f , p , t )
c ,f ,p,t
∀ f ∈ F,
+
∑ (CVssp,p,t + CVssn,p,t ) + ∑ (CVccp,p,t + CVccn,p,t ) s ,p,t
c ,p,t
(2.11)
∀f∈F
∀f∈F
(2.7)
If the new NLP model cannot be solved, which means the current solution cannot meet the specified component quality, the integer cut is applied in the process to remove the current solution. The MILP and NLP model will be solved iteratively until we get the feasible result. The detailed reformulated equations will be illustrated in the Supporting Information and the flowchart of this method is presented in Figure 6.
(2.8)
Tfc , f , t − H(1 − wvc , f , t ) ≤ Tdfvf ≤ Tsc , f , t + H(1 − wvc , f , t ), t∈T
t
H (2NT + 1)
Vss , t − 1 + 2Vssinit
t ∈ T , t < Tn
(2.6)
f ∈ F,
∑ Cinvss(∑ Vss ,t + ∑ t ,i
∀ f ∈ F,
t∈T
∀ c ∈ C,
c ,f ,t
s
c∈C
Tds f + svs pf − svs nf = Tdsvf ,
f
i ,s,t
(2.5)
Tdf f + svf fp − svf fn = Tdfvf ,
∑ FLc ,f ,t ) c∈C
+ Co ∑ FLi , s , t + Cset ∑ wvc , f , t +
c∈C
t∈T
c∈C
− Tsc , f , t ) −
f
Logical Constraints. The logistic constraints are used to control the transportation in the tank area and assign the specific time point for every shipment operation. These constraints can be referred to the Supporting Information. Time Constraints. The model also includes the time constraints in the continuous formulation, which is used to decide the exact operation sequence and time duration for the inventory management. However, as we have discussed before, the incomplete time information could lead to an infeasible status of the inventory management problem. To overcome this situation and provide the related information for the next stage’s rectification process, we use some slack variables to allow the variation in the demand flow rate or change the starting and ending time of different event orders. This constraint is presented in eqs 2.5−2.9.
∑ FLc ,f ,t ≤ DMf ∑ (Tfc ,f ,t
∑ (Tfc ,f ,t c∈C
+ MF ∑ (svf fp − svf fn ) + MS ∑ (svs fp − svs nf )
(2.4)
− Tsc , f , t ),
MPf (DM f
f ∈F ,t∈T
×
∑ FLc ,f ,t ≥ DMLf ∑ (Tfc ,f ,t
∑
+
c∈C
c∈C
i
∀ f ∈ F,
t∈T
∑ Cseai(Tvsi − Tarri)
i
(2.3)
(2.9) 6485
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
between the two sites. The basic logical flowchart of this subcooperative model and framework is represented in Figure 7. As shown in the Figure 7, the initial information for the subcooperative event information comes from the monthly result of the furnace cracking operation and actual possible supply feed flow comes from naphtha inventory management. First, according to the provided information, if the demand order from the furnace cracking operation cannot be satisfied, the naphtha inventory management problem will be solved with additional slack variables. Then these slack variables will be max u l transformed into the deviation level of flvmin e , flve , tvfe , tvfe, which representing the maximum or minimum flow rate and time point deviation level in the new event-based model. When the current month’s optimization problem is feasible on both sides and meets the terminal criterion as the deviation level is within an acceptable level, the naphtha inventory management problem for the following month is calculated in a rolling horizon framework. Based on the transformed variables, two methods are usually adopted to match the flow rate supply and demand information from both sides. If the given demand order only has a slight deviation, we use a heuristic approach to select the suitable furnace and adjust its original operation batch time. Otherwise, a subcooperative model with serial cracking operation is introduced to cope with this situation. In this model, the new naphtha oil demand order is generated based on the criteria that minimize the deviation level and optimize the monthly cracking profit in the furnace groups. These event demand orders of naphtha oil will be returned to the upstream inventory management problem again to test its feasibility. In the next step, once the current month’s supply meets the demand orders, some of the decision variables are returned to the original cyclic scheduling mode with fixed binary decision variables of the former and current month schedules to optimize the cyclic model again under new circumstances. The new order from the next following month will be sent to the feedstock area for the inventory management. Finally, this framework is
Figure 6. Flowchart of the solution process.
5. BALANCING OF THE UPSTREAM AND DOWNSTREAM PROBLEMS 5.1. Synchronized Scheduling Framework. In the beginning of the paper, we indicated that due to asymmetric information of the full space long-term scheduling model, the inventory management problem often generate unsatisfied solutions of the demand orders. The situations like shipment delay or insufficient supply naphtha may cause a serious impact on the furnace cracking process when the running capacity is always under variation. In order to overcome these difficulties, a subcooperative problem is proposed to incorporate the former two scheduling models into a synchronized optimization framework and balance the inconsistent feeding and supply flow
Figure 7. Flow logic of the detailed step in the synchronized framework. 6486
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
calculated iteratively until the optimization results are all aligning with the proper supply chain information and scheduling orders. 5.2. Subcooperative Problem. With the scheduling results of the upstream and downstream problem as the initial parameter, the subcooperative model will adjust the furnace cracking operation with new starting time and ending time and generate new demands for the naphtha feeds in the current month. The event demands are calculated as we have suggested in section 3.3. Besides, considering the maintenance of the robustness of the original cyclic scheduling orders, these time point of the cracking time line should not be deviated from the original point too much. Serial cracking usually occurs under the situation like insufficient or over supply, and the furnace will substitute its current cracking material with another feed while maintaining the same cracking condition. Nevertheless, this situation could lead to an increasing coking thickness and change the coefficient of coking velocity. Thus, the scheduling result should be reoptimized. In the subcooperative model, this serial-cracking operation mode of furnace operation is used to overcome the aforementioned infeasible situation within the permitted time duration. The mathematical model is given as follows: First we assume there are several batch slot n for each furnace j under naphtha cracking mode i. These batches are generated from the scheduling results within a month. Each serial cracking operation can be added into this batch slot under specific situation. The duration of serial cracking s in slot n of furnace j is represented in eq 3.1. Equations 3.2 through 3.4 illustrate the time constraints of the serial cracking operation. It is important to note that the serial cracking operation only occurs in the s batches cracking single naphtha feed INP rather than the c cocracking mode INP to keep the current cracking operation stable. TSHs , j , n = Tsfs , j , n − Tsss , j , n ,
∀ s ∈ S,
mml , s , j , nTLSi , j , n ≤ TSHMs , j , n ≤ mml , s , j , nTUSi , j , n , sc ∀ l ∈ L NP ,
∑
s ∀ i ∈ INP ,
Tsfs , j , n ≤ Tfvi , j , n ,
s i ∈ INP ,
Tsss + 1, j , n ≥ Tsfs , j , n ,
j ∈ J, j ∈ J,
∀ s ∈ S,
n∈N n∈N
j ∈ J,
n∈N
i∈
s INP ,
s ∈ S,
∑ sms ,j ,n ≤ NSj ,n ,
∀ s ∈ S,
i∈
s ∈ S,
s ∈ S,
sc l ∈ L NP
s i ∈ INP ,
∀ e ∈ E,
n∈N
(3.10)
j ∈ J,
j ∈ J,
∀ e ∈ E,
n∈N
(3.11)
∀ e ∈ E,
n∈N
(3.12)
(y′s , e , j , n − 1)M ≤ Tefve − Tsfs , j , n ≤ y′s , e , j , n M , s ∈ S,
j ∈ J,
∀ e ∈ E,
n∈N
(3.13)
It can be deduced that the exact time duration during the event slot e of serial cracking s is usually equal to the minimum value of the three terms Tef ve − Tsss,j,n, Tsfs,j,n − Tesve and Tsfs,j,n − Tsss,j,n. Equations 3.14 through 3.22 are used to calculate the amount of decreased naphtha feeding flow FLs,e,j,n during the current event order e after introducing the serial cracking operation s. Then average flow rate of demand event e is adjusted according to eq 3.23 and eq 3.24.
(3.2) (3.3) (3.4)
FLs , e , j , n ≤ FIl , i , j(Tsfs , j , n − Tsss , j , n) + (4 − xs , e , j , n − x′s , e , j , n − ys , e , j , n − y′s , e , j , n ), s ∈ S,
j ∈ J,
∀ l ∈ L NP ,
i ∈ INP,
n∈N
e ∈ E, (3.14)
FLs , e , j , n ≥ FIl , i , j(Tsfs , j , n − Tsss , j , n) − (4 − xs , e , j , n − x′s , e , j , n − ys , e , j , n − y′s , e , j , n ), s ∈ S,
j ∈ J,
∀ l ∈ L NP ,
i ∈ INP ,
n∈N
e ∈ E, (3.15)
n∈N FLs , e , j , n ≤ FIl , i , j(Tefve − Tsss , j , n) + (3 − xs , e , j , n − x′s , e , j , n
(3.6)
mml , s , j , n = sms , j , n ,
j ∈ J,
(ys , e , j , n − 1)M ≤ Tsfs , j , n − Tesve ≤ ys , e , j , n M ,
s
∑
(3.9)
(x′s , e , j , n − 1)M ≤ Tefve − Tsss , j , n ≤ x′s , e , j , n M ,
(3.5) s INP ,
s i ∈ INP ,
(xs , e , j , n − 1)M ≤ Tsss , j , n − Tesve ≤ xs , e , j , nM ,
∀ s ∈ S,
n∈N
s ∈ S,
(3.8)
For the upstream scheduling model, the exact information we need is based on the demand events of the furnace system. Four auxiliary binary variables xs,e,j,n, xs,e,j,n ′ , ys,e,j,n, ys,e,j,n ′ is introduced to indicate the relationship of the starting time and ending time between the serial cracking operation and the original naphtha demand event e. For example, if the starting time of serial s in furnace j in batch slot n is smaller or equal to the starting time of event e, then the xs,e,j,n = 1. The related mathematical formulation of this kind of constraint is presented in eqs 3.10−3.13.
In order to keep the robustness of the original plan, eq 3.5 suggests that the time duration of serial cracking should be within a specific time range. Besides, in consideration of the adverse affection the serial cracking operation could bring, eq 3.6 gives the maximum number of this operation during one batch slot. In eq 3.7, the material l ∈ LscNP represents the feeding material that can be successively cracked in the same furnace condition with naphtha oil. Equations 3.8 and 3.9 give the time duration for each cracking material l in serial cracking s. sms , j , nTLSi , j , n ≤ TSHs , j , n ≤ sms , j , nTUSi , j , n ,
TSHMl , s , j , n = TSHs , j , n ,
n∈N
n∈N
(3.1)
Tsss , j , n ≥ Tsvi , j , n ,
s i ∈ INP ,
sc l ∈ L NP
j ∈ J,
n∈N
s ∈ S,
− ys , e , j , n + y′s , e , j , n )M ,
n∈N
e ∈ E,
(3.7) 6487
s ∈ S,
∀ l ∈ L NP ,
j ∈ J,
n∈N
i ∈ INP , (3.16)
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research
Article
∑ ∑ ∑ FLs ,e ,j ,n ≤ (Fleoe − Fleel)Tehe
FLs , e , j , n ≥ FIl , i , j(Tefve − Tsss , j , n) − (3 − xs , e , j , n − x′s , e , j , n − ys , e , j , n + y′s , e , j , n )M , e ∈ E,
s ∈ S,
∀ l ∈ L NP ,
j ∈ J,
s∈S j∈J n∈N
i ∈ INP ,
n∈N
+ Tehel(Fleoe − Fleve) − Tehel(Flee − Fleel),
(3.27)
∑ ∑ ∑ FLs ,e ,j ,n ≤ (Fleoe − Fleeu)Tehe
FLs , e , j , n ≤ FIl , i , j(Tsfs , j , n − Tesve) + (3 + xs , e , j , n − x′s , e , j , n − ys , e , j , n − y′s , e , j , n )M , e ∈ E,
s ∈ S,
∀ l ∈ L NP ,
j ∈ J,
s∈S j∈J n∈N
+ Teheu(Fleoe − Fleve) − Teheu(Flee − Fleeu),
i ∈ INP ,
n∈N
e ∈ E,
s ∈ S,
∀ l ∈ L NP ,
j ∈ J,
In addition, the time point of demand orders can be adjusted according to the variation level provided by the scheduling result of naphtha oil inventory management. In this respect, the starting time and ending time is reformulated in eq 3.29 and eq 3.30 with the auxiliary variables Tefpe , Tefne , Tespe , Tesne . Tefpe and Tespe represent the positive part of the variation for the original finishing time Tefoe and starting time Tesoe, Tefne , and Tesne stand for the negative part of the variation. Besides, the finishing time for the last demand event e and the starting time for the first event are fixed. However, if the demand event e is generated from the feed flow rate variation of cleanup operation in furnace j of batch slot n, which is presented by Eclj,n, the duration time is equal to the cleanup time Tsvi′,j,n+1 − Tfvi,j,n. This constraint is illustrated in eq 3.31.
FLs , e , j , n ≤ FIl , i , j(Tsfs , j , n − Tsss , j , n), i ∈ INP ,
e ∈ E,
i ∈ INP ,
n∈N
s ∈ S,
(3.19)
∀ l ∈ L NP , j ∈ J,
n∈N
(3.20)
−(2 + xs , e , j , n − x′s , e , j , n + ys , e , j , n − y′s , e , j , n )M ≤ FLs , e , j , n ≤ (2 + xs , e , j , n − x′s , e , j , n + ys , e , j , n − y′s , e , j , n )M , ∀ l ∈ L NP ,
i ∈ INP ,
e ∈ E,
s ∈ S,
∀e∈E (3.28)
(3.18)
FLs , e , j , n ≥ FIl , i , j(Tsfs , j , n − Tesve) − (3 + xs , e , j , n − x′s , e , j , n − ys , e , j , n − y′s , e , j , n )M ,
∀e∈E
(3.17)
Tefv e = Tefoe + Tefep − Tefen ,
∀ e ∈ E,
e < En (3.29)
j ∈ J,
Tesve = Tesoe + Tesep − Tesen ,
(3.21)
n∈N
∀ e ∈ E,
e ≠ E1 (3.30)
−(2 − xs , e , j , n + x′s , e , j , n − ys , e , j , n + y′s , e , j , n )M ≤ FLs , e , j , n
Tefve − Tesve = Tsvi
≤ (2 − xs , e , j , n + x′s , e , j , n − ys , e , j , n + y′s , e , j , n )M , ∀ l ∈ L NP ,
i ∈ INP ,
e ∈ E,
s ∈ S,
n ∈ N,
j ∈ J,
Fleoe(Tefve − Tesve) −
∑ ∑ ∑ FLs ,e ,j ,n s∈S j∈J n∈N
= Fleve(Tefve − Tesve),
Fleel ≤ Fleve ≤ Fleeu ,
∀e∈E
(3.23)
∀e∈E
Tesve = Tfvi , j , n ,
(3.24)
Tesve = Tsvi , j , n , Tesve = TSOi , j , n , n = 1,
s∈S j∈J n∈N
+
− Fleve) −
−
∀e∈E
e = En ,
(3.25)
∀ e ∈ Ej , n ,s ,
n ∈ N,
∀ e ∈ Ej , n ,s , n ∈ N ,
∀ e ∈ Ej , n ,s ,
i ∈ INP
Tefve = Tesve + 1 ,
∑ ∑ ∑ FLs ,e ,j ,n ≥ (Flevoe − Fleel)Tehe
i ∈ INP i ∈ INP e = 1,
∀ e < En
(3.34)
n ∈ N,
n = N, (3.35) (3.36)
The rest furnaces which process the other feeding materials will also have to reschedule their batch starting and ending time in the current month to avoid simultaneous cracking. These constraints are given in eqs 3.37 and 3.38.
s∈S j∈J n∈N
+ Teheu(Fleoe − Fleve) − Teheu(Flee − Fleel),
n ∈ N,
i ∈ INP
Tefve = TFOi , j , n ,
Fleeu),
∀ e ∈ Ej , n , f ,
(3.33)
∑ ∑ ∑ FLs ,e ,j ,n ≥ (Flee − Fleeu)Tehe Tehel(Flee
(3.31)
(3.32)
The term Fleve (Tef ve−Tesve) in eq 3.23 is a bilinear constraints, it can be linearized through the convex and concave envelopes proposed by McCormick,30 First we rewrite the eq 3.23 using a new variable auxiliary Tehe to represent the term Tef ve−Tesve of the time duration in event e. The linearized formulation is given in eqs 3.25−3.28.
Tehel(Fleoe
∀ e ∈ Ejcl, n ,
Equations 3.32−3.33 are used to present the relationship between the starting/ending time demand order e with the starting/ending time in furnaces that are cracking naphtha oil. This relationship is deduced from the set Ej,n,f and Ej,n,s, which stand for the set of demand e correlating with the ending time and starting time of furnace j in current month’s batch slot n. Equation 3.36 give the sequential constraints for the demand event orders.
(3.22)
n∈N
− Tfvi , j , n , ′,j,n+1 i ∈ INP , i ∈ I
∀e∈E (3.26) 6488
dx.doi.org/10.1021/ie500079w | Ind. Eng. Chem. Res. 2014, 53, 6477−6499
Industrial & Engineering Chemistry Research Tsvi , j , n = TSOi , j , n + Tsip, j , n − Tsin, j , n , n ≠ 1,
Article
∀ j ∈ J,
n ∈ N,
i ∉ INP
SCl , i , j , n = FIl , i , jTsss +
(3.37)
′,j,n
+ FIl , i , j(Tfvi , j , n − Tsfs
) ″,j,n ∀ s , s′ , s″ ∈ S ,
∑ FIl ,i ,j(Tsss + 1,j ,n − Tsfs ,j ,n ) s
Tfvi , j , n = TFOi , j , n + Tfi p, j , n − Tf in, j , n , n ≠ N,
∀ j ∈ J,
s′ = s1,
n ∈ N,
i ∉ INP
l∈
(3.38)
j ∈ J,
l ∈ Lic ,
j ∈ J,
s″ = sn ,
n ∈ N,
Tsss
∫Tsv
+
∑∫
FIl
Tsss + 1, j , n + Tsvi , j , n + SHj , n
Tsfs , j , n + Tsvi , j , n + SHj , n
Tfvi , j , n + SHj , n
∫Tsf
s
+
+ Tsvi , j , n + SHj , n
i , j , n + SHj , n
s
+
′,j,n
maximize
″
+ Tsvi , j , n + SHj , n ,j,n
∑∑∫ s
TSHLl
″,s,j,n
0
l″
FIl
′
FIl
′
+ al
′
−
l′ ∈
Lic ,
SPl , i , j , n =
l″ ∈
∑∫ l′
∀i∈
sc , LNP
I jc
− −
Tsvj , n + SHj , n
−
s i ≠ INP ,
−
∑ cxe(Tefep
+ Tef en ) −
∑ cxe(Tesep + Tesen) e∈E
∑ ∑∑
exj , n(Tsvip, j , n
+ Tsvin, j , n)
∑ ∑ ∑ fx(Tfvip,j , n + Tfvin,j , n) ∑ ∑ ∑ ∑ gx(Tssl , s + 1, j , n − Tsfl , s , j , n ) ∑ ∑ ∑ ∑ kx·TSHLl. s , j , n
−
∑∑∑ s sc i ∈ INP j ∈ J l ∈ LNP
−
∑∑∑ s sc i ∈ INP j ∈ J l ∈ LNP
hx(Tfvi , j , n − Tssl , s , j , n ) ′ ′ ′ ix(Tsfl , s
∑∑∑ ∑ s sc 1
