Integration of Scheduling and Dynamic ... - ACS Publications

Article pubs.acs.org/IECR

Integration of Scheduling and Dynamic Optimization of Batch Processes under Uncertainty: Two-Stage Stochastic Programming Approach and Enhanced Generalized Benders Decomposition Algorithm Yunfei Chu and Fengqi You* Department of Chemical & Biological Engineering, Northwestern University, Evanston, Illinois 60208, United States ABSTRACT: Integration of scheduling and dynamic optimization significantly improves the overall performance of a production process compared to the traditional sequential method. However, most integrated methods focus on solving deterministic problems without explicitly taking process uncertainty into account. We propose a novel integrated method for sequential batch processes under uncertainty. The integrated problem is formulated into a two-stage stochastic program. The first-stage decisions are modeled with binary variables for assignment and sequencing while the second-stage decisions are the remaining ones. To solve the resulting complicated integrated problem, we develop two efficient algorithms based on the framework of generalized Benders decomposition. The first algorithm decomposes the integrated problem according to the scenarios so that the subproblems can be optimized independently over each scenario. Besides the scenario decomposition, the second algorithm further decomposes dynamic models from the scheduling model, resulting in a nested decomposition structure. For a complicated case study with more than 3 million variables/equations under 100 scenarios, the direct solution approach does not find a feasible solution while the two decomposition algorithms return the optimal solution. The computational time for the first algorithm is 23.9 h, and that for the second algorithm is only 3.3 h. Furthermore, the integrated method returns a higher average profit than the sequential method by 17.6%.

1. INTRODUCTION With the advance in enterprisewide optimization,1−4 an increasing number of integrated methods for scheduling and dynamic optimization have emerged in the recent decade. They are applied to various continuous processes5−10 or batch processes.11−15 In an integrated method, scheduling decisions are optimized collaboratively with task operational recipes characterized by dynamic models. Compared to the traditional methods which sequentially solve the dynamic optimization problems and the scheduling problem, the integrated methods significantly improve the overall performance of the entire production process.16−18 However, the integrated method needs to solve a complicated mixed-integer dynamic optimization (MIDO) problem19,20 or a large-scale mixed-integer nonlinear program (MINLP) after discretizing the dynamic models21,22 that could be computationally intractable when the integrated problem is complicated. Owing to the computational complexity, most integrated methods concentrate on deterministic problems where all parameters are assumed to be known with certainty. The assumption is, however, not realistic considering that uncertainty is inevitable in a real-world process, which can soon deviate the production process from what is planned. Under uncertainty, the predetermined optimal solution can quickly become suboptimal or even infeasible.23−26 Though important, the issue of uncertainty is only considered by relatively few literature reports.5,9,14 All these methods adopt a reactive approach where the knowledge of uncertainty is not explicitly taken into account in the model formulation and the solution algorithm. Uncertainty is coped with by resolving the integrated problem © 2013 American Chemical Society

online after observing its realizations. To ensure an instant response, the reactive methods require a very fast online algorithm which is much more computational demanding than the already time-consuming offline algorithm. In addition, the reactive methods may considerably alter the entire offline solution in the response to just a minor uncertain outcome. The dramatic change in the scheduling decisions can cause problems such as shop floor nervousness,27−30 which disrupts personnel assignment, material delivery, and coordination with subsequent processes. After uncertainty realizations are observed, the reactive approach actually solves a deterministic problem with parameters updated according to the uncertainty realizations. By contrast, the uncertainty information can be explicitly formulated in an optimization problem. Robust optimization describes uncertainty by a range of uncertain parameters, and the optimization problem is solved over the entire uncertainty range to ensure that the solution quality in the worst case is still acceptable.31,32 Chance constraint optimization formulates uncertainty into probabilistic constraints such that the probability of violating these constraints under uncertainty is less than a given threshold.33,34 Two-stage stochastic programming describes uncertain parameters by random variables, and the mean objective function value is optimized.35,36 Received: Revised: Accepted: Published: 16851

August 9, 2013 October 3, 2013 October 18, 2013 October 18, 2013 dx.doi.org/10.1021/ie402621t | Ind. Eng. Chem. Res. 2013, 52, 16851−16869

Industrial & Engineering Chemistry Research

Article

Figure 1. Two decomposition algorithms for solving the complicated integrated problem under uncertainty.

two-stage stochastic programming is suitable for problems with a hierarchical structure, e.g. integrated process design and dynamic optimization,37 design and planning,38 and planning and scheduling.39 As the integrated scheduling and dynamic optimization problem also has a hierarchical structure, we formulate it into a two-stage stochastic nonlinear program. For simplicity, we restrict this work to the sequential batch processes40−42 where material splitting and mixing are not allowed. The integrated model is formulated by combining the general precedence

A feature of two-stage stochastic programming is its proactivereactive decision making to address uncertainty. All variables are clustered into two groups: the first-stage variables, which are determined offline before uncertainty realizations are observed, and the second-stage variables, which are determined online after the uncertain realizations are observed. The first-stage decision variables are proactive in hedging against uncertainty and ensuring that the solution performs well under all uncertainty realizations. The second-stage variables are reactive in responding to the observed realization of uncertainty. The 16852

dx.doi.org/10.1021/ie402621t | Ind. Eng. Chem. Res. 2013, 52, 16851−16869


Article

scheduling model43−45 containing the order value function14,46 and dynamic models discretized by the collocation method.21,22 The knowledge of uncertainty is characterized by a limited number of scenarios, and the objective is to maximize the production profit on average over the scenarios. Unlike the purely reactive methods, the proposed method explicitly considers uncertainty in the model formulation. The first-stage decisions are chosen as the binary assignment and sequencing decisions in the scheduling problem while the second-stage variables are the remaining ones in the scheduling model and all dynamic models. Formulating the integrated problem into a two-stage stochastic program has some advantages. (1) The two-stage formulation is consistent with the hierarchical structure of scheduling and dynamic optimization, facilitating implementation of the integrated method. The major scheduling decisions of task assignment and sequencing are determined offline. Thus, only a simple reduced problem needs to be solved online, significantly reducing the online computational burden. (2) The two-stage formulation can prevent a dramatic change in the schedule, as the task assignment and sequencing variables are not changed online under different uncertainty realizations. This is a desirable property for determining a schedule in practice, avoiding the problem of shop floor nervousness.27−30 Though being solved offline, the formulated two-stage stochastic nonlinear program turns out to be a large-scale MINLP which might be intractable for general-purpose MINLP solvers. We propose two decomposition algorithms based on the framework of generalized Benders decomposition (GBD),47,48 which are shown in Figure 1. The first algorithm decomposes the integrated problem according to the scenarios. The master problem is a mixed-integer linear program (MILP) which is solved to determine the binary assignment and sequencing variables. The primal problem is a set of separable restricted problems over every scenario. Each restricted problem is still a large-scale MINLP which contains binary variables in the piecewise linear functions describing order values. All dynamic models have to be optimized simultaneously with the restricted scheduling model. To simplify the primal problem of the first algorithm, the second decomposition algorithm is developed which further decomposes the dynamic models from the restricted scheduling model. The dynamic models are optimized independently over every operational task. A nested decomposition structure is derived by the second decomposition algorithm, where all subproblems are simple and can be solved efficiently. The proposed model formulation under uncertainty and the decomposition algorithms are demonstrated by a complicated case study containing more than 3 million equations/variables with 100 scenarios from sampling over the probability distributions of major process uncertain parameters. The direct solution approach by the MINLP solver fails to find a feasible solution in 50 h while the two decomposition algorithms return the optimal solution. The first single decomposition algorithm solves the problem in 23.9 h, and the second nested algorithm only requires 3.3 h. The average production profit returned by the integrated method is 17.6% higher than that of the sequential method. The remainder of the paper is organized as follows. The problem statement is given in section 2. The formulation of the integrated problem under uncertainty is presented in section 3. The two decomposition algorithms based on the GBD framework are developed in section 4. The integrated method

