Decomposition Algorithm for the Scheduling of Typical Polyvinyl

Article pubs.acs.org/IECR

Decomposition Algorithm for the Scheduling of Typical Polyvinyl Chloride Production by Calcium Carbide Method Miaomiao Tian,† Xiaoyong Gao,† Yongheng Jiang,† Dexian Huang,*,†,‡ and Ling Wang†,‡ †

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

‡

S Supporting Information *

ABSTRACT: In our previous work (Tian et al. Ind. Eng. Chem. Res. 2016, 55(21), 6161−6174), a plantwide scheduling model was presented, which was difficult to solve for industrial scale instances in acceptable time. Because the vinyl chloride monomer (VCM) buffer links the continuous process with the batch process, the whole scheduling problem can be decomposed into the upstream VCM production processes and the downstream polymerization processes. Thus, a decomposition algorithm is presented in this paper to accelerate the computation progress. Using the decomposition algorithm, the polymerization scheduling optimization problem is first conducted and thus the detailed VCM demand schedule is obtained. Then, with an off-line model formulated in advance, the operating states (i.e., start/stop operations) of arc furnaces are optimized in the second step, which would be the hard-to-solve binary variables in the plantwide scheduling model. In the off-line work, the furnaces operating range (i.e., the production rate) is discretized into multiple operational levels and the corresponding optimal furnaces selection scheme is then obtained based on the energy consumption model of the arc furnaces. Finally, the determined binary variables are embedded into the plantwide scheduling model and thus a reduced scale scheduling optimization is executed. Computational results show that the proposed algorithm can accelerate the computation greatly and the scheduling results are close to or even better than those given in our previous work.

1. INTRODUCTION The PVC production scheduling was formulated as an integrated MILP model by the calcium carbide method,1 in which both of the VCM production process (continuous process) and polymerization process (batch process) were considered. The nonlinear electricity consumption in terms of production rates for the energy intense process (i.e., VCM production process, including furnaces and electric bathes) was considered and formulated as a piecewise linear model. The simulation results showed that the energy consumption of the plantwide model was reduced by 1.7%. However, the combinatorial complexity of the proposed integrated model makes it difficult to solve when confronted with actual industry scale cases: (1) Multiple facilities for parallel production, that is, multiple electrolytic baths for chlorine production, multiple arc furnaces for calcium carbide, and multiple pots for polymerization, lead to a large quantity of binary variables. (2) Because of the distinct operating characteristics for a continuous process, adopting a smaller time interval for the continuous processes to describe them will result in a large scale continuous process model. (3) There exist excessive complex constraints on production rates and energy consumption on start/stop operations of arc furnaces. As a result, the model is difficult to solve efficiently for industrial application. © XXXX American Chemical Society

There are many research reports and solvers for solving an MILP problem.2 One of the difficulties is the large integer search space caused by binary variables. Commercial MILP solvers require exponential computation time and even fail to find a feasible solution for the large-scale problems.3−6 The decomposition approach is commonly used for solving a large-scale MILP problem. In some literature, the whole problem was decomposed into several subproblems according to the mechanism of the process.2,7−11 Grossmann et al. presented a decomposition strategy to divide the whole scheduling optimization problem into an assignment subproblem and a sequencing subproblem so as to reduce the solving complexity.2 In another work of Grossmann et al., the special features of a steel making process was considered, and the original problem was split into several smaller problems.8 Mahalec et al. analyzed the gasoline blending scheduling problem and decomposed it into recipe optimization and volumes optimization.12,13 Meyr and Mann decomposed the multiline problem into a series of single-line problems, and Received: Revised: Accepted: Published: A

September 1, 2016 October 20, 2016 October 28, 2016 October 28, 2016 DOI: 10.1021/acs.iecr.6b03375 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. PVC production process.

proposed different decomposition approaches for the multiline master problem and the single-line subproblems.9 Khmelnitsky et al. proposed a time-decomposition approach for large-scale scheduling problems to obtain a near-optimal solution by dividing the original problem into three subproblems, namely sequencing, loading, and timing problems.14 The decomposition technology such as Lagrangian decomposition15,16 and Benders decomposition17−19 are widely used. Niknam et al. proposed a Benders decomposition-based approach to solve the thermal unit commitment problem by decomposing the original problem into an integer master optimization problem and a nonlinear optimization subproblem. In the master problem, the on/off states of the generating units were solved.20 In the PVC production process, VCM is stored in buffer tanks as the product of the synthesis reaction of HCL and C2H2 and the material of further polymerization. The whole process can be naturally divided into the VCM production processes and the polymerization process. As pointed out in our previous work,1 the major energy consumption was due to the continuous process. Compared to the continuous process, the batch process is of little potential margin for improvement. Moreover, the major computation burden was costed in determining the operating condition (i.e., operating states and load) of parallel units in upstream VCM production processes. For the sake of accelerating the solving process, a decomposition strategy for the plantwide model is proposed in this paper. First, the polymerization process scheduling that only considers the product demand is optimized and the detailed VCM demand is obtained. An off-line model based on the energy consumption features of the arc furnaces is constructed and assists the optimization of the operating states (i.e., binary variables of start/stop operations) of multiple arc furnaces. Finally, the plantwide scheduling model with the determined binary variables as the known is performed to determine the operating capacities. This decomposition strategy is able to solve the large-scale industrial cases in acceptable time. The rest of the paper is organized as follows. In section 2, the process description and the problem statement are provided. Then, the decomposition approach is presented in details in section 3. In section 4, several well designed cases are provided

to show the effectiveness of the proposed method. Finally, some conclusions are drawn in section 5.

