Novel Optimization Model and Efficient Solution Method for

Article pubs.acs.org/IECR

Novel Optimization Model and Efficient Solution Method for Integrating Dynamic Optimization with Process Operations of Continuous Manufacturing Processes Hanyu Shi, Yunfei Chu, and Fengqi You* Department of Chemical and Biological Engineering, Northwestern University, Evanston, Illinois 60208, United States ABSTRACT: In this paper, we propose a novel mathematical model for the integrated optimization for production planning, scheduling, and dynamic optimization of continuous manufacturing processes. With a novel approach to estimate the inventory cost, the integrated problem is first formulated as a mixed-integer dynamic optimization (MIDO) problem, which is then reformulated into a large scale mixed-integer nonlinear program (MINLP). By using metamodeling to characterize the detailed linking information between different decision layers, we develop an efficient solution method to decompose the integrated MINLP into a mixed-integer linear program (MILP) for an integrated planning and scheduling problem with flexible recipes, and a set of dynamic optimization problems for generating metamodels. We further apply a bilevel decomposition algorithm to improve the computational efficiency of solving the integrated planning and scheduling problem with flexible recipes. The proposed models and algorithms are demonstrated through case studies of methyl methacrylate (MMA) polymerization processes. The results show that the proposed methods reduce the computational time by more than 2 orders of magnitude compared with the approach of solving the entire integrated optimization problem directly.

1. INTRODUCTION Planning, scheduling, and dynamic optimization of chemical processes are highly interconnected.1−4 The aim of planning is to determine the production assignment of manufacturing processes over a long time horizon. In each planning period, scheduling is used to determine the detailed production sequence, task assignment, processing time, resource allocation, etc. Dynamic optimization determines the optimal trajectories of process inputs and state variables during the transitions between different productions in continuous manufacturing processes.5−7 Though these three problems can be solved in a sequential way, their integration usually results in a better overall result. Thus, the integration of planning, scheduling, and dynamic optimization has been an active research topic recently.8−15 Most previous studies concentrate only on parts of the integrated problem. Some work focuses on the integrated planning and scheduling, where transition times and costs are assumed to be known parameters.16−19 Another group of literature concentrates on the integration of scheduling and dynamic optimization considering variable transition times and costs, determined by time-dependent trajectories.20−28 The full space integrated scheduling and dynamic optimization problem accounts for all of the detailed information on scheduling and dynamic optimization, but it could be computationally intractable due to its scale and complexity.29−34 Recently, a few advances have been made on integrating the process operations and dynamic optimization problems, both for batch processes and for the continuous processes.35,36 The full space problem integrating planning, scheduling, and dynamic optimization could be computationally challenging due to the presence of integer variables and dynamic models. To the best of our knowledge, there are no computationally efficient algorithms to solve the integrated problem for continuous © XXXX American Chemical Society

processes. The aim of this work is to propose efficient solution approaches to solve the integrated process operations and dynamic optimization problems for continuous processes. In this work, we first propose a novel production planning model based on a new approach to estimate the inventory cost. Next, we formulate a mixed-integer dynamic optimization (MIDO) problem for the integrated planning, scheduling, and dynamic optimization. The MIDO problem is then reformulated as a full space mixed-integer nonlinear program (MINLP) problem by discretizing the differential equations. To enhance the computational efficiency, we propose an efficient flexible recipe method based on metamodeling. The flexible recipe method decomposes the full space MINLP problem into a set of inner dynamic optimization problems and an outer planning and scheduling problem. The inner problems are solved offline to generate a number of metamodels that establish the relationship between transition times and costs. The key step of this method is to replace the inner problem with the candidate transition times and costs provided by metamodels, resulting in an integrated planning and scheduling problem with flexible recipes. The integrated planning and scheduling problem with flexible recipes is a mixed-integer linear programming (MILP) problem that is easier to solve than the original full space MINLP problem. To further improve the computational efficiency, a bilevel decomposition algorithm is then applied to the MILP problem for integrated planning and scheduling with flexible recipes. The bilevel decomposition algorithm decomposes the MILP problem into an upper level planning Received: September 30, 2014 Revised: December 30, 2014 Accepted: January 29, 2015

A

DOI: 10.1021/ie503857r Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Scope of the integrated problem.

of differential equations. The results demonstrate that the proposed methods reduce the computational time by more than 2 orders of magnitude compared to the approach of solving the full space integrated problem directly. The major novelties of this work are summarized as follows: • With a new planning model, a novel model formulation for integrated planning, scheduling, and dynamic optimization of multiproduct continuous processes is proposed. • Novel efficient solution methods, including the flexible recipe method and the bilevel decomposition algorithm, are proposed, which can reduce the computational complexity significantly. The remainder of the paper is structured as follows: the problem statement is given in section 2; section 3 introduces the formulation of the integrated planning, scheduling, and dynamic optimization problem; the flexible recipe method is provided in section 4; section 5 describes the bilevel decomposition algorithm; two case studies are used to demonstrate the proposed method in section 6; section 7 concludes the paper.

problem and a lower level detailed scheduling problem. The algorithm solves the upper and lower level problems iteratively. At each iteration, it adds integer and logic cuts to the upper level problem, until the predetermined stopping condition is met. Figure 1 demonstrates the scope of the integrated problem. The planning problem is on the top. The whole planning horizon is divided into a set of planning periods. Order demands are given at the end of each period. Each planning period has a scheduling problem. The scheduling problem is in the middle level. Each scheduling horizon is divided into several time slots. The number of time slots is the same as the number of total products. The scheduling problem determines the production sequence and the production time. Dynamic optimization is at the bottom level. The dynamic optimization problems focus on the transition processes from one product to another. In this work, we use the approximation approach to estimate the inventory cost that is fully explained in the following section. In general, the integrated problem is complex and difficult to solve. In this work, we take advantage of the metamodeling to deal with this difficulty. We replace the dynamic optimization with a set of candidate transition times and transition costs, which enhance the computational efficiency greatly. To demonstrate the applicability of the proposed integration framework and the solution methods, we consider a case study for a methyl methacrylate (MMA) polymerization process with an azobisiso-butyronitrile initiator and a toluene solvent.37,38 Multiple products with different molecular weights are produced during the steady-state production periods. The dynamic model for the transition processes is described by a set

2. PROBLEM STATEMENT We study a multiproduct continuous manufacturing process that manufactures a set of products over the planning horizon. The planning horizon is divided into several planning periods. As shown in Figure 2, the total number of the planning periods is denoted by np. In each period there are ns time slots. Each time slot consists of a production period and a transition period.39,40 In the production period the production process is B


Article


then presented, which is built by incorporating the planning model, the scheduling model, and the dynamic optimization model into a complex MIDO model.14,35,41−43 A list of indices, sets, parameters, and variables is given in the Nomenclature section, where the variables are denoted with upper case symbols, and the parameters are denoted in lowercase symbols or Greek letters. 3.1. Dynamic Model. In this section, we present the general mathematical formulation for the dynamic model and the discretization method to transform the dynamic model into an algebraic model. The dynamic model in slot l of period t is given in a general form as follows. Equation 1 is the compact form of the differential equations that represents the features of the dynamic processes, such as the heat balance and the mass balance. Equation 2 is the compact form of the algebraic equations that denote the phase equilibriums. Equation 3 is used to set the initial conditions of the state variables at the beginning of the transition process. Equation 4 is used to set the final values of the outputs at the end of the dynamic time in slot l of period t. Constraints 5−8 are used to confine all the variables between their corresponding upper and lower bounds. Zlt(Tlt) is the state variables in slot l of period t. Xlt(Tlt) denotes the algebraic variables in slot l of period t. Ult(Tlt) is the manipulate variables in slot l of period t. Ylt(Tlt) is the outputs of the dynamic model in slot l of period t. They are all time-dependent variables. qlt is the time-independent parameters in slot l of period t. Tlt is the time variable, and TLlt is the end time of the dynamic time in slot l of period t. ZIni lt is the initial conditions in slot l of period t. YSet is the output set points in slot l of period t. lt

Figure 2. Planning horizon, planning periods, and time slots.

operated in the steady state, while in the transition period the production process changes from one steady state to another. We note that the transition period can be zero if there is no transition. In general the same product is manufactured in both the last scheduling slot of one planning period and the first scheduling slot of the next planning period. So usually there is no transition process in the last scheduling slot of the planning period. Production assignments of each planning period are determined by the planning model according to order demands. In this work, we assume that all the extra products can be sold out after the order demands are met. The production sequence and the duration of each time slot in a planning period are determined by the scheduling model. The dynamic optimization problem is solved to determine the transition time and the transition cost between two products. Specifically, the integrated problem of planning, scheduling, and dynamic optimization is stated as follows: Process: Multiproduct continuous manufacturing process Given: A set of products Planning horizon, planning periods, and time slots Demand for each product at the end of each period Steady-state operating conditions Unitary operating costs and inventory costs for different products Dynamic models Determine: Production amounts in planning periods Production sequences Production times and production costs Transition times and transition costs Trajectories of dynamic models Objective: To maximize the total profit, that is the revenue minus the inventory cost, the production cost, and the transition cost

dZlt(Tlt ) = flt (Zlt(Tlt ), Xlt (Tlt ), Ult(Tlt ), qlt ), ∀ l , t dTlt

(1)

glt (Zlt(Tlt ), Xlt (Tlt ), Ult(Tlt ), qlt ) = 0, ∀ l , t

(2)

Zlt(Tlt = 0) = ZltIni , ∀ l , t

(3)

Ylt(TLlt ) = Y ltSet ,

(4)

∀ l, t

ZltL ≤ Zlt(Tlt ) ≤ ZltU ,

∀ l, t

(5)

XltL ≤ Xlt (Tlt ) ≤ XltU ,

∀ l, t

(6)

UltL ≤ Ult(Tlt ) ≤ UltU

, ∀ l, t

(7)

YltL ≤ Ylt(Tlt ) ≤ Y ltU ,

∀ l, t

(8)

Equation 9 defines the transition cost CDlt as a function of the process variables in slot l of period t in the dynamic model. And the transition time in slot l of period t in the dynamic model is TLlt. CDlt = ϕlt (Zlt(Tlt ),Xlt (Tlt ),Ult(Tlt ),qlt ),

∀ l, t

(9)

To solve the dynamic optimization problem, the collocation method is usually applied to discretize the differential and algebraic equations.44−48 First, in slot l of period t the time variables are divided into a set of equal-length finite elements, which means the lengths of finite elements in a specific slot are the same. The original continuous dependent variables are then sampled at each collocation point. After this discretization process, the differential equations are approximated by a set of algebraic equations. The discretized dynamic model is shown as follows.

3. MODEL FORMULATION A novel linear approach to estimate the inventory cost is proposed in this section. This new estimation approach reduces the estimation error significantly. With this novel estimation approach, the full space MINLP model for the integrated planning, scheduling, and dynamic optimization problem is C


