Discrete Time Formulation for the Integration of Scheduling and

Article pubs.acs.org/IECR

Discrete Time Formulation for the Integration of Scheduling and Dynamic Optimization Yisu Nie,† Lorenz T. Biegler,*,† Carlos M. Villa,‡ and John M. Wassick§ †

Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States The Dow Chemical Company, Freeport, Texas 77541, United States § The Dow Chemical Company, Midland, Michigan 48674, United States ‡

S Supporting Information *

ABSTRACT: We propose a model-based optimization approach for the integration of production scheduling and dynamic process operation for general continuous/batch processes. The method introduces a discrete time formulation for simultaneous optimization of scheduling and operating decisions. The process is described by the resource task network (RTN) representation coupled with detailed first-principles process dynamic models. General complications in scheduling and control can be fully represented in this modeling framework, such as customer orders, transfer policies, and requirements on product quality and process safety. The scheduling and operation layers are linked with the task history state variables in the state space RTN model. A tailored generalized Benders decomposition (GBD) algorithm is applied to efficiently solve the resulting large nonconvex mixed-integer nonlinear program by exploring the particular model structure. We apply the integrated optimization approach to a polymerization process with two parallel semibatch reactors and continuous storage and purification units. The two polymerization reactors share cooling utility from the same source, and the utility price is dependent on the consumption rate. The optimization objective is to design the process schedule and reactor control policies simultaneously to maximize the overall process profit. The case study results suggest improvements in plant profitability for the integrated approach, in contrast to the typical sequential approach, where recipes of the polymerization tasks are individually optimized but the interactions among process units are overlooked.

■

first to incorporate rigorous dynamic models into batch scheduling problems. Although the scheduling problem was restricted to flowshop plants with a special class of material transfer polices, this work suggested great potential of the integrated optimization approach to improve plant profitability in general. Recently, integration of batch processes starts to gain more attention, building on the advances on batch scheduling methods8−10 and dynamic optimization algorithms.11,12 Mishra et al.13 carried out a comparative study between two scheduling methods. The first was the standard approach with predetermined unit control recipes, while the second contained dynamic models of process units and recipes that were generated together with schedules. The latter approach results in optimization problems with significantly larger sizes but higher overall profits. This indicates potential economic benefits in embedding rigorous process dynamic models into scheduling formulations but also challenges in solving the resulting optimization problem. Research work along this line includes Romero et al.,14 Chu and You,15,16 and Capón-Garcı ́a et al.17 We have also addressed the batch integration problem in our previous work by employing the State Equipment Network to represent the switching dynamic

INTRODUCTION There are very large economic incentives for integrating decision-making functions in process industries, such as the integration of process scheduling and dynamic optimization.1 Typically, batch and semibatch processes require the determination of optimal operating recipes within individual units. These recipes then provide the data to generate schedules of the overall batch process. However, the sequential approach for dynamic optimization and scheduling problems does not lead to strong integration between them and therefore leads to suboptimal performance, as we will see later. Instead, simultaneous optimization of production schedules and unit control strategies has been shown as a promising approach to improve plant profitability in the literature. In the Process Systems Engineering community, early research efforts were mostly focused on continuous processes, especially commercial polymerization reactors. For example, Nyström et al.2 studied a grade sequencing and transition optimization problem of a polymerization process, where the transition trajectories, operating points, and sequencing of grades are determined all-at-once with a mixed integer dynamic optimization formulation. Cost minimization is achieved by reducing raw material use and off-spec products. Integrated scheduling and dynamic optimization problems for continuous processes have been also studied by Flores-Tlacuahuac and Grossmann,3 Terrazas-Moreno et al.,4 Prata et al.,5 and Zhuge and Ierapetritou.6 Integration of batch processes was less often studied in the past, partially due to its modeling and solution complexity. The work by Bhatia and Biegler7 was among the © XXXX American Chemical Society

Special Issue: Scott Fogler Festschrift Received: July 24, 2014 Revised: September 16, 2014 Accepted: September 19, 2014

A

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

Industrial & Engineering Chemistry Research

Article

behavior of batch units, and the dynamics are described with first-principles models. In the case study, we obtain more profitable production schedules together with unit control profiles. However, the computational complexity limits our ability to solve large-scale problems.18 The aforementioned studies have adopted a top-down integration strategy, where the problem is viewed as enriching the higher level scheduling optimization formulation by replacing fixed recipe parameters with detailed dynamic models. Therefore, in these formulations, the continuous time representation is favored as it is borrowed from the original scheduling models. The scheduling time horizon is divided into a number of time periods with variable lengths that are determined through optimization. In a different vein, the integration problem can also be formulated from a bottom-up perspective, which is categorized as optimal control of hybrid discrete/continuous systems in the control literature.19 The guiding methodology is model predictive control (MPC), where task recipes are determined by optimizing the continuous control variables, while production schedules are represented by the discrete control variables indicating task switching signals. For this, Gallestey et al.20 used discrete time grids in these studies as they lend themselves to model predictive control problems. Here, the time horizon consists of multiple time slots of equal lengths, and task lengths are predefined parameters in the optimization formulation as integer multiples of the unit slot duration. A related study that uses a continuous time frame is reported in de Prada et al.21 Continuous and discrete time representations are applicable to the same model instances in general, with their major pros and cons summarized in Table 1.

