A Decomposition Approach for the Solution of ... - ACS Publications

Article pubs.acs.org/IECR

A Decomposition Approach for the Solution of Scheduling Including Process Dynamics of Continuous Processes Jinjun Zhuge and Marianthi G. Ierapetritou* Department of Chemical and Biochemical Engineering, Rutgers University, 98 Brett Road, Piscataway, New Jersey 08854, United States ABSTRACT: Scheduling is often obtained without the consideration of process dynamics that affect the transition between steady states where production takes place. In this work we first formulate the scheduling optimization problem including process dynamics and then propose a decomposition approach that results in the efficient solution of the integrated problem. Optimality Analysis is utilized to prove that the production sequence and transition times are independent of products’ demands. The proof leads to the decomposition of the integrated problem into two subproblems that can be solved separately without the need for iterations. Results of case studies verify the feasibility and effectiveness of the proposed approach in reducing the computational complexity of the integrated problem but also obtaining the optimal solution.

1. INTRODUCTION Cyclic production arises in manufacturing in which production stages are processed sequentially and the production stages are repeated within the production cycle (production wheel). Usually one production stage is dedicated to the production of one product and therefore multiple products can be produced in a predetermined sequence by sharing one unit. Time slots are usually used for scheduling a problem in cyclic production1−3 while in sequential production, processing stages are utilized.4−6 Note that sequential production is essentially a particular case of cyclic production since it can be considered as one cycle. A cyclic production optimization problem typically involves decisions at both scheduling and dynamic process optimization such as production sequence, production time, and manipulations during transitions between production stages. For continuous processes, units are operated at steady state during production stages and the optimal control strategies are applied to reduce off-spec products during transitions between production stages. The overall objective is to maximize profit per cycle time, or profit in a fixed time horizon. A traditional approach to handle production scheduling and dynamic optimization problems in chemical processes is to consider them separately. However, there are variables that link the scheduling and dynamic optimization problems, for example, the states at the beginning and end of the transient periods, duration of the transient periods and steady state periods, and amount of material in production stages. Scheduling and dynamic optimization are naturally connected by the linking variables and thus cannot be considered completely separately. Using an integrated approach, information can be shared between scheduling and the dynamic optimization level, leading to the integrated decision making which improves the profitability of process operations.7−10 Recently a number of publications appeared in the literature to address the integration of scheduling and dynamic optimization problems for continuous processes. Two approaches can be identified in this area: simultaneous modeling and decoupled modeling. With simultaneous © XXXX American Chemical Society

modeling, the process dynamics is treated as constraints and incorporated into the scheduling problems, forming an integrated problem which is a mixed integer dynamic optimization (MIDO) problem.11 The MIDO can be discretized using collocation point methods, and a mixed integer nonlinear programming (MINLP) is obtained.1,12,13 With the use of decoupled modeling, the scheduling problem (master problem) is modeled as mixed integer linear programming for which process dynamics are considered separately in a dynamic optimization. This approach is implemented by iterating between the master problem and the primal problem until convergence is achieved.4,5 In our previous work we proposed a closed loop implementation of the integrated approach which considers process disturbances.14 We also developed an integration using multiparametric model predictive control to handle a disturbance efficiently.15 Chu and You16 developed a fast computation strategy for real time implementation of integrated scheduling and control. Nie et al.17 employed a state equipment network (SEN) to model batch process scheduling and built mixed-logic dynamic optimization (MLDO) for the integration of scheduling and dynamic optimization. They verified the advantages of an integrated approach compared with a conventional scheduling method through a number of case studies. Integration of scheduling and dynamic optimization results in a large size MINLP which is computationally very expensive. To achieve better computational performance, TerrazasMoreno et al.18 proposed Lagrangean decomposition for the integrated model and solved it with iteration between a primal control problem and a master scheduling problem. Chu and You19,20 explored the structure of the integrated problem and applied Benders Decomposition to decrease the computational complexity. Their approach effectively lowered the computaReceived: May 26, 2015 Revised: December 23, 2015 Accepted: January 15, 2016

A

DOI: 10.1021/acs.iecr.5b01916 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research tion time and achieved optimality of the integrated problem. They also model the integrated problem using game theory and efficiently solve the resulting bilevel optimization problem.21 Recently Nie et al. developed an integrated model based on resource task network (RTN) representation.22 They applied a tailored generalized Benders decomposition (GBD) algorithm to efficiently solve the resulting large nonconvex MINLP. Chu and You23 investigated the integrated planning, scheduling, and dynamic optimization problem and used surrogate models to represent the links between subproblems. The decomposition approach built based on the surrogate modeling significantly reduced the computation time. Du et al.24 proposed an explicit and low order model for the closed loop dynamics. The solution approach of the resulting integrated problem demonstrated computational advantages and economic performance comparable to the existing approach that uses perfect process models. In our study, we explore the structure of the integrated optimization problem for continuous processes and establish an efficient decomposition scheme based on the mathematical structure of the corresponding model. We decompose the scheduling problem into an optimization model for the production time and a separate one for production sequence. In this way, two subproblems are formed. One is an MINLP dealing with the production sequence and dynamic profile during transitions, and the other is an NLP dealing with the production time for each product considering the market demands. We prove the separability of the resulting two subproblems analytically. Unlike the Lagrangean decomposition approach, no iterations are needed in the proposed approach. Moreover, the computational complexity is significantly reduced compared to the simultaneous scheduling and dynamic optimization problem under varying demand. The main contribution of this work is that there is no need to solve the entire integrated problem for different demand specifications because the only decision variables that need to be updated are the production times. The rest of the paper is organized as follows. The description of the general integrated model for scheduling and dynamic optimization for continuous processes are presented in section 2, followed by the proposed decomposition methods and analytical proofs in section 3. On that basis, two case studies are presented in section 4 to show the feasibility and advantages of the proposed approach, and finally conclusions are presented in section 5. To keep the flow of the paper and make it easier for the reader to understand, some steps of the analytical proofs in section 3 are shown in the Appendices.

Figure 1. Cyclic production of continuous processes.

which are obtained from the solution of problem (eq 1) in which state variables xk govern the system dynamics at slot k. max Jk (xk , uk , qk ) uk

⎧ xk̇ = f (xk , uk) ⎪ ⎪ s.t.⎨ qk = h(xk , uk) ⎪ ⎪(xk , uk , q ) ∈ Ωxuq ⎩ k

(1)

In problem 1 Jk is the local profit, and uk and qk are the manipulated variables and outputs, respectively, whereas Ωxuq is the domain for xk, uk, qk. With simultaneous modeling, the process dynamics (differential equations derived from mass and energy balance) are considered as constraints and incorporated into the scheduling problem, resulting in the integrated model in formula 2. More specifically, ẋk = f(xk, uk) are the differential equations describing the process dynamics which can be solved under the initial value of state and manipulated variables xk0 and uk0, leading to xk = f f( f, uk, xk0, uk0). Here f f is the solution for xk. It is determined by the form of f and the manipulated variable uk and the initial conditions xk0 and uk0. f f is not necessarily in an explicit form. Note that xk0 and uk0 are equal to the steady state values of the prior time slot, which are relevant to the scheduling solution y. In continuous processes because the transition between products is followed by the steady state period, the sequence of products determines the set point of the transition in a certain slot. Thus, we obtain xk = f f( f, uk, xk0(y), uk0(y)) which is integrated with the scheduling model as shown in problem 2. max J(w , y , x , u , q) w,y,u

⎧ g (w , y) ≥ demand(i) ∀ i ∈ I ⎪ i ⎪ g (w , y) ≤ g ∀i∈I max ⎪ i ⎪ s.t.⎨ xk = ff (f , uk , xk 0(y), uk 0(y)) ∀ k ∈ K ⎪ ⎪ qk = h(xk , uk) ∀ k ∈ K ⎪ ⎪(xk , q , uk , w) ∈ Ωxquw ∀ k ∈ K ⎩ k

2. GENERAL MODEL In this section we present the general scheduling and dynamic optimization models for continuous cyclic productions, and explain how the scheduling and dynamic optimization problems are integrated. 2.1. Simultaneous Scheduling and Dynamic Optimization for Continuous Processes Cyclic Production. In continuous processes cyclic production, a slot-based formulation is adopted1 and in each the slot transition period (dynamic profile) is followed by a production period (steady state profile) (Figure 1). The optimization problem at the scheduling level determines the production time and production sequence in order to satisfy demands. During transitions between steady state production stages, optimal dynamic profiles are needed