Article

Industrial & Engineering Chemistry Research In eq 10, Zrslt is the discretized state variables at collocation point s of finite element r in slot l of period t. Zrlt is the discretized state variable at the beginning of finite element r in slot l of period t. Xrslt denotes the discretized algebraic variables at collocation point s of finite element r. Urslt is the discretized manipulate variables at collocation point s of finite element r. Llt is the length of the finite element at slot l of period t of finite element r. css′ and ds represent the constant matrices that are calculated according to the collocation method. Equation 12 represents the mass and heat balance, etc. Equations 13 and 14 are to set the initial values of the state variables and the final values of the outputs, where nf lt is the number of the finite elements in slot l of period t. Constraints 15−18 confine all variables between their upper and lower bounds.

in slot l of period t is zero if product i is not assigned to that slot. 0 ≤ PTilt ≤ ht ·Wilt ,

Git = ri·∑ PTilt ,

TElt = TSlt + (10)

qlt ) = 0,

Zltr= 1, s = ZltIni , r = nf , s = 3 Ylt lt

∀ l, t , r, s

∀ l, t , s

Y ltSet ,

ZltL ≤ Zltrs ≤ ZltU ,

∀ l, t , r, s

(15)

XltL

XltU ,

∀ l, t , r, s

(16)

UltL ≤ Ultrs ≤ UltU ,

∀ l, t , r, s

(17)

≤

≤

YltL ≤ Yltrs ≤ Y ltU ,

∀ l, t , r, s

TLlt = nflt ·Llt ,

∀ l, t

i

∀ l < ns , t

TEns , t = htt ,

(26)

∀t

(27)

The symmetry breaking constraints 28−32 simplify the solution space when more than one slot is assigned to the same product in a planning period,12,19,36 where the variable NYit is the total number of slots assigned to produce product i in period t, NPit is the binary variable that denotes if product i is produced in period t, npro is the total number of all the products, and bm is a big positive number (the big-M method). bm should be no less than the number of all the products, for example, npro + 1 could be assigned to bm. NYit =

(19)

∑ Wilt ,

∀ i, t (28)

l

NPit ≥ Wilt ,

(20)

We should note that this section only presents a general dynamic model. This general model can be customized for many specific applications. For this work, a methyl methacrylate (MMA) polymerization dynamic process is used in section 6. 3.2. Scheduling Model. The scheduling model is presented in this section. We should note that there are several scheduling problems in the integrated model since there is one scheduling problem in each planning period. The number of planning periods is given as a fixed parameter. Constraint 21 means that exactly one product is assigned to each slot. The binary variable Wilt = 1 implies product i is assigned to slot l of period t.

∑ Wilt = 1,

(25)

Equation 27 ensures the ending time of the last slot in period t is equal to the ending time for period t.

(18)

∀ l, t

∀ l = ns , t < np

TElt = TSl + 1, t ,

The transition cost and the transition time can be formulated in eqs 19 and 20, where CDlt is the transition cost in slot l of period t in the dynamic model and TLlt is the transition time in slot l of period t in the dynamic model. CDlt = ϕlt (Zltrs , Xltrs , Ultrs , qlt ),

(24)

Equation 26 indicates that the ending time of one slot is the same as the starting time of the next one.

(13) (14)

Xltrs

TElt = TS1, t + 1 ,

(12)

∀ l, t

=

∀ l, t

Equation 25 guarantees that the ending time of the last slot in each period is the same as the starting time of the first slot in the next period. (11)

Ultrs ,

∑ PTilt + TTlt , i

s

Xltrs ,

(23)

In constraint 24, the ending time TElt of slot l of period t equals the sum of the starting time TSlt, the production time in that slot, and the transition time TTlt.

Zltr = Zltr− 1 + Llt ·∑ ds·flt (Zltrs , Xltrs , Ultrs , qlt ),