The discrete-time RTN model has a concise structure of a few basic equations, but it is versatile enough to accommodate a wide array of different scheduling constraints without excessive customization.23 In the RTN model, process units and materials are represented as resources that can be generated and/or consumed by operation tasks. Optimization is performed to allocate tasks in time while ensuring the resources levels stay within feasible limits. One of the main reasons for using the discrete time RTN approach is due to its versatility and consistency for different process applications. Moreover, it often provides a low integrability gap in the resulting MILP that compensates for the larger number of binaries used. This study extends this approach by integrating scheduling and dynamic optimization with the discrete time RTN. On the other hand, dynamic models bring in differential-algebraic equations (DAEs) into the integrated optimization formulation that cannot be solved directly. A widely used solution strategy is to discretize the DAEs with the simultaneous collocation method.25 The continuous time domain is replaced by many discrete finite elements with collocation points inside. Differential variables are approximated with orthogonal polynomials defined at these collocation points. The collocation method results in a large number of continuous variables in the optimization problem. The most straightforward approach to tackle the integrated scheduling and dynamic optimization problem is to formulate monolithic models, where both the scheduling and process dynamics are fully represented.26 However, great challenges exist in solving the resulting mixed-integer nonlinear programs (MINLPs). Typical solution methods for convex MINLPs include the nonlinear branch and bound,27 generalized Benders decomposition (GBD),28 outer approximation,29 extended cutting plane,30 and LP/NLP-based branch and bound.31 These methods can be applied to find good feasible solutions of the integrated formulation, though the problem is typically nonconvex and may cause difficulties in convergence. Global MINLP methods such as the branch and reduce algorithm32 guarantee to find the global optimal solution in principle, but the computation load is often too heavy for integrated problems. Solving the integrated problem directly is difficult, and thus several decomposition methods have been tried to improve the solution speed and reliability. Terrazas-Moreno et al.33 applied a Lagrangean heuristic approach for simultaneous cyclic scheduling and optimal control of multigrade polymerization reactors. The original integrated problem was separated into two subproblems for scheduling and control, respectively. The linking variables (variables that appear in both subproblems) were duplicated by adding new equality constraints, which were further dualized using the Lagrangian relaxation. The sum of the objectives of the two subproblems represented the upper bound of the original problem, and the lower bound was calculated through a heuristic procedure where the original model was solved as a nonlinear problem (NLP) by fixing the binary variables to the values obtained in the scheduling subproblem. Their case study results suggested improvements both in solution quality and speed. A generalized Benders decomposition method was studied by Chu and You16 to solve the integrated problem for batch processes of network structures. The integrated problem was decoupled to a linear master scheduling problem and dynamic optimization subproblems for individual operation tasks, which had a separable structure and can be solved in parallel. Benders cuts were derived from the optimal solutions of the subproblems and

Table 1. Pros and Cons for the Two Time Representation Methods model accuracy model size formulation tightness consistency

time-dependent event

continuous: high, continuous timing variables discrete: low, round-up errors in task lengths continuous: small, sparse time grids discrete: large, dense time grids continuous: low, big-M terms in timing constraints discrete: high, low integrality gap continuous: low, customization for different applications discrete: high, concise and consistent structure continuous: indirect, additional constraints needed discrete: direct, known time grid point positions

Many scheduling approaches have been proposed in the literature, especially those based on mixed integer linear programming (MILP) methods.9 For problems with general network structures, the state task network (STN) and resource task network35 (RTN) are widely applied as the modeling framework. Moreover, two time representations, discrete and continuous, give the two types of scheduling models. The discrete-time approach divides the scheduling time horizon into a number of uniform time intervals of fixed lengths such that the beginning and end of a discrete event lie on the interval boundaries, while the continuous-time approach introduces continuous timing variables that can place events at arbitrary time points. The Dow Chemical Company has successfully applied discrete time scheduling formulations in optimizing many of its production facilities, especially with the RTN approach.22−24 B



■

accumulated in the master problem to close the overall optimality gap. The authors highlighted the advantages of the GBD method in its solution speed and ability to deal with largescale problems. The major motivation for these decomposition methods is to separate the discrete scheduling decisions and dynamic process models and therefore lower the complexity of the resulting math programs. The synergy between scheduling and dynamic optimization is explored as solving the separated problems in certain iterative procedures. In this work, we develop a general-purpose formulation for the integration of scheduling and dynamic optimization for batch/continuous processes. The discrete time representation is adopted to enable detailed modeling of process characteristics in both the scheduling and operation layers. The resulting integrated model has a decomposable structure such that the GBD method can be used to solve the model efficiently. We also demonstrate the integration method with a case study on a polymerization process. The process consists of two semibatch polymerization reactors in parallel, a continuous purification unit, and a buffer tank in between. The integrated approach is applied to maximize the overall profit of the process within a given scheduling horizon. The paper is organized as follows: the definition of the integrated problem is given in the next section. Then the integrated formulation is discussed, including the RTN based scheduling problem and process dynamic optimization with the simultaneous collocation method. We show the GBD algorithm by defining its master and primal problems and also the iteration procedure. Next, a case study for a polymerization process is presented and the integrated approach is demonstrated and compared to the typical sequential approach. Lastly, we conclude with the major findings and future extensions of this work.

Article

MODEL FORMULATION

We first present the key equations for process scheduling based on the RTN representation and then deal with dynamic process models. We use the state space RTN model,34 which is a reformulation of the conventional discrete time RTN35,36 with extended variables and equations.23 The state space RTN model is able to handle common scheduling complications such as various transfer policies, product changovers, order fulfillments, etc. The model is also well suited for online applications as it incorporates lifted variables to record task history, and rescheduling can be performed periodically with a rolling horizon scheme.37 For process dynamics, first-principles models strengthen model fidelity, that is, they can provide accurate predictions of plant behavior over a wide operation range.38 However, computational complexity needs to be carefully managed. The counterpart of first-principles models are datadriven models, which may have narrower prediction ranges but are less expensive to develop and solve. One particular advantage of first-principles models for the integrated problem is that process operation flexibility can be explored in an extensive yet reliable manner, and therefore we focus on firstprinciples models in this work. Scheduling with RTN. As noted above, the RTN model represents process units and materials as resources that can be generated and/or consumed by operation tasks. Optimization is performed to allocate tasks in time while ensuring that resource levels stay within feasible limits. Discrete time RTN models have considerable versatility and consistency for different process applications. In this section we review the key elements of the linear state space RTN model that was developed in our earlier paper.34 A detailed derivation of this state space model along with its comparative advantages can be found in that study. One of the most important equations in any RTN models is the resource balance equation, which accounts for the variation of excess resource levels over time, as shown below:

■

PROBLEM DESCRIPTION We show a general approach for integrated production scheduling and dynamic optimization. The problem can be stated as follows: Given: Plant conf iguration ★ a set of process units and their uses ★ a set of process materials ★ production sequences for products Scheduling preknowledge ★ scheduling time horizon ★ target performance index, e.g., maximum profitability ★ transfer policies for materials ★ market price, demand of final products ★ market price of energy and utilities Process dynamics ★ dynamic models of process operations ★ process constraints regarding safety and quality issues Determine: ★ optimal production schedule ★ optimal control profiles of process operations ★ optimal performance index value

τi , m

R r ,t = R r ,t−1 +

∑ ∑ ∑ μi ,m,r ,θ wi̅ ,m,t ,θ i∈0 m∈4 θ=0 τi , m

+

∑ ∑ ∑ ∑ νi ,m,n,r ,θξi̅ ,m,n,t ,θ + Πr ,t , i∈0 m∈4 n∈5 θ=0

∀ r ∈ 9, t ∈ ;

(1)