is demonstrated by a case study in section 5. The conclusion is given in section 6. The detailed constraints of the integrated problem are listed in the Appendix, and the parameters of the case study are listed in Tables 1−3. Table 1. Parameters of Dynamic Models for Reaction Tasks Executed in Reactor RI or RII symbol k1

description

value for job 1, 2, 3

value for job 4, 5, 6

unit

1E7

5E5

h−1

1E10

1E7

h−1

5E3

4E3

K

E1R

pre-exponential factor pre-exponential factor activation energy

E2R

activation energy

8E3

6E3

K

ΔH1 ΔH2 VR ρR

heat of reaction heat of reaction volume of reactor density of reactor

−1E4 −1E3 5 1E3

−5E3 −1E3 5 1E3

kJ/kmol kJ/kmol m3 kg/m3

cR

heat capacity of reactor volume of jacket density of jacket heat capacity of jacket heat transfer coefficient temp of heating fluid temp of cooling fluid unit price of heating fluid unit price of cooling fluid initial conc of species A initial conc of species B initial reactor temp initial jacket temp final conc of species B final reactor temp max flow rate of heating fluid max flow rate of heating fluid

2.5

2.5

kJ/kg·K

1.5 1E3 4

1.5 1E3 4

m3 kg/m3 kJ/kg·K

8E4

8E4

kJ/m2·K·h

350 294 5

350 294 5

K K m.u./m3

1

1

m.u./m3

1

1

kmol/m3

0

0

kmol/m3

294 294 0.7

294 294 0.7

K K kmol/m3

300 15

300 15

K m3/h

15

15

m3/h

k2

VJ ρJ cJ UAJ Thot Tcool PRhot PRcool C0A C0B T0R T0J CFB TFR Fmax hot Fmax cool

Table 2. Parameters of Dynamic Models for Reaction Tasks Executed in Reactor RIII or RIV symbol k3 k4 PRin V0R C0B C0D C0E C0F VCFE VFF Fmax in 16853

description rate constant rate constant unit price of inlet materials initial reactor volume initial conc of species B initial conc of species D initial concentration of species E initial conc of species F final amount of species E final reactor volume max inlet flow rate

value for job 1, 2, 3 2 1 10

value for job 4, 5, 6 1.5 0.5 10

unit 3

m /kmol·h m3/kmol·h m.u./m3

5 0.7 0 0

5 0.7 0 0

m3 kmol/m3 kmol/m3 kmol/m3

0 2.5

0 2.5

kmol/m3 kmol

10 5

10 5

m3 m3/h



Article

general differential and algebraic equations with uncertain parameters so that other types of uncertainties can also be solved by the proposed method. Expressions of uncertainty in an optimization problem also vary. In stochastic programming, uncertain parameters are regarded as random variables. To facilitate the solution algorithm, the uncertain parameters are often assumed to have discrete probability distributions such that they have a finite set of possible values, called scenarios.49 The integrated problem is solved according to the given scenarios representing the uncertainty realization. The objective is to maximize the production profit on average over all scenarios. Generating scenarios is important to a two-stage stochastic program. As the number of uncertain parameters grows, the number of scenarios by enumerating all possible values of the uncertain parameters often increases exponentially. It is difficult to consider all scenarios in a stochastic program. However, various scenario generation and reduction approaches have been developed, 50−52 which can return a small number of representative scenarios. These approaches can be applied to generate scenarios for the integrated problem. We assume the generated scenarios are given as problem parameters. Specifically, the integrated scheduling and dynamic optimization problem under uncertainty is stated as follows:

Table 3. Scheduling Parameters symbol

description

value

unit

PTijkfix

fixed processing time for filtration tasks

0.7

h

PTijkfix

fixed processing time for separation tasks

1.0

h

STi

P jR

unit setup time constant order value

Cjmax

max job completion time

Dj

job due date

CF

fixed cost

0.1 200 10 5 300

h m.u. h h m.u.

2. PROBLEM STATEMENT The model formulation and the solution algorithm for the problems of integrated scheduling and dynamic optimization largely depend on the process structure. This work focuses on the common sequential batch process where operations on material splitting and mixing are not allowed.41,42 In a sequential batch process, the batch integrity is preserved. All materials processed in a unit will be delivered, as a whole, to the subsequent unit. Therefore, the order demand can be quantified by the number of batches. The entire procedure for producing one batch of products starting from raw materials is referred to as a job. An operational stage of a job executed in a unit is called a task. The task execution in a unit can then be described by a dynamic model. The dynamic model determines the task recipe of the processing time and the processing cost, which can be optimized simultaneously with scheduling decisions by solving the integrated problem. The variable recipes enable the integrated method to optimize the entire production process collaboratively. The achieved overall performance is superior to the traditional sequential method where task recipes are fixed parameters for the scheduling problem. The variable task recipes, however, substantially complicate the integrated problem. A batch process is subject to various uncertainties. This work studies process uncertainty occurring in dynamic systems, e.g. uncertain kinetic parameters and initial concentrations for reaction tasks. We formulate process uncertainty by a set of

3. MODEL FORMULATION Basically, the integrated model is formulated by combining the scheduling model and all dynamic models. An illustration of the integrated problem is shown in Figure 2. Jobs are indexed by j or j′, and each job consists of a set of operational stages indexed by k or k′. The combination of a job and a stage denotes a task indicated by (j,k) or (j′,k′). A task is executed in a capable unit indexed by i. The main scheduling decisions are to assign a task to a unit and to determine the task execution sequence when multiple tasks are assigned to the same unit. Binary variables ξijk are introduced to determine the task assignment. When ξijk = 1, task (j,k) is assigned to unit i. The task execution sequence in a unit is represented by binary variables βijkj′k′. When βijkj′k′ = 1, task (j,k) precedes task (j′,k′) in unit i.

16854



Article

Figure 2. Formulation of the integrated scheduling and dynamic optimization problem.

uncertain. It is worthwhile to note that Figure 2 only illustrates one dynamic model. There are actually a great number of dynamic models in the integrated problem, which correspond to every combination of tasks and their capable units. The objective of the integrated problem is to maximize the production profit. It is an increasing function of the order values by completing the jobs and a decreasing function of the processing costs for executing the tasks. The order value, exhibited in Figure 3, is a nonincreasing function with job completion times, dtj, since a penalty term is added to tardy jobs.14,46 When job j is completed before the due date Dj, the fulfilled order has a constant value PRj . However, when the completion time exceeds the due date, the order value decreases, penalizing a tardy job. The upper bound of the completion time is Cmax j . To express the piecewise linear function of the order

Besides the binary assignment and sequencing variables, the scheduling problem needs to determine the starting times tsjk and the ending times tejk of each task. The duration between a starting time and an ending time is the task processing time, denoted by ptijk, which is linked to the dynamic model describing the task execution. The dynamic model is also linked to the scheduling model via the processing cost, denoted by pcijk. The processing time and the processing cost represent the operational recipe of the task, which can be adjusted in the integrated problem. To determine the flexible task recipe, a set of differential and algebraic equations characterizing the dynamic model should be included into the integrated problem as constraints. The task processing time and the processing cost are adjusted by manipulating the process input trajectory, denoted by uijk(tijk), which governs the state trajectory, denoted by xijk(tijk). In the system equations, PDyn ijk denotes model parameters, which can be 16855



Article

