Integrated Scheduling and Dynamic Optimization of Complex Batch

Article pubs.acs.org/IECR

Integrated Scheduling and Dynamic Optimization of Complex Batch Processes with General Network Structure Using a Generalized Benders Decomposition Approach Yunfei Chu and Fengqi You* Department of Chemical and Biological Engineering, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208, United States ABSTRACT: We address the integration of scheduling and dynamic optimization for batch chemical processes. The processes can have complex network structures, allowing material splitting and mixing. The integrated problem is formulated as a mixedinteger dynamic optimization problem where a continuous-time scheduling model is linked to the dynamic models via processing times, processing costs, and batch sizes. To reduce the computational complexity, we develop a tailored and efficient decomposition method based on the framework of generalized Benders decomposition by exploiting the special structure of the integrated problem. The decomposed master problem is a scheduling problem with variable processing times and processing costs, as well as the Benders cuts. The primal problem comprises a set of separable dynamic optimization problems for the processing units. By collaboratively optimizing the process scheduling and the process dynamics, the proposed method substantially improve the overall economic performance of the batch production compared with the conventional sequential method which solves the scheduling problem and the dynamic optimization problems separately. In comparison with the simultaneous method which solves the integrated problem by a general-purpose MINLP solver directly, the proposed method can reduce computational times by orders of magnitude.

1. INTRODUCTION

Even for the batch process, the complexity of the integrated problem depends on the process structure. Generally, there are two types of batch processes: the sequential process where material splitting and mixing are not allowed, and the network process where these operations are allowed.22 In a sequential process, the batch integrity is preserved and the batch size is fixed. Taking advantage of these special features, an efficient integrated method can be developed for the sequential process.19 However, a network batch process is intrinsically much more complex than a sequential batch scheduling problem. Consequently, methods applicable to a sequential process generally cannot be extended to a network process.22,23 Though there are previous studies on integrated problems of a batch process with general network structures,20,21 they did not provide an effective method to solve the integrated problem. A consequence is that the returned solution could be only suboptimal after a reasonable amount of time. These studies solved complicated MIDO problems using the simultaneous method which reformulated an MIDO problem into an MINLP problem by discretizing the differential equations.24,25 A generalpurpose MINLP solver was then invoked to solve the reformulated MINLP problems. However, the MINLPs were largescale and challenging to solve. Except for very simple examples, the simultaneous method encountered difficulties in solving the MINLP problems to optimality with limited computational resources.20,21 After a long computational time, the optimality gap could be still unsatisfactorily large.

In traditional production hierarchy, scheduling and dynamic optimization belong to two different layers. Thus the two problems are usually solved separately in a sequential way. The dynamic optimization problem is first solved for each processing unit to determine the recipe data for the scheduling problem. The recipe data are treated as fixed parameters after the scheduling decisions are optimized. Though widely applied, the sequential method only returns a suboptimal solution for the entire batch process. To improve the performance of the entire process, the scheduling problem and the dynamic optimization problem should be solved in an integrated way such that the recipe data can be optimized in terms of the economic criteria.1−3 However, solving the integrated problem is much more challenging than solving the subproblems sequentially. The complexity of the integrated problem arises from the coupling of binary variables in the scheduling problem with differential equations characterizing the dynamic models. The integration results in a complicated mixed-integer dynamic optimization (MIDO)4,5 or mixed-logic dynamic optimization (MLDO)6,7 problem. An MLDO problem can be reformatted into an MIDO by big-M constraints and a convex-hull formulation.8 Because the integrated problem needs to optimize the entire batch plant with multiple products and units, a number of dynamic models can be included in the integrated problem according to the number of processing units, demanded products, and batch repetitions. Due to the complexity, most integrated scheduling and dynamic optimization approaches are restricted to continuous processes,9−18 while relatively few methods are developed for batch processes.19−21 Compared with the continuous process, a batch chemical process can have a more complicated network structure, leading to a much more complex integrated problem. © XXXX American Chemical Society

Received: February 11, 2013 Revised: April 15, 2013 Accepted: May 13, 2013

A

dx.doi.org/10.1021/ie400475s | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

The remainder of the paper is organized as follows. Section 2 provides the formulation of the integrated problem. The GBD decomposition method is presented in section 3. Three case studies are given in section 4 to demonstrate the proposed methods. The conclusion is provided in section 5.

The objective of this work is to propose an efficient and effective method to solve the challenging integrated problem for batch processes with a general network structure. We first formulate the integrated problem as an MIDO by incorporating a continuous-time scheduling model with differential equations of the dynamic models. Next we discretize the differential equations into algebraic equations by using the orthogonal collocation method.24 The MIDO problem is transformed into a mixed-integer nonlinear programming (MINLP) problem. Finally, we develop a tailored solution method based on the framework of generalized Benders decomposition (GBD)26 to solve the formulated MINLP problem. The proposed GBD method is distinguished from existing GBD methods4,27 for solving a general MIDO problem. These methods decompose the binary variables from nonlinear equations. The resulting master problem is a mixed-integer linear program (MILP) including Benders cut constraints. Nevertheless, the primal problem is still a very complicated dynamic optimization where all dynamic models are coupled with the entire scheduling model. We observe that the integrated scheduling and dynamic optimization has a special structure besides the coupling of binary variables with nonlinear equations. Though the scheduling model and the dynamic models can be very complicated, equations that link the two kinds of problems are relatively simple, with linking variables of processing times, processing costs, and batch sizes. Taking advantage of this special structure, the proposed GBD method separates the entire scheduling problem from all dynamic models. The proposed GBD method has the following features: • Ef fective decomposition. The master problem is the scheduling problem with variable processing times and processing costs, and the primal problem is a set of separable dynamic optimization problems on each processing unit. • Fast computation. In the primal problem, the dynamic optimization problems can be solved separately. Moreover, the infeasibility of the primal problem is straightforward to cope with. Putting a lower bound on the processing time can guarantee the feasibility of each dynamic optimization problem. The tedious procedure to generate the Benders feasibility constraints can be avoided. The decomposition method enables us to solve a complex integrated problem to optimality within a short computational time when the direct solution methods fail. Three case studies with different complexities demonstrate advantages of the proposed method. In the simple case, the proposed method returns the same optimal solution as the simultaneous method. The computational time of the proposed method is, however, reduced by 2 orders of magnitude. In the two complex cases, the proposed method finds the optimal solution while the simultaneous method only finds a suboptimal one within 24 or 50 h. We should note that an integrated problem generally contains non-convex terms from the dynamic models. The non-convexity could incur difficulties for both the simultaneous method, using a local MINLP solver like SBB, and the GBD method. Neither a direct nor decomposed solution strategy can guarantee the global optimality. However, the decomposition method can take advantage of the problem structure to reduce the computational time. Because even a local solver can be time-consuming to solve a complex integrated problem, it is impractical to use a global MINLP solver, though it can guarantee the global optimality in theory.

2. FORMULATION OF INTEGRATED PROBLEM This section provides the model formulation of the integrated scheduling and dynamic optimization problem. The problem statement is given in section 2.1. The integrated problem is formulated as an MINLP problem after discretizing the differential equations in the MIDO problem. To formulate the integrated problem, we first present a continuous-time scheduling model in section 2.2 and then provide dynamic models including the discretization approach in section 2.3. The integrated model is finally presented in section 2.4. 2.1. Problem Statement. The integrated problem determines scheduling decisions and unit dynamics simultaneously. In a common scheduling problem, only the total amounts of materials are monitored. For example, the batch size of a task denotes the total mass of the materials processed in a unit. However, a dynamic model requires more-detailed information. Besides the total amount of all species, a reaction task, for example, should know the concentration of each species. Thus, we include the concentration information in the state task network (STN),28 which is a common graph representation for a batch process. An illustrative example of STN is shown in Figure 1, which is simple but contains all essential components. (A more complicated

Figure 1. Illustrative example of state task network with concentration of species.

example can be seen in Figure 4 in the case studies.) There is one reaction task node, represented by the rectangle, and there are three state nodes, represented by the circles. The coefficients β1 and β2 indicate the fraction values of the input states S1 and S2 to the task. Because only one output state exists, the unit coefficient for the output state S3 is omitted. In a pure scheduling problem, the coefficients represent the mass fractions. To formulate an integrated problem, the concentrations at each state node are required. From the concentrations in the input states, the initial condition in the reaction task can be calculated using the densities ρ1A, ρ2B, ρAB, and the coefficients β1, β2. Similarly, the final value in the reaction task can be expressed from the concentrations in the output state. The initial condition and the final value are also linked by the stoichiometric coefficients and the reaction conversion. In practice, the final concentrations in a task are often predetermined to satisfy the quality requirement. B



Article