For variable names, R is the resource level, w̅ and ξ̅ are the discrete and continuous task history states, τ is the task length, μ and ν are the corresponding task resource interaction parameters, and Π is used to represent external transfer events. For the indices, r is resource, t is time slot, i is task, m is task mode, n is task extent, and θ is a dummy index used for time shift. Two important aspects of the above resource balance equation need to be highlighted. First, the task mode index m is added to enable different operating modes or strategies for the same task. For example, different modes allow the task length to vary and product quality constraints to be set to different threshold values. This provides opportunities for dynamic optimization to have influence on the schedule design by selecting different execution modes. On the other hand, from the perspective of the RTN model (eq 1), resource transformations are still related in the same way for all modes within a given task.

This is a general description of the integrated formulation. For model constraints, other specific complications, either on the scheduling or operation considerations, may be encountered in applications. The formulation should be flexible enough to be extended to accommodate the additional complications. C



Article

time-dependent level limits. For example, it can reduce the number of tasks in modeling certain material transfer operations.23 Inequality constraints are used to ensure that the limits on the resource level are not violated for all resources during the scheduling time horizon:

Another index that is encountered less often is the task extent n. This allows for a single task to have multiple sizing variables, where an extent refers to the amount of material processed in an intermediate task step or resource interaction route. The major advantage of the task extent is to reduce the number of tasks needed to model operations that take place jointly.23,34 The second modification is the representation of the task history states, which are defined by lifting the task assignment and batch sizing variables:37

max R rmin ,t ≤ R r ,t ≤ R r ,t ,

max Vimin , m , nwi̅ , m , t , θ ≤ ξi̅ , m , n , t , θ ≤ Vi , m , nwi̅ , m , t , θ ,

(2a)

∀ i ∈ 0, m ∈ 4, n ∈ 5, t ∈ ; , 0 ≤ θ ≤ τi , m

ξi̅ , m , n , t , θ Δξi , m , n , t − θ , ∀ i ∈ 0, m ∈ 4, n ∈ 5, t ∈ ; , 0 ≤ θ ≤ τi , m

z ̇ = f (z(s), y(s), u(s)),

(3a)

(3b)

Note that the discrete task states can be relaxed as continuous variables except when θ = 0, as eq 3a ensures their integrality. In eq 3, we assume that tasks execute smoothly without any interruptions. However, it is possible to extend the state evolution equations to model disturbance events such as task delays and unit breakdowns with additional parameters.34,37 For resource balance relations, balance equations can be defined with respect to the upper and lower limits of excess resource levels, in addition to the resource balance equation (eq 1): =

+

∑ ∑ ∑ αi ,m,r ,θwi̅ ,m,t ,θ

∑ ∑ ∑ ∑ βi ,m,n,r ,θ ξi̅ ,m,n,t ,θ ,

K

z(s) = z j0 + hj ∑ Ωk(τ )zj̇ , k ,

∀ r ∈ 9, t ∈ ;

k=1

i∈0 m∈4 n∈5 θ=0

τi , m

∑ ∑ ∑ αi′,m,r ,θwi̅ ,m,t ,θ i∈0 m∈4 θ=0 τi , m

+

∑ ∑ ∑ ∑ βi′,m,n,r ,θ ξi̅ ,m,n,t ,θ ,

∀ r ∈ 9, t ∈ ;

i∈0 m∈4 n∈5 θ=0

(4b) max

j∈1 (9)

Here, j ∈ 1 is the index of finite elements with 1 = {1, 2, ...J }, and k ∈{1,2,...,K} is the index for collocation points, hj is the length of finite element j, τ ∈ [0,1] is the normalized time coordinate within an element and s = s0j + hj τ (s0j denotes the starting time of element j), z0j is the value of the differential state at the beginning of element j, żj,k is the value of the first derivatives at collocation points, and Ωk is a polynomial of τ of order K, defined as

(4a) min R rmin ,t = R r ,t−1 +

(8)

The simultaneous collocation method is used to convert the DAE system and constraints from the continuous form to discretized equations.39 The collocation equation for the differential of the state variables in finite element j can be stated as below:

i∈0 m∈4 θ=0 τi , m

+

(7b)

h(z(s), y(s), u(s)) ≤ 0

τi , m

R rmax ,t−1

(7a)

The differential equations are defined for the differential states z in semiexplicit form, and the algebraic states are y. We assume that the DAE system is index 1. Under this assumption, the initial conditions of the differential states are given by z0, while the initial conditions for the algebraic states are obtained from eq 7b. Process control variables are u. We use s to denote the continuous time coordinate in contrast to t, which is used for the discrete time grid. For dynamic optimization, additional constraints on the state variables are enforced to deal with process requirements on product quality and process safety. Also, control variables are bounded within acceptable ranges. Without loss of generality, these constraints can be expressed as

ξi̅ , m , n , t , θ = ξi̅ , m , n , t − 1, θ − 1 ,

R rmax ,t

z(0) = z 0

g (z(s), y(s), u(s)) = 0

wi̅ , m , t , θ = wi̅ , m , t − 1, θ − 1 ,

∀ i ∈ 0, m ∈ 4, n ∈ 5, t ∈ ; , 1 ≤ θ ≤ τi , m

(6)

Operation Optimization with Dynamic Models. Mechanism models for dynamic systems are represented as DAEs in the continuous time domain, often developed from the firstprinciples such as mass/energy balances and thermodynamics. The input control variables provide degrees of freedom for operation optimization. A DAE system can be written in a generic form as follows:

(2b)

The task assignment variables w and ξ are variables used in the conventional RTN models to decide the occurrence and size of a task and thus only effective when a task starts. In contrast, the lifted states w̅ and ξ̅ have an extra dimension of the task length, so that they are able to record the entire history of a task. For example, w̅ i,m,t,θ = 1 denotes an active task i of mode m at time t that has been executed for θ time periods. Also note that w̅ i,m,t,θ and wi,m,t are equivalent when θ = 0. The change of the task history states over time is described by the state evolution equations:

∀ i ∈ 0, m ∈ 4, t ∈ ; , 1 ≤ θ ≤ τi , m

(5)

Finally, the task extent limit constraint is incorporated to specify the upper and lower limits of task extent sizes:

∀ i ∈ 0, m ∈ 4, t ∈ ; , 0 ≤ θ ≤ τi , m

wi̅ , m , t , θ Δwi , m , t − θ ,

∀ r ∈ 9, t ∈ ;

min

Here, R and R are the upper and lower resource limits, respectively. The interaction parameters in eq 4 (α, β, α′, β′) are defined in analogy to the parameters (μ, ν) used in the resource balance equation. The resource limit balance equations are able to efficiently deal with resources that have