2. PROCESS DESCRIPTION AND PROBLEM STATEMENT This paper studies the same PVC process as our previous work.1 A typical system in Figure 1 is used to illustrate the plantwide PVC production model. In the continuous VCM processes, the chlorine and hydrogen are produced in electrolytic baths. The calcium carbide is produced from arc furnaces and stored in the storage tank after being cooled. The cooled calcium carbide synthesizes with water to obtain acetylene (C2H2), and the obtained acetylene reacts with hydrogen chloride from chlorine and hydrogen to produce VCM. In the batch polymerization process, the unreacted VCM will be recycled from the polymerization pots and stored in the buffer tank. From Figure 1, it can be seen that the continuous process supplies the VCM monomer for the batch process as the material of polymerization. The whole process can be naturally decomposed into the continuous process and the batch process. Moreover, the major solving difficulty is to determine the operating conditions of multiple parallel units in the upstream VCM production part. In addition, the energy consumption and scheduling margin of the arc furnaces are much larger than that of the electrolytic baths. Meanwhile, a large number of binary variables are needed to determine the start/stop operations of the arc furnaces, while a determination of the start/stop operation for the electrolytic baths does not consider the safety risk caused by frequent start/stop operations. That is, excessive binary variables involved in the plantwide scheduling model for determining the arc furnaces operating states are hard to solve, and the number of binary variables increases as the arc furnaces number increases. In our previous work,1 decision variables for the arc furnace operating state (i.e., start/ stop) and operating load (i.e., production rates) are lumped together and determined simultaneously, which makes the industrial scale problem difficult to solve. To accelerate the solving process, a natural idea is to decompose the whole problem into several subproblems that are convenient to handle. In this paper, a decomposition method will be proposed, and the detailed procedures will be introduced in the next section. B

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

Article


3. DECOMPOSITION ALGORITHM The proposed decomposition algorithm involves an off-line furnace combination model and three optimization steps (three MILP subproblems), as illustrated in Figure 2 with the following steps:

Figure 2. Framework of the decomposition algorithm.

Figure 3. Flowchart of the decomposition algorithm.

Step 1: Downstream VCM polymerization process optimization, take the demand order satisfying as the target, and obtain the VCM demand schedule that is expected to be as steady as possible. Step 2: Optimize the start/stop operations of the arc furnaces of the continuous process. Step 3: Plantwide optimization with partially determined information from Step 1 and Step 2. The above algorithm starts at the optimization of the batch process, which determines the demand of VCM and further calcium carbide for the continuous process production. In this step, the operating capacities in the continuous process, that is, upper and lower bounds of the production rates of VCM, calcium carbide, and chlorine, are considered in the model. The penalty of the VCM production fluctuations is introduced to the objective function to reduce the energy consumption caused by the VCM production fluctuation, which is consistent with the objective of the integrated optimization of the whole process. Then the operating states of the arc furnaces are optimized to meet the VCM demand from the downstream polymerization process. The start/stop operations decision of multiple arc furnaces leads to a large quantity of binary variables, which is difficult to solve. So, an off-line furnace combination algorithm is proposed to simplify the start/stop decision of the arc furnaces. Finally, the plantwide model is optimized with the fixed binary variables from the former two steps. In this step, the operation load of the arc furnaces and the electrolytic baths are determined, that is, the calcium carbide and chlorine production rates. The flowchart of the proposed decomposition algorithm is shown as Figure 3, where the objective and the constraints of each step are summarized.

Next, the decomposition algorithm and the model formulations of each step are described. The batch process optimization of step 1 is given in section 3.1. The operating state optimization of step 2 is given in section 3.2, including the off-line furnace combination model proposed to simplify the complexity of the optimization. The plantwide model optimization of step 3 is given in the section 3.3. 3.1. Step 1: the Batch Process Optimization. The batch process optimization is to obtain the arrangement of the polymerization pots, satisfying the constraints around the capacity of continuous process and the PVC product order demand. The VCM monomer as the material of polymerization in batch process is expected to be supplied as stable as possible. So, the production of calcium carbide should be steady and the state change operations of both downstream polymerization and upstream continuous process production should be minimized, including the arc furnaces production. The objective of the batch process optimization is to minimize Cb in eq 1, including the changeover cost (the first term of the right side), stockout cost (the second term of the right side), electricity cost of polymerization (the third term of the right side), material cost of calcium carbide (the fourth term of the right side), and chloride (the fifth term of the right side), and the storage cost of VCM (the sixth term of the right side) and PVC (the seventh term of the right side). The detailed definitions are referred to eqs A66−A70 in the Appendix in the Supporting Information. Moreover, to obtain a steady VCM demand plan, a penalty cost penl(t) is introduced in the last term of the right side to calculate the difference between the VCM production rates and the average VCM production rate apV, defined as eq 2. C


Article

Industrial & Engineering Chemistry Research Cb =

arc furnaces should guarantee the constraints of the inventory of the calcium carbide. The inventory of the calcium carbide at time period t is equal to the inventory at time period t − 1 adding the difference between the calcium carbide production and its consumption during period t. The constraints of inventory are as follows,

∑ ∑ ∑ ch(j , k , t )co t∈T j∈J k∈K

+

∑ ∑ δ(k)(R(i , k) − d(i , k)) i∈I k∈K

+

∑ μ ∑ ∑ Y (j , k , t )E P (j , k ) t∈T

0 ToI ICa(t ) = ICa + pCa (t )σ dl − uCa(t ),

j∈J k∈K

+ ωCa ∑

∑ pCa (g , t )

t∈T g∈G

+ ωCl ∑

t ϵT kϵK

∑ penl(t ) (2)

where Tm is the length of the scheduling horizon, R(i,k) is the demand of order i for product k, I0V is the initial inventory of VCM, and IfV is the final inventory of VCM. The VCM inventory should be larger than the final inventory when the scheduling horizon ends. To avoid nonlinearity, eq 2 is reformulated as eqs 4−9,

∀t∈T

τ

(6)

∀t∈T (7)

If δ(t) = 0, then penl(t) = −(pV(t) − apV)*pec,

penl(t ) + (pV (t ) − apV ) ·pec ≥ Lδ(t ),

∀t∈T ∀t∈T

pToI Ca (t)

∀ t ∈ T, τ ∈ Γ

τ

(14)

penl(t ) − (pV (t ) − apV ) ·pec ≤ U ′(1 − δ(t )),

penl(t ) + (pV (t ) − apV ) ·pec ≤ U ′δ(t ),

∀ t ∈ T, g ∈ G