glt (Zltrs ,

∀ i, t

l

s′

∀ l, t , r > 1

(22)

where PTilt is the production time of product i produced in slot l of period t, and ht is the duration of period t. Constraint 23 shows that the production amount of product i in period t is the product of the production rate ri and the total production time in period t.

Zltrs = Zltr− 1 + Llt · ∑ css ′·flt (Zltrs , Xltrs , Ultrs , qlt ), ∀ l , t , r > 1, s

∀ i, l, t

∀ i, l, t

NPit ≤ NYit ≤ npro ·NPit ,

(29)

∀ i, t

NYit ≥ npro − (∑ NPit − 1) − bm ·(1 − Wi ,1, t ),

(30)

∀ i, t

i

(31)

NYit ≤ npro − (∑ NPit − 1) + bm ·(1 − Wi ,1, t ),

∀ i, t

i

(32)

3.3. Planning Model. We present a novel planning model in this section based on previous work.12,13,17,19 The novelty lies in the newly proposed linear approach to estimate the inventory cost. The reason why the estimated inventory cost is used instead of the real inventory cost is that the estimated approach provides a linear way to calculate the inventory cost while the approach to calculate the real inventory cost is nonlinear and nonconvex, as shown in the Appendix. Since a

∀ l, t (21)

Constraint 22 shows the production time is bounded by the length of the planning period. The production time of product i D


Article


Figure 3. Inventory cost. (I) Real inventory cost in a planning period. (II) Trapezoid-shaped inventory cost in a planning period.

inventory cost rate for all products at the beginning of planning period t. INRt is the inventory cost rate for all products at the end of planning period t before the demand is met, as defined in eq 38. cinv i is the inventory cost for unit product t.

linear problem is usually much easier to solve compared to a nonlinear problem, most of the current work uses the linear approach to estimate the inventory cost.12,16−19,35,36 However, the existing approaches to estimate the inventory cost lead to a big estimation error. The estimation error is the difference between the real inventory cost and the estimated inventory cost. Unfortunately, most previous work just ignores this estimation error directly. In this work, we propose a novel linear approach to estimate the inventory cost that reduces the estimation error significantly. Figure 3 illustrates the difference between the real inventory cost and the estimated inventory cost approximated by our newly proposed approach. As shown in Figure 3 (I), the real inventory cost in one planning period is shaped by the shadow area that is a polygon. Existing literature overestimates the inventory cost using the product of the time length of the planning period and the final total inventory cost rate of the planning period before the demand is met,12,16−19,35,36 which is the area of rectangle (ADEB) in Figure 3 (II). However, there is a large approximation error between the estimated inventory cost that is calculated by the rectangle method in previous work and the real inventory cost. To reduce the large approximation error, a novel approach is proposed to estimate the inventory cost. As shown in Figure 3 (II), the trapezoid (CDEB) is used to approximate the inventory cost. This approach allows for a more accurate inventory cost estimation while keeping the linearity of the planning and scheduling model formulation. Equation 33 shows that the inventory level INVit of product i in the first period is the sum of the initial inventory level invInitial i and the production amount Git in the first planning period. INVit = inviInitial + Git ,

∀ i, t = 1

INROt =

INVOit = INVit − Sit ,

∀ i, t > 1 ∀ i, t

∀t=1 (36)

i

INROt =

∑ (INVOi ,t− 1·ciinv),

∀t>1 (37)

i

INR t =

∑ (INVi ,t ·cinv i ),

∀t>1 (38)

i

Equation 39 indicates that the estimated inventory cost is the area of the trapezoid, which can be seen in Figure 3 (II). INVCt is the estimated inventory cost of all products in period t. INVCt = (INROt + INR t ) ·

ht , 2

∀t

(39)

For clarity of presentation, the nonlinear approach to calculate the real inventory cost and the estimated approach used in previous works12,18,19 to estimate the inventory cost are presented in the Appendix. In this work, we use the newly proposed linear approach to estimate the inventory in the optimization model. The newly proposed approach has a few advantages over those presented in the Appendix. First, compared with the nonlinear approach to calculate the real inventory cost, the newly proposed approach allows the integrated planning and scheduling problem to be a mixedinteger linear optimization problem. Second, compared with the estimated approach in previous works,1,8,20,21,24,26,27 the newly proposed approach has a much smaller estimation error and is very close to the real inventory cost. Consequently, we use the newly proposed approach to estimate the inventory cost. Equation 40 ensures that in each time period the sale of one product is greater than or equal to the demand of that product, where dit denotes the order demand for product i in period t, shown as the lower bound for the sale at the end of each period. The inequality also means that all the extra products can be sold out after the order demands are met. We should note that the production process cannot manufacture the amount of product demands that are beyond the capacity of the reaction system. In this problem, the order demands are given parameters that are smaller than the production capacity. Otherwise, the problem will be infeasible.

(33)

Equations 34 and 35 calculate the inventory levels for each period; Sit is the sale amount of product i in period t. INVit = INVOi , t − 1 + Git ,

∑ (inviInitial ·cinv i ),

(34) (35)

where INVit denotes the inventory level at the end of each period before the products are sold, and INVOit denotes the final inventory level at the end of each period after the products are sold, which is also the initial inventory level at the beginning of the next period. The new approach to estimate the inventory cost is introduced as follows. In eqs 36 and 37, INROt is the initial E


Article

Industrial & Engineering Chemistry Research Sit ≥ dit ,

∀ i, t

Equations 49 and 50 indicate that the set point YSet lt of the output of slot l in period t is determined by the product manufactured in the subsequent time slot, where ySetpoint is the i output set point of product i.

(40)

Equation 41 defines the revenue RVN, where pit is the selling price of product i in period t.

RVN =

∑ ∑ pit ·Sit

(41)

Y ltSet =

In eq 42, the total inventory cost is the sum of those over all the planning periods.

Y ltSet =

i

t

Equation 43 defines the total operation cost COP, where coper it is the operation cost for producing unit product i in period t, Git is the amount of product i produced in period t.

∑ ∑ citoper ·Git i

(43)

t

Equation 44 defines the total transition cost CTR over the planning horizon, where CTlt is the transition cost in time slot l of period t.

CTR =

∑ ∑ CTlt l

(44)

t

max PROFIT = RVN − CIV − COP − CTR

(45)

CTlt = CDlt ,

∀ l, t

(46)

TTlt = TLlt ,

∀ l, t

(47)

(51)

(Integrated) max PROFIT given in eq (45) s.t. Dynamic optimization model (10)−(20) Scheduling model (21)−(32) Planning model (33)−(39), (40)−(44) Linking equations (46)−(50)

The integrated problem includes one planning problem over the planning horizon, a set of scheduling problems over each planning period, and a number of dynamic optimization problems relating to the transition processes. Existing publications suggest that the direct solution approach for the integrated MINLP problem could be computationally intractable.14,35 Since in general the formulated MINLP problem is very complex, it is difficult to find the optimal solution to the MINLP problem even for offline calculation.21,25,36,49 Efficient solution methods are needed to reduce the computational time.20,50 In this work, we apply the flexible recipe method and a bilevel decomposition method to increase the computational efficiency.18

where CTlt is the transition cost used in the planning and scheduling model in slot l of period t, TTlt is the transition time used in the planning and scheduling model in slot l of period t. Equation 48 determines the initial value of the state variables Initial ZIni is the lt according to the scheduling assignment, where zi steady-state value of the state variables for manufacturing product i. i

(50)

By integrating the planning, scheduling, and dynamic optimization problems discussed in the previous subsections, the integrated problem can be formulated into a large scale MINLP problem as follows:

where RVN denotes the total revenue for all products, CIV is the total inventory cost, COP is the operation cost for all products, and CTR is the total transition cost throughout the planning horizon. 3.4. Integrated Model. In this section, we present the integrated model by combining the planning and scheduling model with the dynamic model. The dynamic optimization problem mainly deals with the transition process in the transition period. The transition times and costs are major variables in the transition process, which will be provided to the integrated planning and scheduling problem. The way of linking the planning and scheduling part to the dynamic optimization part is based on the model formulation. In this model formulation, the transition times and transition costs are used to link the planning and scheduling part and the dynamic optimization part together. Equations 46 and 47 ensure that the planning and scheduling model is linked to the dynamic optimization models via transition times and transition costs.

∑ Wilt ·ziInitial ,

∀ l = Ns , t < Np

yiSetpoint ·(1 − b) ≤ Yi ≤ yiSetpoint ·(1 + b)

The total profit in eq 45 is divided into 4 parts: the revenue of selling the products, the inventory costs, the operation costs, and the transition costs. The objective function is to maximize PROFIT given below:

ZltIni =

∑ Wi ,1,t + 1·yiSetpoint ,

(49)

During the transition period, the production process transfers from one production period to another production period. We introduce a quality bound around the steady-state value of the output to define the boundary between the production period and the transition period. In this way, we make the products meet the production requirements during the production period. The production period begins when the output falls into the bound and ends when the next transition period begins. During the production period, the output needs to stay between the bounds continuously, as shown in constraint 51. The width of the bound is defined by the value of b. The value of b is typically chosen as 2%.20 Yi is the output of product i.

(42)

t

COP =

∀ l < Ns , t

i

i

∑ INVCt

CIV =

∑ Wi ,l+ 1,t ·yiSetpoint ,

4. FLEXIBLE RECIPE METHOD The structure of the full space problem is shown in Figure 4. In the integrated problem, problems at different levels are solved simultaneously. The planning problem provides the scheduling problem with production assignments, and the scheduling problem returns production costs and production times to the planning problem. The scheduling problem provides the dynamic optimization problem with the production sequence, and the dynamic optimization problem returns transition costs

∀ l, t (48) F


Article


decision determines the product assignment W ilt, the production time Pilt, and the operation cost COP. The dynamic decision determines the transition costs and times. If the planning and scheduling variables can be determined, the total transition cost CTR is decided solely by the dynamic decision variables. Based on the above observation, the original integrated planning, scheduling, and dynamic optimization problem is decomposed into two subproblems, as shown in the following “Decomposed_Integrated” problem. The outer problem is the integrated planning and scheduling problem. The inner problem is the dynamic optimization problem. The outer problem and the inner problem are linked by the linking equations. (Decomposed_Integrated)

Figure 4. Structure of the integrated planning, scheduling, and dynamic optimization problem.

(Outer) PROFIT =

(53)

and transition times to the scheduling problem. Due to the complexity of the full space MINLP model, the direct solution approach is time-consuming, although the stand-alone planning, scheduling, and dynamic optimization subproblems are relatively easier to solve.30 To overcome the computational challenge of solving the full space problem, we propose a decomposition method based on the structure of the full space MINLP problem, as shown in Figure 4. Studying complex phenomena and real world systems or solving challenging enterprise-wide optimization problems is needed by a large number of scientific and engineering fields. The metamodeling techniques have already been widely used in those fields.51 The decomposition method takes advantage of a set of metamodels that features the information on dynamic optimization problems.36 The structure of the integrated problem is shown in Figure 4. The main idea of the decomposed approach is that the dynamic optimization problem in each individual slot can be solved independently if we fix the production sequencing decisions determined by the planning and scheduling model. If the initial Set condition ZIni it and the set point Ylt are known, a dynamic optimization problem implicitly defines the transition cost CTlt as a function of the transition time TTlt. However, the function characterizing the relation between the transition cost and the transition time is a black-box function without closed-form expression.36 Metamodeling is then introduced to replace this complex function in a simple approach.51,52 The metamodel is constructed based on the response of the exact functions to a set of sampling points. The number of the sampling points should be sufficiently large to generate a relative accurate metamodel.20 Those sampling points are the candidate flexible recipes for the integrated planning and scheduling problem with a flexible recipe.21 4.1. Decomposition of the Integrated Problem. To be specific, we rewrite the objective function 45 with decision variables, as shown in function 52: max

{DV P , DV S , DV D}

max {RVN − CIV − COP − CTR }

{DV P , DV S}

s.t. Scheduling model (21)−(32) Planning model (33)−(39), (40)−(44) (Link) Linking equations (46)−(50) (Inner) maxD { −CTR } = minD {CTR }

(54)

{DV }

{DV }

s.t. Dynamic optimization model (10)−(20)

According to eq 44, we rewrite the inner problem as follows: (Decomposed_Inner) min CTR = minD

{DV D}

{DV }

∑ ∑ CTlt , l

∀ l, t (55)

t


Since in the inner problem the planning and the scheduling decision variables are treated as given parameters, the transition sequence is determined. Consequently, the dynamic optimization problem in one time slot is solved independently from that in another time slot. The inner problem is then further simplified as follows: (Decomposed_Inner) min

{DV D}

CTlt } ∑ ∑ CTlt = ∑ ∑ {{min DV } D

l

t

l

(56)

t


The inner problem implicitly defines the transition cost as a function of transition time. Consequently, in eq 56 the inner dynamic problem in each individual time slot can be replaced by an optimal value function φlt(TTlt) that characterizes the relationship between the transition cost and the transition time. The optimal value function solves the dynamic optimization problem to yield the minimum transition cost at a certain transition time, as shown in eq 57. When an inner dynamic optimization problem is solved, the decision variables from the outer problem are assumed as given parameters.

PROFIT = RVN − CIV − COP − CTR (52)

where DVP denotes the decision variables in the planning level problem, DVS denotes the decision variables in the scheduling level problem, and DVD is the decision variables from the dynamic optimization level problem. According to the model formulation in section 3, the planning decision determines the sale of each product Sit, the inventory cost CIV, and the revenue RVN. The scheduling

(Optimal_Value) φlt (TTlt ) = minD CTlt , {DV }

G

∀ l, t

(57)


Article

Industrial & Engineering Chemistry Research s.t. Dynamic optimization model (10)−(20)

(Min_Time) ttikmin = min TLFik ,

The inner dynamic optimization problems are replaced by the corresponding optimal value functions. The decomposed integrated problem with the optimal value functions is reformulated as follows. (Decomposed_Integrated) (Outer) PROFIT =

s.t. dZFik(Tik) = fik (ZFik(TFik),XFik(TFik),UFik(TFik),qfik ) dTFik (60)

max {RVN − CIV − COP − CTR }

{DV P , DV S}

(58)

s.t. Scheduling model (21)−(32) Planning model (33)−(39), (40)−(44) Linking equations (46)−(50) CTR =

∑ ∑ φlt (TTlt ) l

t

(Inner) φlt (TTlt ) = minD CTlt , {DV }

∀ l, t

∀ i, k

(59)

RCik = ϕik (ZFik(TFik),XFik(TFik),UFik(TFik),qfik )

(61)

gik (ZFik(TFik),XFik(TFik),UFik(TFik),qfik ) = 0

(62)

ZFik(TFik = 0) = ziInitial

(63)

YFik(TLFik) = ykSetpoint

(64)

ZFikL ≤ ZFik(TFik) ≤ ZFikU

(65)

XFikL ≤ XFik(TFik) ≤ XFikU

(66)

UFikL ≤ UFik(TFik) ≤ UFikU

(67)

YFikL ≤ YFik(TFik) ≤ YFikU


(68)

ttmin ik

After the minimum transition time is obtained, we then need to determine other sampling points to build the metamodel. One approach to design the sampling points is to set the interval step size between two adjacent sampling points and the total number of the sampling points. This approach is show in eqs 69 and 70:

In the “Decomposed_Integrated” problem, the integrated planning, scheduling, and dynamic optimization problem is decomposed into the outer integrated planning and scheduling problem and a number of inner individual dynamic optimization problems. The dynamic optimization problems are represented by a set of optimal value functions. The outer integrated planning and scheduling problem is an MILP problem, which is easy to solve. Nevertheless, the optimal value functions have no closed-form expressions. The computational complexity remains at the optimal value functions, which are black-box functions. We use metamodels to represent the optimal value functions to overcome this difficulty and improve the computational efficiency. 4.2. Integrated Planning and Scheduling Model Based on Flexible Recipes. The first step to build the metamodels is determining the feasible range of each optimal value function. We first find the minimum transition time in each individual slot, which is determined by solving the dynamic optimization problem with the objective to minimize the transition time. The initial condition ZIni lt depends on the product assigned to the current slot, while the set point YSet lt depends on the product assigned to the following slot. To build a general metamodel, we consider the transition from product i to product k. In this transition, the initial condition zInitial and the set point ySetpoint are regarded as given i k parameters. The transition cost from product i to product k is denoted by RCik to distinguish the transition cost CTlt in slot l of period t. The minimum transition time ttmin ik is obtained from solving the following “Min_Time” problem. The constraints and the variables in the “Min_Time” problem are similar to the constraints and the variables described in subsection 3.1. For example, both Tlt and TFik are time variables and both Zlt and ZFik are state variables, etc. The variables in subsection 3.1 are for the dynamic problem in slot l of period t, while the variables in the “Min_Time” problem are for the dynamic problem from product i to product k.

rtikg = ttikmin + ISLg , ISLg = (g − 1) ·α ,

∀ g ∈ {1, 2, ..., nc} ∀ g ∈ {1, 2, ..., nc}

(69) (70)

where a set of the sampling transition times are given by rtikg and g is the total number of sampling points to choose. The interval step size between two adjacent sampling points is denoted by α. ISLg is the interval length between the minimum transition time and the g-th sampling point. The metamodel is built offline so that a trial-and-error procedure can be used to determine the interval step size and the total number of the sampling points. At each value of the transition time rtikg given in eq 69, the minimum transition cost rcikg is then calculated by solving the following optimization problem (Min_Cost) rcikg = min RCik ,

∀ i, k, g

s.t. Dynamic optimization model (60) − (68) TLFik = rtikg

In this way, we obtain a set of discrete transition times and transition costs as {rtikg , rcikg }g = 1, ···, Nc ,

∀ i, k

(71)

With those discrete pairs of transition time and transition cost, a general metamodel is built. This general approach to build a metamodel is then applied to all potential transitions between different products to build a set of metamodels. For example, if there are 3 different products (A, B, and C) in total, H


Article


transition time in the corresponding metamodel. Since the integrated planning and scheduling problem with the flexible recipes is an MILP problem, it can be solved more efficiently than the full-space MINLP problem. By replacing the dynamic models with a set of metamodels, the flexible recipe method substantially reduces the computational complexity for solving the integrated MINLP problem. In our work, we do not use interpolation between sample points. However, the metamodeling method does allow interpolation between sample points if the model formulation is modified accordingly. For example, linear piece-wise interpolation can be used between sample points when metamodels are built.21 The linear piece-wise interpolation increases the accuracy of metamodeling while it also makes the formulation of the metamodeling more complex. We should note that all the sampling points are feasible and locally optimal for the original integrated problem, so there is no sacrifice of feasibility in this approach. Specifically, the integrated planning and scheduling problem with flexible recipes is listed as follows:

there are 6 metamodels to be built for the potential transitions of A to B, A to C, B to C, B to A, C to A, and C to B. The next step is to link the outer integrated planning and scheduling model with the inner metamodels. The linking binary variables are used to link the outer problems and the inner problems together as follows. In constraints 72 and 73, the linking binary variable ZBikltg is defined. After that, the discrete values in the collection 71 are linked to the transition times and the transition costs in time slot l of period t in constraints 74 and 75, respectively. ZBikltg is a binary variable to select recipes. ZBikltg is equal to one if product i is assigned to slot l of period t, product k is assigned to the subsequent slot, and the g-th flexible recipe is selected. Only a pair of the transition time and the transition cost is selected for each transition.

∑ ZBikltg

≥ Wilt + Wk , l + 1, t − 1,

∀ i , k ≠ i , l < Ns , t

g

(72)

∑ ZBikltg

≥ Wilt + Wk , ll , t + 1 − 1,

(Integrated_SM) max PROFIT given in equation (45)

g

s.t. ∀ i , k ≠ i , t , l = Ns , ll = 1

CTlt =

∑ ∑ ∑ ZBikltg ·rcikg , i

TTlt =

k

k

g

Planning model (33)−(39), (40)−(44) Scheduling model (21)−(32)

∀ l, t (74)

g

∑ ∑ ∑ ZBikltg ·rtikg , i

(73)

Sampling point collection of

∀ l, t

metamodeling (71) (75)

Selection of transition

Since the inner dynamic optimization problems are replaced by the metamodels that represent variable task recipes, this method is referred to as the flexible recipe method. The flexible recipe method reformulates the full space problem in Figure 4 to the flexible recipe problem in Figure 5.

times and transition costs (72)−(75)

5. BILEVEL DECOMPOSITION ALGORITHM In this section, we apply a bilevel decomposition algorithm18,53−56 to improve the computational efficiency of solving the MILP problem for integrated planning and scheduling with flexible recipes. This algorithm decomposes the MILP problem into an upper level planning problem and a lower level scheduling problem with flexible recipes. The production assignment in each time period, the estimated production amount, and the estimated inventory level of each product are determined by the upper level problem. It generates an upper bound of the total profit, since it is a relaxation from the original MILP model. The lower level problem exploits the solution from the upper level problem. In the lower level, the original MILP problem is solved by only searching the subset of the production assignment determined by the upper level problem, i.e. making production assignment decisions only based on the products that are chosen by the upper level problem. The size of the lower level problems is the same as the original integrated planning and scheduling problem with flexible recipes, whereas the feasible region of the lower level problem at each iteration is much smaller compared to the original problem. Since the lower level problem only gets optimal solution within a confined feasible region at each iteration, the solution to the lower level problem at each iteration provides a lower bound of the total profit. The upper level and the lower level problems are solved iteratively until the gap between the upper bound and the lower bound drops under a predefined optimality tolerance. To hasten the convergence rate of upper and lower bounds, integer cuts and logic cuts are added into the upper level problem at

Figure 5. Diagram of the flexible recipe method for the integrated problem.

The diagram of the flexible recipe method is summarized in Figure 5. This method consists of the integrated planning and scheduling problem with flexible recipes and a set of dynamic optimization problems. First, the dynamic optimization problems are solved offline to build a set of metamodels. The metamodels are then incorporated into the integrated planning and scheduling problem as flexible recipes. Only the integrated planning and scheduling problem with flexible recipes is then solved directly. The integrated planning and scheduling problem only chooses the flexible recipes that are located within the range between the upper and lower bounds of the I


Article


level problem can be formulated into a new form shown in eq 76. The revenue, inventory cost, and operation cost terms have the same meaning as those in eq 41 to eq 43. The difference lies in the last term for the transition costs, where Yit is the binary variable to denote if product i is assigned to period t, MITCi is the minimum transition cost from product i to other products, and MATCt is the maximum among the minimum transition costs in period t, which is defined by constraints 77 and 78. We note that, in constraints 78, there an upper bound for transition cost MATCt is posed. However, eq 76 shows that the objective is to maximize the profit and drives MATCt to stay at its upper bound.

each iteration. Previous solutions obtained from the upper level problem are kept out by means of the integer cuts. The logic cuts reject the subsets and supersets of the previous solutions in addition to the ones that cannot satisfy the capacity limitations. The flowchart of the bilevel decomposition algorithm is shown in Figure 6.

max UPP =

∑ ∑ pit ·Sit − ∑ INVCt − ∑ ∑ citoper ·XPit i

−

t

t

i

t

∑ (∑ (MITCi·Yit ) − MATCt ) t

i

MATCt ≥ MITCi·Yit , ∀ i , t MATCt ≤ max{MITCi}, i

(76) (77)

∀t

(78)

Equation 79 shows that the production amount equals the production rate multiplied by the production time. XPit = ri·PTUit ,

∀ i, t

(79)

where XPit is the production amount of product i in period t, and PTUit is the production time of product i in period t in the upper level problem. Constraint 80 sets the length of a planning period as the upper bound of the production time for products assigned to that period.

Figure 6. Flowchart for the bilevel decomposition algorithm.

PTUit ≤ ht ·Yit ,

5.1. Upper Level Problem. The main idea of formulating the upper level problem is to solve the planning problem first without fully considering the scheduling problem and the transition problem. Since the upper level problem concentrates on the production assignment in each period, the specific production sequence in each planning period and the precise transition trajectories from one product to another cannot be obtained by solving the upper level problem. We cannot disregard the transition cost and the transition time completely; otherwise, the relaxation is too loose, which may result in slow convergence of the algorithm. The excessively loose relaxation might lead to inappropriate production assignments, which may cause infeasibility of the lower level problem. To overcome the excessively loose relaxation, we use the minimum transition time and cost for each product assigned to the certain planning period. We note that in general the product manufactured in the last scheduling slot of one planning period is also manufactured in the first scheduling slot of the next planning period. So there is no transition process in the last scheduling slot of the planning period. The transition process between two adjacent planning periods is ignored, which results in a reasonable and tight upper bound for the integrated problem. Since the transition between two adjacent planning periods is ignored, the total number of transitions in one planning period plus one is the same as the total number of products assigned to this period. Thus, we first determine all the minimum transition time for each product assigned to a certain planning period and then neglect the maximum of them in the model. Since we do not consider the detailed production sequence in the upper level problem, the objective function of the upper

∀ i, t

(80)

Constraints 81 to 83 define the time relation in the upper level problem, similar to those for the transition costs. MITTi is the minimum transition time from product i to other products, and MATTt is the maximum among the minimum transition times in period t. Constraint 83 guarantees all the production and transitions are completed within the duration of planning period t. In eq 82, an upper bound for transition time MATTt is posed. However, objective 76, and constraints 81−83 drive MATTt to stay at its upper bound. MATTt ≥ MITTi ·Yit ,

∀ i, t

MATTt ≤ max{MITTi }, i

(81)

∀t

(82)

∑ (PTUit + MITTi ·Yit ) − MATTt ≤ ht , i

∀t (83)

Equations 84 and 85 define that the inventory level of product i before the demands are met for the current period equals the sum of the inventory level of product i at the beginning of this period and the amount of product i produced in this period. INVit = inviInitial + XPit ,

INVit = INVOi , t − 1 + XPit ,

∀ i, t = 1

∀ i, t ≠ 1

(84) (85)

The upper level problem also includes constraints 35 to 40, which have been discussed in section 3.3.57 The upper level problem is summarized as follows: J


Article


index pair (i,t) if product i is not assigned to period t at previous iteration b in the upper level problem. The second cut 89 is used to exclude the subset of the previous feasible solutions from the upper-level problem at iteration b. At previous iteration b, all the subsets of the feasible solutions from the upper-level problem are searched in the lower level of iteration b. Hence, it is not necessary to search the subset of the previous feasible solutions again.

(Upper_Level) max UPP given in equation (76) s.t. Constraints (35)−(39), (40), (77)−(85)

5.2. Lower Level Problem. The lower level problem is used to solve the original integrated planning and scheduling problem with flexible recipes by excluding those products which are not assigned by the upper level problem at the current iteration. In this way, the feasible region of the lower level problem at each iteration is reduced to enhance the computation efficiency. Constraint 86 ensures that production assignment of period t in the lower level is only a subset of the production assignment of period t obtained in the upper level. Integer cut 87 is added into the lower level problem at each iteration to exclude previous solutions obtained from the lower level problem at iteration b. YLit ≤ Yit ,

∑

∀ l, t

YLit −

(i , t ) ∈ YLYb

∑

∑

(87)

where YLYb = = 1} contains all the index pair (i,t) if product i is assigned to period t at previous iteration b, and YLNb = {(i,t)|YLbit = 0} contains all the index pair (i,t) if product i is not assigned to period t at previous iteration b. The lower level problem is summarized as follows:

htnp ≥ [

t

Yit + Yi ′ t ′ − |YLYb|],

∀ (i′, t ′) ∈ YLNb (91)

Scheduling model (21 )−(32) Sampling point collection of metamodeling (71) Selection of transition times and transition costs (72 )−(75) Constraints (86 )−(87)

5.3. Integer and Logic Cuts. In each iteration, the stopping condition is checked immediately after the lower level problem is solved. The stopping condition is examined by checking if the difference between the upper bound and the lower bounds is within the predetermined tolerance. If the value of the difference between the upper bound and the lower bound is smaller than the predetermined tolerance, the algorithm stops; otherwise, integer and logic cuts are added into the upper level problem to get new solutions and a new iteration is then carried out.54 The first cut 88 is used to exclude previous feasible solutions from the upper-level problem at iteration b. In this way, a new production assignment is obtained at a new iteration.

∑

∑

l

5.4. Algorithm Steps. The Bilevel Decomposition Algorithm 1 Initialize: the upper bound UB = ∞; the lower bound LB = −∞; the tolerance ε is chosen according to the optimality requirement. 2 Set the iteration count b = 1. 3 while (UB − LB )/|UB| > ε do 4 Solve the “Upper_Level” problem to determine the production assignment Ybit of iteration b as well as the upper bound UB. Set YUYb = {i,t|Ybit = 1} and YUNb = {i,t|Ybit = 0}. 5 Solve the “Lower_Level” problem to determine the production assignment YLitb, production sequencing, production time, and transition as well as the lower bound LBb. Set YLYb = {i,t|YLbit = 1} and YLNb = {i,t|YLbit = 0}, and LB = max{LBb}. 6 end while