Ωk =

∫0

τ

Sk(τ′) dτ′,

k∈2

(10)

Here, Sk is a basis function for a Lagrange interpolation polynomial, given by D



Article

Figure 1. Time representation of the integrated formulation. K

∏

Sk(τ ) =

k ′= 1, ≠ k

τ − τk ′ , τk − τk ′

Task switching events are described by the task assignment/ history variables in the RTN scheduling model. In the lower layer dynamic optimization problem, process dynamic models are influenced by switching tasks. This may correspond to changes in the model parameters or structure as well as different specifications on the process constraints. Therefore, the dynamic system and constraints are also functions of the task history states in the integrated model. In a time slot t, we obtain

k∈2 (11)

In addition, the continuity condition across finite element boundaries is preserved: K

z j0+ 1 = z j0 + hj ∑ Ωk(1)zj̇ , k , k=1

1≤j≤J−1 (12)

K

z f = zJ0 + hJ ∑ Ωk(1)zJ̇ , k k=1

zṫ = f (zt , yt , ut , wi̅ , m , t , θ , ξi̅ , m , n , t , θ), (12b)

The algebraic states y can be treated similarly with Lagrange polynomials: K

y(s) =

∑ Sk(τ)yj ,k , k=1

j∈1 (13)

zt0 = A(zt f− 1 , wi̅ , m , t , θ , ξi̅ , m , n , t , θ)

(14a)

g (zt , yt , ut , wi̅ , m , t , θ , ξi̅ , m , n , t , θ) = 0

(14b)

h(zt , yt , ut , wi̅ , m , t , θ , ξi̅ , m , n , t , θ) ≤ 0

(14c)

Equations 14a are left in the continuous form for notational convenience, and the discretized version can be obtained through the simultaneous collocation method. The states and controls are defined with respect to the entire process such that they do not have indices for task or task mode dependency. Moreover, the initial condition of the differential states are no longer known parameters but a function of their values at the end of the previous slot (zft−1) and task assignment variables (w̅ i,m,t,θ, ξ̅i,m,n,t,θ). If no task switch occurs, z0t and zft−1 are equal to ensure the continuity of the differential states; otherwise, sudden jumps are allowed at the beginning and end of specific operations. This discrete switching behavior is described by A(·), which is a balance equation of zt across the discrete time slot boundaries, similar to the resource balance equation (eq 1). For the objective function, profit maximization is a common choice for the integrated scheduling and dynamic optimization problem, since it accounts for the trade-offs between product sales and manufacturing costs in a natural way. Product sales are decided by the scheduling level decisions such as batch numbers, while manufacturing costs are calculated using the detailed dynamic models. To summarize, the integrated formulation for scheduling and dynamic optimization can be described as follows:

Because the DAE is index 1, the continuity condition is not needed for algebraic states. Control variables also need to be parametrized, where piecewise constant and linear profiles are often assumed in a given finite element. The algebraic equations in the DAE system and process constraints are only enforced at collocation points, using the above substitutions for z(s), y(s), and u(s) in eq 7b. Integrated Formulation. The integrated formulation ties the RTN scheduling model and process dynamic models together. First, time representation needs to be considered, as a means to coordinate scheduling time slots, finite elements, and collocation points. A unified time grid is defined as shown in Figure 1. In this example, the entire scheduling time horizon, indexed by t ∈ {1,2,...,H} consists of H discrete time slots of unit length. A time slot contains two finite elements with three collocation points inside each element. The arrangement is adjustable for specific applications, e.g., the number of finite elements in a time slot and the type and number of collocation points used. Note that the starting point of a scheduling time slot is aligned with a finite element. By this means, task switching events can be simultaneously addressed in both the scheduling and dynamic control equations.

Here, the overall objective function Φ(·) is written as a linear combination of two subobjectives Φsch(·) and Φdyn(·), E



Article

An alternative approach that constructs a more concise Lagrange function is to introduce duplicated complicating variables (ξ̅d, w̅ d) in eq 16:

depending on the scheduling and dynamic optimization variables, respectively. Also, the constraints can be separated into scheduling and dynamic optimization groups.

■

p

min Φdyn(zt , yt , ut , wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) u

SOLUTION STRATEGY The integrated formulation is posed as a nonconvex MINLP problem, and the direct solution of which is time-consuming or even computationally intractable for large-scale problems. However, the model structure (eq 15) can be efficiently explored with decomposition algorithms. The scheduling equations and dynamic models are only linked by the task history state variables, and the remaining parts are isolated from each other. The GBD method is well suited for this type of problem. In this method, the original model is decoupled to a primal problem and a master problem. In the primal problem, the linking variables (also termed as complicating variables) are temporarily fixed to some value such that the remaining parametrized problem has a more tractable size and structure. The master problem is a projection of the original problem to a restricted variable space, i.e., the variables in the primal subproblem (except for the complicating variables w̅ and ξ̅) are not explicitly included in the master problem. Alternatively, the master problem includes cutting planes which are constructed from the solution and dual information on the primal problem. The optimal solution of the original problem is bounded by the best primal and master problem solutions. In this section we discuss the main components of the GBD method in the context of integrated scheduling and dynamic optimization. Primal Problem. The primal problem for the integrated formulation is the dynamic optimization problem (eq 14) with fixed task history states. Because all discrete variables are fixed, the problem is reduced to a continuous NLP:

t

p

s . t . F dyn(zt , yt , ut , wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) = 0 p

Gdyn(zt , yt , ut , wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) ≤ 0 wi̅ ,pm , t , θ = wi̅ d, m , t , θ p

d

ξ i̅ , m , n , t , θ = ξ i̅ , m , n , t , θ

(18)

To derive the Lagrange function, we obtain d

3(wi̅ d, m , t , θ , ξ i̅ , m , n , t , θ) * p = Φdyn + λ*TF dyn(zt*, yt* , ut*, wi̅ ,pn , t , θ , ξ i̅ , m , n , t , θ) p + γ *TGdyn(zt*, yt* , ut*, wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) p d + σ *T(wi̅ ,pm , t , θ − wi̅ d, m , t , θ) + ω*T(ξ i̅ , m , n , t , θ − ξ i̅ , m , n , t , θ) * = Φdyn + σ *T(wi̅ ,pm , t , θ − wi̅ d, m , t , θ) p d + ω*T(ξ i̅ , m , n , t , θ − ξ i̅ , m , n , t , θ)

(19)

Here, σ and ω are the corresponding Lagrange multipliers for the duplication equations. The simplification is due to the fact that the discretized dynamic equations are not a function of the duplicated variables, and