(2)

where I is the product set and K is the slot set; y is a vector of binary variables representing assignment of products in slots; w is the vector of production times which are decision variables in scheduling problem; and Ωxquw is the domain for w, the state variables x, the manipulated variables u, and the output variables q defined in problem 1. In this study infinite amounts of raw material and resources are considered. In the scheduling problem demands should be satisfied by the production B


Article

Industrial & Engineering Chemistry Research

Figure 2. Decomposition for the solutions of integrated problem of continuous processes.

amount gi. gi is confined by an upper bound in order to avoid significant overproduction. The objective J(w, y, x, u, q) represents the profit in unit time or in general the total profit in a predetermined time period. 2.2. Detailed Model of Integration. The integrated model (2) follows the scheduling formulation introduced by Flores−Tlacuahuac and Grossmann.1 Implicit Runge−Kutta is employed for the discretization of the dynamic model.14 To optimize process operations, we utilize the objective of profit maximization per unit time, which can be computed in eqs 3−5.

Φ = Φ1 − Φ2 Np

Φ1 =

∑ i=1

products at the end of the cycle or delivering the product once it is produced in the cycle. Detailed constraints of the integrated model are presented in Appendix B. By combining the objective functions described in eqs 3−5 and the constraints at both scheduling level that are shown in Appendix B (B.1−B.18) and dynamic optimization level (B.19−B.25) along with the linking constraints (B.26−B.34) we obtain the optimization problem for the integrated model as shown in problem (6). max

yik , zipk , xke , uke , θkt , Θi, Wi , tks , tke , Tc

⎧(B.1)−(B.18) constraints at scheduling level ⎪ ⎪(B.19)−(B.25) constraints of dynamic s.t.⎨ ⎪ optimization at control level ⎪ ⎩(B.26)−(B.34) linking constraints

(3)

CipWi Tc

(4)

Φ2 = (COSTtr + COSTss)/Tc Ns

= (∑ ∑ C k=1 e=1 Ns C r(uk̅ 1 k=1

∑

(uke1

+ ··· +

ukem)hk θkt +

+ ··· + uk̅ m)pk )/Tc

(6)

Note that in the integrated problem 6 the production sequence and transitions are connected. For example, if the sequence is BDEAC, then the transitions from B to D, D to E, E to A, A to C, and C to A have to be determined. Thus, the production sequence and transitions are both part of the solution of the integrated solution (i.e., the overall solution) and they are obtained simultaneously once the integrated problem is solved.

Ne r

Φ1 − Φ2

(5)

where Φ1, Φ2 represent the rate of revenue and raw material cost, respectively. Φ2 is composed of the cost during transition periods COSTtr and cost during steady state production periods COSTss. The total revenue is calculated as the amount produced (Wi) times the product price (Cpi ). Raw material cost is composed of two parts, the raw material consumption during production periods and during the transition periods. Unlike the objective in Flores−Tlacuahuac and Grossmann,1 we do not consider state variation in the objective, since we address a pure economic objective in this study. Another difference with Flores−Tlacuahuac and Grossmann1 is that we do not include inventory cost which can be considered at the planning level. We believe that incorporating inventory cost could complicate the decomposition addressed in this study, not only because the inventory cost could greatly add complexity to the objective but also because of the depleting strategy: delivering all the

3. DECOMPOSITION APPROACH For the purpose of optimality analysis, we introduce an alternative form of objective as shown in Equation 7 which is equivalent to the one in 3−5. J=

∑i Pi Θi − f (yik , θkt ) ∑i Θi + ∑k θkt

(7)

θtk

where Θi is production time for product i, and is the transition time at slot k. Note that eqs 4 and 5 share the same denominator Tc which is equal to the sum of production time and transition time, that is, ∑iΘi + ∑kθtk. Combining eqs 4 and 5 leads to eq 7. ∑iPiΘi represents the profit gained from production and f(yik, θtk) is the total transition cost. Note that Pi C


Article

Industrial & Engineering Chemistry Research = (Cpi (production rate) − Cr (material consumption rate)). Since production rate and material consumption rate are known in a steady state, we group them together to make the expression more compact. Pi can be considered as a revenue coefficient. As shown in Figure 2, the original integrated problem can be decomposed into two subproblems. The first problem involves the transition stages including sequences and transition profiles, and the second one consists of the production stages. The corresponding optimization problems for these two parts are as shown in Figure 3.

where (z*ipk, y*ik , θtk*, x*ke, u*ke) ⊂ {solution set of eq 8} Both subproblems correspond to much smaller problems compared to the original integrated problem. The following proof and case study results in section 4.1 verify the proposed decomposition. It should be noticed that the proposed decomposition approach does not lead to the solution of the scheduling and the dynamic optimization problem separately. Theorem 1 provides a proof that the solution of eq 8 and eq 9 is the same as the one from the integrated problem (6). Theorem 1. Assuming that inventory costs are zero, the decomposition in Figure 3 does not affect the optimality of the integrated problem; that is, solving the subproblems 8 and 9 results in the same solution as the integrated problem (6). Three lemmas are proposed to prove Theorem 1. Lemma 1. The positive perturbation on production time does not affect the transition for any specific time slot. Proof: The profit for each time slot is J=

In Figure 3 subproblem 1 deals with the sequence and transition times, and subproblem 2 deals with production. By assigning a fixed value to the production time in the original integrated problem 6, we obtain subproblem 1 (eq 8). max

P Θ − ftc (θ t *) P Θ − ftc (θ t ) ≥ Θ + θt Θ + θ t*

(8)

Solving (eq 8) generates the optimal production sequence and optimal transition time. Those are then assigned in subproblem 2 (eq 9), which generates all the remaining decision variables in the integrated problem. max

Θi, Wi , tks , tke , Tc

(11)

The proposition that needs to be proved can be described as follows.

Φ1 − Φ2

⎧(B.1)−(B.7) constraints at scheduling level ⎪ ⎪(B.11)−(B.18) constraints at scheduling level ⎪ ⎪(B.19)−(B.25) constraints of process dynamics s.t.⎨ ⎪(B.26)−(B.34) linking constraints ⎪ ⎪ Θ = Θ0 assign production time for each ⎪ i product ⎩

(10)

where P is revenue coefficient introduced in eq 7. Θ, θt, and f tc (θt) correspond to production time, transition time, and transition cost, respectively. Assume that θt* is the optimal transition time when the production time is Θ. Then we have

Figure 3. Decomposition for the optimization problem of the integration of scheduling with dynamic optimization for continuous processes.

yik , zipk , xke , uke , θkt , tks , tke , Tc

ftc (θ t ) P Θ − ftc (θ t ) PΘ − = Θ + θt Θ + θt Θ + θt

Figure 4. Perturbation on production does not affect transition within one slot.

Φ1 − Φ2

⎧(B.8)−(B.18) constraints at scheduling level ⎪ ⎪ zipk = zipk * assign production sequence ⎪ ⎪ y = y* assign production sequence ⎪ ik ik s.t.⎨ t ⎪ θk = θkt * assign transition time for each slot ⎪ ⎪ x = x * assign state profile ke ⎪ ke ⎪ u = u * assign control profile ⎩ ke ke

If there is a small perturbation on Θ shown in Figure 4, θt* is still the optimal transition time which can be expressed mathematically as follows: P(Θ + δ Θ) − ftc (θ t ∗) (Θ + δ Θ) + θ t ∗

≥

P(Θ + δ Θ) − ftc (θ t ) (Θ + δ Θ) + θ t

(12)

The specific procedures to establish this proposition are provided in Appendix A. Lemma 2. A property of θtk*:θtk* is obtained as the minimum value of feasible θtk obtained by solving problem 13.

(9) D


Article

Industrial & Engineering Chemistry Research θkt * =

min

uk(t ), t ∈ [0, θkt ]

θkt

⎧ x ̇ (t ) = f (x (t ), u (t )), t ∈ [0, θ t ] k k k ⎪ k ⎪ xin , k = xk(0) = xk − 1 ̅ ⎪ ⎪ t s.t.⎨ xk(θk ) = xk̅ ⎪ ⎪ xk(t ) ∈ Ωx , t ∈ [0, θkt ] ⎪ ⎪ u (t ) ∈ Ω , t ∈ [0, θ t ] ⎩ k u k

Figure 5. A perturbation on production does not affect all the transitions. (13)