Planning model (33 )−(39), (40 )−(44)

(i , t ) ∈ YUNb

∑ ∑ ∑ PTilt + TTTb·

(i , t ) ∈ YLYb

s.t.

Yit −

∀ (i′, t ′) ∈ YUNb (90)

i

(Lower_Level) max PROFIT given in eq 45

∑

Yit + Yi ′ t ′ ≤ |YUYb| ,

The capacity logic cut 91 exploits the solutions of the previous lower level problems to exclude those solutions whose total production time and transition time exceed the length of the planning horizon. TTTb is the total transition time at iteration b.

{(i,t)|YLbit

(i , t ) ∈ YUYb

(89)

(i , t ) ∈ YUYb

YLit ≤ |YLYb| − 1

(i , t ) ∈ YLNb

∀ (i′, t ′) ∈ YUYb

The third cut 90 is used to exclude the superset of the feasible solutions from previous iteration b in the upper-level problem. At previous iteration b, all the supersets of the feasible solutions have already been searched in the upper level problem. Hence, there is no need to search those supersets again.

(86)

∑

Yit + Yi ′ t ′ ≥ 1,

(i , t ) ∈ YUNb

6. CASE STUDIES To demonstrate the proposed modeling and solution frameworks, we consider a case study for a methyl methacrylate (MMA) polymerization process, which is a nonlinear free radical polymerization using azobisiso-butyronitrile as the initiator and toluene as the solvent.37,38 All optimization problems are modeled in GAMS 24.2.158 and solved in a PC with a 3.10 GHz CPU, 8GB RAM, and Window 7 64-bit operating system. 6.1. Mathematic Model. The dynamic optimization model for the MMA polymerization process is described by a set of differential equations below:37,38

Yit ≤ |YUYb| − 1 (88)

where YUYb = {i,t|Ybit = 1} contains all the index pair (i,t) if product i is assigned to period t at previous iteration b in the upper level problem, and YUYb = {i,t|Ybit = 0} contains all the K


Article


collocation points are generated in each finite element and the constant matrices css′ and ds are given as follows. To determine the transition time, we have to divide the entire finite elements into two parts. The first part is used to model the transition periods, which have 45 finite elements. The second part is used to model the production periods, which have 30 finite elements. The transition time is the same as the end value of the 45-th finite element.

F(Cmin + Cm) dCm 2f *·kI = −(kp + k fm) Cm CI + V dt k Td + k Tc (92)

FI ·C Iin − F ·CI dCI = −kI ·CI + V dt

(93)

dD0 2f *·kI 2f *·kI = (0.5k Tc + k Td) CI + k fm Cm CI dt k Td + k Tc k Td + k Tc −

F ·D0 V

dD1 2f *·kI F ·D1 = Mm(kp + k fm) Cm CI − kp dt k Td + k Tc V y=

D1 D0

u = FI

⎡ 88 − 7 6 ⎢ 360 ⎢ ⎢ 296 + 169 6 css = ⎢ ′ ⎢ 1800 ⎢ ⎢ 16 − 6 ⎢⎣ 36

(94)

(95)

⎡ 16 − 6 ds = ⎢ ⎣ 36

(96) (97)

296 − 169 6 1800 88 + 7 6 360 16 + 6 36

16 + 6 36

−2 + 3 6 ⎤ ⎥ 225 ⎥ ⎥ −2 − 3 6 ⎥ ⎥ 225 ⎥ 1 ⎥ ⎥⎦ 9

1⎤ ⎥ 9⎦

According to the general approach to build the metamodel that is presented in section 4, the first step is to find the minimum transition time. The minimum transition time is obtained by solving the “Min_Time” problem given in section 4.2. The results for the dynamic optimization problem that minimizes the transition time in transition A to D are shown in Figure 7.

There are four state variables in this model: Cm is the monomer concentration, CI denotes the initiator concentration, D0 stands for the molar concentration of dead chains, and D1 indicates the mass concentration of dead chains. The output of the dynamic model is the number-average molecular weight for each grade, denoted by y, and u is the manipulated variable that is the flow rate of the initiator. The products are produced in the steady-state production period in a CSTR, and different products have different molecular weights. All kinds of steady-state values can be calculated in advance according to the different products.20 We should note that we assume the quality of the raw materials is constant in this work. However, in the real production process the quality of the raw materials is in general not a constant. If the quality of the raw materials is a given dynamic parameter, our model can also be modified to take the time-dependent quality function into consideration. The quality functions should be inserted into the dynamic models. The dynamic optimization problems that have the quality functions are solved off-line again to obtain the new candidate transition times and transition costs. The integrated planning and scheduling problem with the new flexible recipes is then solved to determine the production strategy accordingly. The planning horizon can be separated into several equallength planning periods, and the order demands for the products have to be satisfied at the end of each planning period. The scheduling problem is to determine the production sequence and the production time within each planning period according to the production assignment and the order demands. The transition from one product to another can vary dramatically if different transition trajectories are adopted. Planning, scheduling, and dynamic optimization are taken into consideration simultaneously to yield the maximum profit. 6.2. Offline Metamodeling. In this section, the dynamic process from product A to product D is chosen to demonstrate the procedure of building a metamodel. We use the collocation method to solve the dynamic optimization problems.44,46 The number of finite elements is set as 75, which is a key parameter for the collocation method and is chosen by a trial and error approach over all transitions. In this work, the Radau IIA method is used, which shows a robust performance for stiff differential equations.59 In the Radau IIA method, three

Figure 7. Output profile for off-line dynamic optimization of minimum transition time from A to D.

The solution and model statistics for the dynamic optimization problem that minimizes the transition time from product A to product D are listed in Table 1. The dynamic optimization problem that minimizes the transition time is a NLP problem. According to eqs 63 and 64, the steady-state conditions of product A are the initial conditions of the state Table 1. Solution and Model Statistics for the Dynamic Optimization Problem That Minimizes the Transition Time in Transition A to D

L

Optimizer

CONOPT 3

BARON 12

Objective (h) CPU (s) GAP (%) Type Con. Var. Constraints

1.36 0.69 0

1.36 36,000 2,561.19 NLP 3,380 3,563 DOI: 10.1021/ie503857r Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research variables and the steady-state conditions of product D are the final values of the outputs. In Table 1, the dynamic optimization problem of determining the minimum transition time between two different products contains 3,380 continuous variables and 3,563 constraints. The dynamic optimization problem is highly nonlinear and nonconvex. The first optimizer used to solve this problem is CONOPT 3, which is a local optimization solver. The objective value returned by the local optimizer is 1.36 s, which cannot be guaranteed to be the global solution. The computational time by CONOPT 3 is 0.69 CPUs. The global optimizer BARON 12 is also used to solve this dynamic optimization problem. However, BARON 12 fails to converge to the global optimal solution after running for 36,000 s and the relative optimality gap is 2,561.61%. BARON 12 returns the same objective value as CONOPT 3 does, which is also unable to be guaranteed to be the global solution because of the big relative gap. The above comparison shows that the computational efficiency of BARON 12 is relatively low though it is a global optimizer. In most cases of solving the dynamic optimization problem, BARON 12 fails to converge to the global optimal solution within 36,000 s. Theoretically, if the computational time is long enough (for example, several weeks or months), the global optimizer BARON 12 can finally return a global optimal solution. Practically, for a single dynamic optimization problem, a local optimizer that is able to return a local optimal solution within seconds is preferred because of its computational efficiency and relatively good solution quality. Consequently, CONOPT 3 is used to solve the dynamic optimization problems to build the metamodels. The “minimum transition time” recipes for all the possible transition processes are listed in Table 2. After the minimum transition time is determined, the next step is to determine the interval step size between two adjacent sampling points and the number of all sampling points. Via a trial-and-error procedure, the interval step size is set as 0.1 h and the number of the sampling points is set as 16. Thus, for each possible transition,

the upper bound of the transition time is 1.6 h longer than the lower bound of the minimum transition time. We note that the interval step size and the number of the sampling points depend on the specific case study, and different case studies should have different values. At each sampling point, the “Min_Cost” problem presented in section 4.2 is solved to obtain the minimum transition cost. A set of discrete transition times and transition costs is then obtained from solving the dynamic models. With those discrete pairs, the corresponding metamodel from product A to product D is plotted in Figure 8.

Figure 8. Metamodel from product A to D.

In this work, the metamodel is a two-dimensional model and represents the relationship between the transition cost and transition time. We note that it is not a “nonlinear regression” in Figure 8, since we do not use regression approaches in this work. In Figure 8 these sampling points are “exact” and “feasible” recipes for the dynamic model, and the flexible recipe method only chooses the task recipe from the sampling points. The above approach to build a metamodel is then applied to all potential transitions to build all metamodels. If the number of all the products is np and a transition process is allowed between any two different products, the number of metamodels we need to build is np × (np − 1). In one metamodel there are g sampling points. At each sampling point, a dynamic optimization problem is solved to obtain the transition time and cost. In the large scale case study of this work (section 6.4), five products are produced during the continuous process. For each metamodel, 16 sampling points are chosen to build one metamodel. The total number of the dynamic optimization problems to be solved is 5 × (5−1) × 16 = 320. If we use the global optimizer BARON 12 to solve all the dynamic optimization problems and set 10 h as the computational time limit for each dynamic optimization problem, it is likely to take more than 4 months to build all the metamodels. Consequently, CONOPT 3 is used to solve the dynamic optimization problems to build the metamodels. 6.3. Small Scale Example. We first consider a small scale problem which includes three products (A, B, C). The entire planning horizon is divided into three planning periods. The time length of each planning period is 24 h. Each planning period is separated into three time slots. The number of the time slots in each period is equal to the total number of the products. The demands of each product are shown in Table 3. The production cost, the sale price, the inventory cost, and the production rate for each product are listed in Table 4. The predetermined tolerance for the difference between the upper and lower levels is set as zero. The problem is easy to solve due to its small scale. Table 5 presents the problem sizes, the computational times, and the objective values for all the three problems. For the bilevel decomposition method, integer and logic cuts are added into

Table 2. “Minimum Transition Time” Recipes for All the Possible Transition Processes From

To

Transition time (h)

Transition cost ($)

A A A A B B B B C C C C D D D D E E E E

B C D E A C D E A B D E A B C E A B C D

0.44 0.71 1.36 0.34 0.38 0.33 0.49 0.32 0.41 0.19 0.26 0.34 0.45 0.19 0.11 0.39 0.33 0.36 0.61 0.76

685.0 667.8 756.4 740.0 800.5 619.5 614.1 735.7 814.7 642.4 605.7 747.6 830.2 646.5 614.6 761.6 774.3 656.3 680.7 641.8 M


Article


space integrated problem only obtain a local optimal solution within a reasonable computational time limit. We should note that since the integrated planning and scheduling problem with flexible recipes and the bilevel decomposition method use the metamodels which contain information obtained by a local optimizer at a set of discrete points, the objectives by those two methods are also local optimal solutions. Table 6 shows the upper bound and lower bound at each iteration of the bilevel decomposition method. The bilevel

Table 3. Demands for Each Product at Each Period for Case Study 1 Demand (m3)

Period 1

Period 2

Period 3

A B C

5 1 4

4 3 4

4 4 1

Table 4. Cost Data and Production Rate for Case Study 1 Product

Operation cost) ($/m3)