where pTo Cl (t) denotes the aggregate production rates of all the baths at time t. 3.2. Step 2: Start/Stop Operation Optimization of Arc Furnaces. Multiple facilities of the continuous process (i.e., arc furnaces and electrolytic baths) lead to many binary variables and constraints, which causes the continuous process model of industrial scale to be still difficult to solve in acceptable time. The decision variables of the continuous process are the operating states of arc furnaces and the operating load of arc furnaces and electrolytic baths. The scheduling margin of the energy consumption of the arc furnaces are much larger than that of the electrolytic baths. The binary variables of the operating states of the arc furnaces are especially difficult to solve, while the operating states of the electrolytic baths are not considered in the scheduling problem. So, the binary variables of the arc furnaces are the key considerations. In our previous work,1 the operating state and the load decisions for furnaces are lumped together in one problem, making it difficult to yield the optimal solution in a reasonable time limit. Hence, a separate strategy is proposed for the operating states and the load decision. To coordinate the operating states decision making, an off-line furnace combination model is proposed to convert multiple arc furnaces to a whole “furnace combination” and the binary variables are reduced. Moreover, an approximation strategy based on the off-line model is proposed to optimize the operating states of arc furnaces. 3.2.1. The Off-Line “Ordered Arc Furnace Combination” Formulation. The off-line work is provided to describe the relationship between the operating states of all arc furnaces and the corresponding optimal operating load levels. The arc furnace combination model is formulated based on the rule “the furnace with high efficiency is preferred”, which helps to

(4)

where binary variables δ(t) are introduced as auxiliary variables, and L, U, and ε are constants. δ(t) equals to 1 if pV(t) − apV > 0 holds. ε is a small positive constant, and L ≤ 0, U ≥ ε. If δ(t) = 1, then penl(t) = (pV(t) − apV)·pec,

penl(t ) − (pV (t ) − apV ) ·pec ≥ L(1 − δ(t )),

pTof Ca (t),

∑ pClmin (τ) ≤ pClTo (t ) ≤ ∑ pClmax (τ),

(5)

∀ t ∈ T, U ′ ≥ 0

∑

where = t > N − dl. is the aggregate initial production rates of all the arc furnaces, pTof Ca (t) is the aggregate final production rates of the arc furnaces in the last period of the scheduling horizon, and dl is the length of cool-down time of calcium carbide. The constraints of the chloride production in eq 13 and eqs A38−A42 in Appendix 1 in Supporting Information are used. The chlorine production rate is within its lower and upper bounds,

(3)

pV (t ) − apV ≤ Uδ(t ),

≤

max pCa (g ),

(13)

i ∈ I, k ∈ K

∀t∈T

≤

To pCa (t )

g

pTo Ca (t)

k

pV (t ) − apV ≥ (L − ε)(1 − δ(t )) + ε ,

min pCa (g )

g

where apV is the average production rate of the VCM production that can be calculated by the initial inventories and the demand of orders, and pec is the following constant. i

t > dl (12)

∑

t∈T

apV = (∑ ∑ R(i , k) − I V0 + IVf )/Tm ,

1 < t ≤ dl

(1)

t ϵT

penl(t ) = pec|pV (t ) − apV | ,

ToI ICa(t ) = pCa (t )σ dl + ICa(t − 1)σ − uCa(t ),

To ICa(t ) = ICa(t − 1)σ + pCa (t − dl)σ dl − uCa(t ),

t ϵT

+ ωP ∑ ∑ IP(k , t ) +

(10)

(11)

∑ pCl (τ , t ) + ω V ∑ IV(t )

t ∈ T τ ∈Γ

t=1

(8) (9)

where U′ is a constant, and U′ ≥ 0. Other constraints include allocation constraints (eqs A51− A52 in Appendix 1 in the Supporting Information), changeover constraints (eqs A53−A57 in Appendix 1), inventory constraints (eqs A58−A60 in Appendix 1), and delivery constraints (eqs A61−A62 in Appendix 1). Besides, the constraints of VCM production are also considered, that is, VCM production rate constraints (eq A42 in Appendix 1), recycle and storage constraints (eqs A44−A50 in Appendix 1). Moreover, to calciulate the calcium carbide demand for the batch process, eqs 10−13 are used. pTo Ca (t) is the aggregate production rates of all the arc furnaces. The production of the D


Article


As shown in Figure 5, the set of four arc furnaces can be seen as an “ordered furnace combination”. Different production rate

obtain a scheduling scheme that all arc furnaces can work with the high efficiency operating loads. In a discrete-time representation model, the operating state optimization of arc furnaces can be regarded to choose a best set of arc furnace for calcium carbide production at each time period. The priorities of arc furnaces are determined by their efficiency, that is, energy consumption per unit production. Hence, the set of arc furnaces can be converted to an ordered set, and arc furnaces are sorted based on the efficiency. The operating capacity is divided into a series operating levels. Each level has its specific arc furnaces selection scheme according to the ordered set. Then, the problem of operating state determination of all arc furnaces is converted to the problem of operating level determination, i.e. one of the “furnace combination” schemes, which is of much less complexity. First, multiple arc furnaces should be sorted by their efficiency. As the unit consumption−production rate relationship is nonlinear, the efficiency of the furnace cannot be evaluated by the value at a specific condition. To well reflect the real average efficiency, a reference operation range around the optimal point is used. The avarage efficiency under the reference range is used for evaluation. Figure 4 shows the unit consumption-production rate curve of the arc furnace. To avoid the nonlinearity, the piecewise liner

Figure 5. Diagram of furnace combination.

levels based on different number of arc furnaces can be chosen to meet different demand. Meanwhile, the best efficiency production rates and corresponding energy consumption of each production rate level of the “ordered furnace combination” can be obtained. The energy consumption-production rate relationship of the “ordered furnace combination” still can be seen as a piecewise linear function. The breakpoints of the model are the sum of the optimal production rates of arc furnaces in each production rate level f n (f n = 1, 2, ···, n), denoted as op( f n). op(fn ) =

∑ pCaop (g ),

g ∈ Gf

n

g

where f n is the sequence number of production rate level, and the number of arc furnaces is f n in the production rate level f n; Gfn denotes the set of arc furnaces for production rate level f n; for example, in Figure 5, G3 contains arc furnaces 1−3; pop Ca(g) denotes the optimal production rate of arc furnace g. And cp(f n) is the energy consumption corresponding to production rate op( f n),

Figure 4. Average efficiency of arc furnaces.

model is introduced in our previous work.1 Suppose the formulation of the unit consumption in terms of production rate is E = f(p), where p denotes the production rate. In Figure 4, point A is the optimal point with the corresponding production rate pA and unit consumption f(pA). [pb1, pb2] is the reference range aforementioned. The reference operation range is determined by Ear, which denotes the difference between the unit consumption of the bounds of the range and the lowest unit consumption f(pA). That is, f (pb1 ) = f (pb2 ) = f (pA ) + Ear