p

min Φdyn(zt , yt , ut , wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) u

p λ*TF dyn(zt*, yt* , ut*, wi̅ ,pn , t , θ , ξ i̅ , m , n , t , θ) = 0

(20a)

p γ *TGdyn(zt*, yt* , ut*, wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) = 0

(20b)

t

p

s . t . F dyn(zt , yt , ut , wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) = 0 G

dyn

(zt , yt , ut , wi̅ ,pm , t , θ ,

p ξ i̅ , m , n , t , θ)

≤0

at the primal optimum. This concise Lagrange function is equivalent to the one in eq 17 and computationally more favorable.40 The Lagrange function is used to form the cutting plane constraints (a.k.a. Benders cuts) in the master problem. Master Problem. The master problem contains the scheduling equations used in the linear state space RTN model (eqs 1, 3, 4, 5, and 6), as well as the cutting plane constraints, which can be stated as

(16)

The superscript p of w̅ and ξ̅ indicates that they are fixed parameters. The objective function considers only the manufacturing cost Φ dyn (·), and it is rewritten as a minimization problem following the convention of the GBD method. Assuming the primal problem is solved to optimality by NLP algorithms with given history state values, the Lagrange function 3(·) can be constructed as follows:

min −Φsch(R r , t , wi̅ , m , t , θ , ξi̅ , m , n , t , θ) + η

*

3(wi̅ , m , t , θ , ξi̅ , m , n , t , θ) = Φdyn

p + λ*TF dyn(zt*, yt* , ut*, wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ)

max s . t . F sch(wi̅ , m , t , θ , ξi̅ , m , n , t , θ , R r , t , R rmin ,t , R r ,t ) = 0

p + γ *TGdyn(zt*, yt* , ut*, wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ);

max Gsch(wi̅ , m , t , θ , ξi̅ , m , n , t , θ , R r , t , R rmin ,t , R r ,t ) ≤ 0

p 0 ≤ γ *⊥Gdyn(zt*, yt* , ut*, wi̅ ,pm , t , θ , ξ i̅ , m , n , t , θ) ≥ 0

3(wi̅ , m , t , θ , ξi̅ , m , n , t , θ) ≤ η

(17)

Here, λ and γ are the vectors of optimal Lagrange multipliers for Fdyn(·) and Gdyn(·) obtained at the optimal solution, and the optimal objective and variable values are noted with *. The complementarity condition of the inequalities is satisfied between γ* and Gdyn(z*t ,y*t ,u*t ,w̅ pi,m,t,θ,ξ̅ip,m,n,t,θ). The structure of the Lagrange function can be complex as the number of equations in Fdyn(·) and Gdyn(·) are often large after collocation.

(21a)

(21b)

There is an infinite number of cutting planes that can be constructed by the Lagrange function. A typical strategy is to replace the master problem with its relaxation, where only a finite collection of the cutting planes are included: min −Φsch(R r , t , wi̅ , m , t , θ , ξi̅ , m , n , t , θ) + η F

(22a)



Article

max s . t . F sch(wi̅ , m , t , θ , ξi̅ , m , n , t , θ , R r , t , R rmin ,t , R r ,t ) = 0 max Gsch(wi̅ , m , t , θ , ξi̅ , m , n , t , θ , R r , t , R rmin ,t , R r ,t ) ≤ 0

3 b(wi̅ , m , t , θ , ξi̅ , m , n , t , θ) ≤ η , b ∈ )

(22b)

Equation 22 is an MILP and the set of valid cutting planes are denoted by b ∈ ) . The cuts can be collected at different primal optimal solutions. GBD Algorithm. The GBD method applies an iterative procedure to obtain the optimal solution of the original problem, where the primal and relaxed master problems are solved in a loop, as depicted in Figure 2. The algorithm often

Figure 3. Process flowsheet of the polymerization process.

details on the reactor model equations and process RTN representation, please refer to the Supporting Information. The polymerization reactors demand large amounts of cooling capacity as the reactions are highly exothermic. An interesting feature of the plant is that the two reactors share cooling utility from the same source, and the utility price p is a function of its consumption rate r: ⎧1, 0 ≤ r ≤ r0 ⎪ p=⎨ 1 ⎪1 + (r − r0), r0 ≤ r ≤ 3r0 r ⎩ 0

Figure 2. GBD iteration scheme.

starts with solving the primal problem with the complicating variables fixed to an initial point, and the master problem is then constructed with the obtained Benders cut at the primal optimum. The master problem is solved to update the complicating variables and triggers the next GBD iteration. The primal solution gives the upper bound of the optimal objective value, while the relaxed master problem calculates the lower bound. As Benders cuts are accumulated in the master problem, the lower bound is nondecreasing through the iterations. The algorithm converges when the two bounds fall in a close neighborhood specified by optimality tolerances.41 It is worth noting that global optimality is only guaranteed when the primal problem is convex. This is usually not satisfied by the dynamic optimization problem. Also, the following integer cuts are added to the relaxed master problem to avoid cyclic integer assignments:

The utility price equals one money unit (MU) when the total utility consumption rate of the reactors is below a basis value r0, and the price rises linearly beyond r0 until the maximum rate 3r0, where the price reaches 3 MU. Other complexities come from both the scheduling and dynamic optimization decisions. For scheduling considerations, product A is made-to-order while products B and C are made-to-stock. The optimized schedule is required to produce the exact amount of product A according to customer orders, while extra products B and C can be made in addition to ordered amounts. Products A and C belong to the same product family, and product B is an orphan product. Equipment cleaning is required for PU during the transition between different product families. Also, PU needs to carry out of f line maintenance after processing a certain amount of raw polymers. In other words, the processing capacity is consumable, and the maintenance can replenish it after it drops below a threshold value. Finally, mixing of different products is not allowed in any units. For the reactors, a 1 h setup time is needed before starting any polymerization operations. The main constraints considered are the number-average MW, byproduct ratio, unreacted monomer concentration, and maximum cooling duty. Other minor constraints are also included, such as limits on the variation of reactor temperature. The decision variables are the task assignments for scheduling and the control strategy of the reactors that include the temperature trajectory and the monomer feed policy. There are two modes for the polymerization of products B and C with different task lengths but the same quality specs. Product A has only one polymerization mode. PU has a maximum processing rate that is product dependent. The goal of the integrated optimization is to maximize the plant profit within a 3-day scheduling time horizon, which is calculated by the product sales minus the total manufacturing cost. The cost consists of the cooling utility cost, fixed cost for polymerization operations and PU maintenance, as well as penalties for monomer flow rate and temperature fluctuations. We solve the integrated optimization problem with the GBD decomposition approach. The primal problem is defined as dynamic optimization of the polymerization operations in the