Price ($/m3)

Production rate (m3/h)

Inventory cost ($/m3·h)

A B C

263 188 163

288.5 232.8 217.8

0.8 0.9 1.1

0.1 0.1 0.1

Table 6. Upper and Lower Bounds at Each Iteration of the Bilevel Decomposition Algorithm for the Small Case Study

the upper level problem in every iteration. Thus, the size of the upper level problem increases during the iteration process. In Table 5, only the problem size of the upper level problem at the last iteration is shown. The computational time of the bilevel decomposition algorithm is regarded as the sum of the computational time of both the upper and the lower level problems. The number of constraints of the full space problem is 2 orders of magnitude larger than that of the other two problems. The number of continuous variables of the full space problem is more than 1 order of magnitude larger than that of the other two problems. The number of binary variables of the integrated planning and scheduling problem with flexible recipes and the lower level problem of the bilevel decomposition algorithm is larger than that of the full space method, because extra binary variables are introduced to select the pairs of the transition times and costs from the flexible recipes. In Table 5, the solution and model statistics for the small scale problem is listed. The full space problem is first solved by a local optimizer, SBB. The objective is $1,325.42, which is returned by SBB after 13,430.76 CPUs. The full space problem is then solved by a global optimizer, BARON 12. However, BARON 12 fails to converge to a global optimal solution within 36,000 CPUs. The relative gap is 53.09%. The objective by BARON 12 is $1,325.42, which is the same as the objective by SBB and cannot be guaranteed to be a global optimal solution because of the existing relative gap. Both the integrated planning and scheduling problem with flexible recipes and the bilevel decomposition method yield the same solution. The objective of those two methods is $1,387.41, which is obtained within 0.5 s for both methods. The optimality gap is set as zero. The objective values of the flexible recipe method and the bilevel decomposition method are a bit larger than that of the full space method. The reason is that the optimizers for the full

a

Iteration

Upper bound

Lower bound

1 2 3

1,562.05 1,524.12 1,222.30

657.15 1,387.42a 1,387.42

The optimal solution.

decomposition algorithm takes 3 iterations to converge. And the solution of the bilevel decomposition algorithm appears at the second iteration. In the bilevel decomposition method, at the beginning of one new iteration, the integer and logic cuts are added to the upper level problem to exclude previous feasible solutions and the subsets and supersets of the previous feasible solutions. If the iteration does not stop, all the feasible solutions will be excluded gradually. Thus, the upper bound always decreases as the iteration goes, while the lower bound of the bilevel decomposition method is the maximum of all the solutions from the lower level problem at each iteration. To compare the three different approaches to calculate the inventory cost, we choose the inventory information returned by the flexible recipe method and the bilevel decomposition method. The inventory cost by the newly proposed estimated approach is $123.84 throughout the planning horizon, while the real inventory cost for all product is $120.11. The estimation error between the real inventory cost and the estimation is 3.11% of the real inventory cost and only 0.26% of the total profit. The estimated inventory cost by the estimated approach in previous work12,16−19,35,36 is $209.28 which has an estimation error of 74.24% of the real inventory cost and 6.43% of the total profit. The comparison result demonstrates the advantage of the proposed estimated approach, which reduces the approximation error significantly compared with the approach in previous work and also gives a good approximation of the real inventory cost.12,16−19,36 The same results are obtained by the flexible recipe method and the bilevel decomposition method, as shown in Figure 9 (I). The result shows that the production sequence is A-C-B-C.

Table 5. Comparison of the Three Methods for Solving the 3-Product and 3-Period Problem Method

Full space method

Flexible recipe method

Type Bin. Var. Con. Var. Constraints Optimizer Obj. ($) CPU (s) Iteration Gap (%)

MINLP 68 27,284 28,874

MILP 900 1,048 261 CPLEX 12 1,387.41 0.35 − 0

SBB 1,325.42 13,430.76 − 0

BARON 12 1,325.42 36,000 − 53.09% N

Bilevel decomposition algorithm MILP 9 (Upper) 66 (Upper) 99 (Upper)

900 (Lower) 1,051 (Lower) 266 (Lower)

CPLEX 12 1,387.41 0.34 = 0.16 (Upper) + 0.18 (Lower) 3 0 DOI: 10.1021/ie503857r Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Figure 9. Result of the small case study of the flexible recipe method and the bilevel decomposition algorithm. (I) The result over the entire planning horizon; (II) The transition trajectory from product A to C; (III) The transition trajectory from product C to B; (IV) The transition trajectory from product B to C.

Figure 10. Result of the small case study from the full space method. (I) The result over the entire planning horizon; (II) The transition trajectory from product A to C; (III) The transition trajectory from product C to B; (IV) The transition trajectory from product B to C.

The detailed trajectories of the transitions in Figure 9 (I) are shown in (II) − (IV). Both the results of the full space method returned by SBB and by BARON 12 are the same, as shown in Figure 10.

Comparing the results in Figure 10 (I) with those in Figure 9 (I), they share the same production sequence from the overall perspective but have slight differences. For example, the transition profile in the first transition from product A to C O


Article


The size of the bilevel decomposition problem shown in Table 9 is from the last iteration. The total computational time of the bilevel decomposition method is the sum of the computational time from both the upper and the lower level problems. The number of constraints of the full space method is 2 orders of magnitude larger than that of the other two methods. The number of continuous variables of the full space method is close to 10 times larger than that of the other two methods. Similar to the small scale case study, the number of binary variables of either the flexible recipe method or the lower level problem of the bilevel decomposition method is larger than that of the full space method due to the task recipe selection. The upper bound and the lower bounds at each iteration of the bilevel decomposition method for the large scale problem are shown in Figure 11. The bilevel decomposition method

in Figure 10 (II) is different from that in Figure 9 (II). Because there is a slight fluctuation at the end of the transition process in Figure 10 (II) compared with the smooth transition process in Figure 9 (II), the transition cost from product A to C returned by the full space method is a bit larger than the transition costs from product A to C returned by the other two methods. The difference comes from the nonconvexity of the MINLP problem, which results in a local solution, as we have discussed. 6.4. Large Scale Example. The second case study is a large scale problem that contains five products. The entire planning horizon is 1 month, which is separated into four planning periods. The time length of each planning period is a week (168 h), which is divided into five slots. The number of the slots in each period is the same as the number of all the products. The demands for each product are listed in Table 7. The operation Table 7. Demands for Each Product at Each Period for Case Study 2 Demand (m3)

Period 1

Period 2

Period 3

Period 4

A B C D E

0 16.7 15 20 22.2

20 8.3 10 10 33.3

20 2 5 2 44.4

0 15 15 20 22.2

cost, the sale price, the inventory cost, and the production rate are predetermined parameters shown in Table 8. The predetermined tolerance for the bilevel decomposition algorithm is set as zero.

Figure 11. Upper and lower bounds at each iteration of the bilevel decomposition algorithm for the large case study.

Table 8. Cost Data and Production Rate for Case Study 2 Product

Operation cost ($/m3)

Price ($/m3)

Production rate (m3/h)

Inventory cost ($/m3·h)

A B C D E

263 188 163 226 220

288.5 232.8 217.8 293 330.5

0.8 0.9 1.1 1.2 1

0.1 0.1 0.1 0.1 0.1

takes 46 iterations to meet the stopping condition. The optimal solution appears at the 8-th iteration in the lower level problem. In Table 9, we can see that the total computational time of the upper level problem is 13.5 s while the total computational time of the lower level is 189.2 s. The reason is because the upper level problem is a relaxed MILP problem while the lower level problem is a detailed integrated planning and scheduling problem with flexible recipes. The relaxed problem is much easier to solve compared with the detailed integrated problem. The detailed production sequence and production time is shown in Figure 12 (I), as well as the first two detailed transition trajectories. The comparison of the three different approaches to calculate the inventory cost is presented as follows. The inventory cost by the newly proposed estimated approach is $6,973.16 throughout the planning horizon, while the real

The results are shown in Table 9. Both the full space methods by SBB and by BARON 12 fail to obtain a feasible solution within a computational time limit of 36,000 CPUs. The flexible recipe method uses 1,020.96 CPUs to obtain the optimal solution. The bilevel decomposition method exhibits its superiority in the large scale problem. It only uses 202.7 s to obtain the same optimal solution as the flexible recipe method does.

Table 9. Comparison of the Three Methods Solving the 5-Product and 4-Period Problem Method

Full space method

Flexible recipe method

Type Bin. Var. Con. Var. Constraints Optimizer Obj. ($) CPU (S) Iteration Gap (%)

MINLP 234 64,874 68,736 SBB/BARON 12 Infeasible 36,000 − −

MILP 6,520 6,928 971 CPLEX 12 51,151.82 1,020.96 − 0 P

Bilevel decomposition algorithm MILP 20(Upper) 135(Upper) 1,532(Upper)

6,520(Lower) 6,928(Lower) 1,016(Lower)

CPLEX 12 51,151.82 202.7 = 13.5(Upper) + 189.2(Lower) 46 0 DOI: 10.1021/ie503857r Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