(16)

cp(fn ) =

∑ ECaop(g ),

g ∈ Gf

n

g

(17)

Eop Ca(g)

where is the optimal energy consumption corresponding to optimal production rate pop Ca(g). Moreover, the sum of the maximum production rates of all the n arc furnaces should be added to an extra production rate level. It represents the upper bound of the operation load of the furnace combination, which is the (n + 1)th production level of the ordered furnace combination. The production rate and energy consumption of level n + 1 is shown in eqs 18 and 19,

(15)

and the value of Ear is given by combining the actual plant operating practice and the recommended operating conditions, for example, 200kwh/t in Figure 4. Then, the operating state decision problem is reduced to n operating levels selection issue after determining the priority of arc furnaces (n is the number of arc furnaces).

op(n + 1) =

∑ pCamax (g ),

g∈G (18)

g

cp(n + 1) =

∑ ECamax (g ), g

E

g∈G (19) DOI: 10.1021/acs.iecr.6b03375 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Figure 6. Piecewise linearized unit consumption-production rate curves of the arc furnaces.

where Emax Ca (g) is the energy consumption of arc furnace g corresponding to maximum production rate pmax Ca (g). The implementation algorithm of the off-line model is as follows: • Determine the “reference production rate range” of the arc furnace. The difference between the unit consumption of the reference range bounds and the lowest unit consumption f(pA) is first determined, and the upper and lower bounds of the reference range can be obtained. • Calculate the average efficiency of arc furnace. • Arc furnaces are sorted according to the order of the average efficiency, and the “ordered furnace combination” is formed. • Calculate the production rates and energy consumption of different levels of the “ordered furnace combination”. An example of the ordered furnace combination model is given as follows to illustrate the furnace combination modeling: Suppose there are two arc furnaces g1 and g2, and the parameters including their capacity and energy consumptionproduction rate relationships are known. The piecewise linearized unit consumption-production rate curves of the two arc furnaces are shown in Figure 6. The breakpoints of the piecewise linear curve q(g1,1), ..., q(g1,6) and q(g2,1), ..., q(g2,6) (the red diamonds in Figure 6) and corresponding energy consumption f(g1,1), ..., f(g1,6) and f(g2,1), ..., f(g2,6) are shown in Table 1. In addition, the boldface statements in the second column of the table are the optimal production rates of the furnaces, that is, the best efficiency work points. The steps of the offline “ordered furnaces combination” modeling are as mentioned above: 1. Determine the “reference production rate range” of the arc furnace: According to the operating practice and the recommended operating conditions, the value of Ear is determined as 300kwh/t, i.e. the unit consumption of the bounds of the “reference range” is 300kwh/t larger than the lowest unit consumption of the arc furnace. For example, the optimal load of arc furnace g1 is 8.04t/h, and its corresponding unit consumption is 2772.3kwh/t. So the unit consumption of the bounds of the “reference range” is 3072.3kwh/t. With the known energy consumption-production rate curve, the “reference range” can be determined: [4.5, 11.5](t/h). The

Table 1. Breakpoints and Corresponding Energy Consumption of the Arc Furnacesa arc furnace

breakpoints of production rate/t/h

1

q(g1 1) = 4.5 q(g1 2) = 6.24 q(g1 3) = 8.04 q(g1 4) = 9.76 q(g1 5) = 11.42 q(g1 6) = 13.0 q(g2 1) = 4.5 q(g2 2) = 5.95 q(g2 3) = 7.47 q(g2 4) = 9.10 q(g2 5) = 10.98 q(g2 6) = 13.0

2

corresponding energy consumption/kwh f(g1 1) f(g1 2) f(g1 3) f(g1 4) f(g1 5) f(g1 6) f(g2 1) f(g2 2) f(g2 3) f(g2 4) f(g2 5) f(g2 6)

= = = = = = = = = = = =

13422.87 17575.52 22289.04 27892.45 34779.01 43282.02 15668.24 18991.46 23163.10 29518.57 40825.41 59833.19

a

Boldface denotes the optimal production rates of the furnaces, that is, the best efficiency work points.

reference rages are shown in the green shadow areas in Figure 6. 2. Calculate the average unit consumption of the arc furnace: The average unit consumption in the “reference range” can be calculated according to the energy consumption- production rate curve. In this case, the average unit consumptions of furnaces g1 and g2 are 2878.74kwh/t, 3218.75 kwh/t. 3. The arc furnaces are sorted according to the order of the average efficiency, and the “ordered furnace combination” is formed. Sort the arc furnaces according to the average unit consumption obtained in the second step (low to high): g1, g2. Also, the priority of the furnace set is determined. 4. Calculate the production rates and energy consumption of different levels of the “ordered furnace combination”. The production rate levels 1−2 of the furnace combination are load of furnace g1, sum of loads of furnaces g1 and g2. The optimal production rates of the two levels are the sum of the optimal production rates of the arc furnaces in each level: op(1) = 8.04t/h, op(2) = 15.51t/h. In addition, a full-load level is added considering the upper bound of the loads of the two arc F


Article

Industrial & Engineering Chemistry Research furnaces, and op(3) = 26t/h. The corresponding energy consumption can be computed. After that, the individual furnace energy consumption model is replaced with the ordered furnaces combination-based production levels formulation. The combination is determined for the known arc furnaces, and the combination model can be prepared before the scheduling optimization. That is, it needs not to be rebuilt unless the features of facilities (e.g., the energy consumption representation) are changed. Correspondingly, the furnaces operating states decision problem is replaced by the production level optimal selection problem. Clearly, the decision-making complexity is reduced greatly. How to realize the production level-based furnace operating states decision making will be introduced in the following section. An MILP model is formulated to optimize the operating conditions of arc furnaces, which gives a simplified strategy for the problem. 3.2.2. The Operating State (Start/Stop Operation) Optimization Model. On the basis of the proposed off-line formulation, the operating state (start/stop operation) arrangement of the arc furnaces is easy to determine. The objective is to minimize the electricity consumption and the arc furnaces operating states switch cost, as shown in eq 20, Cc =