∑ wi̅ b, m, t ,0 − ∑ wi̅ b, m, t ,0 ≤ |>1b| − 1 >1b

(24)

> b0

{wi̅ b, m , t ,0 ∈ >1b|wi̅ b, m , t ,0 = 1}, {wi̅ b, m , t ,0 ∈ > b0|wi̅ b, m , t ,0 = 0} (23)

The cuts are defined with the binary task history state variables w̅ bi,n,t,0. An integer cut obtained at iteration b is added to the next master problem at b + 1.

■

CASE STUDY We illustrate the discrete time integrated formulation with a polymerization process. The process flowsheet is shown in Figure 3. The plant has two parallel polymerization reactors (Rxr1, Rxr2) of identical processing capacity, and they both connect to the downstream buffer tank (T) followed by the purification unit (PU). Raw polymer products are made in the reactors, and PU removes the catalyst in the raw products to give final products. The polymerization reactors are semibatch units, and both the tank and PU are continuous units. Three types of products are made from the process (A, B, C) with different specs on molecular weight (MW) and byproduct ratio. Rigorous dynamic models are developed for the polymerization reactors,46,47 while the other continuous units are modeled by the linear resource balance in the RTN representation. For G



Article

Longer tasks require more finite elements, as a 1 h finite element is used in this case study. The scheduling problems are solved within a 1200 CPUs limit and the optimality gap ∼5%. The MIQP problem is difficult to solve not only because of its large size but also due to the degenerate nature of the scheduling problems, where many different schedules (typically due to different task times) give the same objective value.45 Note that the number of continuous variables in these RTN scheduling problems is significantly larger than with previous RTN models. This is due to the inclusion of the lifted task state variables. For the integrated formulation, the GBD algorithm is able to solve the problem with three iterations, and model and solution statistics are shown in Table 3. The algorithm starts with the socalled zeroth primal problem, which is the first dynamic optimization problem with given task states. These initial values of the task states are obtained at the optimal solution of the sequential scheduling problem (MIQP). Then the GBD algorithm iterates between the master MILP problem and the primal NLP problem. The major percentage of CPU time is spent on the primal problems, given that NLPs are large. Solution of these NLPs benefits from a systematic initialization of the dynamic optimization models, which solves the collocation equations in an element-by-element manner. This ensures that primal NLPs can be solved in an efficient and robust way. Moreover, fixed recipes can be used to carry out simulation once the schedule is determined, and this offers good starting points for the primal NLPs. The master problems are solved to zero optimality gap and the solution speed is fast, due to the inclusion of the Benders cuts that helps to eliminate the degeneracy of the scheduling problem. This is in contrast to the MIQP problem in the sequential approach, although both models have similar sizes. The optimizer is able to evaluate schedules with different timings with the assistance of the multiplier values for the task history states (σ and ω in eq 19). As a maximization problem, the upper and lower bounds of the optimal objective are given by the primal and master problems, respectively. The optimality gap is reduced to zero after the third master problem. The maximal plant profit is 92.26 MU for the integrated optimization approach, which is 10.5% higher than that of the sequential approach. The optimized schedules are shown in Figure 4 in comparison. Task rectangles are colored to denote different products. Also, different polymerization task modes have different fill patterns. For the buffer tank, the upper small rectangles represent the material transfer operations that load the reactor products to the buffer tank, while the lower rectangles represent the outgoing flows from the tank to PU. In addition, the reactors can be used as temporary storage units to hold batches when the buffer tank is not immediately available.

two reactors to minimize the manufacturing cost, and the master problem deals with the complete process as a linear scheduling problem without the dynamic reactor models. The complicating variables are the task history states associated with the polymerization reactors. We use 1 h scheduling time slots, and only one finite element is included for a time slot. Three Radau collocation points25 are employed in a finite element. The model is formulated in the General Algebraic Modeling System (GAMS),42 and we use Gurobi43 to solve the master MILP problems and Conopt44 to solve the primal NLP problems. All computation work is performed with a laptop with 2.4 GHz CPU and 8 GB memory, running Mac OS X 10.9.4. We compare the integrated optimization approach with the sequential approach, where polymerization recipes are optimized individually and the schedule is determined for these fixed recipes. This requires solving separate dynamic optimization problems for recipe optimization (as described in previous studies46,47) and a discrete optimization problem for schedule design. The scheduling problem is a mixed-integer quadratic program (MIQP) because it includes a quadratic term for the cooling utility cost in its objective function, while all model constraints are linear. The model and solution statistics for the sequential approach is given in Table 2. Five recipe Table 2. Model and Solution Statistics for the Case Study: The Sequential Approach recipe optimization model solution product (mode)

type

A B(1) B(2) C(1) C(2)

NLP NLP NLP NLP NLP

Obj. (MU)

time (s)

5.04 3.3 7.46 3.7 6.99 2.5 10.86 4.0 8.76 2.4 schedule optimization

model solution

model size Var.

Cons. #

1200 832 924 648 740

1239 859 954 669 764

model size

type

Obj. (MU)

time (s)

var. no. (/dis.)

cons. no.

MIQP

82.56

1200.0

30779(1825)

30823

optimization problems are solved for the five polymerization modes (see Table S1 in the Supporting Information for details), where the objective function is the manufacturing cost. All five cases result in small scale NLPs that are solved within seconds. The reactor model has the same structure for the five cases (see the Supporting Information), but model parameters such as the initial charge condition, kinetic constants, and product specs may vary. The difference in model sizes is mainly due to the number of finite elements used in each problem.

Table 3. Model and Solution Statistics for the Case Study: The Integrated Approach model solution

model size

GBD algorithm

iteration

type

Obj. (MU)

time (s)

Var. no. (/dis.)

Cons. no.

zeroth primal first master first primal second master second primal third master

NLP MILP NLP MILP NLP MILP

68.17 120.94 72.23 104.22 79.10 92.26

53.6 4.8 199.2 10.4 120.4 4.5

14820 30561(1825) 14820 30561(1825) 14820 30561(1825)

16570 30666 16570 30668 16570 30670

H

upper

lower

gap (%)

120.94 120.94 104.22 104.22 92.26

85.16 85.16 89.04 89.04 92.26 92.26

29.58 26.38 14.57 11.48 0.00



Article

Figure 4. Optimal production schedules.

Figure 5. Operating profiles for the downstream units.