Proof: The overall objective of the whole production cycle is J=

On the basis of Lemmas 1, 2, and 3 we prove that the production does not affect the optimality of transitions and production sequence. Therefore, Theorem 1 holds and the decomposition in Figure 3 is verified. Note that Lemma 1 and 3 are essential steps in the proof of the decomposition. In Lemma 1 we prove that the production does not affect the optimality of transition, and in Lemma 3 we prove the production does not affect sequence. Thus, we can separate production variables from the integrated problem and achieve decomposition as shown in Figure 2 and Figure 3. In principle the transition time can be precomputed through problem 13. However, it involves significant computation since we need to solve problem 13 repeatedly, as described in the paragraphs following problem 15. To avoid the huge computation burden, we propose a solution procedure that involves two steps for the integrated problem 6. Note that the proposed solution procedure is supported by Theorem 1 and Lemmas 1−3 in the proof. On the basis of the above analysis for decomposition, we propose the following solution procedure for the integrated problem 6: (1) Treat Θi as constants and assign extremely small positive values for them. Solve subproblem (8) and obtain (z*ipk, yik*, θtk*, xke *, uke *). This guarantees that the optimal production time Θi is greater than the assigned value. The difference between the optimal production time and the assign value could be viewed as a positive perturbation. Note that the obtained optimal transitions and sequence in this step are equal to the solutions of the original integrated problem since it is proved in Lemma 1 and Lemma 3 that the positive perturbation on production does not affect the optimality of transitions and sequence. (2) Substitute (z*ipk, y*ik , θtk*, x*ke, u*ke) into the original problem, which results in subproblem (9). Solve it and obtain Θ*i . Since we have proved that production times do not affect production sequence or transition profiles, this solution approach will lead to the same solution as the integrated optimization problem.

∑i Pi Θi − ftc (yik , θkt ) ∑i Θi + ∑k θkt

(14)

where ftc (yik , θkt )

=

∑∑C i

=

k

r

∫0

∑∑∑C i

p

k

r

θkt

ukyik dt

∫0

θipt

ukzipk dt (15)

r

and C is a coefficient for transition cost, for instance, the unit price for raw materials. The transition cost f tc(yik, θtk) can be evaluated as the integral of the manipulated variable over the transition duration θtk in all time slots. Binary variables yik and zipk indicate assignments of production in slots ((B1−B7) in Appendix B). Taking the partial derivative of J with respect to θtk when Θi, θtK−k, yik(∀i ∈ I, ∀k ∈ K) are fixed, leads to the following property for θtk; that is, θtk* ≤ θtk. The details of the proof are provided in Appendix A. With this property, we can calculate θtip* before hand by solving optimization problem 13 and calculate the corresponding transition cost from product i to product p at slot k. However, it involves significant computation since we need to solve problem 13 many times. For a cyclic production where n products are produced, the number of transitions would be p(n, 2) = n!/(n − 2)!. Note that the transition from product A to B is not necessarily equivalent to the transition from B to A. It means that we have to solve problem 13 for 20 times, so this will be computationally impossible for a large number of products. Note that the optimal transition time computed by problem 13 is independent of the product demands. It should be noted that which transitions exist in the scheduling solution is of course a function of the demand but the optimal transition times are not. For example, the transition from D to E can be found through problem 13 without knowing the demands of D or E. However, the transition from D to E may not exist in the overall solution if for example the optimal sequence is BCDAE. Lemma 3. The optimal production sequence is not affected by positive production perturbation. Proof: Suppose z* is the optimal sequence; that is, J(z*) ≥ J(z). What we need to prove is that this inequality holds in the presence of a perturbation in production. For example, the optimal sequence C−A−E−D−B in Figure 5 is retained when δΘ presents. Following the steps in Appendix A, we obtain a sufficient condition (A.36) for this hypothesis.

4. CASE STUDIES To demonstrate the feasibility of the decomposition approach, we studied two numerical cases, one for a single state continuous stirred-tank reactor (CSTR) and the other for methyl methacrylate (MMA) polymerization which involves multiple states. In each case study we investigate three scenarios, and discuss their results. The first scenario is to solve for the original integrated problem MINLP (problem 6) using both local and global optimization solvers. In scenario 2 we apply the proposed decomposition approach and solve the E


Article

Industrial & Engineering Chemistry Research subproblems 8 and 9 sequentially following the procedure presented in section 3. In scenario 3 we solve for the scheduling and dynamic optimization problem in a conventional approach that does not involve integration of them. Given the estimation of transition time θ̃tpi and transition cost COST trpi , scheduling problem 16 is solved under the constraints at scheduling level and linking constraints. Dynamic optimization problem 17 is then solved based on solution of the scheduling problem e.g. y*ik . The dynamic optimization problem minimizes the cost during transitions under the constraints at dynamic optimization level and linking constraints. Note that the estimation of transition time θ̃tpi is made empirically based on the knowledge of the process. A aggressive estimation unnecessarily extends the transition period and thus undermines the profitability, while a reckless one could make the dynamic optimization problem infeasible, that is, a transition too short makes it impossible to drive the states to the desired steady states within the estimated time period. The transition cost from product p to i COST trpi is calculated as the product of the average of initial and desired manipulated variables and the estimated transition time and the price, that is, COST trpi = Cr·0.5(uin,k + u̅k)·θ̃tpi. This is a simplification of the standard form of transition cost in eq 15. Both θ̃tpi and COST trpi can be obtained beforehand. They present in problem 16 as parameters. max

yik , zipk , θkt , Θi, Wi , tks , tke , Tc

Ns

xke , uke

Table 1. Steady States and Product Information of a Single State CSTR Q (L/hour)

CR (mol/L)

demand rate (kg/hour)

product price ($/kg)

A B C D E

10 100 400 1000 2500

0.0967 0.2 0.3032 0.393 0.5

3 8 10 10 10

200 150 130 125 120

(18)

where C0 is feed stream concentration (a parameter in the model), Q is the feed flow rate (manipulated variables), and CR is the concentration of raw material in the outflow (state variable). The detailed integrated problem including the overall objective and the constraints for both scheduling and dynamic optimization levels are presented in Appendix B. We solve the integrated problem for five production cycles in which demand for each product is varying. Table 2 shows the specific demands that we assign in five cycles and Figure 6 demonstrates the solutions for each cycle. It can be seen that the demands for all the five products are varying and correspondingly the productions are varying in order to satisfy the demands. Nevertheless, the transitions prior to the production stages are consistent but do not respond to the varying demands as shown in Figure 6. Obviously the optimal profit is varying with respect to cycles, due to the fluctuating demands. It is worth mentioning that all the cycles generate the same optimal sequence C−B−A−E−D, regardless of the varying demands. Table 2 demonstrates comparisons between the demand rate and the yielding rate that is calculated as W/ Tc. W and Tc are the production amount and cycle time in the solution of the integrated problem, respectively. Product A, B, C, and D are produced in the amount exactly the same as the demand rate while product E is produced much more than demand. This is due to the fact that product E has the highest value of G*Price, where G is the production rate. In other words, producing product E is the most profitable thus the process dedicates as much time resource as possible to produce E. It is also found in this case study that an upper bound on the production amount confines the production of product E but does not affect the production of the other products if the upper bound is greater than their demands. Table 4 summarizes the problem size, computation time, and the results of scheduling and dynamic optimization under all these three scenarios described above, that is, scenario 1 the original integration, scenario 2 the proposed approach, and scenario 3 the conventional approach. All scenarios correspond to product demands of cycle 1 in Table 2. Problems in scenario 3 are solved on the basis of the transition estimation provided by Table 3. Scenarios 1 and 2 result in the same scheduling and dynamic optimization solutions and the same profit, however scenario 2 (the proposed decomposition approach) takes much less time to compute it. Though the problem size of scenario 2 is only slightly smaller than the original integrated problem in scenario 1, the computational complexity for MINLP is significantly reduced. This is because decomposing the decision

(16)

Ne

∑ ∑ Cr(uke1 + ⋯ + ukem)hkθkt k=1 e=1

⎧(B.19) − (B.25) constraints of dynamic ⎪ ⎪ optimization at control level ⎪ ⎪ s.t.⎨(B.14), (B.26) − (B.34) linking constraints ⎪ t t ⎪ θpi = θpĩ ⎪ ⎪ y = y* ⎩ ik ik

product

dC R Q = (C0 − C R ) + rR dt V

Φ1 − Φ2