The final concentrations in the task determine the concentrations in the output states. These states can also be input states for other downstream tasks. Therefore, we can regard the concentrations in each state as fixed parameters, which imply the initial condition and the final value of each task are parameters as well. The fixed final concentrations are determined by material balance. The determined values should be feasible, which can be attained by the dynamic models. This is a design problem,29,30 which is not the focus of this work. We assume the final concentrations are known and feasible for the integrated scheduling and dynamic optimization problem. It should be noted that, for the sequential batch process, the final concentrations can be regarded as variables and determined by the integrated method.19,21 Each production line in the sequential batch process does not interact with others and can be formulated into a multi-stage dynamic optimization problem. However, the network batch process is generally more complex than the sequential batch process. Because of material splitting and mixing, there are much stronger interactions among the tasks and the states. It will be shown that even if the final concentrations for the network batch process are fixed as parameters, solving the integrated problem directly by a general-purpose local MINLP solver is time-consuming. For the variable final concentrations, a number of nonlinear blending equations will be included in the integrated problem, making an already challenging problem mathematically intractable, except for very simple cases. The problem statement of the integrated scheduling and dynamic optimization is given in Chart 1. 2.2. Continuous-Time Scheduling Model. A great variety of batch scheduling models have been developed, including discrete-time models and continuous-time models.31,32 From the perspective of integrated method, the continuous-time model is advantageous to handle variable processing times. Therefore, we use the continuous-time scheduling model in this work. The model we adopt is a global time point model33 where the time points of all processing units are aligned globally. The adopted model has some advantages because it has fewer big-M terms than other common formulations.33 One should note that the adopted model mainly serves as a demonstration for the proposed solution strategy, which can be easily extended to other continuous-time formulations. The continuous-time model33 introduces a dummy task with index i = 0. When a unit is idle, the dummy task is assigned to unit. The model includes the following constraints: Time Slots. The model divides the scheduling horizon sh into a set of time slots indexed by k. The length of a time slot is denoted by SLk, which is constrained by

∑ SLk ≤ sh

Chart 1

one when task i is continuously executing in unit j at time point k, and YEijk, which is one when task i ends in unit j before time point k. The constraints regarding the timing relationships are YR ijk = YR ij(k − 1) + Yij(k − 1) − YEijk ,

(3) (1)

k

Yijk + YR ijk ≤ 1,

The starting time of the time slot SLk is defined as TSk. The time points of the model are the starting times of all time slots, which will be used to specify the location of a task bar on the Gantt chart. As the time slots are consecutive, the starting times of time slots have the relations of TSk = TSk − 1 + SLk − 1 ,

∀ i ∈ Ij , j , 1 < k < nk

∀k>1

∀ i ∈ Ij , j , 1 < k < nk

YR ijk + YEijk ≤ 1,

∀ i ∈ Ij , j , 1 < k < nk

(4) (5)

The set Ij denotes the set of all tasks including the dummy one, which can be executed in unit j. The variables YRijk and YEijk have the binary values; however, they can be relaxed to be continuous variables in [0,1] because only those with value of 0 or 1 can satisfy eqs 3−5.33 Unit Status. To denote the unit status, a variable Zjk is introduced, which is one when unit j begins a task at time point k. The unit status is determined by the timing variables as

(2)

Timing Relationships. A task can only start at a time point, but it can end at or before another time point. A binary variable Yijk is introduced to determine the unit assignment and production sequence, which is one when task i starts in unit j at time point k. To obtain a complete description of the timing relationships, two more variables are introduced: YRijk, which is

Zjk =

∑ Yijk , i ∈ Ij

C

∀ j , 0 ≤ k < nk (6)


Industrial & Engineering Chemistry Research Zjk =

∑ YEijk ,

Article

STsk ≤ caps ,

∀ j , 0 < k < nk (7)

i ∈ Ij

Zjk = 1 −

∑ YR ijk ,

STp(k = nk) ≥ demp ,

(8)

The variables Zjk are binaries; however, they can be assumed to be continuous variables in [0,1] because only those with value of 0 or 1 can satisfy eqs 6−8.33 Processing Time and Startup Time. Because a task can end before a time point, a positive variable TRjk is introduced to represent the remaining time for unit j to complete an ongoing task at time point k. The relating constraint is TR j(k + 1) ≥ TR jk +

∑ Yijkptij + ∑ Yijkstij − SLk , i ∈ Ij

i ∈ Ij

i ∈ Ij

∀ j , k < nk

SALES =

∑ prp STp(k= n ) k

where prp is the unit price of product p. The total material cost is the sum of consumed materials, TMC =

j

i ∈ I *j ∩ Is+

(12)

∀ i ∈ I *j , j , 1 < k < nk

βisBEijk +

∑ ∑ j

i ∈ I *j ∩ Is−

βisBijk ,

Bijk pcijvb (21)

where is the fixed cost and is the coefficient related to the batch size. The variable processing cost in the second term is determined by the dynamic model, and the expression in eq 21 is actually an approximation. We will discuss the detailed expression of the batch size in section 2.3. 2.3. Dynamic Models. The batch scheduling problem is solved based on the recipe data, some of which are actually determined by the dynamic models of processing units. A dynamic model can be generally expressed by a set of nonlinear differential and algebraic equations,

(14)

max where bmin denote the minimum and maximum batch ij and bij size of task i in unit j, respectively. Material Balance. The inventory level of state s at time point k is denoted by STsk. The material balance is expressed by the equality

∑ ∑

∑ i ∈ I *j , j , k < nk

pcfixed ij

(13)

STsk = STs(k − 1) +

Yijkpcijfixed +

i ∈ I *j , j , k < nk

∀ i ∈ I *j , j , k < nk

∀ i ∈ I *j , j , 1 < k < nk

∑

TPC =

(11)

bijminYEijk ≤ BEijk ≤ bijmax YEijk ,

(20)

where cm is the unit cost of materials indexed by m and is the initial inventory value. The total processing cost is expressed as

where I*j is the index set of all tasks except the dummy one which can be executed in unit j. The batch sizes are also constrained by the lower bounds and the upper bounds, which are

bijminYR ijk ≤ BR ijk ≤ bijmax YR ijk ,

k

s0m

(10)

∀ i ∈ I *j , j , k

∑ cm(sm0 − STm(k= n )) m

Batch Size. Three types of variables, Bijk, BRijk, and BEijk, are introduced corresponding to the three assignment variables Yijk, YRijk, and YEijk, respectively. Similar to eq 3, the batch sizes are related by the equality

bijminYijk ≤ Bijk ≤ bijmax Yijk ,

(19)

p

∀ j , 1 < k < nk

BR ijk = BR ij(k − 1) + Bij(k − 1) − BEijk ,

(18)

The sales are equal to the sum of produced products multiplied by their prices:

where ptij denotes the processing time of task i in unit j and stij denotes the task startup time. In a fixed-recipe scheduling problem, the processing times and the startup times are all constant so the startup times can be absorbed into the processing times. However, in the integrated problem, the processing times are variables so we express the two terms separately in constraint (9). The variables representing the remaining time in a unit is bounded from above by

∑ YR ijkstij ,

(17)

PROFIT = SALES − TMC − TPC

i ∈ Ij

∑ YR ijkptij +

∀p

where the products belong to a subset of states indexed by p. Objective Function. The objective function we consider is the maximization of the profit, which is expressed by the sales minus the total material cost and the total processing cost:

(9)

TR jk ≤

(16)

The inventory level of products needs to meet the order demand demp as

∀ j , 0 < k < nk

i ∈ Ij

∀ s, k

pcvb ij

Xḋ (Td) = Fd(Xd(Td), Ud(Td), Vd)

(22)

Gde(Xd(Td), Ud(Td), Vd) = 0

(23)

Gdi (Xd(Td), Ud(Td), Vd) ≤ 0

(24)

Φd = Hd(Xd(Tdfinal), Ud(Tdfinal), Vd)

(25)

Xd(0) = Xdinitial(Vd)

(26)

0 ≤ Td ≤ Tdfinal

(27)

where the subscript d is the index of the dynamic model, which will be linked to the indices of the scheduling problem when the integrated problem is formulated. For the purpose of compact notations, we use the vector form to express the variables x u and the constraints. Xd ∈ nd contains state variables, Ud ∈ nd

∀ s, k

(15)

v

includes time-varying input variables, Vd ∈ nd consists of time-

where βis denotes the coefficient of task i to state s. If the state provides materials for the task, the coefficient is negative. If the state stores materials produced by the task, the coefficient is positive. The coefficient is zero when the state is not linked to the task. The set I+s indexes all tasks with positive coefficients corresponding to state s and the set I−s indexes all tasks with negative coefficients. The inventory level is bounded by the storage capacity caps as

x

invariant input variables. Differential equations Fd ∈ nd (eq 22) characterize the system dynamics while the state and e input trajectories are constrained by equalities Gde ∈ nd (eq 23) i

and inequalities Gdi ∈ nd (eq 24). The scalar Φd ∈  (eq 25) denotes the processing cost. Without loss of any generality, the processing cost is expressed as the state and input variables at the D



Article

initial final time point Tfinal of the state variables is d . The initial value Xd assigned in eq 26. The time variable Td ∈  is bounded by the inequalities in eq 27. There are various solution methods for the dynamic optimization problem constrained by differential equations.34−36 Because the path constraints, regarding the safety and quality issues, are critical to a batch production, we use the simultaneous method24,25 to solve all dynamic optimization problems in this work. The simultaneous method uses the orthogonal collocation approach to discretize the differential equations. The collocation approach is actually an implicit Runge−Kutta method.37 The collocation points and the discrete state values can be expressed explicitly as

Tdrq = Tdr − 1 + cqΔd , Xdrq = Xdr − 1 + Δd

∀ r > 0, q

models, we introduce new variables VPTijk and VPCijk to substitute the original parameters ptij and pcvb ij . We should note that the parameter ptij only depends on task i and unit j, while the corresponding variable VPTijk additionally depends on time point k. When a task is executed multiple times in the same unit at different time points, the processing time of an execution can be different. The difference in the processing times is the outcome of the varying production environment when the task is executed at different times. For the same reason, the variables of processing costs also depend on time point k. In practice, not all tasks are dynamic systems. For these tasks, the processing times and the processing costs remain as parameters. We introduce an index set Idyn to denote the tasks j represented by dynamic models in unit j. The complement of dyn Idyn is denoted by Idyn have j j̅ . Only the tasks belonging to Ij variable processing times. Substituting the parameters of processing times by the corresponding variables, constraints (9) and (10) become