In fact, polymerization continues in the interim, because the catalyst is still effective as the reactor temperature is maintained within the allowable range. PU processes different products, carries out transition cleaning, and turns offline for main-

tenance to replenish its processing capacity. In the schedule of the sequential approach (Figure 4a), two batches of product A, four batches of product B, and three batches of product C are produced. Also, different task modes are observed for products I



Article

Figure 6. Comparison of the optimal control recipes.

B and C. One additional batch of product C is observed in the schedule of the integrated approach (Figure 4b). This is in contrast to the sequential approach, where running two batches of product C in parallel is undesirable as it leads to very expensive utility prices (see Figure 9). This operating flexibility is explored more efficiently in the integrated approach, such that the utility price can be brought down to a level where the additional batch is profitable. The main disadvantage of the sequential method is that it overlooks the interaction between the two reactors. More specifically, dynamic optimization for individual recipes is performed with respect to a single reactor,

while the actual utility price is determined by the summation of the utility consumption rates of the two reactors. However, this synergy cannot be considered by the sequential approach, as the dynamic optimization is determined independently of the production schedule. The inventory profiles of the intermediate tank are shown in Figure 5a, and the process rates of PU are shown in Figure 5b. To show the details and to compare reactor operations, we construct a “base case”, which uses the optimal reactor recipes from the sequential approach (Figure 6) and superimposes these recipe profiles on the schedule determined from the J



Article

Figure 7. Process constraint on the maximum cooling rate.

Figure 8. Process constraint on the product quality indices.

Figure 9. Comparison of the cooling utility consumption and cost.

integrated approach (Figure 4b). This allows a direct comparison of the operating profiles, tank inventories, product quality indices, and cooling constraints in Figures 5−9 over the entire time horizon. The optimal recipe profiles from the integrated formulation are depicted by red solid lines, and the blue dotted lines sketch the sequential (base case) recipe

profiles. The plotted temperature data are in Celsius and calibrated with respect to a basis level value Tb; and the monomer flow rates are scaled. Since the variation of reactor temperature is penalized, gradual transitions in the temperature curves are observed from one batch to another. The monomer is allowed to enter the reactors 1 h after a polymerization task K



Article

starts, due to the setup requirements. It is worth noting that the feed profiles tend to minimize the overlap between the two reactors, R1 and R2. Particularly for the integrated solution, we obtain flow profiles where the sum of the two reactor feeds is significantly less than in the base case. This results from the fact that the feed rates are directly related to the heat generation rates in the reactors, and therefore, smaller total feed rates lead to lower cooling duty and lower utility prices. The maximum cooling capacity and product quality constraints (specs on the byproduct ratio, number-average MW, and unreacted monomer concentrations) are important constraints that should be satisfied for all polymerization batches (refer to eqs 34−36 in the Supporting Information for the definition of the specs). We show these constraint profiles for Rxr1. The utility consumption profiles are plotted together with the maximum cooling capacity limits for both cases in Figure 7. We measure the rate with r0 as the base unit. The cooling capacity limit depends on the temperature of the reactors so that it varies in time (eq 37 in the Supporting Information). The constraint on cooling capacity is active for most of the operation periods. In Figure 8, the byproduct ratio should stay below its spec limits, while the polymer MW should be above its spec at the end of polymerization. In the plots, the target spec values are marked with crosses. The byproduct ratio and MW have product dependent target values. The unreacted monomer concentration also reaches desired levels for all polymerization runs. These quality constraints are only effective at the grid points when raw polymers are transferred to the buffer tank and otherwise relaxed by using big M terms in the model formulation (A(·) in eq 14a). In Figure 9, the cooling utility price, total consumption, and total cost are plotted over time. In the price plot, high price periods occur in the first half of the time horizon, especially at the beginning when polymerizing product C. The optimized recipe is able to reduce the peak values of the high price segments and also their occurrence frequency. The middle graph shows the accumulated total utility consumption of the two reactors. The amounts of consumption are equal for the two cases at the final time, as the same amount of products is made. However, we observe a lower cost for the integrated solution in the utility cost curve because it takes full advantage of the utility price structure.

has advantages over the sequential approach, where the control recipes are determined separately. The method is currently offline, and online implementation is in the scope of our future work. Here the structure of the integrated scheduling problem is closely aligned to the models used in MPC, so that a moving horizon framework can be considered. Also the model needs to incorporate terms to represent process disturbances. For the example process used in our case study, development of the dynamic models of PU is also an interesting direction that allows exploration of interactions between the upstream and downstream units. Additional topics for future research also include the following. First, reactive scheduling problems can be addressed as direct extension of the formulations given in this study. One advantage is that the state-space RTN formulation with lifted variables for the task history states (see eqs 2 and 3) allows for the interruption and resumption of operating schedules (say, due to an upset or outage). Second, further integration to include planning also needs to be considered in future work. Potentially useful strategies include implementation of largescale MINLP decomposition strategies problems over longer horizons. Finally, there are a number of sources of uncertainty, including unmeasured process disturbances such as uncertain inputs (due to feed composition, ambient conditions, and impurities) as well as model mismatch, which were not considered in this study. Handling these requires a more complex optimization framework that remains an open research question. Potentially useful strategies include multiscenario optimization and the introduction of back-off constraints.

■

ASSOCIATED CONTENT

S Supporting Information *

Detailed descriptions of the RTN representation of the polymerization process, including the polymerization reactors, buffer tank, purification unit, and optimization formulation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

■

The authors declare no competing financial interest.

■

CONCLUSIONS We develop a general framework for the integration of production scheduling and dynamic optimization decisionmaking. The method combines the state space RTN representation with rigorous process dynamic models to formulate the integration problem as math programs. The scheduling and control layers are linked through the lifted task history state variables, which provide opportunities to apply the GBD decomposition method for an efficient problem solution. The GBD algorithm decouples the problem into the MILP master scheduling problem and continuous primal dynamic optimization problem. The two problems are solved in a loop, and Benders cuts are derived to close the optimality gap between the two. We have applied our method to a polymerization process with networked structure and multiple polymer products. The polymerization reactors are described by first-principles models. The objective is to maximize the plant profitability. The cooling utility cost brings in interactions between the two reactors, and we show the integrated approach

REFERENCES

(1) Shobrys, D. E.; White, D. C. Planning, scheduling and control systems: why cannot they work together. Comput. Chem. Eng. 2002, 26, 149−160. (2) Nyström, R.; Franke, R.; Harjunkoski, I.; Kroll, A. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29, 2163−2179. (3) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45, 6698−6712. (4) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and optimal control of polymerization reactors. AIChE J. 2007, 53, 2301−2315. (5) 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. (6) Zhuge, J.; Ierapetritou, M. G. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 2012, 51, 8550−8565.