⎧(B.1) − (B.18) constraints at scheduling level ⎪ ⎪(B.26) − (B.31) linking constraints ⎪ ⎪ t t s.t.⎨ θpi = θpĩ ⎪ Np Np ⎪ tr ⎪COST = ∑ ∑ COSTtrpizipk , ∀ k ∈ K ⎪ i=1 p=1 ⎩

min

(Table 1) and are manufactured in a cyclic mode. Mass balance results in a kinetic model in eq 18.

(17)

The purposes of investigating three scenarios are to demonstrate the significance of integration in terms of boosting profitability (comparing scenario 1 and 3), and the advantage of the proposed decomposition approach in reducing computing time (comparing scenario 1 and 2). 4.1. Single State CSTR Cyclic Production. The data for the first case study involving a CSTR is taken from Flores− k

Tlacuahuac and Grossmann.1 Reaction 3R 3R → P occurs in an isothermal CSTR under reaction rate −rR = kCR3. Products A, B, C, D, and E are differentiated by their concentration CR F


Article

Industrial & Engineering Chemistry Research Table 2. Demands and Yields in a Single State CSTR Cyclic Productiona demands and productions (kg/h) product product product product product a

A B C D E

cycle 1

cycle 2

cycle 3

cycle 4

cycle 5

(3, 3) (8, 8) (10, 10) (10, 10) (10, 75.8)

(3.3, 3.3) (6, 6) (12, 12) (8, 8) (11, 60.7)

(2.5, 2.5) (11, 11) (7, 7) (14, 14) (9, 103.4)

(2, 2) (9, 9) (8, 8) (12, 12) (8, 203.4)

(1, 1) (10, 10) (8, 8) (11, 11) (11, 328.3)

The first number in parentheses corresponds to demand, the second number corresponds to yield.

Figure 6. Solutions for integrated problem at different cycles assigned with different demands.

Table 3. Estimation of Transition Time (hour) and Transition Cost ($) between Products (Single State CSTR)a (θ̃tpi, COST trpi ) product product product product product a

A B C D E

product A

product B

product C

product D

product E

(0, 0) (60, 33000) (62, 127100) (65, 328250) (69, 865950)

(0.4, 220) (0, 0) (4, 10000) (6, 33000) (8, 104000)

(0.5, 1030) (0.4, 10000) (0, 0) (2.2, 15400) (3.5, 50750)

(1.8, 9090) (0.8, 4400) (0.9, 6300) (0, 0) (1.6, 28000)

(2.1, 26360) (1.8, 23400) (1.5, 21750) (1.2, 21000) (0, 0)

For instance, transition time from product C to B is 4 h and the corresponding transition cost is $10000.

Table 4. Results of Single State CSTR Cyclic Production: Comparison between Three Scenarios

problem type and solver in GAMS problem size (constraints, variables) CPU (s) production sequence cycle time (h) production time in slots (h) transition time in slots (h) raw material cost ($/h) revenue ($/h) profit ($/h)

scenario 1:

scenario 2:

scenario 3:

original integration problem (6)

decomposition approach

conventional (no integration)

problem (6), MINLP, SBB 4046, 2986

problem (6), MINLP, BARON 4046, 2986


49.3 C−B−A−E−D 140 5.02, 14.00, 46.49, 8.49, 2.31 1.12, 3.74, 56.57, 1.61, 0.63 2353.67 13450.17 11096.49

527.9 C−B−A−E−D 140 5.02, 14.00, 46.49, 8.49, 2.31 1.12, 3.74, 56.57, 1.61, 0.63 2353.67 13450.17 11096.49

17.8

problem (9), NLP, CONOPT 65, 31

0.32 C−B−A−E−D 140 5.02, 14.00, 46.49, 8.49, 2.31

problem (16), MINLP, SBB 1011, 1116 5.2


2.6 C−A−B−D−E 140 5.02, 46.49, 14.0, 2.30, 4.27

1.12, 3.74, 56.57, 1.61, 0.63

3.5, 62.0, 0.4, 0.8, 1.2

2353.67 13450.17 11096.49

2218.30 8929.74 6711.44

G


Article

Industrial & Engineering Chemistry Research variables effectively changes the structure of the optimization problem without undermining the optimality of the integrated problem, as proved by the theorems in section 3. The comparison between scenario 1 and 3 demonstrates the significance of integration of scheduling and dynamic optimization. Scenario 3 follows the conventional approach where the integrated problem is not obtained. It results in much lower profit than the original integrated problem does. This is because the scheduling solutions obtained based on the estimation of transition time and transition costs are not optimal in terms of the overall profit in the integrated problem. Note that all scenarios generate the same optimal cycle time 140 h while solving the corresponding optimization problems under the constraints for the cycle time 90 ≤ Tc ≤ 140. In this case study we obtain a feasible solution using a conventional approach (scenario 3). However, a not good estimation of transition time could cause infeasibility at the control level when for instance, the transition time is too short for the state to change from one steady state to another. To find a good estimation of transition time and transition cost so that the infeasibility can be avoided, an iterative approach is considered. First, we make an initial estimation of the transition time and transition cost and solve for the scheduling problem based on these estimations. Then we implement the scheduling solution to control level by simulation and check whether it is feasible. If it is infeasible, we update the estimations and solve for corresponding scheduling again until we obtain feasible control solution. 4.2. MMA Production with Multiple States. The proposed decomposition approach is applied to an isothermal free radical polymerization process that produces methyl methacrylate (MMA) using azobis(isobutyronitrile) as initiator and toluene as the solvent. The reaction takes place in an isothermal CSTR at the temperature of 335 K. The kinetic model is describe by Mahadevan et al.25 and summarized in the following equations. Equations 19−22 present the dynamics of four state variables monomer concentration Cm, initiator concentration CI, molar concentration of dead chains D0, and mass concentration of dead chains D1. Equation 23 presents an auxiliary variable P0 and eq 24) calculates the output variable the molecular weight D1/D0. The manipulated variable is the initiator flow rate FI. The associated parameters and state/ manipulated variables are provided in Table 5 and Table 6, respectively. The process is illustrated in Figure 7.

Table 6. Variables in the Kinetic Model of MMA Polymerization in a CSTR Cm (kmol/m3) CI (kmol/m3) D0 (kmol/m3) D1 (kg/m3) FI (m3/h) y = D1/D0 (kg/kmol)

kTc = 1.33 × 1010 m3/(kmol h) CIin = 8.00 kmol/m3

temperature of the reactor monomer flow rate reactor volume initiator efficiency propagation rate constant termination by disproportionation rate constant termination by coupling rate constant inlet initiator concentration

Cmin = 6.00 kmol/m3

inlet monomer concentration

kfm = 2.45 × 103 m3/(kmol h) kI = 1.02 × 10−1 h−1 Mm = 100.12 kg/kmol

chain transfer to monomer rate constant initiation rate constant molecular weight of monomer

F(Cmin − Cm) dCm = −(kp + k f m)P0Cm + dt V

(19)

FC dC I I I in − FC I = −kIC I + dt V

(20)

dD0 FD0 = (0.5k T c + k T d)P02 + k f mP0Cm − dt V

(21)

dD1 FD1 = M m(kp + k f m)P0Cm − dt V

(22)

⎛ 2f *kIC I ⎞0.5 P0 = ⎜ ⎟ ⎝ kT d + kT c ⎠

(23)

y=

D1 D0

(24)

Figure 7. MMA polymerization in a CSTR.

In this case study, five products are produced cyclically at steady state and each product grade corresponds to a value of molecular weight. The steady-state values of the state variables, the manipulated variables, and the output variables are provided in Table 7. Similarly to the first case study we again investigate three scenarios and compare their problem size, CPU time, scheduling, and dynamic optimization solution as well as the overall profit as provided in Table 9. Scenario 1 is solved using normal MINLP solver SBB and global optimal solver BARON. They result in similar overall profit. Scenario 2 follows the decomposition approach presented in section 3.1 and the solution procedure described in section 3.2. Scenario 3 is solved based on the estimated transition time and transition cost in Table 8. It can be observed that scenario 2 (the proposed decomposition approach) produces slightly less profit, 562.43, than the global optima, 564.11, but consumes much less CPU time. However,

Table 5. Parameters in the Kinetic Model of MMA Polymerization in a CSTR T = 335 K F = 10.0 m3/h V = 1.0 m3 f * = 0.58 kp = 2.50 × 106 m3/(kmol h) kTd = 1.09 × 1011 m3/(kmol h)