(28)

nc

∑ aqq ′Fd(Xdrq ′, Udrq ′, Vd),

∀ r > 0, q

q ′= 1

(29)

TR j(k + 1) ≥ TR jk +

nc

Xdr = Xdr − 1 + Δd ∑ bqFd(Xdrq , Udrq , Vd),

∀r>0 (30)

q=1

− SLk ,

The time interval is partitioned by a set of discrete points {Trq d }, where r denotes the finite element while q denotes the collocation point. Similarly, the state trajectories are discretized rq rq as Xd(Td) ≈ {Xrq d |Xd = Xd(Td )}, and the input trajectories are rq rq discretized as Ud(Td) ≈ {Ud |Ud = Ud(Trq d )}. The coefficients aqq′, bq, and cq are given by a specific Runge−Kutta method. In this work, we use the Radau IIA method, which has a robust performance for stiff differential equations.38 The corresponding coefficients are

TR jk ≤

∑ i ∈ I jdyn

YijkVPTijk +

∑ i ∈ I jjdyn

Yijkptij +

∑ Yijkstij i ∈ Ij

∀ j , k < nk

∑ i ∈ I jdyn

∑

YR ijkVPTijk +

i ∈ I jdyn

YR ijkptij +

∑ YR ijkstij , i ∈ Ij

∀ j , 1 < k < nk

Bilinear terms YijkVPTijk and YRijkVPTijk appear when the processing times are variables. Since the bilinear terms are products of a binary variable and a continuous variable, the constraints can be linearized as TR j(k + 1) ≥ TR jk +

∑ i ∈ I jdyn

a YTijk +

∑ i ∈ I jdyn

Yijkptij +

∑ Yijkstij − SLk , i ∈ Ij

∀ j , k < nk

(36)

a b YTijk + YTijk = VPTijk ,

∀ j , i ∈ I jdyn , 1 < k < nk

a 0 ≤ YTijk ≤ ptijmaxYijk ,

b 0 ≤ YTijk ≤ ptijmax(1 − Yijk),

Besides the differential equations, other constraints are transformed into discretized values: Gde(Xdrq ,

Udrq ,

Vd) = 0,

Gdi (Xdrq , Udrq , Vd) ≤ 0,

Φd =

(r = n )(q = nq) Hd(Xd r ,

∀ r > 0, q ∀ r > 0, q (r = n )(q = nq) Ud r ,

Vd)

∀ j , i ∈ I jdyn , 1 < k < nk

(38)

∀ j , i ∈ I jdyn , 1 < k < nk (39)

(31)

TR jk ≤

(32)

∑ i ∈ I jdyn

a YRTijk +

∑ i ∈ I jdyn

YR ijkptij +

∑ YR ijkstij , i ∈ Ij

∀ j , 1 < k < nk

(40)

(33)

a b YRTijk + YRTijk = VPTijk ,

Xd(r = 0) = Xdinitial(Vd)

(34)

a 0 ≤ YRTijk ≤ ptijmaxYR ijk ,

Tdfinal /nr

(35)

b 0 ≤ YRTijk ≤ ptijmax(1 − YR ijk),

Δd =

(37)

Equation 35 defines the length of a finite element, which implies the bounds on the time variable in eq 27. After the discretization procedure, the differential algebraic equations (22)−(27) are transformed into pure algebraic equations (28)−(35). 2.4. Integrated Model. Combining the scheduling model in section 2.2 and the dynamic models in section 2.3, the integrated model is formulated in this subsection. Some recipe parameters for the conventional scheduling problem turn out to be variables in the integrated problem. Because the processing times and the processing costs are determined by the dynamic

∀ j , i ∈ I jdyn , 1 < k < nk ∀ j , i ∈ I jdyn , 1 < k < nk

(41) (42)

∀ j , i ∈ I jdyn , 1 < k < nk (43)

ptmax ij

where is the upper bound of processing times. Similarly, when the processing costs are variables determined by the dynamic model, constraint (21) becomes

∑

TPC =

i ∈ I *j , j , k < nk

+

∑ i ∈ I jdyn ∩ I *j , j , k < nk

E

∑

Yijkpcijfixed +

VPCijk

i ∈ I jdyn , j , k < nk

Bijk pcijvb (44)



Article

The batch sizes Bijk in eq 21 are absorbed into the dynamic model, which will affect the processing cost VPCijk. To calculate the processing times and the processing costs, the discretized dynamic models are appended to the scheduling model through the interface constraints: final Tijk = VPTijk ,

Φijk = VPCijk ,

Vijk = Bijk ,

∀ j , i ∈ I jdyn , k < nk ∀ j , i ∈ I jdyn , k < nk

∀ j , i ∈ I jdyn , k < nk

(2) Project the optimization problem, including the objective function and the constraints, on the space of the complicating variables. (3) Decompose the optimization problem into a primal problem and a master problem according to the projection. (4) Solve the original problem through iterations between the primal problem and the master problem. The primal problem returns a lower bound of the optimal solution (for the maximization problem) while the master problem provides an upper bound. The iteration stops when the gap is sufficiently small comparing with a predefined threshold. According to the GBD framework, we first need to identify the complicating variables. For the integrated problem, the model is an MINLP. An apparent structure is that the model contains both discrete and continuous variables. The difficulty of an MINLP mainly stems from the presence of the discrete variables. When the discrete variables are fixed, the restricted problem is an NLP, which is much easier to solve than the original MINLP. Therefore, for an MINLP problem, a GBD method commonly labels the discrete variables as complicating ones.4,27 However, the integrated problem has a more special structure which we can take advantage to develop a more effective and meaningful GBD method. The integrated problem (Integrated_Problem) primarily consists of two types of submodels: the scheduling model and the dynamic models. The scheduling model is coupled with the dynamic models via the linking variables: the processing times Tfinal and VPTijk in eq d 45, the processing costs Φd and VPCijk in eq 46, and the batch sizes Vd and Bijk in eq 47. Despite the complexity of each model, the linking equations are relatively simple. The observation motivates us to decompose the scheduling model from the dynamic models by breaking the linking equations. To implement the decomposition scheme, we first rewrite the objective function in eq 18 in a more explicit expression in terms of the processing cost. Substituting the total processing cost in eq 44 to the profit in eq 18, we obtain

(45) (46) (47)

The index d of dynamic models is defined as the combination of (i, j, k) for the scheduling model. The number of dynamic models is the product of the numbers for the three indices. Thus, the integrated model can include a great number of dynamic models even for a moderate-size scheduling problem. The processing time is linked to the final time in eq 45, and the processing cost is linked to the cost value in eq 46. The batch size is connected to the dynamic model by the time invariant input variable in eq 47. The integrated model is formulated as

The extended scheduling model includes the variable processing times and processing costs for the tasks related with dynamic models. The dynamic models are discretized by the collocation method where the differential equations are transformed into algebraic equations. The integrated problem is formulated as an MINLP. Owing to the existence of the binary scheduling variables and the presence of a great number of dynamic models, the MINLP is really complicated. The simultaneous method solves the MINLP problem by a generalpurpose MINLP solver. However, because of the complexity, the method can be very inefficient with a large optimality gap after a long computational period.19−21

PROFIT = SALES − TMC −

∑

Yijkpcijfixed −

i ∈ I *j , j , k < nk

∑ i ∈ I jdyn ∩ I *j , j , k < nk

Bijk pcijvb −

∑ i ∈ I jdyn , j , k < nk

Φijk (48)

where the processing costs VPCijk are replaced by the values Φijk from the dynamic models according to eq 46. In the objective function, only the last term explicitly depends on the dynamic models. Thus, we classify the complicating variables as

3. SOLUTION STRATEGY BASED ON GENERALIZED BENDERS DECOMPOSITION To obtain an efficient solution for the integrated problem, we develop a decomposition method based on the GBD framework. Benders decomposition was first developed for linear programs39 and then generalized to nonlinear programs.26 There are variant methods following the GBD framework.40 The main idea of GBD is that when a subset of decision variables, called complicating variables, is temporarily fixed, the restricted problem is much easier to solve than the original problem. Based on the idea, a GBD method decomposes a complicated problem into a master problem determining the complicating variables and a primal problem determining the remaining variables. The two subproblems are solved iteratively to determine the optimal solution. Generally, a GBD method consists of several common steps:26,40−42 (1) Identify the complicating variables according to the special structure of a particular problem so that computational benefits can be achieved from the decomposed problems.

{Complicating variables} = ξ ∪ {VPTijk} ∀ j , i ∈ I jdyn , k < nk ∪ {Bijk } ∀ j , i ∈ I jdyn , k < nk

(49)

The subset of ξ is defined as ξ = {All variables in the scheduling model except VP Tijk and Bijk ,

for ∀ j , i ∈ I jdyn , k < nk }

(50)

The remaining variables are those for the dynamic models, which are represented by final {Remaining variables} = ωijk ∪ {Tijk } ∪ {Vijk},

for ∀ j , i ∈ I jdyn , k < nk F

(51)



Article

The subset ωijk is expressed as

Because the first term in the objective function depends only on the complicating variables, it becomes a fixed value denoted by

ωijk = {Variables in the dynamic model indexed by (i , j , k) final except Tijk , Vijk},