L



Article

(7) Bhatia, T.; Biegler, L. T. Dynamic optimization in the design and scheduling of multiproduct batch plants. Ind. Eng. Chem. Res. 1996, 35, 2234−2246. (8) Méndez, C. A.; Cerdá, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art review of optimization methods for short-term scheduling of batch processes. Comput. Chem. Eng. 2006, 30, 913−946. (9) Maravelias, C. T. General framework and modeling approach classification for chemical production scheduling. AIChE J. 2012, 58, 1812−1828. (10) Harjunkoski, I.; Maravelias, C. T.; Bongers, P.; Castro, P. M.; Engell, S.; Grossmann, I. E.; Hooker, J.; Méndez, C.; Sand, G.; Wassick, J. Scope for industrial applications of production scheduling models and solution methods. Comput. Chem. Eng. 2014, 62, 161−193. (11) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processes I. Characterization of the nominal solution. Comput. Chem. Eng. 2003, 27, 1−26. (12) Kameswaran, S.; Biegler, L. T. Simultaneous dynamic optimization strategies: Recent advances and challenges. Comput. Chem. Eng. 2006, 30, 1560−1575. (13) 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. (14) Romero, J.; Espuña, A.; Friedler, F.; Puigjaner, L. A new framework for batch process optimization using the flexible recipe. Ind. Eng. Chem. Res. 2003, 42, 370−379. (15) Chu, Y.; You, F. Integrated scheduling and dynamic optimization of sequential batch processes with online implementation. AIChE J. 2013, 59, 2379−2406. (16) 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, 7867−7885. (17) Capón-Garcı ́a, E.; Guillén-Gosálbez, G.; Espuña, A. Integrating process dynamics within batch process scheduling via mixed-integer dynamic optimization. Chem. Eng. Sci. 2013, 102, 139−150. (18) Nie, Y.; Biegler, L. T.; Wassick, J. M. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 2012, 58, 3416−3432. (19) Bemporad, A.; Morari, M. Control of systems integrating logic, dynamics, and constraints. Automatica 1999, 35, 407−427. (20) Gallestey, E.; Stothert, A.; Castagnoli, D.; Ferrari-Trecate, G.; Morari, M. Using model predictive control and hybrid systems for optimal scheduling of industrial processes. Automatisierungstechnik 2003, 51, 285−293. (21) de Prada, C.; Grossmann, I. E.; Sarabia, D.; Cristea, S. A strategy for predictive control of a mixed continuous batch process. J. Process Control 2009, 19, 123−137. (22) Wassick, J. M. Enterprise-wide optimization in an integrated chemical complex. Comput. Chem. Eng. 2009, 33, 1950−1963 FOCAPO 2008, Selected Papers from the Fifth International Conference on Foundations of Computer-Aided Process Operations. (23) Wassick, J. M.; Ferrio, J. Extending the resource task network for industrial applications. Comput. Chem. Eng. 2011, 35, 2124−2140. (24) 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. (25) Biegler, L. T. Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes; Society for Industrial and Applied Mathematics: Philadelphia, PA, 2010. (26) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133 FOCAPO 2012. (27) Gupta, O. K.; Ravindran, A. Branch and bound experiments in convex nonlinear integer programming. Manage. Sci. 1985, 1533− 1546. (28) Geoffrion, A. M. Generalized Benders decomposition. J. Optim. Theory Applications 1972, 10, 237−260.

(29) Duran, M. A.; Grossmann, I. E. An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math Program. 1986, 36, 307−339. (30) Westerlund, T.; Pettersson, F. An extended cutting plane method for solving convex MINLP problems. Comput. Chem. Eng. 1995, 19, 131−136. (31) Quesada, I.; Grossmann, I. E. An LP/NLP based branch and bound algorithm for convex MINLP optimization problems. Comput. Chem. Eng. 1992, 16, 937−947. (32) Ryoo, H. S.; Sahinidis, N. V. A branch-and-reduce approach to global optimization. J. Global Optim. 1996, 8, 107−138. (33) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Lagrangean heuristic for the scheduling and control of polymerization reactors. AIChE J. 2008, 54, 163−182. (34) Nie, Y.; Biegler, L. T.; Wassick, J. M.; Villa, C. M. Extended discrete-time resource task network formulation for the reactive scheduling of a mixed batch/continuous process. Ind. Eng. Chem. Res. 2014, DOI: 10.1021/ie500363p. (35) Pantelides, C. C. Unified frameworks for optimal process planning and scheduling. In Proceedings on the Second Conference on Foundations of Computer Aided Operations, Crested Butte, Colorado, July 18−23, 1993; pp 253−274. (36) Schilling, G.; Pantelides, C. C. Optimal periodic scheduling of multipurpose plants. Comput. Chem. Eng. 1999, 23, 635−655. (37) Subramanian, K.; Maravelias, C. T.; Rawlings, J. B. A state-space model for chemical production scheduling. Comput. Chem. Eng. 2012, 47, 97−110. (38) Pantelides, C. C.; Renfro, J. G. The online use of first-principles models in process operations: review, current status and future needs. Comput. Chem. Eng. 2012, 51, 136−148. (39) Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process.: Process Intensif. 2007, 46, 1043−1053. (40) 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. (41) Sahinidis, N. V.; Grossmann, I. E. Convergence properties of generalized benders decomposition. Comput. Chem. Eng. 1991, 15, 481−491. (42) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R.; Rosenthal, R. GAMS A User’s Guide; GAMS Development Corporation: Washington DC, 2006. (43) Gurobi. Gurobi Optimizer Reference Manual. Gurobi Optimization, 2012. (44) Drud, A. S. CONOPT−a large-scale GRG code. INFORMS J. Comput. 1994, 6, 207. (45) Velez, S.; Maravelias, C. T. Mixed-integer programming model and tightening methods for scheduling in general chemical production environments. Ind. Eng. Chem. Res. 2013, 52, 3407−3423. (46) Nie, Y.; Biegler, L. T.; Villa, C. M.; Wassick, J. M. Reactor modeling and recipe optimization of polyether polyol processes: Polypropylene glycol. AIChE J. 2013, 59, 2515−2529. (47) Nie, Y.; Biegler, L. T.; Villa, C. M.; Wassick, J. M. Reactor modeling and recipe optimization of ring-opening polymerization: Block copolymers. Ind. Eng. Chem. Res. 2014, 53, 7434−7446.

M


Discrete Time Formulation for the Integration of Scheduling and

Recommend Documents