∑ consump(t ) + ∑ st(t )·stcost, t

Integer variable st(t) is introduced to calculate the number of restart arc furnaces at time period t. Equations 29−33 are equivalent to the expression “If num(t) − num(t − 1) > 0, st(t) = num(t) − num(t − 1), t > 1, else st(t) = 0”. And δ′(t) is an intermediate variable, which equals to 1 if num(t) − num(t − 1) > 0 holds. num(t ) − num(t − 1) ≥ (L − ε)(1 − δ′(t )) + ε ,

num(t ) − num(t − 1) ≤ Uδ′(t ),

t ′< t

st(t ) − (num(t ) − num(t − 1)) ≥ L(1 − δ′(t )),

st(t ) ≤ Uδ′(t ),

∀t∈T (21)

∀t∈T

fn

p(t ) =

(22)

∑ h(fn , t )op(fn ),

∀t∈T

fn

consump(t ) =

(23)

∑ h(fn , t )cp(fn ),

∀t∈T

fn

(24)

∀t∈T

(25)

num(t ) − fn ≥ L(1 − h(fn , t )),

∀t∈T

(26)

If the binary variable h(n + 1,t) equals to 1, it indicates that all the furnaces obtain the maximum production rates at production level n + 1 and the corresponding number of the arc furnaces is n. num(t ) − n ≤ U (1 − h(n + 1, t )),

∀t∈T

(27)

num(t ) − n ≥ L(1 − h(n + 1, t )),

∀t∈T

(28)

(33)

4. CASE STUDY The model is tested on six cases. To make a fair comparison, cases 1−3 are the same as those in our previous work.1 Moreover, to verify the performance of the proposed algorithm, larger scale cases 4−6 are also given. All cases are tested by GAMS win64 24.2.2 on a Dell PC with Intel Xeon Quad CPU processor, 2.50 GHz, Windows Server 2008 R2 OS and 32.0 GB RAM. The MILP models are solved using CPLEX 12.6.0.0. The flowchart is the same as that in Figure 1. The product demands of five different PVC grades for the cases are shown in Table 2. The scheduling horizons are 72 to 120 h (horizon of 72 h is considered for case 1, 96 h for case 2, and 120 h for cases 3−6). Besides, 4 arc furnaces, 4 polymerization pots, and 14 electrolytic baths are considered for cases 1−3, while 8 arc furnaces, 8 polymerization pots, and 28 electrolytic baths are considered in cases 4−6 with a larger size.

Equations 25−28 indicate that “If h( f n,t) = 1, then num(t) = f n”, where num(t) denote the number of working arc furnace in time period t. num(t ) − fn ≤ U (1 − h(fn , t )),

t>1

The energy cost of electrode heating is ignored in this step when the arc furnace is idle, as the cost is very small compared to the operating state switching cost. Given the start/stop operation obtained by the optimization, the production rate should be adjusted according to the practical rules, which was introduced in our previous paper.1 Since this issue is not considered as constraints in the optimization model of step 2, the results can be improved. If a restart furnace shuts down later, time of shut-down operation of the particular furnace should be pushed back for an hour to compensate the production of the first period after restart. In this way, the optimal production rates in the furnace combination are more likely to be obtained. 3.3. Step 3: the Plantwide Process Optimization. The model of this step is the same as the plantwide model in Appendix 1 in the Supporting Information. The binary variables (i.e., operating states of all the units, including start/stop operations of furnaces and polymerization pots) transferred from the step 1 and step 2 optimization are imported into this step as the whole process scheduling optimization. As a result, the majority of discrete variables in the plantwide scheduling model are taken as known, and hence the optimization search space is reduced significantly, and can be solved rapidly. Via this step, the operating conditions, that is, the production rates of arc furnaces and electrolytic baths are optimized, and the final scheduling scheme is obtained.

where p(t′) is the total calcium carbide production rate of all the arc furnaces, and ∑t′1

(20)

t ′< t

(30)

(31)

t>1

The production quantity of the calcium carbide should be larger than the batch process demand calculated from step 1.

∑ p(t′) ≥ ∑ pCaTo (t′),

t>1

st(t ) − (num(t ) − num(t − 1)) ≤ U (1 − δ′(t )),

∀t∈T

t

(29)

t>1

G


Article

Industrial & Engineering Chemistry Research Table 2. Demand Profiles for All Cases case

order

product grades

demand/t

due time/h

1

1 2 3 4 1 2 3 4 1 2 3 4 5 1 2 3 4 5 6 1 2 3 4 5 6 7 1 2 3 4 5 6 7 8 9

k1 k2 k3 k4 k4 k1 k2 k3 k2 k4 k5 k1 k3 k2 k1 k4 k5 k1 k3 k2 k1 k4 k5 k1 k3 k3 k1 k3 k2 k4 k2 k5 k1 k3 k4

450 350 650 700 480 800 700 1100 400 900 400 1300 500 600 1000 500 1200 2500 1500 600 900 500 1800 2500 400 1500 300 500 800 200 1000 1600 600 700 1100

24 48 72 72 24 48 72 96 24 48 72 96 120 24 48 48 72 96 96 24 48 48 72 96 96 120 24 24 48 48 72 72 96 96 120

2

3

4

5

6

Table 4. Model Statistics of the Batch Process Optimization and the Start/Stop Operations Optimization of the Decomposition Models case 1 2 3 4−6

Table 3. Model Statistics of the Plantwide Models equations number

continuous variables number

binary variables number

nonzeros

1 2 3 4−6

24 374 32 572 41 426 78 898

13 142 18 085 22 859 43 259

11 232 14 976 19 440 38 640

104 707 142 163 185 749 348 497

Step Step Step Step Step Step Step Step

1 2 1 2 1 2 1 2

continuous variables number

binary variables number

nonzeros

2 675 1 611 3 561 2 139 5 127 2 667 7 167 3 627

1 259 362 1 683 482 2 357 602 2 357 602

1 800 576 2 400 768 3 720 960 7 320 1 440

11 187 9 671 14 971 15 191 23 357 21 863 37 701 25 223