dynamic models should be optimized simultaneously with the scheduling model. Process uncertainty is described by a set of scenarios with varying parameter values in the dynamic models. The structure of the integrated problem under uncertainty is displayed in Figure 4. The binary assignment and sequencing variables are the firststage variables. The remaining variables are the second-stage variables, including the flexible recipes expressed by the task processing times and processing costs, other scheduling variables such as task starting times and ending times, and all variables in the dynamic models. The first-stage variables are determined to maximize the average production profit over the scenarios while the second-stage variables are adjusted according to the realization of a specific scenario. According to the integration structure, the variables are clustered into different categories. The binary assignment and sequencing variables are denoted by

Figure 3. Order value function.

ζ = {Binary assignment and sequencing variables}

value, a binary variable γj is introduced to denote if the completion time exceeds the due date. There is a time−cost trade-off in the profit function. The order value for a job is favored by a small job completion time and, in turn, by small task processing times related with the job. A small processing time may, however, cause a large processing cost, which can lower the profit. To make a better trade-off, all

As the first-stage variables, their value does not depend on the realization of uncertainty. The remaining second-stage variables rely on a specific realization of uncertainty, and thus, they are indexed by scenario index s. The flexible recipes are represented by the task processing times and processing costs as

Figure 4. Structure of integrated model under uncertainty. Red horizontal line 1 divides the integrated model according to the classification of the firststage and second-stage variables. This line also denotes the cutting point from which both GBD algorithms decompose the integrated problem. Red horizontal line 2 divides the dynamic models from the scheduling model, which is the cutting point from which the nested GBD algorithm decomposes the integrated problem. 16856



Article

complicating variables, are temporarily fixed, the restricted problem is easy to solve.47,48 The integrated problem under uncertainty shares such special structures. When the first-stage variables are fixed, the restricted problem becomes a set of separable ones over each scenario. Furthermore, after a deeper investigation of the problem structure shown in Figure 4, it is seen that the dynamic models for each task can be optimized separately when the recipe variables linking the dynamic models to the scheduling model are fixed. Specifically, the first algorithm presented in section 4.1 decomposes the integrated problem according to the uncertainty scenarios, and the second algorithm presented in section 4.2 further decomposes the integrated problem so that each dynamic model for a task can be optimized independently. 4.1. Single GBD Algorithm. A GBD algorithm generally consists of several common steps: (1) Fixing the complicating variables to formulate the restricted problem into the primal problem. (2) Solving the primal problem to obtain a lower bound of the original (maximization) problem and to return the dual information for updating the complicating variables. (3) Using the returned dual information to formulate the master problem, which is a relaxation of the original problem. The master problem is solved to provide an upper bound of the original (maximization) problem and to determine a new value at which the complicating variables are fixed. (4) Iterating between the primal problem and the master problem until the gap between the lower bound and the upper bound is smaller than a given tolerance level. The core of a GBD algorithm is to take advantage of the special model structure to simplify the optimization procedure. To explicitly express the special structure which will be exploited to develop the decomposition algorithm, the integrated problem under uncertainty (Integrated_Problem) is reformulated into a bilevel program as follows:

{ptijks} = {Task processing times} {pcijks} = {Task processing costs}

The remaining variables in the scheduling model are cast into the set of ωs = {Remaining variables in the scheduling model}

The state variables and the input variables in the dynamic models are grouped into {xijks} = {State variables of dynamic models} {uijks} = {Input variables of dynamic models}

Using the variable classification, the integrated scheduling and dynamic optimization problem under uncertainty is summarized as (detailed formulation can be found in Appendix A1) (Integrated_Problem) max Cζ +

∑ (Hsωs + ∑ Vijksptijks − ∑ Wijkspcijks) s

i ,j,k

i ,j,k

(1)

s.t. Aζ ≤ B Gsζ + Ssωs +

(2)

∑ Tijksptijks ≤ Es , s (3)

i ,j,k Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0,

(j , k) ∈ JK ,

i ∈ Ijk , s (4)

Dyn pcijks = ηijks(xijks , uijks , ptijks , Pijks ),

i ∈ Ijk , s

(j , k) ∈ JK , (5)

(Bi-level_L1)

The first term in the objective function in eq 1 denotes the contribution from the first-stage variables, where C is a constant coefficient vector. The second term denotes the average contribution from the second-stage variables over the scenarios where Hs, Vijks, and Wijks are weight vectors before corresponding variables. Inequality 2 represents the constraints on the first-stage binary variables, where A is a parameter matrix and B is a parameter vector. Inequality 3 represents the linking constraints between the first-stage variables and the second-stage variables, in which Gs, Ss, Tijks, and Es are all parameters. Inequality 4 encapsulates the dynamic models including discretized system equations and path constraints. The parameters of the dynamic models, including initial conditions and set points, are represented by PDyn ijks . Equation 5 defines the processing cost pcijks as a function of the state xijks, the input uijks, the processing time ptijks, and the model parameter PDyn ijks .

max Cζ +

∑ γs(ζ ) s

s.t. Aζ ≤ B

(Bi-level_L2) γs(ζ ) = max Hsωs +

∑

Vijksptijks −

i ,j,k

∑

Wijkspcijks , s

i ,j,k

s.t. Gsζ + Ssωs +

∑

Tijksptijks ≤ Es , s

i,j,k Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0, (j , k) ∈ JK , i ∈ Ijk , s

4. GBD ALGORITHMS The integrated problem under uncertainty (Integrated_Problem) formulated in the previous section is an extremely largescale MINLP, which is intractable for a general-purpose MINLP solver. Actually, even solving a deterministic integrated problem without uncertainty is time-consuming.12−14 To solve the complicated problem efficiently, we develop two decomposition algorithms in this section. Both decomposition algorithms, shown in Figure 1, are derived based on the GBD framework. GBD is a framework for solving complicated problems with special structures; that is, when a subset of variables, called

Dyn pcijks = ηijks(xijks , uijks , ptijks , Pijks ), (j , k) ∈ JK , i ∈ Ijk , s

In the top level problem (Bi-level_L1), the second-stage variables are encapsulated by the optimal-value functions (or the recourse functions) γs(ζ), which are evaluated by solving the bottom level problems (Bi-level_L2). When the first-stage variables contained in ζ are fixed in the top level problem, the bottom level problems can be solved independently over each scenario s. 16857



Article

Figure 5. Flow diagram of the single GBD algorithm.

Consequently, the primal problem is formulated as the bottom level problem with a consensus constraint

The consensus constraint (eq 6) represents the linking equation between the master problem and the primal problem, where ζ̅(m) is the fixed value determined by the master problem at iteration m. In the primal problem, the variables in ζ are relaxed to be the continuous ones shown in constraint 7. Actually, the consensus constraint (eq 6) guarantees that variables in ζ can only have binary values. The master problem is formulated as a relaxed top level problem with Benders cuts and integer cuts

(Bi-level_Primal) max γs = Hsωs +

∑ i ,j,k

Vijksptijks −

∑

Wijkspcijks , s

i ,j,k

s.t. Ssωs +

∑

Tijksptijks ≤ Es − Gsζ , s

(Bi-level Master)

i ,j,k

max ψ = Cζ +

Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0, (j , k) ∈ JK , i ∈ Ijk , s

s


ζ = ζ̅

(m)

ζ ∈ [0, 1]Nζ

∑ δs

s.t. Aζ ≤ B (6)

δs ≤ ρs(m ′)ζ + νs(m ′) , s , m′ ≤ m

(7) 16858

(8)



Article