algorithm solved the upper and the lower level problems alternately. The integer and logics cuts were added to the upper level problem at each iteration, until the predetermined tolerance was reached. The model formulations proposed in this work are general and valid for all the integrated planning, scheduling, and dynamic optimization problems of continuous processes, and our proposed efficient solution methods are to address the computational complexity of the full space MINLP problem of the integrated planning, scheduling, and dynamic optimization problem. So our proposed solution approaches are feasible for all the integrated planning, scheduling, and dynamic optimization problems. The advantages of the proposed solution methods and the estimated approach were demonstrated via two case studies of the MMA polymerization processes. The proposed methods reduced the computational time by more than 2 orders of magnitude compared with directly solving the full space integrated problem. In the small scale problem of three products and three planning periods, both the flexible recipe method and the bilevel decomposition algorithm solved the integrated problem within half a second, while the full space method took several hours to obtain a solution. In a large scale problem of five products and four planning periods, the bilevel decomposition algorithm solved the problem with 202.7 s and the flexible recipe method solved the problem with 1,020.96 s, while the full space method failed to obtain a feasible solution within 10 h. The result also demonstrated that the novel linear approach to approximate the inventory cost reduces the estimation error significantly compared with the existing approach.12,16−19,35,36 Disturbance is one important and general issue in continuous production processes.23 The focus of our work is on the dynamic optimization problem, which is solved to determine the optimal transition trajectory for the lower-level advanced controllers. Process disturbances would be handled by closeloop or feedback controllers,38,60,61 which track the optimal transition trajectory determined by the proposed framework. Additionally, when disturbance occurs and the transition time changes, the integrated planning and scheduling problem with a flexible recipe can be resolved with updated parameters to determine the new production strategy, which will be passed to the process controller to control the continuous manufacturing process. For online application, efficient computational efficiency and fast computational time are critical. As we can see in the case study section, for the small scale production process, our proposed approach is fast enough for the online application since the production strategy can be determined within seconds, while, for the medium scale production process, our proposed approach takes hundreds of seconds to determine the production strategy. For larger scale production processes with 30−40 different products, more efficient solution methods might be needed for online application, and this would be a direction of future research.

Figure 12. Result of the large case study of the flexible recipe and bilevel decomposition method. (I) The result over the entire planning horizon; (II) The transition trajectory from product D to C; (III) The transition trajectory from product C to B.

inventory cost for all products is $6,905.85. The estimation error between the real inventory cost and the estimation is 0.97% of the real inventory cost and only 0.13% of the total profit. The inventory cost by the estimated approach in previous work is $12,580.49. The estimation error is 82.17% of the real inventory cost and 11.09% of the total profit. In this case study, we can see the approximation error of the proposed estimated approach goes below 1%. The result demonstrates the proposed estimated approach is a better approach of the inventory cost compared to the approach in previous work. The same results are obtained by the flexible recipe method and the bilevel decomposition method, as shown in Figure 12 (I). The result shows that the production sequence is E-B-D-CA-E-B-D. The detailed trajectories of the first two transitions in Figure 12 (I) are shown in Figure 12 (II) and (III).

7. CONCLUSION Integration of planning, scheduling, and dynamic optimization for continuous manufacturing processes is of great importance, because it could lead to a better overall performance. In this work we proposed a novel planning model, in which a new linear approach to approximate the inventory cost reduced the estimation error significantly. With this new planning model, a full space MINLP model for the integrated planning, scheduling, and dynamic optimization problem of the multiproduct continuous manufacturing process was then presented. To address the computational complexity of the full space MINLP problem, we took advantage of a set of metamodels and proposed an efficient flexible recipe method, which decomposed the full space model into an integrated planning and scheduling problem with flexible recipes, and a number of dynamic optimization problems for generating the metamodels. The bilevel decomposition algorithm was then used to further improve the computational efficiency for solving the integrated planning and scheduling problem with flexible recipes. The

■

APPENDIX The mathematic approach to calculate the real inventory cost is presented first as follows. Equations A1 and A2 define a 0−1 binary variable ZWiklt, which is equal to 1 if product i is assigned to slot l of period t and a different product k is assigned to the subsequent slot. Q


Article

Industrial & Engineering Chemistry Research ZWiklt ≥ Wilt + Wk , l + 1, t − 1,

∀ i , k ≠ i , l < Ns , t

t planning period b iteration

(A1)

Sets

ZWiklt ≥ Wilt + Wk , ll , t + 1 − 1, ∀ i , k ≠ i , t , l = Ns , ll = 1

I set of all the products indexed by i YLYb index collection that contains all the index pair (i,t) if product i is assigned to period t at previous iteration b in the lower level problem YLNb index collection that contains all the index pair (i,t) if product i is not assigned to period t at previous iteration b in the lower level problem YUYb index set that contains all the index pair (i,t) if product i is assigned to period b at previous iteration k in the upper level problem YUNb index collection that contains all the index pair (i,t) if product i is not assigned to period t at previous iteration b in the upper level problem

(A2)

TEIilt,

Equation A3 defines the end time of one time slot l of period t to which product i is assigned. TEiltI = Wilt ·TElt

(A3)

TEIIit ,

Equation A4 defines the end time of the production process of product i in period t. If product i is not assigned to period t, TEIIit would be 0. TEitII =

TEiltI ·ZWiklt + TEiI, l = ns , t

∑

(A4)

k ≠ i , l < ns

Parameters

Equations A5 and A6 define the real inventory area of product i in period t, INAit. Since both production time PTilt and production amount Git are continuous variables, each PTilt· Git is a bilinear term and eqs A5 and A6 are nonlinear equations. INAit = inviInitial i ·ht +

∑l PTilt ·Git 2

bss′ cinv i coper it dit ht htt invInitial i bm nc

+ (htt − TEitII ) ·Git ,

∀ i, t = 1

(A5)

INAit = INVOit ·ht +

∑l PTilt ·Git 2

+ (htt − TEitII ) ·Git ,

∀ i, t ≠ 1

np npro ns nf lt pit qlt ri rcikg

(A6)

In eq A7, the real inventory cost CIVreal is calculated by the sum of the inventory cost of each product in each period. Because nonlinear eqs A5 and A6 are involved in eq A7, the approach to calculate the real inventory cost is a nonlinear approach. CIV real =

∑ ∑ (INAit ·ciinv) t

rtikg

(A7)

i 12,16−19,35,36

The approach in previous work to approximate the inventory area is a linear method; however, it leads to a large estimation error between the estimated inventory cost and the real inventory cost. The existing estimated approach is shown in eq A8: CIV previous =

■

∑ (INR t ·ht ) t

ttmin ik zInitial i ySetpoint i α ε

Continuous variables (A8)

CDlt CTlt

AUTHOR INFORMATION

Corresponding Author

CIV

*Tel: +1 847 467 2943; fax: +1 847 491 3728; E-mail: you@ northwestern.edu.

CIVreal CIVprevious

Notes

The authors declare no competing financial interest.

■

COP CTR INAit INRt

NOMENCLATURE

Indices

i g l r s

constant matrix according to specific collocation method inventory cost for unit product i operation cost for producing unit product i in period t demand of the product i in period t duration of period t ending time for period t initial inventory level of product i a big positive number used in the big-M method total number of candidate pairs from one product to another total number of the planning period total number of all the products total number of the slot in each period total number of finite element in slot l of period t price of product i in period t time-independent parameters in slot l of period t production rate for product i g-th transition cost candidate from product i to product k g-th transition time candidate from product i to product k minimum transition time from product i to product k steady-state value for manufacturing product i output set point of product i time step used in candidate transition time tolerance between upper bound and lower bound

product candidate transition time and cost slot finite element collocation point

INROt INVit R

transition cost in slot l of period t in dynamic modeling transition cost in time slot l of period t in planning modeling inventory cost estimated by the newly proposed approach real inventory cost inventory cost estimated by the approach in previous work total operation cost total transition cost real inventory area of product i in period t inventory cost rate for all products at the end of planning period t before the demand is met inventory cost rate for all products at the beginning of planning period t inventory level of product i DOI: 10.1021/ie503857r Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research XLik

INVCt INVOit INRt

inventory cost for all products in period t inventory level of product i after sale inventory cost rate for all products at the end of planning period t before the demand is met INROt initial inventory cost rate for all products at the beginning of planning period t Llt length of the finite element at slot l and period t LB lower bound for bilevel decomposition method LBb lower bound for bilevel decomposition method at iteration b MITCi minimum transition cost of product i MATCt maximum of minimum transition costs of products assigned to period t NYit number of slots used to produce product i in period t Git production amount of product i produced in period t PTilt production time of product i produced in slot l of period t PTUit production time of product i in period t in the upper level problem PROFIT total profit RCik transition cost from product i to product k RVN revenue for all products Sit sale of product i in period t Tlt time variable in the dynamic optimization model of slot l in period t TLlt final value of the time variable in the transition of slot l in period t TLFik final value of the time variable in transition from product i to product k TElt ending time of slot l in period t TEIilt end time of one time slot l of period t to which product i is assigned TEIIit end time of the production process of product i in period t TFik time variable in the dynamic optimization model from product i to product k TTTk total transition time at iteration k TSlt starting time of slot l in period t TTlt transition time of slot l in period t Ult manipulate variables in slot l of period t ULit lower bound of manipulate variables in slot l of period t UUlt upper bound of manipulate variables in slot l of period t ULik lower bound of manipulate variables from product i to product k UUik upper bound of manipulate variables from product i to product k Urslt discretized input variable vector at the s-th collocation point in the r-th finite element UFik input vector in the dynamic optimization model from product i to product k UMt maximum of minimum transition times of products assigned to period t UPP total profit obtained from the upper level problem UB upper bound for the bilevel decomposition method UBb upper bound for the bilevel decomposition method at iteration b Xlt algebraic variables in slot l of period t XLlt lower bound of algebraic variables in slot l of period t XUlt upper bound of algebraic variables in slot l of period t

lower bound of algebraic variables from product i to product k upper bound of algebraic variables from product i to product k discretized state variable vector at the s-th collocation point in the r-th finite element state variable vector in dynamic optimization model from product i to product k production amount of product i in period t output of product i output in slot l of period t lower bound of output set points in slot l of period t upper bound of output set points in slot l of period t lower bound of output set points from product i to product k upper bound of output set points from product i to product k discretized output set points at the s-th collocation point in the r-th finite element output set point of slot l in period t output vector in dynamic optimization model from product i to product k state variables in slot l of period t lower bound of state variables in slot l of period t upper bound of state variables in slot l of period t lower bound of state variables from product i to product k upper bound of state variables from product i to product k initial conditions in slot l of period t discretized state variables at the beginning of finite element r in slot l of period t discretized state variables at collocation point s of finite element r in slot l of period t state variables from product i to product k

XUik Xrslt XFik XPit Yi Ylt YLlt YUlt YLik YUik Yrslt YSet lt YFik Zlt ZLlt ZUlt ZLik ZUik ZIni lt Zrlt Zrslt ZFik

Binary variables

NPit Wilt Yit Ybit YLit YLbit ZBikltg ZWiklt

■

binary variable to denote if product i is produced in period t binary variable to denote if product i is assigned to slot l of period t binary variable to denote if product i is produced in period t binary variable to denote if product i is produced in period t in iteration b binary variable to denote if product i is produced in period t binary variable to denote if product i is produced in period t in iteration b binary variable to denote product i is assigned to slot l of period t, product k is assigned to the following slot, and the g-th transition candidate is selected binary variable to denote product i is assigned to slot l of period t and product k is assigned to the following slot

REFERENCES

(1) Harjunkoski, I.; Nyström, R.; Horch, A. Integration of scheduling and controlTheory or practice? Comput. Chem. Eng. 2009, 33 (12), 1909−1918. (2) Chu, Y.; You, F. Model-based integration of control and operations: Overview, challenges, advances, and opportunities. Comput. Chem. Eng. 2015, Submitted. S