state: monomer concentration state: initiator concentration state: molar concentration of dead chains state: mass concentration of dead chains manipulate: initiator flow rate output: molecular weight

H


Article

Industrial & Engineering Chemistry Research Table 7. Steady State Information for Each Product (MMA Case Study) product

Cm (kmol/m3)

CI (kmol/m3)

D0 (kmol/m3)

D1 (kg/m3)

FI (m3/h)

D1/D0 (kg/kmol)

demand (m3/h)

price ($/m3)

A B C D E

5.73 5.63 5.57 5.46 5.41

0.0366 0.0713 0.0984 0.1615 0.1963

0.0007 0.0012 0.0015 0.0023 0.0028

27.0324 37.0444 43.0516 54.0648 59.0708

0.0463 0.0900 0.1242 0.2039 0.2479

40084 31938 28293 23153 21294

0.4 0.7 1.1 0.9 0.8

170 150 130 125 110

Table 8. Estimation of Transition Time (hour) and Transition Cost ($) between Products (MMA Case Study)a (θ̃tpi, COST trpi ) product product product product product a

A B C D E

product A

product B

product C

product D

product E

(0, 0) (0.12, 40.89) (0.18, 76.72) (0.28, 175.14) (0.29, 213.29)

(0.14, 47.70) (0, 0) (0.16, 85.68) (0.25, 183.68) (0.27, 228.08)

(0.16, 68.20) (0.14, 74.97) (0, 0) (0.17, 139.44) (0.22, 204.65)

(0.23, 143.86) (0.19, 139.60) (0.15, 123.03) (0, 0) (0.15, 169.42)

(0.28, 205.94) (0.22, 185.84) (0.19, 176.74) (0.16, 180.72) (0, 0)

For instance, transition time from product C to B is 0.16 h and the corresponding transition cost is $85.68.

Table 9. Results of MMA Case Study, Comparison between Three Scenarios

problem type and solver in GAMS problem size (constraints, variables) CPU (s) production sequence cycle time (h) production time in slots (h) transition time in slots (h) raw material cost ($/h) revenue ($/h) profit ($/h)

scenario 1:

scenario 2:

scenario 3:

original integration problem (6)

decomposition approach

sequential (no integration)


problem (6), MINLP, BARON 5430, 4035


275.4 C−B−E−D−A 70 7.7, 9.7, 5.6, 6.3, 40.0 0.10, 0.10, 0.12, 0.07, 0.27 960.97 1523.41 562.43

1154.9 C−E−D−A−B 70 7.7, 5.6, 6.3, 40.0, 9.80

129.1


0.43 C−B−E−D−A 70 7.7, 9.7, 5.6, 6.3, 40.0

0.06, 0.10, 0.07, 0.27, 0.08 960.85 1524.96 564.11

problem (16), MINLP, SBB 592, 390 16.9


24.4 C−A−D−E−B 70 7.7, 40.0, 6.3, 5.6, 9.4

0.10, 0.10, 0.12, 0.07, 0.27

0.14, 0.18, 0.23, 0.16, 0.27

960.97 1523.41 562.43

974.22 1516.78 542.56

In this study the integrated problem is solved offline and implemented in an open-loop manner. What we present in this paper is an solution approach that reduces the computation complexity of the integrated problem. It should be emphasized that the decomposition approach proposed in this study is not simply decomposing the integrated problem into the scheduling and dynamic optimization problems, though these two problems are integrated as one in problem 6. In subproblem 1(8) the sequence variables are associated with the scheduling level while the transition variables are associated with the dynamic optimization level. In other words, solving the integrated problem in a decomposition scheme does not change the nature of integration. The problem that is addressed is indeed the integrated scheduling and dynamic optimization problem, and the decomposition is proposed as a way to improve the computational efficiency of the solution approach. The fact that the transition periods are proven to be independent of the production in continuous cyclic production decreases the computation complexity of the integrated scheduling and dynamic optimization problem. This facilitates the integration with a planning level, since the only linking variables between planning and integrated scheduling and dynamic optimization are the production variables.

scenario 3 produces the least profit, 542.56, owing to the lack of integration.

5. CONCLUSIONS In this work, we propose a decomposition approach for the integrated scheduling and dynamic optimization for continuous cyclic production. Although the main focus is in reaction processes since this is the example used for demonstration, the approach is general to handle other processes that have the same basic behavior. The decomposition approach of a continuous process was applied to the CSTR cyclic production with general kinetic models (B.19). Note that the decomposition for the continuous case is feasible when a certain condition (A.36) is satisfied and inventory cost is not included in the objective function (eq 3). Through optimality analysis we prove that for continuous processes the production sequence is retained when demands are varying. This property is of particular importance in real applications where the demands are varying for different production cycles, since there is no need to update production sequence for the future cycles and the production can be updated by solving the subproblem regarding production periods only. Results of case studies verify this property and also show that the decomposition approaches generate the same results as the integrated approach but using less computation time. I


Article

Industrial & Engineering Chemistry Research

■

Note that u is a function of time, and u(t) is valid when t ∈ [0, θtk]. u does not depend on θtk, instead, u can be any value in the range of [umin, umax]. u(θtk) represents the value of u at t = θtk. It does not correspond to a function of θtk. It is assumed that ∑iPiΘi − f tc(yik, θtk) ≥ 0, that is, the process is profitable. Thus, the numerator of eq A.8 is negative, and we have

APPENDIX A (SUPPORTING MATERIAL FOR THE DECOMPOSITION PROOF IN SECTION 3.1)

Lemma 1

The perturbation on production time does not affect the transition for any specific time slot. Continued proof: Rearranging inequality 11, we obtain

∂J ∂θkt

P Θθ t − Θftc (θ t *) − θ tftc (θ t *) + Θftc (θ t ) − P Θθ t * + θ t *ftc (θ t ) ≥ 0

(A.1)

≤0 Θi, θ −t k , yik

(A.10)

θtk*

which implies that should be the minimum value of the feasible region of θtk; that is, θtk* ≤ θtk. Figure 8 illustrates a transition from slot k − 1to slot k, where two dynamic profiles (xk*1(t), uk*2(t)), t ∈ [0, θtk*] and (xk*2(t),

Rearranging inequality 12 we obtain [P Θθ t − Θftc (θ t *) − θ tftc (θ t *) + Θftc (θ t ) − P Θθ t * + θ t *ftc (θ t )] + [Pδ Θθ t − δ Θftc (θ t *) − Pδ Θθ t * + δ Θftc (θ t )] ≥ 0

(A.2)

Comparing eqs A.1 and A.2, we can establish a sufficient condition so that eq A.2 holds. Pθ t − ftc (θ t *) − Pθ t * + ftc (θ t ) ≥ 0

(A.3)

which can be rearranged as follows ftc (θ t ) − ftc (θ t *) θ t − θ t∗

≥ −P

(A.4) t θmax ]

Assume f tc(•) is continuous on interval [0, and differentiable on interval (0, θtmax). Applying Lagrange’s Mean Value Theorem, there exists a point θt′ in (θt, θt*) if θt < θt* or (θt*, θt) if θt* < θt such that ftc ′(θ t ′) =

ftc (θ t ) − ftc (θ t ∗) θ t − θ t∗

(A.5)

Figure 8. Two dynamic profiles correspond to one optimal transition duration.

Based on eq A.5 we obtain a sufficient condition to establish eq A.4. min

f ′(θ t ) ≥ −P

t θ t ∈ (0, θmax ) tc

(A.6)

uk*2(t)), t ∈ [0, θtk*] correspond to one optimal transition

t

Since f tc(θt) = Cr∫ θ0 u dt (This is a simplified version of eq 15 addressing only one transition), we have

duration θkt*. If a one-to-one mapping between optimal

t

dftc (θ ) dθ t

transition time and optimal transition profile can be proven

= C ru(θ t ) > 0

(A.7)

(i.e. xk*1(t), xk*2(t) and uk*1(t) = uk*2(t), ∀t ∈ [0, θtk*]), we only

Apparently it verifies eq A.6, and as a result, eqs A.2−A.5 all hold. Thus, the proposal in eq 12 is proven. Q.E.D.

need to investigate the optimal transition time rather than the

Lemma 2

whole transition profile, which effectively simplifies the proof.

A property of θtk*:θtk* is obtained as the minimum value of feasible θtk. Continued proof: − ∂J ∂θkt