(l)

(l) (l) φ̅ (l) = φ(ξ ̅ , {VPTijk } ∀ j , i ∈ I jdyn , k < nk , {Bijk ̅ } ∀ j , i ∈ I jdyn , k < nk )

(52)

The fixed value can be moved outside the maximization procedure, so we only need to maximize the second term. There is a minus sign before the total processing cost. Thus, the maximization operation is changed to minimization when the minus sign is moved out. The constraints for the dynamic models in eq 55 are separable, so the dynamic models are actually independent of each other. Consequently, we can optimize them separately; i.e., minimizing the sum of the processing costs is equivalent to the sum of the minimum processing costs. (3) Primal Problem and Master Problem. The primal problem is derived from eq 58, which is a set of separable dynamic optimization problems minimizing the processing costs for fixed processing times and batch sizes. We express the dynamic optimization problem, which is indexed by (i,j,k) such that ∀j, i∈ Idyn j , k < nk, as

The complicating variables belong to the scheduling model, while the remaining variables belong to the dynamic models. To emphasize the linking variables, we separate them from the others in the above notation. Using the notation of the complicating variables and the remaining variables, we can express the integrated problem into a general but compact form as max{φ(ξ , {VPTijk} ∀ j , i ∈ I jdyn , k < nk , {Bijk } ∀ j , i ∈ I jdyn , k < nk ) −

∑

final Φijk (ωijk , Tijk , Vijk)}

i ∈ I jdyn , j , k < nk

(53)

s.t.

g Sch(ξ , {VPTijk} ∀ j , i ∈ I jdyn , k < nk , {Bijk } ∀ j , i ∈ I jdyn , k < nk ) ≤ 0

(54)

Dyn final g ijk (ωijk , Tijk , Vijk) ≤ 0,

(55)

final Tijk = VPTijk ,

Vijk = Bijk ,



(56)


(57)

In the objective function, φ denotes the sum of the terms in eq 48, except for the last one. The last term of the total processing cost is the sum of those over the dynamic models. The constraints in the vector gSch represent the scheduling model, while the constraints in the vector gDyn ijk represent the dynamic model indexed by (i,j,k). Equations 56 and 57 are linking constraints. Because the processing costs have been substituted into the objective function, the original linking constraint regarding the processing costs is discarded. Based on the compact formulation of the integrated problem, we develop a GBD method in the following steps: (1) Complicating Variables. We have already identified the complicating variables and the remaining variables in eqs 49 and 51. (2) Projected Problem. The projected problem is the restricted one with fixed complicating variables. In this work, we use the variable with the bar to represent the fixed value of the corresponding variable. The fixed value is dependent on the iteration number, which is denoted by l. When we fix the complicating variables at given values, i.e., ξ = ξ̅(l), VPTijk = (l) VPT (l) ijk , and Bijk = B̅ ijk , the projected objective function turns out to be

The optimal value of the primal problem is dependent on the (l) fixed complicating variables VPT (l) ijk and B̅ ijk , and we introduce ηijk to refer to the optimal value function. When the primal problem is solved, a feasible solution for the integrated problem is obtained, providing a lower bound of the optimal value: ν ̅ (l) = φ̅ (l) −


(l) (l) ηijk (VPTijk , Bijk ̅ )

v* = max{φ(ξ , {VPTijk} ∀ j , i ∈ I jdyn , k < nk , {Bijk } ∀ j , i ∈ I jdyn , k < nk ) −


ηijk (VPTijk , Bijk )} (65)

Because ηijk is evaluated by solving the primal problem, there is no explicit functional expression of ηijk using VPTijk and Bijk. To address the difficulty, the optimal value of the primal problem is expressed by the dual form

(l)

k

= max{φ̅ (l) −

∑

j

ηijk =

k



∑ ∑ dyn

i∈I j

= φ̅ (l) −

, j , k < nk

final [Φijk (ωijk , Tijk , Vijk)

(66)

where σijk, λijk, and γijk are Lagrangean multipliers corresponding to the constraints. The outer maximization problem with the Lagrangean multipliers does not need to be solved explicitly. We can express the maximum value by a set of equivalent inequalities as


, j , k < nk

∑ dyn

i∈I j

min

−γijk(Bijk − Vijk)]}

dyn i ∈ I j , j , k < nk