Article

Industrial & Engineering Chemistry Research (3) Capon-Garcia, E.; Moreno-Benito, M.; Espuna, A. Improved Short-Term Batch Scheduling Flexibility Using Variable Recipes. Ind. Eng. Chem. Res. 2011, 50 (9), 4983−4992. (4) Wassick, J. M.; Agarwal, A.; Akiya, N.; Ferrio, J.; Bury, S.; You, F. Q. Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products. Comput. Chem. Eng. 2012, 47, 157−169. (5) Takeda, M.; Ray, W. H. Optimal-grade transition strategies for multistage polyolefin reactors. AIChE J. 1999, 45 (8), 1776−1793. (6) Seborg, D.; Edgar, T. F.; Mellichamp, D. Process dynamics & control; John Wiley & Sons: 2006. (7) Debling, J. A.; Han, G. C.; Kuijpers, F.; Verburg, J.; Zacca, J.; Ray, W. H. Dynamic modeling of product grade transitions for olefin polymerization processes. AIChE J. 1994, 40 (3), 506−520. (8) Grossmann, I. Enterprise-wide optimization: A new frontier in process systems engineering. AIChE J. 2005, 51 (7), 1846−1857. (9) Chu, Y. F.; You, F. Q. Integrated Scheduling and Dynamic Optimization of Complex Batch Processes with General Network Structure Using a Generalized Benders Decomposition Approach. Ind. Eng. Chem. Res. 2013, 52 (23), 7867−7885. (10) Zhuge, J. J.; Ierapetritou, M. G. Integration of Scheduling and Control with Closed Loop Implementation. Ind. Eng. Chem. Res. 2012, 51 (25), 8550−8565. (11) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Lagrangean heuristic for the scheduling and control of polymerization reactors. AIChE J. 2008, 54 (1), 163−182. (12) Erdirik-Dogan, M.; Grossmann, I. E. A decomposition method for the simultaneous planning and scheduling of single-stage continuous multiproduct plants. Ind. Eng. Chem. Res. 2006, 45 (1), 299−315. (13) Birewar, D. B.; Grossmann, I. E. Simultaneous Production Planning and Scheduling in Multiproduct Batch Plants. Ind. Eng. Chem. Res. 1990, 29 (4), 570−580. (14) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45 (20), 6698−6712. (15) Chu, Y.; You, F.; Wassick, J. M.; Agarwal, A. Integrated planning and scheduling under production uncertainties: Bi-level model formulation and hybrid solution method. Comput. Chem. Eng. 2015, 72 (0), 255−272. (16) Kallrath, J. Planning and scheduling in the process industry. OR Spectrum 2002, 24 (3), 219−250. (17) Li, Z. K.; Ierapetritou, M. G. Integrated production planning and scheduling using a decomposition framework. Chem. Eng. Sci. 2009, 64 (16), 3585−3597. (18) Erdirik-Dogan, M.; Grossmann, I. E. Simultaneous planning and scheduling of single-stage multi-product continuous plants with parallel lines. Comput. Chem. Eng. 2008, 32 (11), 2664−2683. (19) Yue, D. J.; You, F. Q. Planning and Scheduling of Flexible Process Networks Under Uncertainty with Stochastic Inventory: MINLP Models and Algorithm. AIChE J. 2013, 59 (5), 1511−1532. (20) Chu, Y. F.; You, F. Q. 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. (21) Chu, Y. F.; You, F. Q. Integrated Scheduling and Dynamic Optimization of Sequential Batch Processes with Online Implementation. AIChE J. 2013, 59 (7), 2379−2406. (22) Chu, Y. F.; You, F. Q. Integration of production scheduling and dynamic optimization for multi-product CSTRs: Generalized Benders decomposition coupled with global mixed-integer fractional programming. Comput. Chem. Eng. 2013, 58, 315−333. (23) Chu, Y. F.; You, F. Q. Integration of Scheduling and Dynamic Optimization of Batch Processes under Uncertainty: Two-Stage Stochastic Programming Approach and Enhanced Generalized Benders Decomposition Algorithm. Ind. Eng. Chem. Res. 2013, 52 (47), 16851−16869.

(24) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous Scheduling and Control of Multiproduct Continuous Parallel Lines. Ind. Eng. Chem. Res. 2010, 49 (17), 7909−7921. (25) 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 (11), 3416−3432. (26) 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. (27) Muñoz, E.; Capón-García, E.; Moreno-Benito, M.; Espuña, A.; Puigjaner, L. Scheduling and control decision-making under an integrated information environment. Comput. Chem. Eng. 2011, 35 (5), 774−786. (28) Chu, Y.; You, F. Integrated scheduling and dynamic optimization by Stackelberg game: Bi-level model formulation and efficient solution algorithm. Ind. Eng. Chem. Res. 2014, 53 (13), 5564− 5581. (29) Bansal, V.; Sakizlis, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. New algorithms for mixed-integer dynamic optimization. Comput. Chem. Eng. 2003, 27 (5), 647−668. (30) Allgor, R. J.; Barton, P. I. Mixed-integer dynamic optimization I: problem formulation. Comput. Chem. Eng. 1999, 23 (4−5), 567−584. (31) Allgor, R.; Barton, P. Mixed-integer dynamic optimization I: problem formulation. Comput. Chem. Eng. 1999, 23 (4), 567−584. (32) Chachuat, B.; Singer, A. B.; Barton, P. I. Global mixed-integer dynamic optimization. AIChE J. 2005, 51 (8), 2235−2253. (33) Biegler, L. T. Nonlinear programming: concepts, algorithms, and applications to chemical processes; SIAM: 2010; Vol. 10. (34) Li, X.; Tomasgard, A.; Barton, P. I. Nonconvex generalized Benders decomposition for stochastic separable mixed-integer nonlinear programs. J. Opt. Theory Appl. 2011, 151 (3), 425−454. (35) Tlacuahuac, A. F.; Gutierrez-Limon, M.; Grossmann, I. E. A MINLP formulation for the Simultaneous Planning, Scheduling and Control of Short-Period Single Unit Processing Systems. Ind. Eng. Chem. Res. 2014, 53, 14679. (36) Chu, Y.; You, F. Integrated Planning, Scheduling, and Dynamic Optimization for Batch Processes: MINLP Model Formulation and Efficient Solution Methods via Surrogate Modeling. Ind. Eng. Chem. Res. 2014, 53 (34), 13391−13411. (37) Congalidis, J. P.; Richards, J. R.; Ray, W. H. Feedforward and Feedback-Control of a Solution Copolymerization Reactor. AIChE J. 1989, 35 (6), 891−907. (38) Daoutidis, P.; Soroush, M.; Kravaris, C. Feedforward FeedbackControl of Multivariable Nonlinear Processes. AIChE J. 1990, 36 (10), 1471−1484. (39) Floudas, C. A.; Lin, X. X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Comput. Chem. Eng. 2004, 28 (11), 2109−2129. (40) Schilling, G.; Pantelides, C. C. A simple continuous-time process scheduling formulation and a novel solution algorithm. Comput. Chem. Eng. 1996, 20, S1221−S1226. (41) Mitra, K.; Gudi, R. D.; Patwardhan, S. C.; Sardar, G. Resiliency Issues in Integration of Scheduling and Control. Ind. Eng. Chem. Res. 2010, 49 (1), 222−235. (42) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and optimal control of polymerization reactors. AIChE J. 2007, 53 (9), 2301−2315. (43) Chu, Y. F.; You, F. Q. Moving horizon approach of integrating scheduling and control for sequential batch processes. AIChE J. 2014, 60 (5), 1654−1671. (44) Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process. 2007, 46 (11), 1043−1053. (45) Ogunnaike, B. A.; Ray, W. H. Process dynamics, modeling, and control; Oxford University Press: New York, 1994; Vol. 9. (46) Cuthrell, J. E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33 (8), 1257−1270. (47) Nystrom, R. H.; Franke, R.; Harjunkoski, I.; Kroll, A. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29 (10), 2163−2179. T


Article

Industrial & Engineering Chemistry Research (48) Nystrom, R. H.; Harjunkoski, I.; Kroll, A. Production optimization for continuously operated processes with optimal operation and scheduling of multiple units. Comput. Chem. Eng. 2006, 30 (3), 392−406. (49) 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 (11), 4022−4034. (50) Nie, Y.; Biegler, L. T.; Villa, C. M.; Wassick, J. M. Discrete Time Formulation for the Integration of Scheduling and Dynamic Optimization. Ind. Eng. Chem. Res. 2014. (51) Gorissen, D.; Couckuyt, I.; Demeester, P.; Dhaene, T.; Crombecq, K. A surrogate modeling and adaptive sampling toolbox for computer based design. J. Machine Learning Res. 2010, 11, 2051− 2055. (52) Ong, Y. S.; Nair, P. B.; Keane, A. J. Evolutionary optimization of computationally expensive problems via surrogate modeling. AIAA J. 2003, 41 (4), 687−696. (53) Erdirik-Dogan, M.; Grossmann, I. E.; Wassick, J.; Bi-level, A. Decomposition Scheme for the Integration of Planning and Scheduling in Parallel Multi-Product Batch Reactors. 17th Eur. Symp. Comput. Aided Process Eng. 2007, 24, 625−630. (54) Iyer, R. R.; Grossmann, I. E. A bilevel decomposition algorithm for long-range planning of process networks. Ind. Eng. Chem. Res. 1998, 37 (2), 474−481. (55) You, F. Q.; Grossmann, I. E.; Wassick, J. M. Multisite Capacity, Production, and Distribution Planning with Reactor Modifications: MILP Model, Bilevel Decomposition Algorithm versus Lagrangean Decomposition Scheme. Ind. Eng. Chem. Res. 2011, 50 (9), 4831− 4849. (56) Lima, R. M.; Grossmann, I. E.; Jiao, Y. Long-term scheduling of a single-unit multi-product continuous process to manufacture high performance glass. Comput. Chem. Eng. 2011, 35 (3), 554−574. (57) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldivar-Guerra, E. Dynamic optimization of HIPS open-loop unstable polymerization reactors. Ind. Eng. Chem. Res. 2005, 44 (8), 2659−2674. (58) Richard, E. GAMS-A User’s Guide; GAMS Development Corporation: Washington, DC, USA, 2013. (59) Hairer, E.; Wanner, G. Stiff and Differential-algebraic Problems; Hairer, E.; Wanner, G. Springer: 1996. (60) Doyle, J. C.; Stein, G. Multivariable feedback design: concepts for a classical/modern synthesis; IEEE Trans. on Auto. Control; 1981; Citeseer: 1981. (61) Yue, D.; Han, Q.-L.; Peng, C.State feedback controller design of networked control systems. Proceedings of the 2004 IEEE International Conference on Control Applications; IEEE: 2004; pp 242−247.

U


Novel Optimization Model and Efficient Solution Method for

Recommend Documents