∂ftc (yik , θkt ) ∂θkt

In the following, through eqs A.11−A.27 we prove the one-toone mapping between optimal transition time and optimal

(∑i Θi + ∑k θkt ) − (∑i Pi Θi − ftc (yik , θkt ))

transition profile; that is, the case shown in Figure 8 does not

yik , θ t

−k

=

exist. We need to prove that

(∑i Θi + ∑k θkt )2

Θi, θ−t k , yik

(A.8)

As has been proven before ∂ftc (yik , θkt ) ∂θkt

θkt * ⇔ (uk*(t ), xk*(t )), t ∈ [0, θkt *]

= C ru(θkt ) ≥ 0 yik , θ t

−k

(A.11)

A general formula for the optimal control problem in formula (A.9)

13 is J


Article

Industrial & Engineering Chemistry Research max

u(t ), t ∈ [0, τ ]

J = L(x(τ ), u(τ )) +

∫0

⎧ ⎫ ∂V (x , t ) ∂V (x , t ) f (x , u)⎬ = 0 + min⎨Φ(x , u) + ⎭ u ⎩ ∂x ∂t

τ

Φ(x(t ), u(t )) dt

⎧ x(̇ t ) = f (x(t ), u(t )), t ∈ [0, τ ] ⎪ ⎪ x(0) = x 0 ⎪ ⎪ s.t.⎨ x(τ ) = x1 ⎪ ⎪ u(t ) ∈ Ωu , t ∈ [0, τ ] ⎪ ⎪ x(t ) ∈ Ω , t ∈ [0, τ ] ⎩ x

(A.21)

where V (x , t ) =

H[x(t ), u(t ), λ(t ), t ] (A.13)

Assume u*(t), t ∈ [0, τ*] is the optimal solution for eq A.12) and x*(t) is the corresponding trajectory. According to the Pontryagin Maximum Principle,26 the necessary conditions for u*(t) to be the optimum are eqs A.14−A.16: ∂ H(x∗(t ), u∗(t ), λ ∗(t )) ∂λ

λ ∗̇ (t ) = −

∂ H(x∗(t ), u∗(t ), λ ∗(t )) ∂x

⎧ ⎫ ∂V (x , t ) ∂V (x , t ) f (x , u)⎬ = 0, + min⎨1 + ⎭ u ⎩ ∂x ∂t

(A.14)

∀ t ∈ [0, θkt ]

(A.15)

V (x , t ) = (A.16)

∂f ∂x

∂f ∂H = λ*(t ) =0 ∂u ∂u

(A.17)

(A.18)

min

u(t ′), t ′∈ [t , θkt ]

∫t

θkt

dt ′

(A.24) 1

Propose that there are two optimal control profiles u* (t′) and u*2(t′) satisfying

(A.19)

u∗1(t ′) ≠ u∗2(t ′), ∃ t ′ ∈ [t , θkt ]

1 + λ*(t )f (x*(t ), u*(t ), t ) ≥ 1 + λ*(t )f (x*(t ), u(t ), t ) ∀ u(t ) ∈ Ωu , 0 ≤ t ≤ θkt *

(A.23)

⎧ x(̇ t ′) = f (x(t ′), u(t ′)), t ′ ∈ [t , θ t ]x ̇ k ⎪ ⎪ x(t ) = x ⎪ ⎪ s.t.⎨ x(θkt ) = x ̅ ⎪ ⎪ u(t ′) ∈ Ωu , t ′ ∈ [t , θkt ] ⎪ ⎪ x(t ′) ∈ Ω , t ′ ∈ [t , θ t ] ⎩ x k

and the necessary conditions A.18−A.20. λ*̇ (t ) = −λ*(t )

Φ(x(t ′), u(t ′), t ′)

(A.22)

Applying eqs A.13−A.16 to problem 13 we obtain Hamiltonian H[x(t ), u(t ), λ(t ), t ] = 1 + λ(t )f (x(t ), u(t ), t )

τ

Note that in formula A.22 and A.24 we use t′ in order to differentiate it from t. Both t′ and t represent time. In the case of problem 13 we have HJB eqs A.23 and A.24.

H[x∗(t ), u∗(t ), λ ∗(t )] ≥ H[x∗(t ), u(t ), λ ∗(t )], ∀ u(t ) ∈ Ωu , 0 ≤ t ≤ τ ∗

∫t

⎧ x(̇ t ′) = f (x(t ′), u(t ′)), t ′ ∈ [t , τ ] ⎪ ⎪ x(t ) = x ⎪ s.t.⎨ x(τ ) = x1 ⎪ ⎪ u(t ′) ∈ Ωu , t ′ ∈ [t , τ ] ⎪ ⎩ x(t ′) ∈ Ωx , t ′ ∈ [t , τ ]

Define Hamiltonian H as

x∗̇ (t ) =

u(t ′), t ′∈ [t , τ ]

L(x(τ ), τ ) +

dt ′

(A.12)

= Φ(x(t ), u(t ), t ) + λ T(t )f (x(t ), u(t ), t )

min

θtk]

(A.25)

If there exist t′ ∈ [t, such that u* (t′) ≠ u* (t′), we say that the profiles of u*1(t′) and u*2(t′) are different over region [t, θtk]. Plugging u*1(t′) and u*2(t′) into eq A.23 yields

(A.20)

Here we assume λ*(t) ≠ 0, otherwise eqs A.18−A.20 are always satisfied and u*(t) is free in Ωu, which is abnormal. If u* could attain any value in the feasible region, this becomes a trivial solution of eqs A.18−A.20, which is not complete nor correct in the real case, that is, u must be governed and be able to drive the state to the following production stage. This point can also be verified by checking the Hamilton−Jacobi−Bellman equation in the following. In addition, we assume that ∂f/∂u ≠ 0, that is, the dynamic is affected by control input. Therefore, condition eq A.19 is not satisfied, which indicates that there does not exist a stationary trajectory for u; that is u*(t) could locate at the boundary or non-differentiable point. Since umin ≤ u(t) ≤ umax, 0 ≤ t ≤ θtk we infer that the optimal u must be obtained at the boundary, u(t) = umin, 0 ≤ t ≤ θtk or u(t) = umax, 0 ≤ t ≤ θtk. Here we prove the one-to-one mapping (A.11) through checking the Hamilton−Jacobi−Bellman equations27 (A.21 and A.22).

1

2

⎫ ∂V (x , t ) ⎧ ∂V (x , t ) f (x , u*1)⎬ = 0, + ⎨1 + ⎩ ⎭ ∂t ∂x ∀ t ∈ [0, θkt ]

(A.26)

⎫ ∂V (x , t ) ∂V (x , t ) ⎧ + ⎨1 + f (x , u*2)⎬ = 0, ⎩ ⎭ ∂x ∂t ∀ t ∈ [0, θkt ]

(A.27)

Note that both x and u are trajectories with respect to time. f is a compound function and therefore f is also a trajectory with respect to time. If there exists t ∈ [0, θtk] such that u*1(t) ≠ u*2(t), the resulting f trajectories would be different, given the same state trajectories. In other words, there exists t ∈ [0, θtk] such that f(x(t), u*1(t)) ≠ f(x(t), u*2(t)). Note that when we K


Article

Industrial & Engineering Chemistry Research compare two functions f(x(t), u*1(t)) and f(x(t), u*2(t)) we are comparing their profiles over the whole region of t. Even though their profiles could overlap at some points, the profiles are considered different over the whole region. Therefore, the only way to make both eqs A.26 and A.27 hold is that ∂V(x, t)/ ∂x = 0; that is, the state at t does not affect V(x, t). However, this is not true. Suppose x(t) = x,̅ then obviously V(x, t) = 0 given by A.24. But if x(t) ≠ x,̅ we have V(x, t) > 0 according to A.24. This fact contradicts ∂V(x, t)/∂x = 0. Thus, the proposition in inequality A.25 is not established and there is unique optimal control solution u*. Therefore, θtk* ⇔ (uk*, xk*) is proven. Through eqs A.21−A.27 we introduce a contradiction to prove that eq A.25 does not hold. In other words, the optimal control profile is unique. On the basis of the above derivation, we infer that optimal control input for problem 13 is unique and the corresponding optimal state path is unique as well. In other words, there is a one-to-one mapping between optimal transition time and optimal dynamic profile. Thus, it’s reasonable to investigate the optimal transition time rather than dynamic profile, for the sake of simplification. Q.E.D.

Tc =

i

∑i Pi Θi −

ftc =

ftc ≥

p

k

*) ftc* = ftc (θipt* , zipk

zipk),