Next, the decomposition algorithm is compared to the original full-space formulation in our previous work.1 Table 5 shows the objective and solution time of the proposed decomposition algorithm (denoted as “algorithm” in the table) and those required by commercial solver for the corresponding MILP original full-space model (denoted as “solver (1h)” and “solver (10h)” in the table). For the decomposition algorithm, the stopping criterion of Step 1 optimization is that the optimality gap is smaller than 0.5%, and the stopping criteria of Step 2 and Step 3 are that the optimality gap is 0. The stopping criterion of the “solver (1h)” is that the optimality gap is smaller than 1% or the maximum computation time is no more than 1 h. The stopping criterion of the “solver (10h)” is that the optimality gap is smaller than 1% or the maximum computation time is no more than 10 h. Relative and absolute differences between the decomposition algorithm results and the results of the full-space model by commercial solver are given in Table 5. The absolute difference is determined by the results of the decomposition algorithm minus the results of the commercial solver. So, the value is negative if the result of the decomposition algorithm is smaller than that of the commercial solver. The relative difference is equal to the absolute difference divided by the commercial solver results. Because of the real-time response in practical use, the solution time of 1 h is relatively more acceptable. From the comparative results of the computation time, it is clear that the CPU time of the proposed algorithm is less than 10% of that of the commercial solver. Moreover, the objectives of cases 3−6 by the proposed algorithm are even smaller than that of the commercial solver, which implies that the solution of the decomposition algorithm is better when the model size is larger in a limited and industry-acceptable computation time to satisfy the demand of the real-time scheduling. The results of “solver (10h)” is used to understand the accuracy of results of the proposed decomposition algorithm. In all cases the obtained objective of the proposed algorithm is very close to that of the “solver (10h)” (for most cases the difference is smaller than 0.1%), and for case 6 the objective of the decomposition algorithm is even better. As for solution time, the former costs less than 1% CPU time of the latter for all the cases. Table 6 shows the computation time of each subproblem of the decomposition algorithm. In step 1 (batch process scheduling), the priority of the changeover variables ch(j,k,t) is set to the highest as the most important variables can be given the highest priority during a branch and bound search. The solution time is reduced by roughly half compared to that without the priority setting. But the computation time of step 1

Table 3 shows the model statistics of the plantwide models of all the cases. Table 4 shows the model sizes of the subproblems

case

subproblems

equations number

of the decomposition models, including the batch process optimization model (called “Step 1” in the table) and the start/ stop operations optimization model (called “Step 2” in the table). In the last step the plantwide model is optimized with variables assigned by the former 2 subproblems. It is worth mentioning that the model size of start/stop operation subproblem is reduced with the off-line strategy in section 3.2. The parameters used can be found in Table S1 in Supporting Information of our previous work,1 including production rate limits of calcium carbide, chlorine, and VCM production, capacities of polymerization pots, and conversion rates of all grades of PVC in a polymerization process, which are originated from the actual production. H


Article

Industrial & Engineering Chemistry Research Table 5. Comparisons between Results of the Proposed Algorithm and the CPLEX Solver difference case

objective (¥)

1

algorithm solver (1h) solver (10h) algorithm solver (1h) solver (10h) algorithm solver (1h) solver (10h) algorithm solver (1h) solver (10h) algorithm solver (1h) solver (10h) algorithm solver (1h) solver (10h)

2

3

4

5

6

absolute value (¥)

4127755 4126639 4126639 6092326 6087843 6087842 6958933 6964432 6950297 13556120 13610087 13554970 15422212 15431788 15397423 12456538 12504851 12463342

difference

relative value (%)

1116 1116

0.03 0.03

4483 4484

0.07 0.07

−5498 8636

−0.08 0.12

−53967 1150

−0.40 0.01

−9576 24789

−0.06 0.16

−48313 −6804

−0.39 −0.05

CPU time (s) 42.7 3345.1 3345.1 14.0 3600 6997.7 45.3 3600 9780.0 39.2 3600 36000 89.3 3600 36000 343.3 3600 36000

absolute value (s)

relative value (%)

GAP (%)

−3302.4 −3302.4

−98.7 −98.7

1 1

−3586.0 −6983.7

−99.6 −99.7

1 1

−3554.7 −9734.7

−98.7 −99.5

1.8 1

−3560.8 −35960.8

−98.9 −99.9

3.1 2.6

−3510.7 −35910.7

−97.5 −99.8

2.6 2.3

−3256.7 −35656.7

−90.5 −99.0

4.8 4.6

Table 6. Execution Times of Each Step of the Decomposition Algorithm CPU time (s) case

step 1

step 2

step 3

1 2 3 4 5 6

39.2 10.4 37.8 29.9 67.9 330.6

0.5 0.6 3.8 4.2 14.9 5.8

3.0 3.0 3.7 5.2 6.5 6.8

is still much larger than that of step 2 (optimization state optimization) and step 3 (plantwide process optimization with fixed binary variables), because the size of the batch process model is larger than that of the others that are not approximated or simplified. So, further research on the batch process model optimization can be carried out to reduce the solution time. As for the production rate of calcium carbide of case 5, Figure 7 shows the result of decomposition algorithm, and Figures 8 and 9 show the results of the commercial solver. Furthermore, Figure 8 shows the result of 1 h computation

Figure 8. Production rate of calcium carbideresult of the full-space model by the commercial solver (1 h computation time).

Figure 9. Production rate of calcium carbideresult of the full-space model by the commercial solver (10 h computation time).

time, and Figure 9 shows the results of 10 h computation time. The horizon lines in the figures denote the production rates of the production levels of the ordered furnace combination, compared to the situation that the start/stop operations of

Figure 7. Production rate of calcium carbideresult of the decomposition algorithm. I


Article


there exists buffer storage between the subprocesses. The proposed off-line and online combined strategy can be extended to those systems with parallel multiple facilities. In addition, the workload or other characteristics of the systems can be analyzed in the off-line work as the assistance to reduce the complexity of scheduling optimization.