= φ̅ (l) − min{

max {

final , Vijk σijk ≥ 0, λijk , γijk ωijk , Tijk

Dyn final final − σijk g ijk (ωijk , Tijk , Vijk) − λijk (VPTijk − Tijk )

dyn i ∈ I j , j , k < nk

= φ̅ (l) + max{−

(64)

The master problem maximizes the expression in eq 58 by freeing the complicating variables:

(l) (l) } ∀ j , i ∈ I dyn , k < n , {Bijk v ̅ (l) = ν(ξ ̅ , {VPTijk ̅ } ∀ j , i ∈ I dyn , k < n ) j

(59)

final min Φijk (ωijk , Tijk , Vijk)

(58) G


Industrial & Engineering Chemistry Research ηijk ≥

min

final ωijk , Tijk , Vijk

Article

final {Φijk (ωijk , Tijk , Vijk)

Dyn final final (ωijk , Tijk , Vijk) − λijk (VPTijk − Tijk ) − σijk g ijk

− γijk(Bijk − Vijk)},

∀ σijk ≥ 0, λijk , γijk

(67)

Because there are infinite inequalities for any value of the Lagrangean multipliers, the problem with constraint (67) cannot be solved in practice. So we need to select a finite number of them. When the Lagrangean multipliers are the optimal dual variables, σijk * , λijk * , and γijk * , from the primal problem, we have final σijk * gDyn ijk (ωijk, Tijk , Vijk) = 0. Then a relaxed expression for constraint (67) is obtained as final * * (VPTijk − Tijk *) ηijk ≥ Φ*ijk − λijk ) − γijk*(Bijk − V ijk

(68)

where Φ*ijk is the optimal value of the primal problem, and final Tijk * and Vijk * are determined by the optimal solution. Constraint (68) is known as the Benders cut. Using the relaxed expression of the optimal value function ηijk, the master problem is formulated as

The constraints include the scheduling model (54). The Benders cut constraint (68) are appended sequentially to the master problem as expressed in eq 70, where l′ and l denote the iteration steps. The coefficients in the constraints are calculated as final * * Tijk * + γijk*V ijk Φ̅ (ijkl′) = Φ*ijk + λijk (l ) ′ λijk

Figure 2. Flow chart of the GBD method for solving integrated scheduling and dynamic optimization.

* = λijk

γijk(l ′) = γijk*

optimization problem where each dynamic model can be optimized independently. The solution of the primal problem is used to update the lower bound. The best solution found till the current iteration is recorded. If the gap between the upper bound and the lower bound falls below the specified value ε, the iteration stops, returning the best solution. Otherwise, the optimal dual variables are recorded. A new Benders cut constraint is generated and appended to the master problem. The MILP master problem (GBD_Master) then is solved with the Benders cuts. The optimal solution is used as the fixed value for the complicating variables in the next iteration, while the optimal function value is the new upper bound. The proposed GBD method is different from a common one where only the discrete variables are complicating ones and determined in the master problem.4,27 In the proposed GBD method, the master problem is an extended scheduling problem containing both discrete and continuous variables. Thus, the integer cut constraints in a common GBD are not included in the proposed method.

(72)

where Φ*ijk, Tfinal ijk *, V* ijk, λ* ijk, and γ* ijk are the optimal function value, solution, and dual variables from solving the primal problem (GBD_Primal) in the l′th iteration. Constraint (71) is introduced to guarantee the feasibility of the primal problem. The constraint of the lower bound, ptmin ij , for the processing time is required because in practice the process inputs are always constrained by the device and the processing time cannot be arbitrarily short. The lower bound of the processing time can be determined from the dynamic optimization of minimizing the processing time. As the processing cost terms are relaxed, the optimal value v* of the master problem (GBD_Master) is an upper bound for the integrated problem. (4) Iteration between the Two Subproblems. The master problem and the primal problem only provide the upper and lower bounds for the optimal value. We need to iterate between the two subproblems to reduce the gap between the upper bound and the lower bound. The flow chart of the GBD method is shown in Figure 2. We initialize the method by fixing the complicating variables at some value. Since they are variables in the scheduling model, we can simply fix them at the optimal solution of the scheduling problem solved by the sequential method. The primal problem (GBD_Primal) is then solved with the fixed complicating variables. It is a separable dynamic

4. CASE STUDIES To demonstrate the proposed GBD method, we solve an integrated problem for a batch plant in this section. To have an extensive investigation of the proposed method, we solve the integrated problem in three cases with different scheduling H



Article

FC in the jacket. The model parameters depend on the task and unit combination, which are listed in Appendix. Apart from the reaction tasks, other tasks executed in the mixing tank, the separation unit, and the packing unit have fixed recipe data. Just as those for a conventional scheduling problem, the processing times and the processing costs of the tasks are predefined parameters, given in the Appendix. The batch plant manufactures two products denoted by P1 and P2. The STN network is displayed in Figure 4. There are 13 states and 7 tasks. The dilution task is executed in the mixing tank. The three reaction tasks can be executed in either reactor. The separation task is executed in the separator. The two packing tasks are performed in the packing unit. The STN has a network structure including merging nodes, e.g. the three task nodes, and splitting nodes, e.g. the task node of Reaction 2. The batch size of each task is a variable which should be determined simultaneously with other scheduling decisions. Allowing material mixing and splitting is the main difference between a network batch process and a sequential one. In this example, we assume the volume of materials after mixing in an equipment unit or processing by a task is not changed. The assumption is made to simplify the calculation of the initial concentration and the final concentration for a reaction task. For example, the initial concentration of species B in Reaction 1 can be calculated in a straightforward way as (1/2 × 2)/(1/2 + 1/2) = 1 kmol/m3. We should note that this assumption can be removed by calculating the initial concentration with the densities of the mixer and the component species as shown in Figure 1. However, to prevent a reader being distracted by the tedious unit conversion procedure for each dynamic model, we adopt the assumption of the unchanged volume after mixing in the case study. Because the 100% pure products are typically costly to obtain in practice, the products of each task are impure. The purity is formulated as the quality constraints according to the user’s requirement. The parameters of the scheduling model are shown in the Appendix. After building the dynamic models and the scheduling problem, we can formulate the integrated problem through the linking constraints (45)−(47). For the particular differential equations (73), we have the following expressions for the processing times, the batch sizes, and the processing costs:

horizons and order demands. In each case, we also make a comparison with the sequential method and the simultaneous method to demonstrate the advantage of the proposed method. In this work, all optimization problems are modeled in GAMS 24.0.1 and solved using a PC with Intel Core i5-2400 CPU @ 3.10 GHz, 8 GB RAM, and Window 7 64-bit operating system. 4.1. Problem. The batch plant consists of five processing units displayed in Figure 3, including one mixing tank, two reactors (one large and one small), one separation unit, and one packing unit. Among the five units, we model the dynamics of the two reactors. The tasks performed in the reactors are described by a set of differential equations. We can change the recipe data of these tasks, e.g. the processing times and the processing costs, by manipulating the input variables of the reactors. The reactors are surrounded by heating/cooling jackets. The input variables are the flow rates of the heating/cooling fluid of the jackets. Manipulating these flow rates can change the jacket temperature, then the reactor temperature through the heat exchanger, and finally the chemical kinetics by altering the reaction rates. The dynamic model of a reactor is described by a set of differential and algebraic equations as ⎧d ⎪ C1(t ) ⎪ dt ⎪ ⎪ ⎪ ⎪ d Cc(t ) ⎪ dt ⎪ ⎪ ⎪ ⎪d ⎨ Cnc(t ) ⎪ dt ⎪ ⎪d ⎪ dt TR (t ) ⎪ ⎪ ⎪ d T (t ) ⎪ dt J ⎪ ⎪ ⎪ ⎩

nb

=

∑ νb1R b b=1

⋮ nb

=

∑ νbcR b b=1

⋮ nb

=∑ νbncR b b=1 n

=

=

∑bb= 1 R b( −Δhb) ρR cR

+

ua(TJ − TR ) VR ρR c R

FH(TH − TJ) + FC(TC − TJ)

+

vJ ua(TR − TJ) vJρJ c J

Tdf : = t f

(73)

Vd : = VR nc

R b = Kb∏ (Cc)νbc

Φd : = c hot

c=1

Kb = zb e−erb/ TR

∫0

tf

FH(t ) dt + c cool

∫0

tf

FC(t ) dt

The definition of the final time Tfd is straightforward. The batch size Vd is related to the reactor volume. The processing cost Φd is defined as the cost of the consumed heating/cooling fluid, which is equal to the unit prices chot and ccool times the consumed fluid. The consumed values are calculated by the integrals over the time interval. The primary goal of the decomposition method is to cope with the complexity stemming from integrating the scheduling problem with the dynamic optimization problem. Therefore, an individual dynamic model is set to be a simple one. Such simple dynamic models are often used in the previous studies.19−21 Also, following these studies, some tasks are not described by differential equations, which have the fixed processing times. However, even if the constituent dynamic models are not complex, the integrated problem is really complicated. First, the

For clarity, the index of the dynamic model given by the scheduling problem is removed. The state variables are the species concentrations, Cc for c = 1, ..., nc, the reactor temperature TR, the jacket temperature TJ. The differential equations describe the chemical kinetics, the change in the reactor temperature and the jacket temperature. The chemical kinetics follows the mass rate law where the reaction order of a species is determined by its stoichiometric coefficient, vbc, in a participant reaction. The rate coefficients Kb depend on the reactor temperature, expressed according to the Arrhenius equation where the pre-exponential factor zb and the activation energy erb are parameters. The control variables are the flow rate of the heating fluid FH and the flow rate of the cooling fluid I



Article

Figure 3. Diagram of the batch plant in the case study.

Figure 4. STN of the batch plant. A task has the same color as the processing units which can execute the task in Figure 2.

integrated problem includes a number of dynamic models owing to the repetition of batch productiona large-scale problem results. Second, the large-scale dynamic optimization problem is substantially complicated by the coupling with the binary variables in the scheduling model. Third, the batch process has a network structure which is more complicated than the sequential batch process. The complexity will be demonstrated by long computational times for the simultaneous methods. In the scope of integrated scheduling and dynamic optimization, the case studies given in this section are more complicated than those in the literature.19−21 Therefore, we do not include complex dynamic models because they can also complicate the integrated problem. 4.2. Results. The integrated problem is solved in three cases with different parameters listed in Table 1. The model complexity increases from case I to case III. To avoid unnecessary repetitions, we only report detailed results for case II. Results for case I and case III are summarized. For each case, the three solution strategies are implemented and compared.

Table 1. Scheduling Horizons and Order Demands in Three Cases value parameter

description

case I

case II

case III

sh dem1 dem2

scheduling horizon demand for product P1 demand for product P2

14 h 1 kmol 1 kmol

24 h 3 kmol 2 kmol

30 h 4 kmol 3 kmol

In this example, we solve the dynamic optimization problem by the full discretization approach presented in section 2.3. This approach is an essential component for the simultaneous method, by which the MIDO problem is reformulated as an MINLP problem. For the sequential method and the decomposition method, other approaches can also be used to solve the dynamic optimization problem. However, using the same approach for the dynamic optimization provides a fair comparison because the settings of the discretization method can also affect the J



Article

Table 2. Results for Different Numbers of Time Points in Case II no. of time points objective (m.u.) CPU (s)

7

8

9

10

11

12

infeasible 0.1

893.5 0.28

1428.1 0.31

1720.6 2.2

1720.6 48.9

1720.6 790.8

Table 3. Model and Solution Statistics for the Three Methods To Solve the Integrated Problem in Case II method sequential

simultaneous d

proposed

total

obj. (m.u.) CPU (s) gap (%) iterations

1720.6 4.1 1.0a −

2191.7 86 400 27.5e −

2650.5 385.6 1.0 6

dynamic models/primal problem

finite elements equations variables CPU (s) solver

30 1449 × 6b 1377 × 6b 1.9c CONOPT 3

30 − − − −

30 1509 × 54f 1436 × 54f 37.5g CONOPT 3

scheduling model/integrated model/master problem

type time points equations

MILP 10 1397

MINLP 10 88338

variables binary var. CPU (s) gap (%) solver

930 145 2.2 1.0 CPLEX 12

83772 145 86400 27.5 SBB with CONOPT 3

MILP 10 2471 (1st iter.), 2795 (6th iter.) 2030 145 348.1h 1.0 CPLEX 12

a

The gap for the sequential method is the same as that for the scheduling model. bThe equation and variable numbers are those for a single dynamic model times 6. The sequential method solves 6 dynamic optimization problems separately. cThe computational time is the total one for solving the 6 dynamic optimization problems. dThis is the best value found when the resource limit of 24 h is reached. The solution is found at the 2254th node as the computational time is equal to 66786 s. eThis is the gap when the solver terminates. The upper bound is 3023.5. fThe equation and variable numbers are those for a single dynamic model times 54. The primal problem of the proposed method consists of 54 separate dynamic optimization problems. gThe computational time is the total one for the primal problems over all iterations. hThe computational time is the total one for the master problems over all iterations.

To determine the number of time points required, we solve the scheduling problem by a trial-and-error procedure. The results are listed in Table 2. The scheduling problem is infeasible when the number of time points is equal to 7. With the number starting from 8, an optimal solution is obtained. The objective function increases as the number grows until the number is 10. Further increasing the number to 11 or 12 does not improve the optimal value even much more computational time is spent. Based on the observations, we determine the number of time points as 10. The model and solution statistics for the sequential method are listed in Table 3. The dynamic optimization problems are solved by the NLP solver CONOPT 3 after the discretization procedure. The scheduling problem is solved by the MILP solver CPLEX 12. The advantage of the sequential method is its short computational time. When the dynamic optimization problems are separated from the scheduling problem, each subproblem is straightforward to solve. The total computational time for solving all subproblems is merely 4.1 s. Though computationally efficient, the sequential method suffers from the poor solution quality because it is unable to optimally coordinate the two types of subproblems. As determined by the local criterion of minimizing the processing times, the recipe data are not necessarily the optimal ones for the scheduling problem.

computational performance. A critical user-tuning parameter in the full discretization approach is the number of finite element. It should be determined according to a particular problem. Using the trial-and-error procedure, we determine the number of the finite elements as 30. Considering this number is used for the three methods, the specific value is not a main factor determining the difference in their computational performance. We also note the number is independent of the scheduling horizon and the order demands in the scheduling problem. Thus the number of 30 is used for all cases. To set up the scheduling problem, we also have a critical user-tuning parameter, the number of time points. The number is dependent on the scheduling parameters, so we need to tune it for each case. We determine the number using the sequential method for its efficiency. 4.2.1. Results for Case II. We first apply the sequential method to solve the problem with the parameters of case II. The sequential method solves the dynamic optimization problems of minimizing the processing times prior to the scheduling problem. The dynamic optimization problems are solved for each dynamic model. Since the batch plant consists of 3 reaction tasks and 2 reactors, there are 6 dynamic models. After the dynamic optimization problems are solved, the returned recipe data serve as the fixed parameters for the scheduling problem. K



Article

The objective value returned by the sequential method is only 1720.6 m.u. The value can be improved by applying the simultaneous method to solve the integrated problem. The model and solution statistics for the simultaneous method are summarized in Table 3. By using the integrated optimization approach, the objective value increases significantly to 2191.7 m.u., resulting in a 27.4% improvement comparing with the sequential method. The improvement demonstrates the advantage of the simultaneous method. Meanwhile, the problem complexity also increases dramatically. A complicated MINLP problem has to be solved, which includes 88 338 equations and 83 772 variables (145 binary ones). The complexity leads to a very long computational time. The MINLP problem is solved by SBB which is the common solver for the integrated problems.11,13,20,21 We initialize the problem by the solution of the sequential method. After 24 h, the optimality gap still remains as large as 27.5%. We should note that the computational inefficiency of such a long computational time and a large gap is also reported by other studies.19−21 For such a large gap, one can naturally doubt if a potentially better solution exits. To reduce the complexity of the integrated MINLP problem, we apply the proposed decomposition method. The model and solution statistics are summarized in Table 3. The iteration starts from fixing the complicating variables at the solution of the scheduling problem returned by the sequential method. With the gap of 1% set as the same as the simultaneous method, the iteration between the primal problem and the master problem stops at the sixth step. The evolution of the upper/lower bound is shown in Figure 5a, and the computational times are shown in Figure 5b. The master problem accounts for the most portion

Figure 6. Gantt charts returned by (a) the sequential method, (b) the simultaneous method, and (c) the proposed method. The notations of the processing units are as follow: M, mixing unit; RI, large reactor; RII, small reactor; S, separator; and P, packing unit. The notations of the tasks are as follow: D, dilution; R1, Reaction 1; R2, Reaction 2; S, separation; P1, packing 1; and P2, packing 2. The vertical dashed lines along with the beginning and ending solid lines indicate the locations of the time points. There are task transition times in a unit except the first task. The number on a task bar shows the batch size of the task.

of the total computational time. The master problem is the scheduling problem with variable processing times and processing costs, as well asthe Benders cut constraints. From Table 3, we can see that the master problem contains much more equations than the scheduling model for the sequential method. The increase in the problem size is mainly attributed to the linearization of the bilinear terms in eqs 37−43. The contribution of the Benders cut constraints is marginal because only 54 equations are added in each iteration: 54 = 3 (#reaction tasks) × 2 (#reactors) × 9 (#time points − 1). Since the varying processing times and processing costs increase the problem complexity,

Figure 5. Iteration results of the decomposition method in case II: (a) bounds and (b) computational times. L



Article

Figure 7. Dynamic trajectories of the dynamic model which describes Reaction 2 starting from time slot 2 in the large reactor (RI).

the computational time of the master problem is larger than the scheduling problem with fixed processing times and processing costs. The proposed decomposition method significantly improves the profit to 2650.5 m.u., 54.1% higher than the sequential method. The computational time is 385.6 s, which is much shorter by 2 orders of magnitude than the simultaneous method. More importantly, the proposed method solves the integrated problem to the 1.0% optimality gap. As a result, it returns a much higher profit by 21.0% than the simultaneous method. Due to the large optimality gap, the best solution (2191.7 m.u.) searched by the simultaneous method is still a suboptimal one. To get more insights into the results returned by the three methods, the Gantt charts for the three methods are shown in Figure 6. Comparing the Gantt charts, we see that the simultaneous method and the decomposition method increase the processing times of the reaction tasks. The processing times are minimized by the sequential method because short processing times can have a positive impact on the profit. When the processing times are reduced, the reactors can perform more reaction tasks within a given scheduling horizon. Thus, the short processing times increase the production and in turn boost the sales. Despite the positive effect, reduction in the processing times also entails a negative impact by increasing the processing cost. Since the processing time and the processing cost are determined by the dynamic model of a reaction task, their relation

Table 4. Components in the Objective Function Returned by the Three Methods method profit (m.u.) sales (m.u.) material cost (m.u.) processing cost (m.u.)

sequential

simultaneous

proposed

1720.6 6840.0 2575.0 2544.4

2191.7 5880 2338.3 1350

2650.5 6600.0 2508.4 1441.1

can be revealed by showing the trajectories of the state and input variables. For an instance, the trajectories of Reaction 2 starting from the second time point in reactor I (the large one) are shown in Figure 7. From the Gantt charts in Figure 6, we can see this is the second reaction task executed in reactor I in all schedules determined by the three methods. The concentration trajectories have very similar patterns for the three methods. The concentrations of two reactants A and D gradually decrease while that of the product E constantly increases. The reaction task finishes when the product concentration reaches the predetermined quality value (1.8 kmol/m3). The temperature trajectories are also similar for the three methods. The reactor temperature is initially increased to accelerate the reaction rate. It then stays at some value for a period in which the temperature-dependent rate coefficient keeps the largest value. Finally, the reactor temperature is decreased so that the outlet materials can be cooled enough for a downstream processing M



Article

Table 5. Model and Solution Statistics for the Three Methods To Solve the Integrated Problem in Case I method sequential

simultaneous

proposed

total

obj. (m.u.) CPU (s) gap (%) iteration

444.8 2.0 1.0 −

866.8 1591.5 1.0 −

865.7 16.8 1.0 4



30 1449 × 6 1377 × 6 1.9 CONOPT 3

30 − − − −

30 1485 × 30 1421 × 30 16.3 CONOPT 3



MILP 6 745

MINLP 6 52910


538 85 0.1 1.0 CPLEX 12

50240 85 1591.5 1.0 SBB with CONOPT 3

MILP 6 1387 (1st iter.), 1477 (4th iter.) 1206 85 0.5 1.0 CPLEX 12

Table 6. Model and Solution Statistics for the Three Methods To Solve the Integrated Problem in Case III method sequential

simultaneous

proposed

total

obj. (m.u.) CPU (s) gap (%) iteration

2342.1 23.9 1.0 −

3184.3 180000.0 22.5 −

3684.9 26658.3 1.0 6



30 1449 × 6 1377 × 6 1.9 CONOPT 3

30 − − − −

30 1485 × 66 1421 × 66 43.1 CONOPT 3



MILP 12 1723

MINLP 12 106052


1126 175 22.0 1.0 CPLEX 12

100538 175 86400 22.5 SBB with CONOPT 3

MILP 12 3013 (1st iter.), 3343 (6th iter.) 2442 175 26615.2 1.0 CPLEX 12

to make a trade-off between these conflicting factors results in a poor performance of the sequential method. By contrast, both integrated methods increase the processing time to reduce the processing cost. By making a better balance of the conflicting factors, the overall performance of the batch production is improved. The trade-offs made by the integrated method can also be investigated from the components in the objective function in Table 4. The sequential method generates the highest sales and largest margin between the sales and the material cost. The results demonstrate the positive effect of the minimum processing times. However, the resulting processing cost is also the largest one. The simultaneous method and the proposed method both reduce the processing cost by sacrificing some

task. Such a three-stage pattern of the reactor temperature is a commonplace in process control.43 Owing to the heat exchange, the reactor temperature is controlled by the flow rate of the heating/cooling fluid through the jacket. To generate the temperature pattern, the heating fluid is used to increase the reactor temperature at the beginning of the reaction while the cooling fluid is used to reduce the temperature at the final stage. The temperature and flow rate trajectories demonstrate that a short processing time leads to a large processing cost by consuming more heating/cooling fluid. The processing times and the processing costs are actually conflicting factors in the integrated problem. The minimum processing time maximizes the sales, meanwhile brings out the large processing cost. Inability N



Article

Table A.1. Models of Reaction Kinetics a

task Reaction 1

reaction kinetics K1

B + C → 2D K1 = z1 e−er1/ TR z1 = 1 × 107 m 3/(kmol ·h) er1 = 5 × 103 K

Reaction 2

K2

−er2 / TR

z 2 = 1.2 × 103 m 3/(kmol ·h) er2 = 2 × 103 K

Reaction 3

differential equations

reactor

symbol

description

value

unit

dC B(t ) = − K1C BCC dt

RI (large)

vIR,max vIJ ρIJ cIJ uaI TIH TIC TIR,max

capacity of reactor volume of jacket density of jacket heat capacity of jacket heat transfer coefficient temperature of heating water temperature of cooling water maximum temperature of reactor maximum temperature of jacket maximum flow rate of heating and cooling water

1.5 0.5 1 × 103 4.2 3 × 104 370 300 370

m3 m3 kg/m3 kJ/(kg·K) kJ/(h·K) K K K

370

K

10

m3/h

vIIR,max

volume of reactor

1

m3

vIIJ ρIIJ cIIJ uaII TIIH TIIC TIIR,max

volume of jacket density of jacket heat capacity of jacket heat transfer coefficient temperature of heating water temperature of cooling water maximum temperature of reactor maximum temperature of jacket maximum flow rate of heating water

0.3 1 × 103 4.2 2 × 104 370 300 370

m3 kg/m3 kJ/(kg·K) kJ/(h·K) K K K

370

K

5

m3/h

dCC(t ) = − K1C BCC dt dC D(t ) = 2K1C BCC dt

A + D → 2E K2 = z2 e

Table A.4. Parameters of Reactors

dCA(t ) = − K 2CAC D dt

TIJ,max

dC D(t ) = − K 2CAC D dt

FImax

dC E(t ) = 2K 2CAC D dt

K3

RII (small)

dCC(t ) = − K3CCC E dt

C+E→F K3 = z 3 e−er3/ TR z 3 = 2 × 104 m 3/(kmol ·h) er3 = 3 × 103 K

dC E(t ) = − K3CCC E dt dC F(t ) = K3CCC E dt

a Parameters z1, z2, and z3 are pre-exponential factors. Parameter er1, er2, and er3 are normalized activation energies which absorb the gas constant.

TIIR,max FIImax

Table A.2. Other Parameters of Reaction Tasks task

symbol

description

Reaction 1

Δh1 ρ1R c1R

heat of reaction density of reactor heat capacity of reactor

1 × 103 1 × 103 3

value

kJ/kmol kg/m3 kJ/(kg·K)

Reaction 2

Δh2 ρ2R c2R


−2 × 103 1 × 103 3.5


Reaction 3

Δh3 ρ3R c3R


5 × 103 1 × 103 4


Table A.5. Price Value

unit

state

initial

final

unit

Reaction 1

CB CC CD TR TJ

1.0 1.0 0 300 300

0.05 0.05 1.9 ≤320 −

kmol/m3 kmol/m3 kmol/m3 K K

Reaction 2

CA CD CE TR TJ

1.0 0.95 0 300 300

0.1 0.05 1.8 ≤320 −


Reaction 3

CC CE CF TR TJ

1.0 0.9 0 300 300

0.2 0.1 0.8 ≤320 −


description

value

unit

pr1 pr2 cA cB cC chot ccool

sale price of product P1 sale price of product P2 unit cost of raw material A unit cost of raw material B unit cost of raw material C unit cost of heating fluid unit cost of cooling fluid

700 1200 100 150 200 10 1

m.u./kmol m.u./kmol m.u./kmol m.u./kmol m.u./kmol m.u./m3 m.u./m3

sales. For example, the proposed method substantially reduces the processing cost from 2544.4 m.u. of the sequential method to 1441.1 m.u. Meanwhile, the margin between the sales and the material cost only slightly decreases from 4265.0 m.u. of the sequential method to 4091.6 m.u. Obviously, a much better trade-off is made by the proposed method, yielding a much higher profit. 4.2.2. Results for Cases I and III. We have made a detailed comparison of the proposed method with the sequential method and the simultaneous method according to the parameters of case II in Table 1. Since the scheduling horizon and the order demand strongly affect the problem complexity, we make more comparisons using parameters of cases I and III in Table 1. In case I, the scheduling horizon and the order demand are decreased. A smaller problem is obtained than that in case II. Similar to what we do in case II, we use the sequential method to determine the number of time points. For the small problem in case I, 6 time points are enough. The model and solution statistics of the three methods are listed in Table 5. For the sequential method, the dynamic optimization problem is solved independently with the scheduling problem. So the results of the dynamic optimization is not changed when

Table A.3. Initial Conditions and Final Values of Reaction Tasks task

symbol

O



Article

Table A.6. Fixed Batch Cost (m.u./batch) dilution

Reaction 1

Reaction 2

Reaction 3

separation

packing 1

packing 2

10 − − − −

− 30 20 − −

− 30 20 − −

− 30 20 − −

− − − 100 −

− − − − 50

− − − − 50

mixing tank large reactor small reactor separator packing unit

than the sequential method. The reduction in the problem size also benefits the proposed method. The computational times for both subproblems are shorter than case II, especially for the master problem. The primal problem accounts for 97.0% of the total computational time which is equal to 16.8 s. In terms of the 1.0% optimality gap, the proposed method returns the same optimal solution with the simultaneous method. However, the computational time is reduced by nearly 2 orders of magnitude. When the problem complexity reduces, the advantage of the proposed method in computational efficiency is still evident. To investigate the performance of the proposed method in a more complex problem than case II, we solve the integrated problem in case III. The scheduling horizon increases from 24 to 30 h. More time points are required and a larger problem is obtained. The number of the time points is determined as 12 by applying the sequential method. The model and solution statistics of the three methods are listed in Table 6. Because the number of time points determines the number of binary variables, increase in the number of time points can substantially increase the problem complexity due to the combinatorial nature. This is reflected on the computational time for the sequential method to solve the scheduling problem. Though only 2 more time points are added compared with case II, the computational time increases from 2.2 to 22.0 s. Because the dynamic optimization problems in the sequential method are not changed with the scheduling horizon, the computational time of the dynamic optimization is still 1.9 s. The total computational time of the sequential method is 23.9 s in case III. The increase in the problem complexity has a more significant effect on the integrated methods. The simultaneous method returns a solution with 22.5% optimality gap when the resource limit of 50 h is reached. Though suboptimal owing to the large gap, the profit returned by the simultaneous method is still much larger than that of the sequential method (3184.3 vs 2342.1 m.u.). The suboptimal solution can be improved by the proposed method. The decomposition algorithm stops after 6 iterations. The returned profit is 3684.9 m.u., 15.7% larger than that of the simultaneous method. The increase in the time points makes the master problem more difficult to solve. The computational time of the master problem is 26615.2 s (7.4 h). It accounts for 99.7% of the total computational time while that of the primal problem is only 74.4 s.

Table A.7. Unit of Variable Batch Cost (m.u./m3) mixing tank separator packing unit

dilution

separation

packing 1

packing 2

30 − −

− 100 −

− − 50

− − 50

Table A.8. Upper and Lower Bounds of Batch Size symbol

description

value

unit

bmax M bmax RI bmax RII bmax S bmax P bmin M bmin RI bmin RII bmin S bmin P

maximum batch size in mixing unit maximum batch size in large reactor maximum batch size in small reactor maximum batch size in separator maximum batch size in packing unit minimum batch size in mixing unit minimum batch size in large reactor minimum batch size in small reactor minimum batch size in separator minimum batch size in packing unit

2 1.5 1 2 1 0.2 0.15 0.1 0.2 0.1

m3 m3 m3 m3 m3 m3 m3 m3 m3 m3

Table A.9. Storage Capacity state

type

capacity

M1, M2, M3 S1 W1 I1, I2, I3, I4 I5, I6

raw materials solvent waste intermediate products intermediate products

unlimited unlimited unlimited 2 m3 5 m3

Table A.10. Fixed Processing Time (h) mixing tank separator packing unit

dilution

separation

packing 1

packing 2

1.5 − −

− 3.0 −

− − 1.5

− − 1.5

the scheduling horizon varies. However, the scheduling problem is substantially reduced and the computational time is only 0.1 s. The total computational time of the sequential method is 2.0 s. For the simultaneous method, the size of the integrated problem becomes small. Therefore, the formulated MINLP problem can be solved to the 1% optimality gap. The simultaneous method increases the profit value to 866.8 m.u., much larger than 444.8 m.u. of the sequential method. Meanwhile, the computational time is 1591.5 s, much longer

5. CONCLUSION The integrated scheduling and dynamic optimization on a batch process is formulated into a large-scale MINLP, which is challenging

Table A.11. Task Starting Time (h) mixing tank large reactor small reactor separator packing unit

dilution

Reaction 1

Reaction 2

Reaction 3

separation

packing 1

packing 2

0 − − − −

− 0.2 0.2 − −

− 0.5 0.5 − −

− 0.4 0.4 − −

− − − 3.0 −

− − − − 0

− − − − 0

P



Article

p product q, q′ collocation point r finite element s state

to solve. The continuous-time scheduling model is linked to dynamic models through processing times, processing costs, and batch sizes. Such a complicated MINLP is difficult to directly solve to optimality within limited computational resources. To overcome the computational complexity, we propose an efficient and effective decomposition method based on the GBD framework. The method decomposes the dynamic models from the scheduling model by breaking the linking equations on processing times, processing costs, and batch sizes. The decomposed primal problem includes a series of separable dynamic optimization problems while the master problem is the scheduling problem with Benders cuts. The proposed method solves the integrated problem for the network batch process where material splitting and mixing are allowed. The initial and final values for the dynamic models are assumed to be known and feasible. The performance of the proposed decomposition method is demonstrated by a batch plant. In case I with the shortest scheduling horizon, the simultaneous method and the proposed method return the same optimal solution. However, the computational time of the proposed method is much less than the simultaneous method by nearly 2 orders of magnitude. In case II with a medium scheduling horizon, the simultaneous method only returns a suboptimal solution within the resource limit set as 24 h. The remaining gap is as large as 27.5%. The suboptimal solution is further improved by 21.0% using the proposed method, which solves the problem to 1.0% optimality gap in about 6.5 min. In case III with the large scheduling horizon, the proposed method solves the integrated problem to the optimal solution with 1.0% gap in 7.5 h. However, the gap of the simultaneous method is 22.5% after 50 h, leading to a profit 15.7% less than the proposed method. In all case studies, the decomposition method solves the integrated problem to optimality with 1% gap while the direct simultaneous method only returns a suboptimal solution when the integrated problem is complex.

■ ■

Set

Ij Ij* Idyn j Idyn j̅ I+s I−s

Binary Variable

Yijk equal to 1 if task i begins in unit j at time point k Greek Letter

Φd Δd ε ηiijk γiijk λiijk ωijk βis σiijk ξ Bijk BRijk

batch size of task i starting in unit j at k batch size of task i continuously executed in unit j at k BEijk batch size of task i ending in unit j at k PROFIT profit SALES sales SLk length of time slot k STsk inventory level of material s at time point k TMC total material cost TPC total processing cost Td time variable in dynamic model d Tfinal final time in dynamic model d d Trd grid time point for finite element r in dynamic model d Trq time point at collocation point q for finite element r d in dynamic model d TRjk remaining time in unit j to complete the ongoing task at time point k TSk starting time of slot k Ud time-varying input variables in dynamic model d Urq input values at collocation point q for finite element r d in dynamic model d Vd time invariant input variables in dynamic model d VPCijk processing cost of task i in unit j at time point k VPTijk processing time of task i in unit j at time point k YEijk equal to 1 if task i ends in unit j at time point k YRijk equal to 1 if task i is continuously executing in unit j at time point k Xd state variables in dynamic model d Xinitial initial state values in dynamic model d d Xrd state values at grid point Trd rq Xd state values at collocation point q for finite element r in dynamic model d Zjk equal to 1 if a task begins in unit j at time point k

AUTHOR INFORMATION

Corresponding Author

*Phone: (847) 467-2943. Fax: (847) 491-3728. E-mail: you@ northwestern.edu. Notes

The authors declare no competing financial interest.

NOMENCLATURE

Abbreviation

GBD MIDO MILP MINLP MLDO

generalized Benders decomposition mixed-integer dynamic optimization mixed-integer linear programming mixed-integer nonlinear programming mixed-logic dynamic optimization

Index

b c d i j k l, l′ m

processing cost in dynamic model d length of finite element in dynamic model d gap for stopping GBD iteration optimal value function of primal problem Lagrangean multipliers Lagrangean multipliers collection of non-complicating variables for dynamic system indexed by (i,j,k) except Tfinal ijk and Vijk coefficient of task i to state s Lagrangean multipliers collection of complicating variables except VPTijk and Bijk

Variable

APPENDIX Tables A.1−A.11 provide the parameters for the dynamic models and the scheduling model in the case study.

■

indices of all tasks which can be executed in unit j indices of non-dummy tasks which can be executed in unit j indices of tasks described as dynamic models which can be executed in unit j complement set of Idyn in Ij j indices of non-dummy tasks with ρis > 0 indices of non-dummy tasks with ρis < 0

reaction in dynamic model concentration in dynamic model dynamic model task unit time point GBD iteration raw material Q



Article

Parameter

aqq′ bmin ij bmix ij bq chot ccool cm cq caps concp demp nc nk nr pcfixed ij pcvb ij pcmax ij prp ptij ptmax ij ptmin ij s0m sh stij

■

(14) 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, 463−476. (15) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of tubular reactors: single production lines. Ind. Eng. Chem. Res. 2010, 49, 11453−11463. (16) Yue, D.; You, F. Planning and scheduling of flexible process networks under uncertainty with stochastic inventory: MINLP models and algorithm. AIChE J. 2013, 59, 1511−1532. (17) Mitra, K.; Gudi, R. D.; Patwardhan, S. C.; Sardar, G. Resiliency Issues in Integration of Scheduling and Control. Ind. Eng. Chem. Res. 2010, 49, 222−235. (18) Chu, Y.; You, F. Integration of cyclic scheduling and dynamic optimization using generalized Benders decomposition methods coupled with global mixed-integer fractional programming approaches. Comput. Chem. Eng. 2013, in press. (19) Chu, Y.; You, F. Integrated scheduling and dynamic optimization of sequential batch processes with online implementation. AIChE J. 2013, DOI: 10.1002/aic.14022. (20) Mishra, B. V.; Mayer, E.; Raisch, J.; Kienle, A. Short-term scheduling of batch processes. A comparative study of different approaches. Ind. Eng. Chem. Res. 2005, 44, 4022−4034. (21) Nie, Y. S.; Biegler, L. T.; Wassick, J. M. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 2012, 58, 3416−3432. (22) Mendez, C. A.; Cerda, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art review of optimization methods for shortterm scheduling of batch processes. Comput. Chem. Eng. 2006, 30, 913−946. (23) Reklaitis, G. V. Overview of scheduling and planning of batch process operations. Batch processing systems engineering; Springer: Berlin/Heidelberg, 1996; pp 660−705. (24) Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process. 2007, 46, 1043−1053. (25) Cuthrell, J. E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33, 1257−1270. (26) Geoffrion, A. M. Generalized benders decomposition. J. Optim. Theory Appl. 1972, 10, 237−260. (27) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal design of dynamic systems under uncertainty. AIChE J. 1996, 42, 2251−2272. (28) Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A general algorithm for short-term scheduling of batch operationsI. MILP formulation. Comput. Chem. Eng. 1993, 17, 211−227. (29) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Integration of design and control for chemical processes: A review of the literature and some recent results. Annu. Rev. Control 2009, 33, 158−171. (30) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. N. Recent advances in optimization-based simultaneous process and control design. Comput. Chem. Eng. 2004, 28, 2069−2086. (31) Floudas, C. A.; Lin, X. X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Comput. Chem. Eng. 2004, 28, 2109−2129. (32) Yue, D.; You, F. Sustainable scheduling of batch processes under economic and environmental criteria with MINLP models and algorithms. Comput. Chem. Eng. 2013, 54, 44−59. (33) Sundaramoorthy, A.; Karimi, I. A. A simpler better slot-based continuous-time formulation for short-term scheduling in multipurpose batch plants. Chem. Eng. Sci. 2005, 60, 2679−2702. (34) Bryson, A. E. Applied optimal control: optimization, estimation, and control; Taylor & Francis: Oxford, 1975. (35) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processesI. Characterization of the nominal solution. Comput. Chem. Eng. 2003, 27, 1−26. (36) Voudouris, V. T.; Grossmann, I. E. Optimal synthesis of multiproduct batch plants with cyclic scheduling and inventory considerations. Ind. Eng. Chem. Res. 1993, 32, 1962−1980.

collocation coefficient lower bound of batch size for task i in unit j upper bound of batch size for task i in unit j collocation coefficient unit cost of heating fluid unit cost of cooling fluid unit cost of raw material m collocation coefficient inventory capacity of material s concentration of product p order demand for product p number of collocation points in a finite element number of time points number of finite elements fixed processing cost of task i in unit j coefficient of variable cost of task i in unit j upper bound of processing cost of task i in unit j sale price of product p processing time of task i in unit j upper bound of processing time of task i in unit j lower bound of processing time of task i in unit j initial inventory of raw material m scheduling horizon starting time of task i in unit j

REFERENCES

(1) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (2) Grossmann, I. Enterprise-wide optimization: A new frontier in process systems engineering. AIChE J. 2005, 51, 1846−1857. (3) Wassick, J. M.; Agarwal, A.; Akiya, N.; Ferrio, J.; Bury, S.; You, F. Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products. Comput. Chem. Eng. 2012, 47, 157−169. (4) Bansal, V.; Sakizlis, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. New algorithms for mixed-integer dynamic optimization. Comput. Chem. Eng. 2003, 27, 647−668. (5) Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S. Dynamic optimization in a discontinuous world. Ind. Eng. Chem. Res. 1998, 37, 966−981. (6) Bemporad, A.; Morari, M. Control of systems integrating logic, dynamics, and constraints. Automatica 1999, 35, 407−427. (7) Oldenburg, J.; Marquardt, W.; Heinz, D.; Leineweber, D. B. Mixed-logic dynamic optimization applied to batch distillation process design. AIChE J. 2003, 49, 2900−2917. (8) Turkay, M.; Grossmann, I. E. Logic-based MINLP algorithms for the optimal synthesis of process networks. Comput. Chem. Eng. 1996, 20, 959−978. (9) Capon-Garcia, E.; Moreno-Benito, M.; Espuna, A. Improved short-term batch scheduling flexibility using variable recipes. Ind. Eng. Chem. Res. 2011, 50, 4983−4992. (10) 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. (11) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45, 6698−6712. (12) Nystrom, R. H.; Franke, R.; Harjunkoski, I.; Kroll, A. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29, 2163−2179. (13) Zhuge, J. J.; Ierapetritou, M. G. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 2012, 51, 8550−8565. R



Article

(37) Iserles, A. A first course in the numerical analysis of differential equations; Cambridge University Press: Cambridge, 2008; Vol. 44. (38) Hairer, E.; Wanner, G. Solving ordinary differential equations II: Stiff and differential-algebraic problems; Springer: Berlin, 2004; Vol. 2. (39) Benders, J. F. Partitioning procedures for solving mixedvariables programming problems. Numerische Mathematik 1962, 4, 238−252. (40) Floudas, C. A. Nonlinear and Mixed-Integer Optimization: Fundamentals and Applications; Oxford University Press: New York, 1995. (41) Sahinidis, N. V.; Grossmann, I. E. Convergence properties of generalized benders decomposition. Comput. Chem. Eng. 1991, 15, 481−491. (42) Li, X.; Chen, Y.; Barton, P. I. Nonconvex Generalized Benders Decomposition with Piecewise Convex Relaxations for Global Optimization of Integrated Process Design and Operation Problems. Ind. Eng. Chem. Res. 2012, 51, 7287−7299. (43) Luyben, W. L. Chemical reactor design and control. Environ. Eng. Manage. J. 2007, 6, 597−599.

S


Integrated Scheduling and Dynamic Optimization of Complex Batch

Recommend Documents