ftc ∗Tc − ∑i Pi Θi(Tc − Tc*) Tc*

∀ z ∈ Ωz

ftc ≥ ftc * − Pi(Tc − Tc*)

∀ z ∈ Ωz

(A.34)

Now the problem is transformed into the following: given eq A.33, prove eq A.34. Obviously a sufficient condition to establish eq A.34 is ftc *Tc − ∑i Pi Θi(Tc − Tc*) ≥ ftc * − Pi(Tc − Tc*) Tc* ∀ z ∈ Ωz

(A.35)

Rearranging eq A.35, we obtain * (ftc * + PT i c −

∑ Pi Θi)(Tc − Tc*) ≥ 0

∀ z ∈ Ωz

i

(A.36)

Thus, eq A.36 is a sufficient condition to verify the hypothesis A.30 and A.31; that is, the optimality is preserved in the presence of production variation. Equation A.36 consists of two terms. The first term (the one in the first parentheses) is calculated as the revenue gained by dedicating the whole production cycle to one product plus all the transition cost minus the nominal revenue. This term will be greater than zero in most common cases but there can be some extreme cases (price Pi is extremely low) where this term can be negative. The second term is positive since we proved that the optimal transition time is obtained as the minimum value given any production sequence. Therefore, eq A.36 is satisfied in most cases. In our case studies, f tc* + piT* − ∑iPiΘi ≥ 0, T* ≤ T, ∀z ∈ Ωz, therefore eq A.36 is satisfied. Q.E.D.

(A.28)