Then the primal problem (Bi-level_Primal) is solved to return the optimal solution, the optimal function value, and the dual value of the consensus constraint. Based on the returned values, the objective function of the integrated problem is evaluated and compared to the lower bound. If the objective function value is greater than the lower bound, the lower bound is updated and the best solution is recorded. The algorithm terminates if the gap between the lower bound and the upper bound is less than the given threshold. Otherwise, the parameters of the master problem are updated. The master problem (Bi-level_Master) is solved to find a new value for the first-stage variables and an upper bound. By fixing the first-stage variables at the newly determined value, the next iteration begins. This single GBD algorithm is consistent with the L-shaped method and its variations49,53−55 that are widely used for solving two-stage stochastic programming problems. Due to the nonconvexity of the MIDO problem, this algorithm cannot guarantee global optimality of the solution, but it can solve the problem to a good feasible solution, as will be illustrated in the case study. 4.2. Nested GBD Algorithm. The GBD algorithm in the previous subsection decomposes the integrated problem according to the scenarios. However, the resulting primal problem is still a large-scale MINLP, where all dynamic models have to be optimized simultaneously. Thus, we further decompose the MINLP problem so that each dynamic optimization problem can be solved independently. For this purpose, a nested GBD algorithm is developed in this subsection. To reveal the special problem structure that will be exploited in the nested decomposition algorithm, the integrated problem (Integrated_Problem) is reformulated into a trilevel program as

∑ ζ − ∑ ζ ≤ |I1(m′)| − 1, m′ ≤ m I1(m ′)

I0(m ′)

(9)