furnaces can be observed. Taking Figure 7 as an example, the working furnace number increases to 8 from 7 at the first time period, and the eighth furnace shuts down at the 28th time period. Then, the other 7 working arc furnaces keep on working until the 83rd time period, and one more arc furnace shuts down that is the lowest efficiency furnace of the 7 arc furnace. At the 97th time period, another 2 furnaces shut down, and the number of working furnace turns to 4. Similarly, the 10 h computation time plantwide model solution by the commercial solver is shown in Figure 9, where the eighth furnace starts at the sixth time period, and the furnace shutdown operations occur at the 78th and the 97th time periods. The production curve shown in Figure 8 with the 1 h computation time solution of the plantwide model is more fluctuant than those in Figure 7 and Figure 9. That is, start-up and shut-down operations occur frequently in the time horizon. Figure 10 shows the production curve of VCM by the decomposition algorithm. It can be seen that the production of

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b03375. Equations of the plantwide model of our previous work (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +86 10 62784964. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the National High-tech 863 Program of China (No. 2013AA 040702), the National Natural Science Foundation of China (No. 61273039), and the National Science Fund for Distinguished Young Scholars of China (No. 61525304).

■ Figure 10. Production of VCM.

■

VCM is smooth and steady. The production rate is almost constant in the time horizon.

5. CONCLUSIONS The PVC plantwide scheduling model is a large-scale MILP problem. It is an open challenge to solve the problem efficiently, especially for the practical application. In this paper, a decomposition algorithm is proposed with three sequential steps. Step 1 optimizes the batch process considering the stable production of the VCM. The demand of calcium carbide is determined for step 2 (the optimization of operation states of arc furnaces). The arrangement of polymerization pots is determined for step 3 (the plantwide process optimization). The off-line model is proposed based on the priority of the efficiency of arc furnaces, which can be preformulated. The operating states of arc furnace are determined by optimization using the off-line model, and the plantwide optimization can be done with the fixed binary variables. Computational results show that the objectives of the proposed decomposition algorithm are close to or even better than those required by the mature solver for the corresponding full-space model. The decomposition algorithm can obtain a good enough solution in less than 10% of the time required by a MILP solver. The proposed decomposition strategy can also be considered for other classes of scheduling problems with the following characteristics: (1) the whole process is composed by several separable and sequential subprocesses; (2) some component plays a dominant role for the optimization margin; (3) and

ABBREVIATIONS PVC = polyvinyl chloride VCM = vinyl chloride monomer MILP = mixed integer linear programming NOMENCLATURE g = arc furnace i, i′ = demand order j = polymerization pot k, k′ = PVC production grade lCa(g) = sequence number of partitioning point in piecewise linearization of calcium carbide energy consumption representation of arc furnace g lCl(τ) = sequence number of partitioning point in piecewise linearization of chlorine energy consumption representation of electrolytic bath τ t, t′, t1, t2, t3 = time period τ = electrolytic bath

Sets

G = polymerization pots I = orders J = arc furnaces K = PVC grades T = time periods T1 = time periods of batch process Γ = electrolytic baths Subscripts

C = C2H2 Ca = calcium carbide H = HCL P = PVC V = VCM J


Article

Industrial & Engineering Chemistry Research Parameters

max pmin V , pV

the average production rate of the VCM production energy consumption corresponding to production rate op(f n) cy the polymerization cycle dl length of cool-down time of calcium carbide dt(i) due time of order i e the limit of electricity supply in each time period Ear the difference between the unit consumption of the bounds of the range and the lowest unit consumption Emax the upper bound of the energy consumption Ca (g) of arc furnace g Eop (g) the optimal energy consumption correCa sponding to optimal production rate pop Ca(g) EP(j, k) the electricity cost of polymerization in pot j for product k f Ca(qCa(lCa(g), g) energy consumption at the production rate qCa(lCa(g), g) f Cl(qCl(lCl(τ), τ)) energy consumption at the production rate qCl(lCl(τ), τ) hp the cost of heating power when the furnace is idle I0Ca the initial inventory of calcium carbide IfCa the final inventory of calcium carbide IlCa, IuCa the lower and upper bounds of calcium carbide inventory IuP upper bound of inventory of PVC I0V the initial inventory of VCM IfV the final inventory of VCM IlV, IuV the lower and upper bounds of VCM inventory Ini_lenCa(g) Ini_lenCa(g) = n denotes that the furnace has been idle for n hours before the scheduling horizon begins Ini_onCa(g) Ini_onCa(g) = n denotes that the furnace has been working for n hours before the scheduling horizon begins L a small negative number N number of time periods in scheduling horizon nCa(g) number of partitioning points in piecewise linearization of calcium carbide energy consumption representation of arc furnace g nCl(τ) number of partitioning points in piecewise linearization of chlorine energy consumption representation of electrolytic bath τ op( f n) sum of the optimal production rates of the arc furnaces in each production rate level f n ( f n = 1,2,···, n) max pmin the upper and lower bounds of the Ca (g), pCa (g) production rate of arc furnace g pop optimal production rate of arc furnace g Ca(g) max pmin the upper and lower bounds of the Cl (τ), pCl (τ) production rate of electrolytic bath τ pICa(g, t) the production rate of furnace at dl time periods before t pfCa(g, t) the production rate of furnace at the last dl time periods of the scheduling horizon pec a constant of penalty

qCa(lCa(g), g)

apV cp(f n)

qCl(lCl(τ), τ)

R(i, k) stcost U U′ V(j) α β ε γ η ρ(k) λ1 λ2 σ μ ξ ωCa ωCl ωP ωV

the upper and lower bounds of VCM production rate partitioning point of number lCa(g) in piecewise linearization of calcium carbide energy consumption representation of arc furnace g partitioning point of number lCl(τ) in piecewise linearization of chorine energy consumption representation of electrolytic bath τ demand of order i for product k the energy cost of a start-up operation a positive number, U ≥ ε a large positive number the charging capacity of polymerization pot j the mass ratio of chlorine and HCL in theory the mass ratio of C2H2 and consumption of calcium carbide in theory a small positive constant the mass ratio of HCL and VCM in theory the mass ratio of HCL and C2H2 in theory conversion rate of product k production rate after restart growing velocity of production rate the weathering loss ratio of per unit inventory at per unit time the price of electricity increasing time to the maximum production rate the raw material cost of unit calcium carbide the raw material cost of unit chlorine the storage cost of unit PVC per unit time the storage cost of unit VCM per unit time

Continuous Variables

consump(t) = electricity cost d(i, k) = delivery of production k for order i ECa(g, t) = energy consumption of furnace g at time t ECl(τ, t) = energy consumption of bath τ at time t f αCa(g, t) = energy consumption at the production rate qαCa(g, t) ele(t) = electricity consumption at time t ICa(t) = inventory of calcium carbide at time t IP(k, t) = inventory of product k at time t IV(t) = inventory of VCM pC(t) = production rate of C2H2 pCa(g, t) = production rate of calcium carbide of arc furnace g at time t pTo Ca (t) = the aggregate production rates of all the arc furnaces pToI Ca (t) = the aggregate initial production rates of all the arc furnaces pTof Ca (t) = the aggregate final production rates of the arc furnaces in the last period of the scheduling horizon pCl(τ, t) = production rate of chlorine of electrolytic bath τ at time t pTo Cl (t) = the aggregate production rates of all the baths at time t penl(t) = a penalty cost to compute the difference between the VCM production rates and the average VCM production rate pH(t) = production rate of HCL pV(t) = production rate of VCM qαCa(g, t) = a positive production rate of arc furnace g K


Article


(11) Khmelnitsky, E.; Kogan, K.; Maimon, O. A time-decomposition method for sequence-dependent setup scheduling under pressing demand conditions. IEEE Trans. Autom. Control 2000, 45 (4), 638− 652. (12) Castillo, P.; Mahalec, V.; Kelly, J. Inventory pinch algorithm for gasoline blend planning. AIChE J. 2013, 59 (10), 3748−3766. (13) Castillo, P.; Mahalec, V. Inventory pinch based, multiscale models for integrated planning and scheduling-part I: Gasoline blend planning. AIChE J. 2014, 60 (6), 2158−2178. (14) Gupta, A.; Sivakumar, A. Job shop scheduling techniques in semiconductor manufacturing. International Journal of Advanced Manufacturing Technology 2006, 27 (11−12), 1163−1169. (15) Mouret, S.; Grossmann, I.; Pestiaux, P. A new Lagrangian decomposition approach applied to the integration of refinery planning and crude-oil scheduling. Comput. Chem. Eng. 2011, 35 (12), 2750−2766. (16) Cordeau, J.; Soumis, F.; Desrosiers, J. A Benders decomposition approach for the locomotive and car assignment problem. Transportation Science 2000, 34 (2), 133−149. (17) Mercier, A.; Cordeau, J.; Soumis, F. A computational study of Benders decomposition for the integrated aircraft routing and crew scheduling problem. Computers & Operations Research 2005, 32 (6), 1451−1476. (18) Shi, L.; Jiang, Y.; Wang, L.; et al. Refinery production scheduling involving operational transitions of mode switching under predictive control system. Ind. Eng. Chem. Res. 2014, 53 (19), 8155−8170. (19) Niknam, T.; Khodaei, A.; Fallahi, F. A new decomposition approach for the thermal unit commitment problem. Appl. Energy 2009, 86 (9), 1667−1674.

Rec(t) = the quantity of recycled VCM at time t tr(k, t) = delivery of product k at time t uCa(t) = consumption of calcium carbide at time t uV(t) = consumption of VCM at time t αCl(lCl(τ), τ, t) = continuous variable in piecewise linearization of chlorine energy consumption representation of electrolytic bath τ αCa(lCa(g), g, t) = continuous variable in piecewise linearization of calcium carbide energy consumption representation of arc furnace g Integer Variables

num(t) = the number of working arc furnace in time period t Binary Variables

ch(j, k, t) = ch(j, k, t) = 1 denotes a changeover from production k to a different product at/after time period t in pot j h(f n, t) = intermediate variables hCl(lCl(τ), τ, t) = binary variable in piecewise linearization of chlorine energy consumption representation of electrolytic bath τ hCa(lCa(g), g, t) = binary variable in piecewise linearization of calcium carbide energy consumption representation of arc furnace g X(j, k, t) = X(j, k, t) = 1 if production k is the first production task in pot j in the time period t, t + 1, ··· N XCa(g, t) = XCa(g, t) = 1 if furnace g is producing at time t Y(j,k, t) = Y(j, k, t) = 1 denotes that a polymerization process for product k begins at time period t in polymerization pot j YCa(g, t) = YCa(g, t) = 1 denotes that start-up of furnace g occurs at time period t δ(t),δ′(t) = intermediate variables δ1(g, t), δ2(g, t) = intermediate variables

■

REFERENCES

(1) Tian, M.; Gao, X.; Jiang, Y.; Wang, L.; Huang, D. A plant-wide scheduling model for the typical PVC production by calcium carbide method. Ind. Eng. Chem. Res. 2016, 55 (21), 6161−6174. (2) Harjunkoski, I.; Grossmann, I. Decomposition techniques for multistage scheduling problems using mixed-integer and constraint programming methods. Comput. Chem. Eng. 2002, 26 (11), 1533− 1552. (3) Li, J.; Karimi, I. Scheduling gasoline blending operations from recipe determination to shipping using unit slots. Ind. Eng. Chem. Res. 2011, 50 (15), 9156−9174. (4) Shah, N.; Ierapetritou, M. Short-term scheduling of a large-scale oil-refinery operations: Incorporating logistics details. AIChE J. 2011, 57 (6), 1570−1584. (5) Kondili, E.; Pantelides, C.; Sargent, R. A general algorithm for short-term scheduling of batch operationsI. MILP formulation. Comput. Chem. Eng. 1993, 17 (2), 211−227. (6) Shin, Y.; Choi, K.; Sakurai, T. Power optimization of real-time embedded systems on variable speed processors. Proc. 2000 IEEE/ ACM Int. Conf. Comput.-Aided Design. 2000, 365−368. (7) Castro, P.; Méndez, C.; Grossmann, I.; et al. Efficient MILPbased solution strategies for large-scale industrial batch scheduling problems. Comput.-Aided Chem. Eng. 2006, 21, 2231−2236. (8) Harjunkoski, I.; Grossmann, I. A decomposition approach for the scheduling of a steel plant production. Comput. Chem. Eng. 2001, 25 (11), 1647−1660. (9) Meyr, H.; Mann, M. A decomposition approach for the General Lotsizing and Scheduling Problem for Parallel production Lines. European Journal of Operational Research 2013, 229 (3), 718−731. (10) Wu, D.; Ierapetritou, M. Decomposition approaches for the efficient solution of short-term scheduling problems. Comput. Chem. Eng. 2003, 27 (8), 1261−1276. L


Decomposition Algorithm for the Scheduling of Typical Polyvinyl

Recommend Documents