*) ∑i Pi Θi − ftc (θipt* , zipk * ∑i Θi + ∑i ∑p ∑k θipt*zipk ∑i Pi Θi − ftc (θipt* , zipk) ∑i Θi + ∑i ∑p ∑k θipt*zipk

(A.29)

We need to prove that * )|Θ + δ ≥ J(θipt* , zipk)|Θ + δ J *(θipt* , zipk i Θ i Θ i

i

■

(A.30)

Substituting eq A.28 into A.30 yields *) ∑i Pi Θi + Piδ Θi − ftc (θipt* , zipk

APPENDIX B (DETAILED INTEGRATED PROBLEM FOR CONTINUOUS PROCESSES)

Scheduling Constraints

* ∑i Θi + δ Θi + ∑i ∑p ∑k θipt*zipk

Product Assignment. Ns

∑i Pi Θi + Piδ Θi − ftc (θipt* , zipk) ∑i Θi + δ Θi + ∑i ∑p ∑k θipt*zipk

∑ yik

= 1,

∀i∈I (B.1)

k=1

(A.31)

Np

After fraction manipulations, we obtain a sufficient condition for eq A.31.

∑ yik

= 1,

∀k∈K (B.2)

i=1

* Pi ∑ ∑ ∑ θipt*zipk − Pi ∑ ∑ ∑ θipt*zipk i

+

p

k

ftc (θipt* ,

i

zipk) −

ftc (θipt* ,

p

(A.33)

Rearranging eq A.32, we obtain

≥ J |Θi

≥

i

ftc (θipt* ,

Rearranging eq A.29 we obtain

Substituting eq A.28) into J*(z*) ≥ J(z) yields

=

k

*, ∑ Θi + ∑ ∑ ∑ θipt*zipk

*) zipk

* ∑i Θi + ∑i ∑p ∑k θipt*zipk

J *|Θi =

p

i

The optimal production sequence is not affected by production perturbation Continued proof: Given that J∗ =

i

Tc* =

Lemma 3

ftc (θipt* ,

∑ Θi + ∑ ∑ ∑ θipt*zipk ,

yik′ = yi , k − 1 ,

k

*) ≥ 0 zipk

yi′1 = yi , N ,

(A.32)

s

Prove eq A.32: Let

∀ i ∈ I, k ∈ K, k ≠ 1

∀i∈I

zipk ≥ yik′ + ypk − 1, L

(B.3) (B.4)

∀ i ∈ I, p ∈ I, k ∈ K

(B.5)


Article

Industrial & Engineering Chemistry Research zipk ≤ yik′ ,

∀ i ∈ I, p ∈ I, k ∈ K

(B.6)

zipk ≤ ypk ,

∀ i ∈ I, p ∈ I, k ∈ K

(B.7)

xkė n = f n (tke , xke1 , ···, xken , uke1 , ···, ukem), ∀ n ∈ Sx , k ∈ K , e ∈ E n K1nke = f n (tke + 0.2113hk , xken + hk (0.25K1ke

where binary variables yik indicate the assignment of products in slots, yik = 1 if product i is produced in slot k, otherwise yik = 0; yi1′ are auxiliary variables. Constraints B.1 enforces that each product can only be produced in one slot and B.2 implies that only one product is produced in one slot. In constraints B.5−B.7 if binary zipk = 1, there is a transition from product i to product p within slot k; otherwise zipk = 0. In the case studies of this work it is assumed that Ns = Np. Demand Constraints. Wi ≥ DiTc ,

∀i∈I

(B.8)

Wi = Gi Θi ,

∀i∈I

(B.9)

0

Gi = F (1 − Xi),

∀i∈I

− 0.0387K 2nke), uke1 ···, ukem),

K 2nke = f n (tke + 0.7887hk , xken + hk (0.5387K1nke + 0.25K 2nke), uke1 ··· , ukem),

∑ θik ,

(B.12)

∀k∈K

Inequality B.11 sets up an upper bound for the time allowed for producing product i in slot k. In eqs B.12 and B.13, Θi, pk define the production time of product i and slot k, respectively. Timing Constraints. Np

θkt

=

∑∑

(B.24)

m m umin ≤ ukem ≤ umax ,

∀ m ∈ Su , k ∈ K , e ∈ E

(B.25)

Initial Values and Final Values of States and Manipulated Variable at Each Slot. Np

xinn , k

∀k∈K

=

∑ xss,n iyi ,k ,

∀ n ∈ Sx , k ∈ K

i=1

(B.14)

i=1 p=1

t1s

∀ n ∈ Sx , k ∈ K , e ∈ E

Constraints that Link Scheduling and Control Level

Np

θpitzipk ,

n n xmin ≤ xken ≤ xmax ,

Inequality B.24 and B.25 provide lower and upper bounds for the state and manipulated variables at each sample point.

(B.13)

i=1

(B.23)

The derivatives of the state variables at each sample step are calculated using eq B.19 where f represents the detailed dynamic model; and subindex e represent elements in each slot. Through the computing of intermediate variables K1, K2 in eqs B.20 and B.21, the state of the next sample step is obtained by eq B.23. Equation B.22 defines the sample size (length of each element) in slot k. Lower and Upper Bounds of States and Manipulated Variables.

∀i∈I

∑ θik ,

(B.22)

∀ n ∈ Sx , k ∈ K , e ∈ E

Np

pk =

∀k∈K

xkn, e + 1 = xken + hk (0.5K1nke + 0.5K 2nke),

(B.11)

k=1

θkt , Ne

hk =

Ns

Θi =

∀ n ∈ Sx , k ∈ K , e ∈ E (B.21)

(B.10)

∀ i ∈ I, k ∈ K

∀ n ∈ Sx , k ∈ K , e ∈ E (B.20)

Inequality B.8 implies that the total production amount of each product Wi, which is calculated by eq B.9, should be greater than the demand that is a product of demand rate Di and cycle time. Equation B.10 computers the production rate as a product of feed flow rate F0 and conversion Xi. Processing Times. θik ≤ θ maxyik ,

(B.19)

=0

Np

xk̅ n

(B.15)

=

∑ xss,n iyi ,k+ 1 ,

∀ n ∈ Sx , k ∈ K , k ≠ Ns

i=1

Np

tke = tks + pk +

∑ ∑ θpitzipk ,

Np

∀k∈K

x ̅Nns =

(B.16)

∑ xss,n iyi ,1 ,

∀ n ∈ Sx

(B.28)

i=1

tke

=

tke− 1 ,

≤ Tc ,

(B.27)

Np

i=1 p=1

tks

(B.26)

∀ k ∈ K, k ≠ 1 ∀k∈K

(B.17)

Np

uinm, k =

(B.18)

∑ uss,m iyi ,k ,

∀ m ∈ Su , k ∈ K

i=1

Equation B.14 computes the time of transition θtk from product i to p in slot k. The starting tsk and the ending time point tek for each slot are presented by eqs B.16 and B.17. The starting point is equal to the ending point of the prior slot. The ending point is calculated as the sum of the starting point, the transition time, and the production time. As shown in inequality B.18 all slots should finish in the cyclic time.

(B.29)

Np

uk̅ m =

∑ uss,m iyi ,k+ 1 ,

∀ m ∈ Su , k ∈ K , k ≠ Ns − 1

i=1

(B.30) Np

u ̅Nms =

Dynamic Optimization (Constraints at Control Level)

∑ uss,m iyi ,1 ,

∀ m ∈ Su

i=1

Discretization of the Dynamic Model Using Implicit RK-4 Method.

xkn,1 = xinn , k , M

∀k∈K

(B.31) (B.32) DOI: 10.1021/acs.iecr.5b01916 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research ukm,1 = uinm, k ,

∀k∈K

(B.33)

xkn, Ne = xk̅ n ,

∀k∈K

(B.34)

COSTtr = cost during transition period COSTtrpi = estimation of transition cost from product p to i COSTss = cost during steady state production Θ = production time of a product θt = transition time between products pk = process time at slot k tke = time point for each element in each slot umke = mth manipulated variable at element e of slot k xnke = nth state variable at element e of slot k xnk̅ = desired value of nth state at slot k u̅mk = desired value of mth manipulated variable at slot k xnss,i = nth steady state value of product i umss,i = mth steady manipulated value of product i xnin,k = nth state’s initial value in slot k umin,k = mth manipulated variable’s initial value in slot k hk = length of element in slot k θik = process time of product iat slot k θtk = transition time at slot k θtpi = transition time from product p to i θ̃tpi estimation of transition time from product p to i f tc = transition cost Θi = total processing time of product i K1, K2 = intermediate variables in implicit Runge−Kutta methods Gi = production rate for product i gi = production amount of product i gmax = maximum production amount of product i Tc = total production cycle time Wi = amount produced for product i P = revenue coefficient

The steady state values for each slot xnss,i, umss,i are computed beforehand by simulating the process at steady state. Equations B.27 B.28, and B.30 B.31 calculate the desired state and manipulated variables at each slot. Equations B.26 and B.29 provide the initial values of state and manipulated variables in slot k. Equations B.32 and B.33 link the initial values with the state and control value of the first element in the transitions; B.34 enforces that the state values at the end of transition are equal to the desired values (steady state values) in the current slot. Constraints B.26−B.33 demonstrate how scheduling and control are linked. For example, sequence y in B.26 and dynamic in B.23 are linked through B.32. Nomenclature: Indices

products = i slots = k elements = e states = n manipulated variables = m Sets

product set = I slot set = K element set = E state set = Sx manipulated variable set = Su

■

Constants

Np = number of products Ns = number of slots Ne = number of elements within each slot Nx = number of states Nu = number of manipulated variables Nspl = number of samples in the production cycle Di = demand rate for product i Cpi = price of product i Cr = unit cost of raw material θmax = upper bound on processing time Xi = conversion in steady state of continuous processes Fo = feed stream volumetric flow rate xnmin, xnmax = minimum and maximum value of nth state variable n n umin , umax = minimum and maximum value of mth manipulated variable

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge financial support from the NSF under Grant CBET 1159244.

■

REFERENCES

(1) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous Cyclic Scheduling and Control of a Multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45 (20), 6698−6712. (2) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous Cyclic Scheduling and Control of Tubular Reactors: Single Production Lines. Ind. Eng. Chem. Res. 2010, 49 (22), 11453−11463. (3) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous Scheduling and Control of Multiproduct Continuous Parallel Lines. Ind. Eng. Chem. Res. 2010, 49 (17), 7909−7921. (4) Nystrom, R. H.; Franke, R.; Harjunkoski, I.; Kroll, A. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29 (10), 2163−2179. (5) Nystrom, R. H.; Harjunkoski, I.; Kroll, A. Production optimization for continuously operated processes with optimal operation and scheduling of multiple units. Comput. Chem. Eng. 2006, 30 (3), 392−406. (6) Prata, A.; Oldenburg, J.; Kroll, A.; Marquardt, W. Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor. Comput. Chem. Eng. 2008, 32 (3), 463−476.

Common Variables

x = state variables u = manipulated variables q = output variables y = binary variable indicating assignment of production in slots z = binary variable indicating transfer between products J = objective function Tc = cycle time W = amount of product tek = ending time at slot k tsk = starting time at slot k Continuous case

Φ = profit Φ1 = revenue Φ2 = cost

N


Article

Industrial & Engineering Chemistry Research (7) Harjunkoski, I.; Nystrom, R.; Horch, A. Integration of scheduling and control–Theory or practice? Comput. Chem. Eng. 2009, 33 (12), 1909−1918. (8) Mitra, K.; Gudi, R. D.; Patwardhan, S. C.; Sardar, G. Resiliency Issues in Integration of Scheduling and Control. Ind. Eng. Chem. Res. 2009, 49 (1), 222−235. (9) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (10) Baldea, M.; Harjunkoski, I. Integrated production scheduling and process control: A systematic review. Comput. Chem. Eng. 2014, 71, 377−390. (11) Allgor, R. J.; Barton, P. I. Mixed-integer dynamic optimization I: problem formulation. Comput. Chem. Eng. 1999, 23 (4−5), 567−584. (12) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and optimal control of polymerization reactors. AIChE J. 2007, 53 (9), 2301−2315. (13) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous design, scheduling, and optimal control of a methylmethacrylate continuous polymerization reactor. AIChE J. 2008, 54 (12), 3160−3170. (14) Zhuge, J.; Ierapetritou, M. G. Integration of Scheduling and Control with Closed Loop Implementation. Ind. Eng. Chem. Res. 2012, 51 (25), 8550−8565. (15) Zhuge, J.; Ierapetritou, M. G. Integration of scheduling and control for batch processes using multi-parametric model predictive control. AIChE J. 2014, 60 (9), 3169−3183. (16) Chu, Y.; You, F. Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and large-scale global optimization algorithm. Comput. Chem. Eng. 2012, 47, 248−268. (17) Nie, Y.; Biegler, L. T.; Wassick, J. M. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 2012, 58 (11), 3416−3432. (18) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Lagrangean heuristic for the scheduling and control of polymerization reactors. AIChE J. 2008, 54 (1), 163−182. (19) Chu, Y.; You, F. Integration of production scheduling and dynamic optimization for multi-product CSTRs: Generalized Benders decomposition coupled with global mixed-integer fractional programming. Comput. Chem. Eng. 2013, 58, 315−333. (20) Chu, Y.; You, F. Integrated Scheduling and Dynamic Optimization of Complex Batch Processes with General Network Structure Using a Generalized Benders Decomposition Approach. Ind. Eng. Chem. Res. 2013, 52 (23), 7867−7885. (21) Chu, Y.; You, F. Integrated scheduling and dynamic optimization by stackelberg game: bilevel model formulation and efficient solution algorithm. Ind. Eng. Chem. Res. 2014, 53 (13), 5564− 5581. (22) Nie, Y.; Biegler, L. T.; Villa, C. M.; Wassick, J. M. Discrete Time Formulation for the Integration of Scheduling and Dynamic Optimization. Ind. Eng. Chem. Res. 2015, 54 (16), 4303−4315. (23) Chu, Y.; You, F. Integrated Planning, Scheduling, and Dynamic Optimization for Batch Processes: MINLP Model Formulation and Efficient Solution Methods via Surrogate Modeling. Ind. Eng. Chem. Res. 2014, 53 (34), 13391−13411. (24) Du, J.; Park, J.; Harjunkoski, I.; Baldea, M. A time scale-bridging approach for integrating production scheduling and process control. Comput. Chem. Eng. 2015, 79, 59−69. (25) Mahadevan, R.; Doyle, F. J.; Allcock, A. C. Control-relevant scheduling of polymer grade transitions. AIChE J. 2002, 48 (8), 1754− 1764. (26) Pontryagin, L. S.; Boltyanskil, V. G.; Gamkrelidge, R. V.; Mishchenko, E. F. The Mathematical Theory of Optimal Processes; Interscience: New York, 1962. (27) Kirk, D. E. Optimal Control Theory: An Introduction. PrenticeHall: London, 1970.

O


A Decomposition Approach for the Solution of ... - ACS Publications

Recommend Documents