The optimal-value functions γs(ζ) in the bilevel program are estimated by the Benders cut (eq 8), where m denotes the current iteration and m′ indicates all previous iterations. The ) coefficient ρ(m s ′ is the dual value of the consensus constraint (eq ) 6), and the constant ν(m s ′ is νs(m ′) = γs*(m ′) − ρs(m ′)ζ ̅

(m ′)

, s , m′ ≤ m

(10)

where γs* ′ is the optimal value of the primal problem at iteration m′ and ζ̅(m′) is the fixed value when the primal problem is solved. The integer cut (eq 9) excludes all binary values ) searched in the previous iterations. The set I(m 1 ′ contains the position indices at which the element variables in ζ̅(m′) equal 1, ) (m ) while the set I(m 0 ′ includes the indices of those in ζ̅ ′ having ) ) values equal to zero. |I(m ′ | calculates the size of the set I(m 1 1 ′. In the master problem (Bi-level_Master), we do not include the infeasibility cuts in a common GBD algorithm. The infeasibility of the primal problem can only occur when a job completes after the limit Cmax shown in the order function value j in Figure 3. The infeasibility can be easily handled by increasing Cmax to a sufficiently large value. According to the piecewise linear j function, if the job completion time is large, a low negative order value can be returned. The information of the poor objective function value will enter the master problem through the optimality cut (eq 8). This means the optimality cut can automatically exclude the values of the first-stage variables in ζ, resulting in a large job completion time. Therefore, there is no need to explicitly include infeasibility cuts in the master problem. Figure 5 is the flow diagram of the single GBD algorithm. To initialize the GBD algorithm, the integrated problem under uncertainty is first solved in a sequential way. The tasks recipes are determined by solving a set of dynamic optimization problems minimizing the processing times (Minimize_PT) (m )

(Tri-level_L1) max ψ = Cζ +

s

s.t. Aζ ≤ B

min ptijks , , (j , k) ∈ JK , i ∈ Ijk , s

(Tri-level_L2)

s.t.

γs(ζ ) = max Hsωs +

Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0, (j , k) ∈ JK , i ∈ Ijk , s

Gsζ + Ssωs +

s

i ,j,k

∑

Wijksφijks(ptijks) , s

i,j,k

∑

Tijksptijks ≤ Es , s

(Tri-level_L3) φijks(ptijks) = min pcijks ,

*) Wijkspcijks

(j , k) ∈ JK ,

i ∈ Ijk , s

i ,j,k

s.t.

s.t.

Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0,

Aζ ≤ B

Gsζ + Ssωs +

∑

i,j,k

(Scheduling_FR) * − Vijksptijks

Vijksptijks −

s.t.

The optimal values of pt*ijks and pc*ijks returned by the dynamic optimization problems (Minimize_PT) are used as fixed recipes to solve the scheduling problem under uncertainty

∑ (Hsωs + ∑

∑ i,j,k


max Cζ +

∑ γs(ζ )

∑


* ≤ Es , s Tijksptijks

i ,j,k

(j , k) ∈ JK ,

i ∈ Ijk , s

(j , k) ∈ JK ,

i ∈ Ijk , s

In the scheduling problem (Scheduling_FR), the processing times and the processing costs are all fixed parameters. The MILP problem is solved to determine the starting value for the first-stage variables ζ0. The upper and lower bounds are initialized at positive infinity and negative infinity, respectively.

The top level problem (Tri-level_L1) is the same as the one (Bilevel_L1) in the bilevel formulation. However, the bottom level problem (Bi-level_L2) in the bilevel formulation is expressed by another bilevel program, resulting in a trilevel program overall. In 16859



Article

Figure 6. Flow diagram of the nested GBD algorithm.

the middle level problem (Tri-level_L2), the dynamic models are encapsulated by the optimal-value function ϕijks(ptijks), which is evaluated by solving the bottom level problem (Tri-level_L3). The trilevel program is solved by a nested GBD algorithm where a GBD algorithm solving the bottom two subproblems is embedded inside another GBD algorithm solving the top two

subproblems. The former embedded algorithm is called the inner GBD algorithm, while the later one is called the outer GBD algorithm. We start from the inner GBD algorithm. The primal problem is derived from the bottom level problem (Tri-level_L3) by adding the consensus constraint as 16860



Article

Figure 7. Diagram of the case study.

(Tri-level_Inner_Primal) min pcijks ,

(j , k)∈ JK ,


i∈ Ijk , s

(11)

(j , k) ∈ JK ,

i ∈ Ijk , s

s.t. Dyn δijks(xijks , uijks , ptijks , Pijks ) ≤ 0,

(j , k) ∈ JK ,

(n) ptijks = ptijks ,

i ∈ Ijk , s 16861

(j , k) ∈ JK , s



Article

Figure 8. Iteration results of the single GBD algorithm: (a) upper and lower bounds and (b) computational times.

Figure 9. Iteration results of the nested GBD algorithm: (a) upper and lower bounds and (b) computational times. The master problem of the inner GBD loop is also the primal problem of the outer GBD loop.

In the consensus constraint (eq 11), the processing time ptijks is fixed at pt (n) ijks determined by the master problem where the iteration is indexed by n. The master problem of the inner GBD algorithm is derived from the middle level problem (Trilevel_L2) as

To ensure that the primal problem (Tri-level_Inner_Primal) has a feasible solution, the task processing times cannot be too small, because variables in a real-world dynamic system cannot change as fast as desired, owing to practical constraints on the system, e.g. the bounds imposed on the input variables. The lower bounds of the processing times, denoted by PTmin ijks , are determined by solving the dynamic optimization problems (Minimize_PT). By adding the feasibility constraint (eq 13) directly, the Benders feasibility cut constraints can be omitted. The master problem of the inner GBD algorithm is linked to the first-stage variables in the outer GBD algorithm, which is formulated by the consensus constraint (eq 14). The first-stage variables are fixed at ζ̅(m), where m indexes the iteration in the outer GBD algorithm. As shown in the constraint in eq 15, the values of ζ can be relaxed to continuous ones. The consensus constraint (eq 14) guarantees the binary values for ζ. The master problem (Tri-level_Inner_Master) in the inner GBD algorithm is also the primal problem in the outer GBD algorithm

(Tri-level_Inner_Master) max γs = Hsωs +

∑

Vijksptijks −

i ,j,k

∑

Wijksϕijks , s

i ,j,k

s.t. Ssωs +

∑

Tijksptijks ≤ Es − Gsζ , s

i ,j,k (n ) (n ) ′, ′ ptijks + μijks ϕijks ≥ πijks

(j , k) ∈ JK ,

s , n′ ≤ n min ptijks ≥ PTijks ,

ζ = ζ̅

i ∈ Ijk , (12)

(j , k) ∈ JK ,

(m)

ζ ∈ [0, 1]Nζ

i ∈ Ijk , s

(13)

(Tri‐level_Outer_Primal) = (Tri‐level_Inner_Master)

(14)

The master problem of the outer GBD algorithm is identical to that for the bilevel formulation

(15)

(Tri‐level_Outer_Master) = (Bi‐level_Master)

In the Benders cut (eq 12), n denotes the current iteration and n′ indicates a previous one. The Benders cut is derived from the ) dual value π(n ijks′ of the consensus constraint (eq 11) when the ) primal problem is solved at iteration n′. The constant μ(n ijks′ is (n ′) (n ′) (n ′) *(n ′) − πijks μijks = pcijks ptijks

) except that the dual value ρ(m s ′ is obtained from the consensus constraint (eq 14) for the last iteration in the inner GBD ) algorithm. The constant ν(m s ′ is calculated in the same way according to eq 10. The flow diagram of the nested GBD algorithm is shown in Figure 6. The outer iteration is initialized in the same way as that in the single GBD algorithm. The processing times in the inner iteration are initialized at the minimum values. The primal

(16)

where pcijks *(n′) is the optimal value of the primal problem and (n ) pt ijks′ is the fixed processing time. 16862



Article

problems (Tri-level_Inner_Primal) in the inner iteration are solved to obtain the optimal solution, the optimal function value, and the dual value. These values are used to evaluate the lower bound of the inner iteration as well as the lower bound of the outer iteration. If the gap of the inner iteration is not smaller than the given threshold, the master problem (Tri-level_Inner_Master) is solved to determine a new upper bound for the inner iteration and new values for fixing the processing times at the next iteration. If the gap of the inner iteration is smaller than the threshold, a similar convergence test is applied to the outer iteration. If the outer iteration does not terminate, the master problem (Tri-level_Outer_Master) is solved based on the dual information from the primal problem (Tri-level_Outer_Primal) which is actually the master problem (Tri-level_Inner_Master) in the inner loop. When the outer iteration terminates, the optimal solution and the optimal function value are returned.

5. CASE STUDY To demonstrate the integrated method, we apply it to a complicated example, shown in Figure 7. There are six jobs

Figure 11. Results returned by (a) the sequential method and (b) the integrated method. The task indexes (j,k) are shown on each task bar. The task assignment and sequencing decisions are made by the firststage variables. The task starting times and ending times are the secondstage variables, dependent on the scenarios. The figures show task starting times and ending times under the first scenario.

Figure 10. Comparison of the objective function values returned by the sequential method and the integrated method.

species F is the byproduct. The processing cost is the amount of the byproduct E produced after the reaction task is executed. As each job has two reaction tasks and each reaction task has two capable parallel units, there are 6 (#jobs) × 2 (#reaction tasks) × 2 (capable units) = 24 dynamic models in the integrated problem. Every dynamic model is discretized by the collocation method with 30 finite elements. The total length of the finite elements is a variable so that the task processing time is adjustable. After the discretization procedure, the integrated problem is reformulated into a large-scale MINLP. Process uncertainty is expressed by the uncertain rate coefficients k1, k2, k3, and k4, and initial reactant concentration C0A in the first reaction task and C0B in the second reaction task. The rate coefficients can have a value of 80%, 100%, or 120% of their nominal value with the equal probability of 1/3. The initial concentrations can have a value of 90%, 100%, or 110% of their nominal value with the equal probability of 1/3. To formulate the two-stage stochastic program, 100 scenarios are generated. The final concentrations of a reaction task specify a quality measure that the dynamic system should meet. Thus, they are the given parameters for optimizing the dynamic system. The set points can be reached by manipulating dynamic trajectories and the processing time according to the uncertainty realization. These variables characterize the operating condition and belong to the second-stage variables. The batch production has two reaction tasks. Ideally, the final concentrations of the first reaction task are equal to the initial concentrations of the second

processed in the batch plant. Each job consists of four operational stages: reaction, filtration, reaction, and separation. Each stage can be processed in one of the two parallel units. The filtration task and the separation task have fixed operational recipes where the processing times and the processing costs are parameters. The reaction tasks are described by dynamic models and their operational recipes of the processing times, and the processing costs can be adjusted by manipulating the input trajectories of the dynamic models. The differential and algebraic equations of the dynamic models are listed in Figure 7, and the model parameters are given in Appendix A2. The first reaction task occurring in reactor RI or RII includes nonisothermal consecutive reactions. Each batch reactor is surrounded by a heating/cooling jacket. The raw material is species A, and the product is species B. Species C is the byproduct. The process inputs are the flow rates of heating fluid and cooling fluid. Manipulating the flow rates changes the jacket temperature and then changes the reactor temperature via the heat exchanger. The reaction rates are finally changed according to the Arrhenius equations. The processing cost is associated with the heating fluid and cooling fluid consumed to process the task. The second reaction task includes competitive isothermal parallel reactions taking place in a semibatch reactor RIII or RIV where reactant D is added gradually. The process input is the inlet flow rate of reactant D. Species E is the product while 16863



Article

Figure 12. Dynamic trajectories of task (6,1) returned by the sequential method (left panel) and the integrated method (right panel). All variables in the dynamic model belong to the second-stage variables dependent on the scenarios. The figures show the trajectories under the first scenario.

reaction task. However, the two reaction tasks are separated by a filtration task, which can change the concentrations. Therefore, the initial concentrations of the second reaction task are regarded as uncertain parameters. The resulting two-stage stochastic nonlinear program for the integrated problem under uncertainty is a large-scale MINLP, which has 3,193,245 equations, 3,187,693 continuous variables, and 204 binary variables. To solve an integrated problem, the MINLP solver SBB is often applied.6,9,10,12,13 However, SBB failed to find a feasible solution within 50 h for the integrated problem in this case study. By contrast, the integrated problem can be solved by the decomposition algorithms. The single GBD algorithm presented in section 4.1 spends 85959.7 s (about 23.2 h) in solving the integrated problem under uncertainty. The iteration procedure is shown in Figure 8. The terminating gap between the upper bound and the lower bound is set as 1%, and the algorithm stops after 75 iterations. The maximum average profit returned is 450.9 m.u. The algorithm decomposes the integrated problem according to the scenarios. The primal problem is a set of separable MINLP problems for every scenario. Each MINLP has 36,607 equations, 36631 continuous variables, and 6 binary variables. The master problem is an MILP containing constraints on the binary assignment and

sequencing variables as well as the Benders cuts and the integer cuts. There are 1,142 equations, 4,901 continuous variables, and 204 binaries in the master problem at the first iteration. Owing to the complexity of the primal problem, most computational time (more than 99%) is spent in solving the primal problem. The nested GBD algorithm presented in section 4.2 spends 11,737.7 s (about 3.3 h) in solving the integrated problem under uncertainty. The iteration procedures are shown in Figure 9. The gaps of both inner iteration and outer iteration are set as 1%. The algorithm stops after 50 iterations in the outer loop. The nested GBD algorithm further decomposes the MINLP primal problem of the single GBD algorithm into an MILP master problem (363 equations, 117 continuous variables, and 6 binaries at the first iteration) and a primal problem containing a set of separable NLP problems (each having 1,362 equations and 1,378 continuous variables). The dynamic models for each task can be optimized independently in the inner GBD loop. The master problem of the inner loop is also the primal problem of the outer loop. The master problem of the outer loop is identical to that of the single GBD algorithm except that the dual information is returned from different primal problems. Similar to the single GBD algorithm, the primal problem involving dynamic optimization in the nested GBD algorithm 16864



Article

Figure 13. Dynamic trajectories of task (6,3) returned by the sequential method (left panel) and the integrated method (right panel). All variables in the dynamic model belong to the second-stage variables dependent on the scenarios. The figures show the trajectories under the first scenario.

accounts for the most computation resources. However, the computational time for the primal problem per iteration is reduced from 1134.9 to 229.2 s, showing the advantage of the nested GBD algorithm. The maximum average profit returned by the nested GBD algorithm is the same as that of the single GBD algorithm. For comparison, we also solve the integrated problem under uncertainty by the traditional sequential method, which solves the dynamic optimization problems separately from the scheduling problem. The dynamic optimization problems determine the task recipes of the processing times and the processing costs, which are fixed parameters for solving the scheduling problem. Because the dynamic models are optimized before the scheduling problem is solved, the objective function value of the profit is not available for the dynamic optimization. Thus, a local criterion has to be chosen. A common criterion is to minimize the task processing times, specifically solving the problems (Minimize_PT). Small task processing times would reduce the job completion times and might increase the production profit. After the task recipes are determined, the scheduling problem (Scheduling_FR) is solved. The sequential method is actually used to initialize the two GBD algorithms. Each dynamic model of a task under a scenario can be optimized

independently in the sequential method, resulting in a fast solution procedure. However, the task recipes determined by optimizing the local criterion may not be the optimal ones for the global criterion of the production profit. A small processing time may imply a large processing cost. The best trade-off can only be made by the integrated method. The average profit returned by the sequential method is only 371.6 m.u., 17.6% less than that returned by the integrated method. The components of the objective function computed by the two methods are displayed in Figure 10. By minimizing the task processing times, the sequential method results in slightly higher sales than the integrated method. However, the minimum processing times accompany large processing costs. The sequential method causes a much larger total processing cost than the integrated method. Because the integrated method is able to make a better trade-off between the sales and the total processing cost, it returns a higher profit than the sequential counterpart. The Gantt charts returned by the integrated method and the sequential method are shown in Figure 11, respectively. As the integrated method can adjust processing times, the processing units can be utilized in a more effective way with fewer idle periods than those scheduled by the sequential method. As an 16865



Article

The constraints in the scheduling models are listed as follows.

example of dynamic optimization results, the dynamic trajectories of task (6,1) and task (6,3) are shown in Figures 12 and 13, respectively. The trade-off between the processing times and the processing costs can also be revealed from the dynamic trajectories. For task (6,1), the flow rates of heating fluid and cooling fluid, determined by the sequential method, nearly reach the maximum values to reduce the processing time, meanwhile leading to a large consumption of utilities. By contrast, a smaller amount of utilities are consumed by the integrated method by allowing a larger processing time. The same conclusion can be drawn from the dynamic trajectories for task (6,3). From the dynamic trajectories, we can see that the sequential method and the integrated method return distinct process inputs. In the first reaction task, the process inputs are the heating and cooling flow rates, which change the jacket and the reactor temperature directly. Thus, the temperature trajectories generated by the two methods are also distinct. However, the concentrations are changed indirectly by the process inputs via the temperatures so the concentrations cannot be changed abruptly. Though the tendencies of the smooth concentration trajectories of the two methods look similar, the processing times are quite different. Similarly, the two methods result in distinct process inputs for the second reaction task but similar tendencies of the concentration trajectories.

Each task is executed exactly once

∑ ξijk = 1,

(A1)

Relation between the assignment variables and the sequencing variables ξijk + ξij

− βijkj k − βij k jk ≤ 1, (j , k) ∈ JK , ′k′ ′′ ′′ (j′, k′) ∈ JK , (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij

′k′

(A2)

≤ ξijk , (j , k) ∈ JK , (j′, k′) ∈ JK , ′k′ (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij k ′′

(A3)

≤ ξij k , (j , k) ∈ JK , (j′, k′) ∈ JK , ′′ ′k′ (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij k ′′

(A4)

βijkj

βijkj

Cuts on infeasible assignment and sequencing ξijk = 0, (j , k) ∉ JK , ori ∉ Ijk

6. CONCLUSION Integration of scheduling and dynamic optimization for batch processes under uncertainty is an important, while often overlooked, issue. In this work, a new formulation of the integrated problem under uncertainty is presented. The integrated problem is formulated into a two-stage stochastic nonlinear program which is then transformed into an MINLP by discretizing the dynamic models. To solve the complicated integrated problem, two decomposition algorithms are developed based on the GBD framework. The single GBD algorithm decomposes the integrated problem so that the subproblems can be solved independently over all the scenarios. The nested GBD algorithm further decomposes the integrated problem so that the dynamic models of each task can be optimized separately. The proposed model formulation and the decomposition algorithms are demonstrated by a complicated case study including more than 3 million variables/equations with 100 scenarios. The two algorithms return identical results. The computational time of the single GBD algorithm is 23.2 h while that of the nested algorithm is only 3.3 h. Compared to the sequential method, the integrated method increases the average profit by 17.6%.

■

(j , k) ∈ JK

i ∈ Ijk

(A5)

= 0, (j , k) ∉ JK , or ′k′ (j′, k′) ∉ JK , or(j , k) = (j′, k′),

βijkj

or

i ∉ Ijk ∩ Ij

′k′

(A6)

Relation between the starting and ending times of the tasks with f ixed processing times PTfix ijk tejks = tsjks +

∑

ξijkPTijkfix ,

(j , k) ∈ JK F , s

i ∈ Ijk

(A7)

Relation between the starting and ending times of the tasks with variable processing times ptijks (xptijks = ξijkptijks) tejks = tsjks +

∑

xptijks ,

(j , k) ∈ JK D , s

i ∈ Ijk

(A8)

0 ≤ xptijks ≤ ptijks ,

(j , k) ∈ JK D ,

i ∈ Ijk , s

(A9)

max xptijks ≤ ξijkPTijks ,

(j , k) ∈ JK D ,

i ∈ Ijk , s

(A10)

max xptijks ≥ ptijks − (1 − ξijk)PTijks ,

APPENDIX

(j , k) ∈ JK D ,

i ∈ Ijk , s

A1. Formulation of Integrated Scheduling and Dynamic Optimization Problem

(A11)

Execution sequence of tasks belonging to the same job

Scheduling Model under Uncertainty. The general precedence model43−45 is adopted to formulate the scheduling problem. The main decision variables are illustrated in Figure 2. In an integrated problem, there is no need to express all task execution procedures by dynamic models. Some tasks, such as packing tasks, typically have fixed operational recipes. Thus, the entire task set JK is partitioned into two mutually exclusive subsets as JK = JKF ∪ JKD, where JKF includes the tasks with fixed recipes and JKD includes the tasks having flexible recipes determined by dynamic models.

tsjks ≥ tej(k − 1)s ,

(j , k) ∈ JK ,

k ≥ 2, s

(A12)

Task sequence in a unit tsjks + PTijkfix + STi − BM(1 − βijkj k ) ≤ tsj k s , ′′ ′′ (j , k) ∈ JK F , (j′, k′) ∈ JK F , (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij 16866

′k′

,s

(A13)



Article

Initial conditions

tsjks + ptijks + STi − BM(1 − βijkj k ) ≤ tsj k s , ′′ ′′ (j , k) ∈ JK D , (j′, k′) ∈ JK D , (j , k) ≠ (j′, k′),

Dyn I (r = 0) xijks = Xijks ∈ Pijks ,

,s (A14) ′k′ Dynamic Models under Uncertainty. The dynamic model for executing a task can be expressed by a set of ordinary differential and algebraic equations

dtijks

s

Final conditions R (r = Nijk )

xijks

Dyn F = Xijks ∈ Pijks ,

(j , k) ∈ JK D ,

i ∈ Ijk ,

s (A19)

Processing times

Dyn ), = fijks (xijks(tijks), uijks(tijks), Pijks

(j , k) ∈ JK D ,

i ∈ Ijk ,

(A18)

i ∈ Ijk ∩ Ij

dxijks(tijks)

(j , k) ∈ JK D ,

R ptijks = lijksNijk ,

∀ (j , k) ∈ JK D ,

i ∈ Ijk ,

s

(A20)

i ∈ Ijk , s

Processing costs 0 ≥ hijks(xijks(tijks), uijks(tijks),

Dyn Pijks ),

D

(j , k) ∈ JK ,

R (r = Nijk )

pcijks = φijks(xijks

∀ (j , k) ∈ JK D ,

),

i ∈ Ijk ,

i ∈ Ijk , s

s (A21)

Objective Function. The objective of the integrated problem is to maximize the average production profit defined by

Each dynamic model is indexed by unit i, task (j,k), and scenario s. Correspondingly, all variables in the differential equations are indexed by ijks, including the time variable tijks. For a compact expression, the states xijks(tijks) and the inputs uijks(tijks) are expressed by the vector forms. The uncertain parameters are packed into PDyn ijks . The vector forms are also used to express other variables and constraints in the dynamic models. To facilitate numerical solution of the dynamic optimization problem, the differential equation is discretized by the collocation method.21,22

profit =

s

where As is the probability weight of scenario s and C denotes the fixed cost term. The sales term is expressed by the sum of the order values prjs, saless =

∑ prjs ,

s (A23)

j 14,46

The order value prjs is a piecewise linear function (shown in Figure 3) of the job completion time dtjs, which is expressed by

rq rq rq Dyn r−1 ′, uijks ′, Pijks xijks = xijks + lijks ∑ Aqq fijks (xijks ), ′ q′

(j , k) ∈ JK ,

(A22)

s F

Discretized dif ferential equations

D

∑ As saless − ∑ As costVs − C F

i ∈ Ijk ,

r > 0, s , r , q

(A15)

rq rq Dyn r r−1 xijks = xijks + lijks ∑ Bq fijks (xijks , uijks , Pijks ),

⎛ ⎞ 1 prjs = P jR ⎜⎜1 − max dt jsII ⎟⎟ , j , s Cj − Dj ⎝ ⎠

(A24)

dt js = dt jsI + dt jsII , j , s

(A25)

γjsDj ≤ dt jsI ≤ Dj , j , s

(A26)

0 ≤ dt jsII ≤ γjs(Cjmax − Dj), j , s

(A27)

dt js = tej(k =|K j|)s , j , s

(A28)

q

(j , k) ∈ JK D ,

i ∈ Ijk ,

r > 0, s , r

(A16)

The coefficients Aqq′, Bq are given by a specific collocation method. In this work, we use the Radau IIA method, which has a robust performance for stiff differential equations.56 The corresponding coefficients are ⎡ 88 − 7 6 ⎢ 360 ⎢ ⎢ 296 + 169 6 [Aqq ] = ⎢ ′ ⎢ 1800 ⎢ ⎢ 16 − 6 ⎢⎣ 36

⎡ 16 − 6 [Bq ] = ⎢ ⎣ 36

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

16 + 6 36

dt js ≤

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

s, r , q

rq uijks ,

Dyn Pijks )

≤ 0,

j, s

(A29)

costVs ,

The variable cost term, is the sum of processing costs, pcijks, over tasks described by dynamic models costVs =

∑

ξijkpcijks , s

(j , k) ∈ JK D , i ∈ Ijk

The bilinear terms can be linearized by the following constraints costVs =

1⎤ ⎥ 9⎦

∑

xpcijks , s

(j , k) ∈ JK D , i ∈ Ijk

0 ≤ xpcijks ≤ pcijks ,

Discretized algebraic equations rq hijks(xijks ,

Cjmax ,

(j , k) ∈ JK D ,

max xpcijks ≤ ξijksPCijk (j , k) ∈ JK D ,

D

(j , k) ∈ JK ,

i ∈ Ijk ,

i ∈ Ijk , s

i ∈ Ijk , s

max xpcijks ≥ pcijks − (1 − ξijk)PCijk (j , k) ∈ JK D ,

(A17)

(A31) (A32)

i ∈ Ijk , s (A33)

16867



Article

xrijks

Integrated Model. The integrated model is summarized as max profit (A22−A33)

xrq ijks

s. t. scheduling model (A1−A14)

xpcijks xptijks

dynamic models (A15−A21)

■

Parameters

AUTHOR INFORMATION

Aqq′ As Bq BM PCmax ijk CF Cmax j Dj Es Gs Hs NRijk

Corresponding Author

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

The authors declare no competing financial interest.

■

NOMENCLATURE

Abbreviations

GBD MILP MINLP NLP

generalized Benders decomposition mixed-integer linear program mixed-integer nonlinear program nonlinear program

PDyn ijks

Indices

i j, j′ k, k′ (j,k) m, m′ n, n′ q, q′ r s

PRj PTfix ijk PTmax ijks PTmin ijks Ss STi Tijks Vijks Wijks XIijks

processing units jobs stages tasks iterations of outer GBD iterations of inner GBD collocation points finite elements scenarios

Sets

Ijk JK JKD JKF

state value at ending point of finite element r for dynamic model describing task (j,k) executed in unit i under scenario s discretized states at finite element r and collocation point q for dynamic model describing task (j,k) executed in unit i under scenario s product of ξijk and pcijks product of ξijk and ptijks

units capable for executing task (j,k) all tasks tasks described by dynamic models tasks having fixed recipes

XFijk

■

Binaries

collocation matrix probability weight of scenario s collocation vector sufficiently large value upper bound of pcijks total fixed cost maximum delivery time for job j parameter in order value function of job j parameter in compact form of integrated problem parameter in compact form of integrated problem weight parameter before ωs number of finite elements in dynamic model describing task (j,k) executed in unit i parameters in dynamic model describing task (j,k) executed in unit i under scenario s constant price for job j fixed processing time of task (j,k) executed in unit i upper bound of ptijks under scenario s lower bound of ptijks under scenario s parameter in compact form of integrated problem setup time of unit i parameter in compact form of integrated problem weight parameter before ptijks weight parameter before pcijks initial value of state variables in dynamic model describing task (j,k) executed in unit i under scenario s final value of state variables in dynamic model describing task (j,k) executed in unit i

REFERENCES

(1) Grossmann, I. Enterprise-wide optimization: A new frontier in process systems engineering. AIChE J. 2005, 51 (7), 1846−1857. (2) Grossmann, I. E. Advances in mathematical programming models for enterprise-wide optimization. Comput. Chem. Eng. 2012, 47, 2−18. (3) Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F. Enterprisewide modeling & optimizationAn overview of emerging research challenges and opportunities. Comput. Chem. Eng. 2007, 31 (5−6), 692−711. (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) Chu, Y.; You, F. Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and largescale global optimization algorithm. Comput. Chem. Eng. 2012, 47, 248− 268. (6) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45 (20), 6698−6712. (7) 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. (8) Prata, A.; Oldenburg, J.; Kroll, A.; Marquardt, W. Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor. Comput. Chem. Eng. 2008, 32 (3), 463−476.

ξijk equal to one when task (j,k) is assigned to unit i βijkj′k′ equal to one when task (j,k) precedes task (j′,k′) in unit i γjs equal to one if completion of job j exceeds the due date under scenario s Variables

costVs dtjs dtIjs dtIIjs lijks

total processing cost under scenario s delivery date (completion time) of job j under scenario s first component of dtjs under scenario s second component of dtjs under scenario s length of finite element for dynamic model describing task (j,k) executed in unit i under scenario s prof it profit prjs order value for executing job j under scenario s pcijks processing cost of task (j,k) executed in unit i under scenario s ptijks processing time of task (j,k) executed in unit i under scenario s saless sales under scenario s tejks ending time of task (j,k) under scenario s tsjks starting time of task (j,k) under scenario s urq discretized inputs at finite element r and collocation point ijks q for dynamic model describing task (j,k) executed in unit i under scenario s 16868



Article

(9) Chu, Y.; You, F. Integration of production scheduling and dynamic optimization for multi-product CSTRs: Generalized Benders decomposition coupled with global mixed-integer fractional programming. Comput. Chem. Eng. 2013, 58, 315−333. (10) Angel Gutierrez-Limon, M.; Flores-Tlacuahuac, A.; Grossmann, I. E. A multiobjective optimization approach for the simultaneous single line scheduling and control of CSTRs. Ind. Eng. Chem. Res. 2012, 51 (17), 5881−5890. (11) 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. (12) 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. (13) 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. (14) 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. (15) Chu, Y.; You, F. Integrated scheduling and dynamic optimization of complex batch processes with general network structure using a generalized Benders decomposition. Ind. Eng. Chem. Res. 2013, 52 (23), 7867−7885. (16) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (17) Harjunkoski, I.; Nystrom, R.; Horch, A. Integration of scheduling and controlTheory or practice? Comput. Chem. Eng. 2009, 33 (12), 1909−1918. (18) Munoz, E.; Capon-Garcia, E.; Moreno-Benito, M.; Espuna, A.; Puigjaner, L. Scheduling and control decision-making under an integrated information environment. Comput. Chem. Eng. 2011, 35 (5), 774−786. (19) 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. (20) Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S. Dynamic optimization in a discontinuous world. Ind. Eng. Chem. Res. 1998, 37 (3), 966−981. (21) Biegler, L. T.; Cervantes, A. M.; Wachter, A. Advances in simultaneous strategies for dynamic process optimization. Chem. Eng. Sci. 2002, 57 (4), 575−593. (22) Cuthrell, J. E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33 (8), 1257−1270. (23) Sahinidis, N. V. Optimization under uncertainty: state-of-the-art and opportunities. Comput. Chem. Eng. 2004, 28 (6−7), 971−983. (24) Chu, Y. F.; Wassick, J. M.; You, F. Q. Efficient scheduling method of complex batch processes with general network structure via agentbased modeling. AIChE J. 2013, 59 (8), 2884−2906. (25) Li, Z. K.; Ierapetritou, M. Process scheduling under uncertainty: Review and challenges. Comput. Chem. Eng. 2008, 32 (4−5), 715−727. (26) 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. (27) Vieira, G. E.; Herrmann, J. W.; Lin, E. Rescheduling manufacturing systems: A framework of strategies, policies, and methods. J. Scheduling 2003, 6 (1), 39−62. (28) Raheja, A. S.; Subramaniam, V. Reactive recovery of job shop schedulesA review. Int. J. Adv. Manuf. Technol. 2002, 19 (10), 756− 763. (29) Herroelen, W.; Leus, R. Robust and reactive project scheduling: a review and classification of procedures. Int. J. Prod. Res. 2004, 42 (8), 1599−1620. (30) Ouelhadj, D.; Petrovic, S. A survey of dynamic scheduling in manufacturing systems. J. Scheduling 2009, 12 (4), 417−431. (31) Ben-Tal, A.; El Ghaoui, L.; Nemirovski, A. Robust optimization; Princeton University Press: 2009.

(32) Mulvey, J. M.; Vanderbei, R. J.; Zenios, S. A. Robust Optimization of Large-Scale Systems. Oper. Res. 1995, 43 (2), 264−281. (33) Wendt, M.; Li, P.; Wozny, G. Nonlinear chance-constrained process optimization under uncertainty. Ind. Eng. Chem. Res. 2002, 41 (15), 3621−3629. (34) Zhang, Y.; Monder, D.; Forbes, J. F. Real-time optimization under parametric uncertainty: a probability constrained approach. J. Process Control 2002, 12 (3), 373−389. (35) Acevedo, J.; Pistikopoulos, E. N. Stochastic optimization based algorithms for process synthesis under uncertainty. Comput. Chem. Eng. 1998, 22 (4−5), 647−671. (36) Balasubramanian, J.; Grossmann, I. E. Approximation to multistage stochastic optimization in multiperiod batch plant scheduling under demand uncertainty. Ind. Eng. Chem. Res. 2004, 43 (14), 3695− 3713. (37) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Ross, R.; van Schijndel, J. M. G. Simultaneous design and control optimization under uncertainty. Comput. Chem. Eng. 2000, 24 (2−7), 261−266. (38) Cheng, L.; Subrahmanian, E.; Westerberg, A. W. Design and planning under uncertainty: issues on problem formulation and solution. Comput. Chem. Eng. 2003, 27 (6), 781−801. (39) Engell, S.; Markert, A.; Sand, G.; Schultz, R. Aggregated scheduling of a multiproduct batch plant by two-stage stochastic integer programming. Optim. Eng. 2004, 5 (3), 335−359. (40) Kallrath, J. Planning and scheduling in the process industry. Or Spectrum 2002, 24 (3), 219−250. (41) Maravelias, C. T. General framework and modeling approach classification for chemical production scheduling. AIChE J. 2012, 58 (6), 1812−1828. (42) Mendez, C. A.; Cerda, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art review of optimization methods for short-term scheduling of batch processes. Comput. Chem. Eng. 2006, 30 (6−7), 913−946. (43) Harjunkoski, I.; Grossmann, I. E. Decomposition techniques for multistage scheduling problems using mixed-integer and constraint programming methods. Comput. Chem. Eng. 2002, 26 (11), 1533−1552. (44) Mendez, C. A.; Cerda, J. Dynamic scheduling in multiproduct batch plants. Comput. Chem. Eng. 2003, 27 (8−9), 1247−1259. (45) Pan, C. H. A study of integer programming formulations for scheduling problems. Int. J. Syst. Sci. 1997, 28 (1), 33−41. (46) Wassick, J. M.; Ferrio, J. Extending the resource task network for industrial applications. Comput. Chem. Eng. 2011, 35 (10), 2124−2140. (47) Geoffrion, A. M. Generalized benders decomposition. J. Optim. Theory Appl. 1972, 10 (4), 237−260. (48) Floudas, C. A. Nonlinear and mixed-integer optimization: fundamentals and applications; Oxford University Press: 1995. (49) Birge, J. R.; Louveaux, F. Introduction to Stochastic Programming, 2nd ed.; Springer: New York, 2011; pp 3−+. (50) Dupacova, J.; Growe-Kuska, N.; Romisch, W. Scenario reduction in stochastic programmingAn approach using probability metrics. Math. Program. 2003, 95 (3), 493−511. (51) Heitsch, H.; Romisch, W. Scenario reduction algorithms in stochastic programming. Comput. Optim. Appl. 2003, 24 (2−3), 187− 206. (52) Heitsch, H.; Römisch, W. A note on scenario reduction for twostage stochastic programs. Oper. Res. Lett. 2007, 35 (6), 731−738. (53) Kall, P.; Mayer, J. Stochastic Linear Programming: Models, Theory, and Computation, 2nd ed.; Springer: New York, 2011; Vol. 156, pp 1−+. (54) You, F. Q.; Wassick, J. M.; Grossmann, I. E. Risk Management for a Global Supply Chain Planning Under Uncertainty: Models and Algorithms. AIChE J. 2009, 55 (4), 931−946. (55) You, F.; Grossmann, I. Multicut Benders decomposition algorithm for process supply chain planning under uncertainty. Ann. Oper. Res. 2013, 210 (1), 191−211. (56) Hairer, E.; Wanner, G. Solving ordinary differential equations II: Stiff and differential-algebraic problems; Springer: 2004; Vol. 2.

16869


Integration of Scheduling and Dynamic ... - ACS Publications

Recommend Documents