Integrated Planning, Scheduling, and Dynamic Optimization for Batch

Subscriber access provided by Penn State | University Libraries

Article

Integrated planning, scheduling, and dynamic optimization for batch processes: MINLP model formulation and efficient solution methods via surrogate modeling Yunfei Chu, and Fengqi You Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie501986d • Publication Date (Web): 01 Aug 2014 Downloaded from http://pubs.acs.org on August 11, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Integrated planning, scheduling, and dynamic optimization for batch processes: MINLP model formulation and efficient solution methods via surrogate modeling Yunfei Chu, Fengqi You* Department of Chemical and Biological Engineering, Northwestern University, Evanston, IL 60208 USA

ABSTRACT: We solve the challenging problem of integrated planning, scheduling, and dynamic optimization for sequential batch processes with fixed batch sizes. The integrated problem is first formulated into a complicated mixedinteger dynamic optimization (MIDO) problem which is then discretized into a large-scale mixed-integer nonlinear programing (MINLP) problem. There are a planning model, multiple scheduling models in planning periods, and a number of dynamic models describing task execution processes. To efficiently solve the complex MINLP problem, we develop two efficient methods which separate the subproblems using surrogate models to represent the linking functions. The first method decomposes the dynamic optimization problems from the integrated planning and scheduling problem where the surrogate models represent task processing costs dependent on the processing times. The second method further decomposes the scheduling problems from the planning problem where the surrogate models represent production costs dependent on production quantities. Compared to the direct solution approach, the proposed methods reduce the computational time by more than four orders of magnitude in the case studies.

*

To whom all correspondence should be addressed. Phone: (847) 467-2943; Fax: (847) 491-3728; E-mail: [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Page 2 of 48

Introduction Process industries are facing more challenges than ever before, such as fierce global competitions, rising

material and labor costs, growing complexity of production processes, extending scales of supply chains, rapidly changing customer markets, and increasingly stringent regulatory oversights.1 Today’s process manufacturing enterprises call for new solutions so as to survive and grow in this highly competitive international marketplace.2 To address the emerging challenges, enterprise-wide optimization (EWO) has attracted much interest.3-5 It has been recognized as a major goal in the process industries.6 However, there are still many challenges emerging in EWO, especially the computational complexity.7 We propose a novel framework to integrate three major decision levels in EWO, i.e. planning, scheduling, and dynamic optimization. Though important, the integration of the three levels is scarcely studied due to the overwhelming computational complexity of the integrated problem. EWO aims to integrate different decision layers in a manufacturing system, such as planning, scheduling, realtime optimization, and control. Conventionally, these problems are solved independently without considering the interactions among different decision layers, even though the interactions play the key role in coordinating the subsystems. As a consequence, when the interactions are taken into account, the optimal solutions for the individual subsystems may become suboptimal or even infeasible.8 To overcome the drawback of the conventional methods, the subsystems should be incorporated into an integrated system and different decision layers should be optimized simultaneously. Under the EWO framework, an increasing number of novel methods are recently presented to integrate planning and scheduling,9-17 scheduling and dynamic optimization,18-25 scheduling and control,26-28 design and control,29, 30 design and scheduling,31, 32 etc. An integrated optimization method proves to achieve a better performance than a conventional method which optimizes the subproblems independently in a sequential way.33, 34 However, the integrated method usually suffers from severe computational complexity. The integrated problem which is formulated by incorporating all subproblems is typically large-scale and much more difficult to solve than its constituent subproblems.35 Consequently, most existing integrated methods only focus on incorporating two decision layers. However, to exploit the full potential of EWO, more decision layers should be integrated and optimized collaboratively. After the integrated problem is augmented by incorporating more decision layers, its complexity increases dramatically, leading to a very complex problem that is challenging to solve. It is a major goal of this work to devise efficient solution methods to solve such a complex problem integrating more than two decision layers. In this work, we propose a novel framework to solve the integrated planning, scheduling, and dynamic optimization problem for batch processes. We consider a multi-period planning problem or a lot-sizing problem.36, 37

The planning horizon is divided into a number of planning periods. According to the order demands assigned to

the end of each planning period, the planning problem is solved to determine the production quantities of the products that should be manufactured in the planning periods. The process is then scheduled in every period to 2 ACS Paragon Plus Environment

Page 3 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fulfill the determined production quantities. In the batch processes, a product is manufactured through a series of tasks. Each task is executed in one of the capable units. The batch size for executing a task is fixed. The scheduling problem is solved to assign tasks to capable processing units, sequence task execution, and determine the task starting times and ending times. The task execution procedure is characterized by a dynamic system. Changing the input trajectories manipulates the process dynamics and, consequently, the task processing time and processing cost. They are variables determined simultaneously with the planning and scheduling decision variables by solving the integrated problem. The integrated problem is formulated as a mixed-integer dynamic optimization problem (MIDO),38, 39 which contains binary variables associated with planning and scheduling decisions, and constraints of differential equations representing the dynamic models. The differential equations are then discretized by the collocation method,40,

41

and they are transformed into a large set of nonlinear algebraic equations. The discretization

procedure reformulates the integrated problem into a large-scale mixed-integer nonlinear program (MINLP), which can be modeled by a commonly used optimization software, such as GAMS (General Algebraic Modeling System). However, the formulated MINLP is too complex to be solved by simply invoking a commercial generalpurpose MINLP solver, even a local solver. To conquer the computational complexity, we develop efficient solution methods based on surrogate modeling.42-44 The proposed methods exploit the special structure of the integrated problem to reduce the computational complexity. The integrated problem has a tree structure shown in Figure 1. It has a planning problem in the top level, a number of scheduling problems in the middle level, and a larger number of dynamic optimization problems in the bottom level. The problems in two adjacent levels are coupled by linking functions. Representing these linking functions by surrogate models breaks the coupling between the upper level problems and the lower level problems. Specifically, the lower level problems are solved to generate data points for building the surrogate models. The surrogate models are then substituted into the upper level problems replacing the detailed lower level problems. Next, the upper level problems are solved with the surrogate models. The determined variables of the upper level problems are passed to the lower level problems which are solved accordingly. We develop two solution methods shown in Figure 1. The first method decomposes the dynamic optimization problems from the planning and scheduling problem. The large set of nonlinear constraints is replaced by the simple surrogate models representing the task operational recipes. The integrated MINLP is decomposed into a mixed-integer linear program (MILP) for the planning and scheduling problem and a set of separable dynamic optimization problems. The resulting MILP is much easier to solve than the complicated MINLP. Based on the first method, the second method further decomposes the planning and scheduling problem. The scheduling problem in a planning period is solved first to generate data points for building the surrogate models, which



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 48

represent the production costs as a function of the production quantities. The generated surrogate models are used to solve the planning problem replacing the detailed scheduling problems.

Figure 1. Outline of proposed methods The proposed methods significantly reduce the computational complexity of the integrated problem. The computational efficiency is demonstrated in two case studies, including a large problem where the direct solution method fails to find a feasible solution after a long time period. The case studies also demonstrate the advantage of the integrated optimization methods in comparison to the conventional sequential method which neglects the interaction between the dynamic optimization problems and the scheduling problems. The remainder of the paper is organized as follows. The integrated problem we investigate is first specified in section 2. Next, the mathematical formulation of the integrated problem is presented in section 3. The efficient methods based on surrogate modeling are then proposed in section 4. Their performances are demonstrated in the case studies in section 5. In the end, the conclusion is given in section 6.


Page 5 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2


Problem Statement There are a great variety of planning problems and scheduling problems as well as the dynamic optimization

problems. A single integration framework cannot cover all of these problems. In this section, we specify the problems we will integrate and solve in the following sections.

Order demand

Order demand

Product A Product B Product C

Product A Product B Product C

Order demand Product A Product B Product C

Planning period Planning horizon

Planning Production quantity Product A Product B Product C

Production quantity Product A Product B Product C

Production quantity Product A Product B Product C

Scheduling horizon

Scheduling Material

Task

Job Intermediate

Task

Product

Unit 1 Unit

Unit

Unit 2 Time

Dynamic optimization

Reactor temperature

Output

Input

Time

Concentration

Time

Figure 2. Diagram of the integrated problem. The diagram of the integrated problem is exhibited in Figure 2. We consider a batch plant that manufactures different products. The planning horizon is partitioned into a set of planning periods. At the end of each period, the order demands for the products are given. According to the order demands, the planning problem is solved to 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 48

determine the production quantities of the products which should be manufactured in a planning period. For an intermediate period, backorders are allowed by incurring backorder costs. However, the total demand over the entire planning horizon is satisfied. Under these constrains, the objective is to minimize the total cost that is comprised of inventory costs, backorder costs, and production costs. When the production quantities in a period are determined, the production process is scheduled to fulfil the production quantities. We consider the batch scheduling problem where material splitting and mixing are not allowed. Thus, the batch sizes are fixed parameters. A batch of a product is manufactured through multiple processing stages. Each processing stage is called a task and the combination of all tasks associated with the batch is called a job. Obviously, the number of batches that should be manufactured for a product is equal to the number of jobs corresponding to the product. The scheduling problem is solved to determine the Gantt chart that records task assignments, execution sequences, task starting times and ending times. Determining a schedule requires information of task processing times and processing costs, which are changed by manipulating the dynamic systems. For example, a reaction task occurring in a batch reactor is given in Figure 2. The process input is the reactor temperature. Manipulating the time-dependent temperature trajectory changes the dynamic profile of the species concentration, and further the reaction time and the reaction cost. The flexible task processing times and processing costs provide degrees of freedom to coordinate the lower level dynamic systems with the upper level operational decisions. The integrated problem investigated in this work is summarized as follows. Assumptions •

The planning problem contains multiple periods and backorders are allowed for an intermediate period.

•

The batch scheduling problems have fixed batch sizes where material splitting and mixing are not allowed.

•

The dynamic models are represented by ordinary differential and algebraic equations.

Given •

Planning horizon and planning periods

•

Order demands at the end of each planning period

•

Unitary inventory cost, backorder cost and production cost

•

Operational tasks for a product

•

Capable units for executing a task

•

Dynamic models along with their parameters, initial conditions, and final conditions

Determine •

Inventory costs, backorder costs, and production costs

•

Production quantities in every planning period

•

Task assignment and execution sequence


Page 7 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


•

Task starting times and ending times

•

Task processing times and processing costs

•

Trajectories of dynamic models representing task execution processes

Objective •

To minimize total cost = inventory cost + backorder cost + production cost In the integrated problem, the initial and final conditions of a dynamic model are known parameters. Thus, a

dynamic model is linked to a scheduling model via the task processing time and processing cost.

3

Formulation of Integrated Problem Based on the problem statement in the previous section, the mathematical model of the integrated problem is

formulated in this section. The integrated model incorporates the planning model, the scheduling models in the planning periods, and the dynamic models for the tasks.

3.1 Planning model The objective is to minimize the total cost expressed by toc = cpr + civ + cbo

(1)

where toc denotes the total cost, cpr is the production cost, civ is the inventory cost, and cbo is the back order cost. The production cost is the sum of those over the planning periods, which is cpr = ∑ pct

(2)

t

where pct denotes the production cost in period t. The production cost in a period is further expressed by the sum over the tasks executed in the period as

pct = ∑ i

∑

( j , k )∈JKIi

prcijkt , ∀t

(3)

where i is the index of processing units, j is the index of jobs, k is the index of operational stages, and thus the combination (j, k) indicates the task at stage k of job j. The set JKIi contains all tasks capable to be processed by unit i. The cost of task (j, k) executed by unit i in period t is represented by prcijkt. The inventory cost is the sum over the products and the planning periods as civ = ∑∑ C pIV iv pt t

(4)

p

where a product is indexed by p, ivpt is the inventory level of product p at the end of planning period t, and C pIV is the unit inventory cost of product p. The backorder cost is expressed by



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

cbo = ∑∑ C pBO bo pt t

Page 8 of 48

(5)

p

where bopt denotes the backorder of product p at the end of planning period t, and C pBO is the unit back order cost of product p. The inventory levels satisfy the mass balance as given by

iv pt = iv p (t −1) + prpt − shpt , ∀p, t

(6)

where prpt denotes the production quantity of product p in period t, and shpt represents the shipment of the product. Similarly, the backorder is modeled by the equation below:

bo pt = bo p(t −1) + Dpt − shpt , ∀p, t

(7)

where the parameter Dpt represents the order demand of product p in period t. The backorders are allowed for an intermediate period but the total demand is satisfied at the end of the planning horizon, i.e. the total shipment is equal to the total order demand as

∑ sh

pt

t

= ∑ D pt , ∀p

(8)

t

In the scheduling problem, the batch sizes are fixed. The batch size of product p is denoted by BSp. Then, the production quantity can be expressed in terms of the number of batches (or jobs) for the product that are processed in a planning period. To count the executed jobs, a binary variable δjt is introduced. If δjt = 1, job j is executed in period t. Using the binary variables and the batch sizes, the production quantities are expressed by

prpt =

∑δ

jt

BS p , ∀p, t

(9)

j∈JPp

where set JPp contains all jobs for manufacturing product p. When a job is completed, one batch of the product is manufactured. As the batch size is fixed, the number of jobs (or batches) required for the product can be calculated from the total order demand as   JPp =  ∑ D pt BS p  , ∀p  t 

where the number of jobs for product p is represented by the size of the set JPp. The number is a ceiling function of the total demand divided by the batch size. The ceiling function returns the smallest integer no less than the inside fractional expression. We should note that the number of jobs for a product is a parameter determined from the given order demands and the batch size so that the ceiling function is not a constraint of the scheduling model. In the planning horizon, each job is executed exactly once, resulting in the following equality constraint:

∑δ

jt

= 1, ∀j

(10)

t


Page 9 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2 Scheduling models There are a number of scheduling problems in the integrated model, one in each planning period. The scheduling problems are represented by the general precedence model.8, 45 The general precedence model is a widely-used formulation for the scheduling problem. Other formulations, such as the intermediate precedence model and the time-slot model, can also be used for represent the scheduling problem.8, 46 The performance of a scheduling model is typically problem-dependent. In the context of integrated scheduling and dynamic optimization, there is no conclusive preference for a specific scheduling model. A binary variable αijkt is introduced to assign a task to a capable processing unit. When α ijkt = 1 , task (j, k) is assigned to unit i in period t. A task can be assigned only when the job that it belongs to is executed. The constraint is expressed as

∑α

ijkt

= δ jt , ∀ j , k ∈ K j , t

(11)

i∈I jk

where set Ijk contains indices of units capable for executing task (j, k), and set Kj indicates the operational stages of job j. When job j is executed in period t (δjt = 1), any of its tasks is executed exactly once in a capable units. Otherwise, δjt = 0 implies αijkt = 0. The starting time of task (j, k) in unit i in period t is denoted by tsijkt. It is zero if the task is not assigned to the unit according to

tsijkt ≤ α ijkt TTt , ∀i ∈ I jk , j , k ∈ K j , t

(12)

where TTt is the ending time of period t. When the task is executed in period t, the task starting time is larger than the ending time of the previous period according to

tsijkt ≥ α ijkt TT(t −1) , ∀i ∈ I jk , j , k ∈ K j , t

(13)

Similar to the starting time, the task ending time, denoted by teijkt, is constrained by the assignment variable as

teijkt ≤ α ijk TTt , ∀i ∈ I jk , j , k ∈ K j , t

(14)

For a job, the task in stage k can start only when its predecessor in stage k – 1 is finished, leading to the inequality

∑ ts

ijkt

i∈I jk

≥

∑

i∈I j( k −1)

teij ( k −1)t , ∀j , k ∈ K j , k ≥ 2, t

(15)

The task ending time is equal to the starting time plus the processing time as

teijkt = tsijkt + prtijkt , ∀i ∈ I jk , j , k ∈ K j , t where prtijkt denotes the processing time of task (j, k) executed by unit i in period t.


(16)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 48

The task execution sequence in a processing unit is modeled by a binary variable β ijkj ′k ′t . If βijkj ′k ′t = 1 , task (j, k) is executed before task ( j ′, k ′ ) in unit i in period t. The starting time of task ( j ′, k ′ ) is no less than the ending time of task (j, k) plus the transition time, which is enforced by tsij ′k ′t ≥ teijkt + UTi − (1 − βijkj ′k ′t ) TTt − ( 2 − α ijkt − α ij ′k ′t ) TTt , ∀i ∈ I jk ∩ I j ′k ′ , j , k ∈ K j , j ′, k ′ ∈ K j ′ ,

(17)

( j, k ) ≠ ( j′, k ′) , t

where UTi is the transition time in unit i. If β ijkj ′k ′t = 0 , task (j, k) is executed after task ( j ′, k ′ ) in unit i in period t, resulting in the inequality tsijkt ≥ teij ′k ′t + UTi − βijkj′k ′tTTt − ( 2 − α ijkt − α ij ′k ′t ) TTt , ∀i ∈ I jk ∩ I j ′k ′ , j , k ∈ K j , j ′, k ′ ∈ K j ′ ,

(18)

( j, k ) ≠ ( j′, k ′) , t

In a scheduling problem, some tasks have dynamic models to represent their processing process. A set JKDyn is introduced to represent the tasks with dynamic models. For those tasks, their processing times and processing costs are variables. The remaining tasks that do not belong to JKDyan have fixed processing times and processing costs. Their fixed processing times are calculated from

( j, k ) ∉ JK Dyn , k ∈ K j , t

prtijkt = α ijkt PTijkF , ∀i ∈ I jk ,

(19)

where the parameter PTijkF denotes the fixed processing time of task (j, k) executed by unit i. Similarly, the processing cost are calculated from

prcijkt = α ijkt PCijkF , ∀i ∈ I jk ,

( j, k ) ∉ JK Dyn , k ∈ K j , t

(20)

where the parameter PCijkF denotes the fixed processing cost of task (j, k) executed by unit i. For the tasks with dynamic models, their processing times are denoted by ptijkD and their processing costs are indicated by pcijkD . Both ptijkD and pcijkD are variables determined by the dynamic models in the following subsection. When ptijkD is known, we have

prtijkt = α ijkt ptijkD , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

(21)

Unlike eq. (19), this is a nonlinear equation as ptijkD becomes a variable. The bilinear term is linearized by the following constraints as

prtijkt ≤ ptijkD , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

prtijkt ≤ α ijkt PTijkMax , ∀i ∈ I jk ,

(22)

( j, k ) ∈ JK Dyn , t

prtijkt ≥ ptijkD − (1 − α ijkt ) PTijkMax , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

where PTijkMax is an upper bound of ptijkD . Similarly, the processing cost is


(23) (24)

Page 11 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


prcijkt = α ijkt pcijkD , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

(25)

The bilinear term is linearized by

prcijkt ≤ pcijkD , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

prcijkt = α ijkt PCijkMax , ∀i ∈ I jk ,

(26)

( j, k ) ∈ JK Dyn , t

prcijkt ≥ pcijkD − (1 − α ijkt ) PCijkmax , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

(27) (28)

where PCijkMax is an upper bound of pcijkD . The solution efficiency of the scheduling problem can be improved by introducing some symmetry breaking constraints.13, 17 The jobs belonging to a product are identical so that their relative sequence can be changed without affecting the schedule quality. The symmetry can be broken by the constraint as

∑δ t ′≤t

jt ′

≥ ∑ δ j ′t ′ , ∀p, t , j ∈ JPp , j ′ ∈ JPp , j < j ′

(29)

t ′≤t

In constraint (29), j and j ′ are indices for two jobs belonging to product p. According to the constraint, the job corresponding to the smaller index j is executed in a previous planning period before or in the same planning period with the job corresponding to the larger index j ′ . When the two jobs are assigned in the same planning period, the symmetry is broken by the constraint on the the sequence variable as

βijkj ′kt ≥ α ijkt + α ij ′kt − 1, ∀i ∈ I jk ∩ I j ′k , p, j ∈ JPp , j ′ ∈ JPp , j < j ′, k ∈ K j , t

(30)

The constraint stipulates that if the tasks in stage k of job j and job j ′ are assigned to unit i in period t ( α ijkt = α ij ′kt = 1 ), then the task of job j precedes that of job j ′ ( βijkj ′kt = 1 ).

3.3 Dynamic models Each task ( j , k ) ∈ JK Dyn has a dynamic model that characterizes the task processing process. The dynamic model is represented by a set of differential and algebraic equations as

d xijk ( t ) = fijk ( xijk ( t ) , uijk ( t ) ) dt

(31)

gijk ( xijk ( t ) , uijk ( t ) ) ≤ 0

(32)

xijk ( t = 0 ) = X ijkIni

(33)

xijk ( t = TijkFinal ) = X ijkFinal

(34)

ptijkD = TijkFinal

(35)

(

pcijkD = hijk xijk ( t = TijkFinal )

)

(36)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 48

The dynamic model is indexed by ijk. For compact expressions, the variables and functions in the dynamic model are represented by the vector forms. Eq. (31) represents the system equations where the state vector and the input vector are expressed by xijk(t) and uijk(t), respectively. Inequality (32) represents the path constraints. The initial condition and the final value of the state variables are specified by eq. (33) and eq. (34), where X ijkIni is the initial condition and X ijkFinal is the final value. The final time is denoted by TijkFinal which is equal to the task processing time according to eq. (35). The task processing cost is defined by eq. (36). To ease numerical solution of a dynamic optimization problem, the dynamic model is often discretized by a collocation method.41 The time interval is partitioned by a set of finite elements and each finite element is divided by collocation points. The state trajectory xijk(t) is then discretized by sampling values at the collocation points, ( ) denoted by xijk where f is the index of the finite elements and c is the index of the collocation points. Similarly, fc

( ) the input trajectory uijk(t) is discretized into uijk . The discretization approach transforms the differential and fc

algebraic equations (31) - (36) to a set of purely algebraic equations as follows:

( ( x(

) )≤0

( fc ) ( fc ) 0 f ijkDis xijk , uijk , xijk , xijkF = 0

gijkDis

fc ) ijk

( fc ) 0 , uijk , xijk , xijkF

(37) (38)

0 xijk = X ijkIni

(39)

xijkF = X ijkFinal

(40)

ptijkD = htijk N F

(41)

pcijkD = hijkDis ( xijkF )

(42)

The constraints (37) - (42) are the counterparts of the original model (31) - (36). The number of finite elements is represented by NF and the number of collocation points is represented by NC. The length of a finite element is a variable denoted by htijk, which has a positive value. The detailed derivation of the discretized dynamic model can be seen in the references.41, 47

3.4 Integrated model By combining the models of the planning, scheduling, and dynamic optimization problems for batch processes presented in the previous subsections, the integrated problem is formulated into a large-scale MINLP as (Integrated_MINLP)

min toc (1)

s.t. Planning model (2) – (10)


Page 13 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Scheduling models (11) – (20), (22) – (24), (26) – (30) Dynamic models (37) – (42) The integrated problem consists of one planning model, a number of scheduling modes in the planning periods, and a number of discretized dynamic models corresponding to the tasks and the processing units. Solving such a complex MINLP by a general-purpose MINLP solver might be intractable. We note that even solving a subproblem of integrated scheduling and dynamic optimization is time consuming. After a long computation time, a commercial solver returns a suboptimal solution with a large optimality gap, or even fails to find a feasible solution.20, 21 When the planning model is incorporated, the integrated problem becomes far more challenging as it includes a number of integrated scheduling and dynamic optimization problems rather than a single one.

4

Solution methods by surrogate modeling To solve the challenging integrated problem efficiently, we propose two efficient methods in this section.

These methods exploit the structure of the integrated problem and decompose the problem by using surrogate models to represent the linking functions among the subproblems. To reveal the structure of the integrated problem and facilitate the derivation of the proposed methods, we express the integrated problem (Integrated_MINLP) by a compact form as

min φ P ( y P ) + ∑ pct

(Integrated_Compact)

(43)

t

s.t.

ψ P ( yP ) ≤ 0

(44)

pct = φtS ( ytS ) +

∑ ∑α

( j , k )∈JK

(

Dyn

ijkt

i∈I jk

ψ tS pr( p =1)t , L , pr( p = N )t , ytS , P

pcijkD , ∀t

{ pt }) ≤ 0, ∀t

ψ ijkD ( ptijkD , pcijkD , yijkD ) ≤ 0, ∀i ∈ I jk ,

D ijk

( j, k ) ∈ JK Dyn

(45) (46) (47)

The variables in the planning problem are stacked into yP. The variables in the scheduling problem in period t are represented by ytS . Similarly, the variables of the dynamic models corresponding to task (j, k) executed in unit i D are presented by yijk . To manifest the linking variables between the subproblems, some original variables, like pct,

are not assimilated into the aggregate variables. Eq. (43) represents the total cost where the first component φ P ( y P ) comprises the backorder cost and the inventory cost while the second component is the sum of the production costs over the planning periods. Eq. (44) 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 48

denotes the planning model. The production cost in a planning period is expressed by eq. (45) where the first term denotes the cost from the tasks without dynamic models and the second term is the sum of costs from the tasks with dynamic models. The scheduling model in a planning period is expressed by eq. (46) where the scheduling decisions depend on the production quantities. For a compact expression, we use the curly brackets to represent a collection of the variables. The notation { ptijkD } is the set of ptijkD over the three indices. Eq. (47) represents a dynamic model for a task, which stipulates the relationship between the task processing time, the processing cost, and other variables represented by the aggregate variable.

4.1 Flexible recipe method The first method is to separate the dynamic models from the planning and scheduling model. The detailed dynamic models are replaced by surrogate models which characterize the task recipes. The integrated planning and scheduling problem is then solved with the flexible task recipes. Thus, this method is referred to as the ‘flexible recipe method’. From the compact expression (Integrated_Compact), we observe that a dynamic model is linked to the upper level model by the processing time ptijkD and the processing cost pcijkD . We also observe that the total cost is an increasing function with the processing cost if the planning and scheduling variables are fixed. Based on the observations, we formulate the integrated problem into a bi-level program as


(Integrated_Bilevel)

t

s.t.

ψ P ( yP ) ≤ 0 pct = φtS ( ytS ) +

(

∑ ∑α

( j , k )∈JK Dyn i∈I jk

ijkt

ϕijkD ( ptijkD ), ∀t

)

ψ tS pr( p =1)t , L , pr( p = N )t , ytS , { ptijkD } ≤ 0, ∀t P

ϕijkD ( ptijkD ) = min pcijkD , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn

s.t.

ψ ijkD ( ptijkD , pcijkD , yijkD ) ≤ 0 The bi-level program has one upper level problem but a number of lower level problems. The upper level problem is linked to each lower level problem by the optimal value function, which is evaluated by solving the lower level problems.48 When a lower level problem is solved, the variables of the upper level problem are regarded as parameters. Thus, the optimal value functions depend on the variables of the upper level problem, specifically the 14 ACS Paragon Plus Environment

Page 15 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


linking variables which also appear in the lower level problems. In the problem (Integrated_Bilevel), an optimal value function is represented as ϕ ijkD ( ptijkD ) , which is evaluated by minimizing the processing cost pcijkD of a dynamic system for the given processing time ptijkD . If the optimal value functions are known, the upper level problem can be solved without the detailed dynamic models. By linearizing the bilinear term α ijktϕijkD ( ptijkD ) using constraints (26)–(28), the upper level problem turns out to be an MILP. Consequently, the computational complexity is significantly reduced. However, the optimalvalue function is a black-box function without a closed-form expression. To address this difficulty, we represent the optimal-value function by a surrogate model. To build the surrogate model, we first determine the feasible range of each optimal value function. A practical dynamic system usually has a minimum transition time when it transits from an initial state to a final state. Therefore, there is a lower bound for the processing time, which is determined by solving the dynamic optimization problem below: (DO_MinPTijk)

PTijkMin = min ptijkD s.t.

ψ ijkD ( ptijkD , pcijkD , yijkD ) ≤ 0 The minimum value of ptijkD is denoted by PTijkMin . In practice, the task processing time is also bounded from above as it is not desirable to let a task occupy a processing unit for a long time. The maximum value of ptijkD is denoted by PTijkMax After determining the feasible interval of a processing time, the surrogate model is built. As this is a onedimensional function, the surrogate model can be represented by a set of sampling points, which discretize the optimal value function. For this purpose, the processing time ptijkD is first discretized by a set of sampling points PTSijks where a sampling point is indexed by s. The sampling points have the relation of

PTijjMin = PTSijk ( s =1) ≤ L ≤ PTSijks ≤ L ≤ PTSijk ( s = NS ) = PTijkMax , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn

(48)

As the sampling points are generated from a single variable, they can be the grid points as PTSijks = PTSijk ( s −1) + DTSijk ,

∀i ∈ I jk ,

( j , k ) ∈ JK Dyn , s ≥ 2

(49)

where DTSijk is the time step between two adjacent grid points. We should note that the surrogate model is built offline so that a trial-and-error procedure can be easily applied to determine the sampling points. At each PTSijks, the processing cost is minimized by solving the dynamic optimization problem given by,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(DO_MinPCijks)

Page 16 of 48

PCSijks = min pcijkD s.t.

ψ ijkD ( ptijkD , pcijkD , yijkD ) ≤ 0 ptijkD = PTSijks The optimization problems are solved over the tasks, the processing units, and the sampling points so as to generate a set of discrete points (PTSijks, PCSijks). These points are samples from the optimal value function D pcijk = ϕijkD ( ptijkD ) .

These discrete points form the candidate processing times and processing costs. The actual processing time and processing cost of a task are selected from the candidates. A binary variable γijkst is introduced for the selection. If γijkst = 1, then prtijkt = PTSijks and prcijkt = PCSijks. Representing the optimal value functions by their sampling points, the upper level problem in the bi-level program (Integrated_Bilevel) becomes (PS_SM)

min φ P ( y P ) + ∑ pct (eq. (1)) t

s.t.

ψ P ( y P ) ≤ 0 (eq. (2) – (10)) pct = φtS ( ytS ) +

(

∑ ∑ prc


ijkt

, ∀t

)

ψ tS pr( p =1)t , L , pr( p = N )t , ytS , { ptijkD } ≤ 0, ∀t (eq. (11) – (20))

∑γ

P

ijkst

= α ijkt , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn , t

s

(50)

prtijkt = ∑ γ ijkst PTSijks , ∀i ∈ I jk ,

( j , k ) ∈ JK Dyn , t

(51)

prcijkt = ∑ γ ijkst PCSijks , ∀i ∈ I jk ,

( j , k ) ∈ JK Dyn , t

(52)

s

s

D The constraints (50) – (52) formulate the surrogate models representing pcijk = ϕijkD ( ptijkD ) . The candidate

selection is applied only when the task is assigned to a processing unit, which is enforced by eq. (50). If the assignment variable αijkt is zero, then γijkst is zero. Otherwise, only one γijkst over s is equal to one, implying exactly one candidate is selected. The selected values are assigned to the actual processing times and processing costs according to eq. (51) and (52).


Page 17 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Solve

Represent

PS _ SM

pcijkD = ϕijkD ( ptijkD )

Return

by

( PTS

ijks

y P* , ytS *

, PCSijks )

ptijkD* , pcijkD*

Solve Retrieve

DO_MinPCijks

yijkD*

Return PCSijks

Solve DO_MinPTijk Return PTijkMin

Figure 3. Flow diagram of the flexible recipe method The surrogate models decompose the dynamic optimization problems from the integrated planning and scheduling problem. As a consequence, the integrated MINLP problem (Integrated_MINLP) is transformed into the MILP problem (PS_SM). The flow diagram of the flexible recipe method is exhibited in Figure 3. The dynamic optimization problems (DO_MinPTijk) are solved first to determine the minimum processing times. Starting from the minimum processing times, the candidate processing times are sampled. Then the dynamic optimization problems (DO_MinPCijks) are solved to calculate the processing costs corresponding to the processing times. Afterwards, the candidate task recipes are generated. These recipes are stored in a database along with the associated dynamic trajectories. Using the database, the integrated planning and scheduling problem (PS_SM) is solved according to the order demands. The optimal planning variables and scheduling variables are returned as well as the optimal task processing times and processing costs. Accordingly, the index of the selected recipe for a task is known. Using the index, the dynamic trajectories are retrieved from the database. We should note that all dynamic optimization problems are solved only once. When uncertainty occurs, the parameters of the integrated problem are updated according to the observation of the uncertain outcomes. Then the integrated problem with flexible task recipes is re-solved. New task recipes are selected from the database and the associated dynamic trajectories are applied. After the decomposition, all 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 48

dynamic optimization problems are solved offline and only the integrated planning and scheduling problem with flexible recipes is solved online.

4.2 Aggregate method The flexible recipe method substantially reduces the complexity of the integrated problem by separating the dynamic optimization problems from the planning and scheduling problem. However, the integrated planning and scheduling problem can be still difficult to solve when many planning periods are present. To further reduce the computational complexity, we decompose the scheduling problems from the planning problem. The decomposition is conducted using surrogate models to replace the detailed scheduling problems, leading to an aggregate planning problem. Thus, this decomposition is called the ‘aggregate method’. The method starts from the problem (PS_SM) of the flexible recipe method, which can be expressed by a bilevel program as (PS_Bilevel)

(

min φ P ( y P ) + ∑ϕ S pr( p =1)t , L, pr( p = NP )t t

)

s.t.

ψ P ( yP ) ≤ 0

(

ϕ S pr( p =1)t , L , pr( p = N

P

)t

) = min pc , ∀t t

s.t.

pct = φtS ( ytS ) +

(

∑ ∑ prc


ijkt

)

ψ tS pr( p =1)t , L , pr( p = N )t , ytS , { ptijkD } ≤ 0

∑γ

P

ijkst


( j , k ) ∈ JK Dyn

s


( j, k ) ∈ JK Dyn


( j , k ) ∈ JK Dyn

s

s

The planning problem in the upper level is linked to a scheduling problem in the lower level by the optimal value

(

)

function ϕ S pr( p =1)t , L , pr( p = N P )t , which is the minimum production cost in a planning period dependent on the production quantities of the period. Following a similar idea of the previous method, we use a surrogate model to represent the optimal value function. To build the surrogate model, we first need to determine the feasible range of the production quantities. According to eq. (9), the production quantity of a product determines the number of jobs that should be executed 18 ACS Paragon Plus Environment

Page 19 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in a planning period. If the production quantity is too large, the required jobs cannot be completed in the period, resulting in an infeasible scheduling problem. Conversely, a small production quantity implies that fewer jobs are required to be executed in the period. In this case, the scheduling problem is more likely to be feasible. Therefore, the feasible range of the production quantities is around the lower-left corner shown in Figure 4 (top). The feasible range can be closely approximated by a right triangle. The sides are the axes of the production quantities and the hypotenuse can be expressed by a straight line as

∑ AP pr p

pt

= 1, ∀t

(53)

p

where APp is a coefficient of the linear expression. Of course, the triangle is an approximate of the exact feasible range. There can be some infeasible points inside the triangle. As the production quantities have discrete values, the infeasible points inside the triangle can be easily excluded once they are detected. The parameters determining the boundary of the feasible range are calculated from solving the following optimization problem: (Scheduling_MaxPrp)

PRpMax = max prp s.t.

(

)

ψ S 0, L , prp , L , 0, y S , { ptijkD } ≤ 0 ptijkD = PTijkMin , ∀i ∈ I jk ,

( j, k ) ∈ JK Dyn

The problem is solved to find the maximum production quantity for product p, denoted by PR pMax . Because the planning periods commonly have identical length, the optimization problem can be solved in any period. For simplicity, the index t is dropped. To produce as many jobs of a product in a planning period as possible, the task processing times are set to the minimum values.


Production quantity of product B Production cost

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 48

Production quantity of product B


Figure 4. Illustration of the surrogate model representing the production cost. The partitions of the feasible range is indexed by the numbers 1 – 4.


Page 21 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The maximum production quantity PR pMax is the interception point of the plane defined by eq. (53) crossing the is the interception of the hypotenuse with the horizontal axis. axis of product p. For example, in Figure 4, PR(Max p = A) Solving the problem (Scheduling_MaxPrp) for any product p, we have Np extreme points like

( 0, L, PR

Max p

, L, 0 ) . Using these points to fit the linear equation in eq. (53), the coefficients are calculated as

APp = 1 PRpMax

(54)

The linear constraint (53) provides a simple way to define the feasible range as the scheduling problem under investigation has fixed batch sizes. After the feasible range is determined, the surrogate model of the production cost can be built. We use the piecewise linear function to form the surrogate model. The triangular feasible range can be partitioned into strips shown in Figure 4 (middle). Each partition is defined by the inequalities of

RP( r −1) < ∑ APp prpt ≤ RPr , ∀t

(55)

p

where r is the index of the partitions and RPr is a grid point dividing two partitions. According to constraint (53), these parameters are constrained by

0 = RP( r =0) ≤ L ≤ RPr ≤ L ≤ RP( r = N R ) = 1

(56)

where the number of partitions is denoted by NR + 1. In each partition r, a number of random production quantities are generated as

( PRS(

p =1) rw

, L , PRS prw , L, PRS( p = N P ) rw

)

(57)

where the sampling points are indexed by w. At each sampling point, the following scheduling problem is solved (Scheduling_MinPCrw) PDSrw = min pc s.t.

pc = φ S ( y S ) +

∑ ∑ prc


(

ijk

)

ψ S pr( p =1) , L , pr( p = N ) , y S , { ptijkD } ≤ 0 P

prp = PRS p , ∀p

∑γ

ijkst


( j , k ) ∈ JK Dyn

s


( j, k ) ∈ JK Dyn


( j , k ) ∈ JK Dyn

s

s



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 48

where index t is dropped for simplicity. Using the calculated PDSrw and PRSprw, the optimal value function of the production cost in partition r is represented by a linear function

pc = ∑ CPpr prp + DPr , ∀r

(58)

p

where the parameters CPpr and DPr are determined by fitting the data via the least squares method. To build the surrogate model, the scheduling problem (Scheduling_MinPCrw) needs to be solved repeatedly to generate data points. The scheduling problem is much more difficult to solve than the dynamic optimization problem (DO_MinPCijks) due to the presence of binary variables. For a complex industrial process, the scheduling problem itself can be very challenging to solve, even if it is solved offline.8 In this work, we focus on tackling the complexity which arises from the integration, rather than that for a single problem. The efficient scheduling methods will be the future work. We should note that the surrogate models encapsulate the scheduling problem so that it can be solved by different methods. Besides the mathematical programing method, heuristic methods, such as the agent-based scheduling methods,49, 50 can also be applied to solve a complex scheduling problem. Based on the surrogate model, the upper level problem in the bi-level program (PS_Bilevel) becomes


(Integrated_Agg)

t

s.t.

ψ P ( yP ) ≤ 0   pct = ∑  ∑ CPpr lprprt + λrt DPr , ∀t r  p 

(59)

lprprt ≤ prpt , ∀p, r ≥ 1, t

(60)

lprprt ≤ λrt PRpMax , ∀p, r ≥ 1, t

(61)

lprprt ≥ prpt − (1 − λrt ) PRpMax , ∀p, r ≥ 1, t

(62)

∑λ

(63)

rt

= 1, ∀t

r ≥1

zpt = ∑ APp prpt , ∀t

(64)

zpt = ∑ zrrt , ∀t

(65)

zrrt > λrt RPp( r −1) , ∀r ≥ 1, t

(66)

zrrt ≤ λrt RPpr , ∀r ≥ 1, t

(67)

p

r

Eq. (59) represents the production cost as a piecewise linear function of the production quantities. The binary variable λrt is introduced to denote which linear function is applied. If λrt = 1, then the production quantities of 22 ACS Paragon Plus Environment

Page 23 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


period t fall in partition r. Inside the sum of eq. (59), the variable lprprt is equal to the product λrtprpt. The bilinear term is then linearized by constraints (60)-(62). Eq. (63) enforces that exactly one component of the piecewise linear function is selected to express the production cost. The first value of λrt is zero, i.e. λ( r =0)t = 0, ∀t . The feasible range is partitioned according to the value zpt defined in eq. (64), which is expressed by the sum of the components zrrt in Eq. (65). These components are determined by the binary variables according to eq. (66) and (67). An illustration of the surrogate model is shown in Figure 4 (bottom).

Represent

(

pc = ϕ S pr( p =1) , L , pr( p = N P )

Solve

)

Integrated _ Agg Return

by

y P* , prpt*

piecewise linear function

Solve Scheduling_MinPC rw with Fixed production Solve

quantities Return

Scheduling_MinPC rw Return

ytS* , pct* PDS rw

ptijkD* , pciDjk*

Calculate cost

φ P ( y P* ) + ∑ pct* t

Solve Scheduling_MaxPrp

Retrieve yijkD*

Return PR pMin

Figure 5. Procedure of the aggregate method. The procedure of the aggregate method is shown in Figure 5. The aggregate method is based on the flexible recipe method. It starts from solving the scheduling problem (Scheduling_MaxPrp) to determine the maximum 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 48

production quantity for each product. The feasible range of the production quantities is then determined and partitioned. In each partition, sampling points for the production quantities are generated. At each sampling point, the scheduling problem (Scheduling_MinPCrw) is solved to return the production cost. The data of the production quantities and the corresponding production cost are used to build a linear response surface model in a partitioned feasible region. These functions are incorporated into a piecewise linear function, which is the surrogate model for the production cost. The piecewise linear function and the feasible range are stored in a database. The aggregate problem (Integrated_Agg) is then solved according to the order demands. The optimal production quantities are returned. They are used to solve the scheduling problem in each period. If the scheduling problem is infeasible, an integer cut is added to the surrogate model, which excludes the infeasible point. The aggregate planning problem is then solved again to determine new production quantities. When the scheduling problem is feasible, the optimal *

solution is returned. The minimum production cost pct is used to update the total cost. Because the aggregate problem uses the surrogate model to calculate the production costs, slight discrepancies can appear between the true production costs and the values estimated by the surrogate model. The discrepancies are eliminated by *

calculating the production cost pct evaluated by solving the detailed scheduling problem. Similar to the flexible recipe method, the optimal dynamic trajectories are retrieved once the selection of the task recipes is known.

5

Case studies To demonstrate the integrated methods, we apply them to a batch plant in this section. The performance of the

two proposed methods are compared with two common methods: the full space method and the fixed recipe method. The four methods compared in the case studies are visualized in Figure 6. The first existing method is the full space method, which solves the integrated MINLP (Integrated_MINLP) directly by a general-purpose solver. In this work, we use SBB to solve the MINLP as the solver is commonly used to solve integrated scheduling and dynamic optimization (or control) problem.21, 23, 26 The MINLP has a very large scale that could be very difficult to solve. The full space method encounters a severe computational burden. The second conventional method shown in Figure 6 is the fixed task recipe method, which solves the dynamic optimization problems independently from the planning and scheduling problem. This method is similar to the flexible recipe method; however, it neglects the interactions among the dynamic optimization problems and the upper level operational problems. The task recipes are fixed parameters which cannot be optimized simultaneously with the planning and scheduling decisions. Though it is easy to implement, the fixed recipe method returns an inferior solution compared to other methods. By contrast, the interactions between the dynamic optimization problems and the upper level problems are represented by the surrogate models in the flexible recipe method. Therefore, the subproblems can be optimized collectively, achieving a better solution than the fixed recipe method. Based on the flexible recipe 24 ACS Paragon Plus Environment

Page 25 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method, the aggregate method further decomposes the scheduling problems from the planning problem. New surrogate models are built to represent the interactions between the planning problem and the scheduling problems.

Figure 6. Diagram of different methods solving the integrated problem

In the case studies, all optimization problems are modeled by GAMS 2451 in a computer with an Intel CPU (3.10 GHz) with 4 cores and 8 GB memory. The operating system is Windows 7 Professional (64bit).

5.1 Problem description The diagram of the batch plant is shown in Figure 7. There are 8 processing units in the batch plant: a mixing tank (TI), four batch reactors (RI – RIV), one filter (FI), and two separators (SI and SII). The batch plant manufactures 3 products, denoted by PI, PII, and PIII, respectively. A job for processing a product from raw materials consists of 5 operational stages. The task of the first stage is mixing, executed in the mixing tank TI. The second stage is a reaction task executed in either reactor RI or reactor RII. The outlet materials of the reaction go through the filtration task executed in the filter FI in the third stage. The fourth stage is another reaction task, 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 48

which can be executed in either reactor RIII or reactor RIV. The following separation task at the fifth stage can be executed in either separator SI or separator SII.

r2 C  →D

r1 A + B  →C

Figure 7. The batch plant for case studies In Figure 7, the chemical reactions are expressed by a general form for different products. The species are actually associated with the products. For example, the final species for product PI is DPI and that for product PII is DPII . However, to simplify the following notations, the subscripts of the products are omitted. All reaction tasks executed in the four reactors are described by dynamic systems. The dynamic model which represents the reaction task taking place in a batch reactor RI or RII is Reaction

r1 A + B  →C

Reaction rate

r1 ( t ) = k1e

Differential equations

dC A ( t ) dt dCB ( t ) dt dCC ( t ) dt

Initial conditions

− ER1 TR1 ( t )

C A ( t ) CB ( t )

= − r1 ( t ) = − r1 ( t ) = r1 ( t )

C A ( t = 0 ) = C A0 CB ( t = 0 ) = CB0 CC ( t = 0 ) = CC0


Page 27 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Final conditions

C A ( t = t F ) = C AF C B ( t = t F ) = CBF CC ( t = t F ) = CCF

Processing time

pt = t F

Processing cost

pc = C P ∫

tF

0

(TR ( t ) − TR ) dt 0

1

The reactants are represented by A and B while the product is denoted by C. The reaction rate r1 follows the Arrhenius equation where the rate constant depends on reactor temperature TR1. The three differential equations delineate the change of concentrations CA, CB, and CC using the reaction rate. The initial conditions of the systems are specified by the parameters CA0 , CB0 , and CC0 . The final concentrations are specified by C AF , CBF , and CCF . The reaction time is denoted by tF, which is also the processing time of the reaction task. The processing cost is calculated by the integral of the reactor temperature.21, 28 The unitary cost is indicated by CP. The trajectory of the reactor temperature TR1(t) is optimized by solving the integrated problem. The optimized trajectory is then imported into a control system as the reference. The objective of the control system is to manipulate process inputs, e.g. the flow rates of the heating fluid and the cooling fluid, so as to track the reference trajectory. To facilitate controller design, the reference trajectory is often required to have a regular profile without abrupt changes and frequent oscillations. For this purpose, the reference trajectory of the reactor temperature is typically designed as a trapezoid.52, 53 A typical example is shown in Figure 8. The reaction period is divided into 3 phases with the ending times of tI, tII, and tF. During the first heating phase in [0, tI], the temperature rises gradually from the initial value TR0 to the maximum value TRM. During the second isothermal phase in [tI, tII], the temperature stays at the maximum value TRM. During the third cooling phase in [tII, tF], the temperature drops from TRM back to TR0. The initial and final temperature TR0 is a parameter equal to the ambient temperature. The values of variables TRM, tI, tII, and tF are determined by solving the integrated problem.

TR M

TR 0 0

tI

t II

Figure 8. Typical trajectory of the reactor temperature.


tF


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 48

Similarly, the dynamic model which represents the reaction task taking place in a batch reactor RIII or RIV is Reaction

r2 C  →D

Reaction rate

r2 ( t ) = k2 e

Differential equations

dCC ( t ) dt dCD ( t ) dt

Initial conditions

− ER2 TR2 ( t )

CC ( t )

= − r2 ( t ) = r2 ( t )

CC ( t = 0 ) = CC0 CD ( t = 0 ) = CD0

Final conditions

CC ( t = t F ) = CCF C D ( t = t F ) = CDF

Processing time

pt = t F

Processing cost

pc = C P ∫

tF

0

(TR ( t ) − TR ) dt 0

2

The reaction rate follows the law of mass reaction and the Arrhenius equation. The processing time and the processing cost are calculated in the same way as those for reactor RI or RII. The parameters of the dynamic models describing the reaction tasks are given in Appendix.

5.2 Case 1: Small problem of 2 products and 2 periods We start with a small-scale problem having 2 products and 2 planning periods. The duration of each planning period is 12 hours. The order demands at the end of the planning periods are listed in Table 1. Because the batch sizes are fixed, the order demands can be quantified directly by the number of batches, which are the original demand values divided by the batch sizes. For generality, we allow the order demands to be fractional values. Backorders are allowed in an intermediate planning period by incurring backorder costs in the objective function. However, the total demand of every product in the entire planning horizon should be satisfied. According to the total demand of a product, the number of batches, which should be executed in the planning horizon, is calculated. It is equal to the number of jobs corresponding to the product. Specifically, the number of jobs for a product is the smallest integer no less than the total demand represented in terms of batches.

Table 1. Order demands of the small problem. (The demand unit is batch.) Product Period 1 Period 2 Total PI 3.2 0.5 3.7 PII 0.2 3.6 3.8 28 ACS Paragon Plus Environment

Job index 1, 2, 3, 4 5, 6, 7, 8

Page 29 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


According to the procedures of the proposed methods, the dynamic optimization problems are solved first. As an example, the results of minimizing the processing times of the two reaction tasks for product A in two reactors are shown in Table 2. The dynamic models are discretized by the collocation method with 60 finite elements. Each finite element includes 3 collocation points. After the discretization, the dynamic optimization problem becomes an NLP which is solved by CONOPT 3. When being solved independently, each dynamic optimization problem is easy to solve in less than a second. For each reaction task, the candidate recipes are generated. Starting from the minimum processing time, 11 processing times are sampled with a step size equal to 0.1 hour. Because the dynamic optimization problems are solved offline, the number of sampling points and the step size for a processing time are determined by a trialand-error approach. In this case, they are determined when a finer discretization procedure does not result in a noticeable improvement. At each processing time, a dynamic optimization problem is solved to minimize the corresponding processing cost. In this way, 11 candidate recipe points of (processing time, processing cost) are generated for the task. The number of the dynamic optimization problems solved to generate all task candidate recipes is 2 (# products) × 2 (# reaction tasks) × 2 (# capable units) × 11 (recipe candidates) = 88. Since the dynamic optimization problems can be solved independently, the total computational time is less than one minute. We should note that these dynamic optimization problems are solved only once. When the order demands are changed, the same candidate recipes are used. Of course, the final selected recipes can be different, depending on the scheduling and planning problem.

Table 2. Results of solving the dynamic optimization problems which minimize the processing times of the reaction tasks for product PI in reactors RI and RIII respectively. Reactor RI RIII Obj. (h) 1.6 1.4 CPU (s) 0.5 0.3 Type NLP NLP Constraints 1,988 1,628 Continuous variables 1,991 1,631 Solver CONOPT 3 CONOPT 3 Based on the generated task recipes, the flexible recipe method can be applied where the task recipes are selected simultaneously with the scheduling and planning decisions. For the aggregate method, the surrogate models representing the production costs dependent on the production quantities are required. They are built from the sampling data. For this simple problem, the total number of all feasible solutions is small so we can exhaustively explore every combination of the production quantities for the two products. The feasible range is shown in Figure 9, which is a triangle at the bottom left corner. The maximum production quantity of either product is 5 batches or 5 jobs in a planning period. At each feasible point, the scheduling problem with the 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 48

flexible task recipes is solved and the minimum production cost is returned and exhibited in Figure 9. It is seen that the production cost increases as the sum of the production quantities increases.

1200 1000 800 600 400 200 0

Figure 9. The function of the production cost evaluated with the flexible task recipes in the small problem. Using the flexible task recipes, the scheduling problem can be solved by simultaneously optimizing the task operations. To reveal the advantage of the flexible task recipes, the scheduling problem with the fixed task recipes is solved to evaluate the production cost. The fixed task recipes are those with the minimum processing times. Compared to the fixed task recipes, the reduction in the production cost evaluated by the flexible recipes is visualized in Figure 10. The average reduction is 19.9%, revealing the advantage of the flexible task recipes.

0.3 0.25 0.2 0.15 0.1 0.05 0

Figure 10. The relative reduction in the procduction costs evaluated by the flexible task recipes against those evaluated by fixed task recipes. 30 ACS Paragon Plus Environment

Page 31 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fitting the production cost data shown in Figure 9, the surrogate model is built. We partition the feasible range into two regions. One region contains the production quantities such that their sum is less than or equal to 4 and the other region consists of those such that the sum of the production quantities is equal to 5. In each region, a linear response surface model is built to fit the data.54 The data calculated by solving the scheduling problem against the data estimated by the fitted surrogate model are displayed in Figure 11. Obviously, the surrogate model accurately represents the functional relationship between the production cost and the production quantities. The relative fitting error is only 0.61%.

1400 Calculated 1200

Estimated

1000 800 600 400 200 0

Figure 11. Prediction of the surrogate model which fits the production cost data in the small problem.

Table 3. Comparison of the 4 methods solving the 2-product and 2-period problem Method Obj. (m.u.) CPU (s) Type Constraints All variables Binary variables Solver Gap (%)

Full space 2,411.2 180,000 MINLP 61,033 59,517 1,040 SBB 47.9

Fixed recipe 2,691.1 2.6 MILP 2,969 1,611 1,040 CPLEX 12 0

Flexible recipe 2,215.7 17.6 MILP 3,033 2,315 1,744 CPLEX 12 0

Aggregate note 1 2,215.7 1.8 MILP 111 + 1,493 × 2 83 + 1,144 × 2 20 + 864 × 2 CPLEX 12 0

Note. 1. The computational time of the aggregate method includes the time for solving the aggregate planning problem and the time for solving the scheduling problems with flexible recipe in the planning periods according to the production quantities determined by the aggregate problem. When the production quantities are determined, the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 48

scheduling problems are solved independently. The model statistics include those for the aggregate problem (before “+”) and those for the two scheduling problems (after “+”). 2. The gap for the full space method is the one for the integrated MINLP problem when the solver terminates. The gaps for the remaining three methods are the ones for the corresponding MILPs which are derived from the MINLP problem. These gaps are not the ones for the integrated problem.

Based on the surrogate models which represent the linking relationship between the production levels, the proposed methods are applied to solve the integrated problem. To verify their performances, the integrated problem is also solved by the conventional methods. Table 3 shows the results. Though the individual planning and scheduling problems are small, the integrated problem including dynamic models turns out to be a complicated MINLP, which has about sixty thousand variables and equations along with more than one thousand binary variables. Such a complicated MINLP fails to be solved efficiently by the commercial general-purpose MINLP solver, SBB, which is widely used to solve integrated scheduling and dynamic optimization (or control) problem.23, 55 Another common MINLP solver DICOPT did not find a feasible solution. Due to the complexity of the integrated model, SBB spent 50 hours in merely finding a suboptimal solution. The gap returned by the solver is as large as 47.9%. As a consequence, the cost returned by the full space method is larger than the proposed methods. However, though suboptimal, the solution returned by the full space method is still smaller (by 10.4%) than that of the fixed recipe method, because the full space method optimizes the task execution recipes simultaneously with the scheduling and planning decisions. To further analyze the solutions searched by the two methods, the Gantt charts returned by the full space method and the fixed recipe method are shown in Figure 12. The two methods result in different schedules in a planning period and even different production quantities assigned to a period. The differences are attributed to the fact that the task recipes are optimized simultaneously with upperlevel operational decisions in the full space method. By contrast, the fixed recipe method lacks such coordination between the dynamic optimization problems and the scheduling problems. Obviously, when the task recipes are changed by the full space method, a substantial reduction in the production cost is obtained. Furthermore, the reduction in the production cost causes the smaller total cost. Comparison of the full space method and the sequential method reveals the advantages of integrated optimization approaches.


Page 33 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12. Gantt charts for the small problem solved by different methods. The number in a task bar is the job index. The tasks for product PI are marked in blue while the tasks for product PII are marked in green.



The performance of integrated optimization approaches can be further improved by the proposed methods. The flexible recipe method first solves the dynamic optimization problems to generate task recipes, which are used as the surrogate models replacing the detailed dynamic models. The integrated problem turns out to be a planning and scheduling problem with flexible task recipes. This is an MILP much easier to solve than the full-space MINLP problem. The MILP is solved in just 17.6 seconds to the zero optimality gap. The cost returned by the flexible recipe method is smaller (by 8.1%) than the full space method while the computational time is reduced by more than four orders of magnitude. The Gantt chart of the flexible recipe method is displayed in Figure 12. The computational time can be further reduced by the aggregate method which uses an additional high-level surrogate model to decompose the scheduling problems from the planning problem. The integrated problem is reduced to a simple MILP problem. The aggregate method determines the production quantities assigned to each planning period. According to the production quantities, the scheduling problem with flexible task recipes in a planning period is solved to determine the detailed schedule and to select the task recipes. The aggregate method takes only 1.8 seconds in total to solve the integrated problem. The efficiency comes from the fact that the scheduling problem in a planning period can be solved independently from those in other periods after the

320

380

315

360

Temperature (K)

Temperature (K)

production quantities are given.

310 305 300

340 320 300

295 0

1

0

2

0.5

1

1.5

Time (hour)

Time (hour) 320

320

315

315

Temperature (K)

Temperature (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 48

310 305 300 295

310 305 300 295

0

1

2

0

Time (hour)

1

2

Time (hour)

Figure 13. The temperature profiles of the first reaction task for job 3, which are returned by the four methods.


Page 35 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It is seen from Figure 12 that the four compared methods determine different task operational recipes. As an example, the temperature profiles of the first reaction task for job 3 are displayed in Figure 13. The fixed recipe method minimizes the task processing time, incurring a large amount of heat exchange and a large processing cost. By contrast, the remaining three methods optimize the dynamic profile collaboratively with the planning and scheduling decisions. The processing time is prolonged so that a small amount of heat exchange is required, resulting in a small processing cost. However, a long processing time along with a small processing cost is not always a good choice. The long processing time can delay the product delivery time and incur a backorder cost. It is also seen from Figure 12 that all four methods return a small processing time for the second reaction task for job 3. This reveals the fact that we cannot predetermine an optimal task operational recipe independently from the planning and scheduling problem. Because the integrated methods, including the full space method, the fixed recipe method, and the aggregate method, can optimize the task recipes with planning and scheduling decisions, they are able to improve the performance of the entire production system compared to the conventional sequential method.

5.3 Case 2: Large problem of 3 products and 5 periods In this problem, we apply and compare the four methods using a large problem, which includes 3 products and 5 periods. In addition, the length of a planning period is increased to 24 hours, implying a larger scheduling problem in a period than the small problem. The order demands for the problem is shown in Table 4. According to the total demand of a product, the number of jobs required to fulfill the order for the products is determined. There are totally 55 jobs created to meet the order demands.

Table 4. Order demands of the large problem. (The demand unit is batch) Product Period 1 Period 2 Period 3 Period 4 Period 5 PI 6.2 4.5 2.6 2.0 2.3 PII 3.2 3.6 3.5 5.3 3.2 PIII 1.8 2.0 4.6 3.0 6.2

Total 17.6 18.8 17.6

Job index 1 ─ 18 19 ─ 37 38 ─ 55

As the problem is augmented without changing the dynamic models, the same candidate task recipes generated in the small problem can be used to solve the large problem. There is no need to solve the dynamic optimization problems again. This is an evident advantage of the proposed methods. By contrasts, the full space method has to re-optimize all the dynamic models according to the changed planning and scheduling configurations, which is cumbersome considering that these configurations, e.g. the order demands, are frequently varied.



4000

Production cost (m.u.)

Calculated Estimated 3000

2000

1000

0 0

5

10

15

20

25

30

35

40

35

40

Index of sampling point

4000 Calculated

Production cost (m.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 48

Estimated 3000

2000

1000

0 0

5

10

15

20

25

30

Index of sampling point

Figure 14. Fitting errors (top) and prediction errors (bottom) for the surrogate model representing the production cost in the large problem. As the planning period is prolonged and an additional product is added, a new surrogate model of the production cost needs to be built for the aggregate method. The largest number of jobs for a product to be processed in the planning period is 15. The feasible range is partitioned into 3 regions. In the first region, the sum of the production quantities is less than 12. The linear response surface model in this region is built from 20 randomly sampled data points. In the second region, the sum of the production quantities is equal to 13 or 14. Ten data points are sampled from the region to build the linear response surface model. In the third region, the sum of the production quantities is equal to 15. The linear response surface model is built by 10 randomly sampled data points. Incorporating the linear response surfaces in the three regions, the piecewise linear surrogate model is formed. The fitting errors of the surrogate model on the 40 data points are shown in Figure 14 (top). These data points belong to the training set as they are used to fit the model. To further validate the surrogate model, a new set of 40 data points are randomly generated, which belong to the validation set. The prediction errors of the 36 ACS Paragon Plus Environment

Page 37 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


surrogate model over the validation set are shown in Figure 14 (bottom) as well. The average of relative fitting errors over the training set is only 0.8%, and the average over the validation set is 1.0%. The small relative errors indicate an accurate representation of the production cost by the surrogate model. The proposed methods solve the integrated problem and the results are compared to the two conventional methods in Table 5. For this large example, the integrated problem is inordinately complicated with more than half a million equations and variables along with more than one hundred thousand binary variables. For such a complex MINLP, the local solver SBB fails to find a feasible solution in 50 hours. The fixed recipe method solves the problems using about 15 minutes. As the dynamic models are replaced by the fixed processing times and processing costs, the integrated problem reduces to an MILP, accelerating the solution procedure. However, the fixed recipe method does not coordinate the task operations with the scheduling and planning decisions, leading to a large total cost. The sub-optimality of the fixed recipe method is evidenced by the proposed methods. The flexible recipe method returns a smaller cost by 19.1% than the fixed recipe method. The computational time of the flexible recipe method is about 25 minutes, which is short considering that a very complicated problem is solved. As the number of jobs in this problem is large, there are a number of binary variables in the scheduling models. However, the number of products is small. The symmetry breaking constraints can significantly improve the computational efficiency because the jobs belonging to a product are identical and their relative sequence is predetermined by the symmetry breaking constraints.

Table 5. Comparisons of methods solving the five-period problem Method Full space Fixed recipe Obj. (m.u.) Infeasible 16,712.5 CPU (s) 180,000 870.4 Type MINLP MILP

Flexible recipe 13,514.8 1,448.5 MILP

Constraints

676,611

274,887

275,987

All variables

528,225

130,241

142,341

Binary variables

121,275

121,275

133,375

Solver Gap (%)

SBB

CPLEX 12 1.0

CPLEX 12 1.0

Aggregate 13,593.6 119.6 MILP 2,682 + 54,704 × 5 486 + 28,451 × 5 290 + 26,620 × 5 CPLEX 12 1.0

Note. 1. The table is formatted in the same way with Table 3.

The shortest computational time is achieved by the aggregate method. In this large problem, as the planning problem and the scheduling problems become complicated, the efficiency of the aggregate method becomes more



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 48

evident than that in the small problem. The computational time is only about two minutes, which is much shorter than that of the flexible recipe method and even that of the fixed recipe method, because those two methods solve the planning problem simultaneously with all scheduling problems. Meanwhile, the objective function value of the aggregate method is slightly larger than that of the flexible recipe method by 0.6%. The difference stems from the approximation inaccuracy of the surrogate model. A finer surrogate model can further reduce the difference.

Figure 15. Gantt charts returned by different methods in the large problem. The jobs corresponding product PI, PII, and PIII are marked in blue, green, and red respectively. 38 ACS Paragon Plus Environment

Page 39 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The Gantt charts of the fixed recipe method, the flexible recipe method, and the aggregate method are shown in Figure 15. The backorder costs returned by the three methods are all zero. The inventory costs are also the same. However, the production costs returned by the fixed recipe method are larger than those returned by the proposed methods since the task operations are not optimized simultaneously with the planning and scheduling decisions. It is seen from the Gantt charts that the fixed recipe method does not make an effectual use of the processing units and there are more idle time periods in a unit than the proposed methods. The proposed methods collectively optimize the dynamic systems with planning and scheduling decisions via the flexible task recipe. The idle periods of a processing unit are reduced by extending the task processing times meanwhile reducing the processing costs.

6

Conclusion Integration of planning, scheduling, and dynamic optimization results in a complex large-scale MINLP, which

includes a planning model, multiple scheduling models, and a number of dynamic models. Because of this production cost complexity, direct solution methods by general-purpose MINLP solvers are time-consuming. We proposed two efficient methods to solve the integrated problem. The integrated problem is decomposed by using surrogate models to represent the linking functions among the subproblems. Specifically, the flexible recipe method uses surrogate models to represent the flexible task recipes which replace the detailed dynamic models. The integrated MINLP problem is transformed into an MILP which is the integrated planning and scheduling with flexible task recipes. The aggregate method further decomposes the scheduling problems from the planning problem by a surrogate model representing the production cost dependent on the production quantities in a planning period. The advantage of the proposed methods was demonstrated in two case studies. In the small problem, the proposed methods returned a lower cost than the full space method which only found a suboptimal solution in 50 hours. Besides, the proposed methods reduced the computational times by more than four orders of magnitude. Compared with the fixed recipe method, the cost returned by the proposed methods is lowered by 17.7% because the proposed methods were able to coordinate the dynamic optimization problems with the planning and scheduling decisions. In the large problem with more than half million constraints and variables, the full space method failed to return a feasible solution in 50 hours. However, the proposed methods returned the solutions in less than 25 minutes. The returned costs were smaller than that of the fixed method by 18.7%.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 40 of 48

Page 41 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nomenclature Abbreviation MIDO

mixed integer dynamic optimization

MILP

mixed integer linear program

MINLP

mixed integer nonlinear program

Index c

collocation point

f

finite element

i

processing unit

j

job

k

operational stage

(j, k)

task

p

product

r

partition for feasible range of production quantities

s

sampling point of processing cost

t

planning period

w

sampling point of production cost

Set units capable for executing task (j, k)

Ijk JK

Dyn

tasks with dynamic models

JKIi

tasks capable to be processed by unit i

JPp

jobs belonging to product p

Kj

operational stages of job j

Parameter APp

coefficient of surrogate model for product p

BSp

batch size of product p

C pBO

unitary backorder cost of product p

C pIV

unitary inventory cost of product p

CPpr

coefficient of surrogate model for product p in partition r

Dpt

order demand of product p at the end of period t 41 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

DPr

parameter of surrogate model partition r

DTSijk

time step for discretizing ptijkD

NC

number of collocation points

NF

number of finite elements

NP

number of products

NR

number of partitions for surrogate model

NS

number of candidate task recipes

PCijkF

fixed processing cost of task ( j , k ) ∉ JK Dyn executed by unit i

PCijkMax

upper bound of pcijkD

PCSijks

D discretized pcijk at sampling point s

PDSrw

discrete production cost at sampling point w in partition r

PR pMax

maximum value of prp

PRSprw

discrete production quantity of product p at sampling point w in partition r

PTijkF

fixed processing time of task ( j , k ) ∉ JK Dyn executed by unit i

PTijkMax

upper bound of ptijkD

PTijkMin

D lower bound of ptijk

PTSijks

D discretized ptijk at sampling point s

RPr

grid value of partition r

TTt

ending time of period t

UTi

transition time of unit i

X ijkFinal

final value for dynamic model representing task (j, k) executed by unit i

X ijkIni

initial condition for dynamic model representing task (j, k) executed by unit i

Binary variable αijkt

=1, if task (j, k) is assigned to unit i in period t

β ijkj ′k ′t

=1, if task (j, k) is executed before task ( j ′, k ′ ) in unit i in period t

λrt

=1, if production quantities in period t belong to partition r

δjt

=1, if job j is executed in period t

Continuous variable 42 ACS Paragon Plus Environment

Page 42 of 48

Page 43 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bopt

backorder of product p at the end of period t

cbo

backorder cost

civ

inventory cost

cpr

production cost

htijk

length of finite elements for dynamic model representing task (j, k) executed by unit i

ivpt

inventory level of product p at the end of period t

lprprt

product of λrt and prpt

pct

production cost in period t

pcijkD

variable processing cost of task ( j , k ) ∈ JK Dyn executed by unit i

prpt

production quantity of product p in period t

prcijkt

processing cost of task (j, k) executed by unit i in period t

prtijkt

processing time of task (j, k) executed by unit i in period t

ptijkD

variable processing time of task ( j , k ) ∈ JK Dyn executed by unit i

shpt

shipment of product p at the end of period t

toc

total cost

teijkt

ending time of task (j, k) by unit i in period t

tsijkt

starting time of task (j, k) by unit i in period t

uijk

input vector of dynamic model representing task (j, k) executed by unit i

( ) uijk

discretized uijt at collocation point c of finite element f

xijk

state vector of dynamic model representing task (j, k) executed by unit i

( ) xijk

discretized xijt at collocation point c of finite element f

yP

aggregate of planning variables

ytS

aggregate of variables for scheduling model in period t

yijkD

aggregate of variables for dynamic model representing task (j, k) executed in unit i

yijkD

aggregate of variables for dynamic model representing task (j, k) executed in unit i

zpt

sum of production quantities in period t

zrrt

component of zpt in partition r

fc

fc



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 48

Appendix The appendix lists the data for the case studies.

Table A1. Parameters of dynamic models describing reaction tasks in reactor RI or RII PI

PII

PIII

k1 (m /kmol h)

0.2

5

1

ER1 (K)

1000

2000

1500

CA0 (kmol/m3)

1000

1000

1000

CB0 (kmol/m3)

1000

1000

1000

CC0 (kmol/m3)

0

0

0

C AF (kmol/m3)

50

50

50

CBF (kmol/m3)

50

50

50

CCF (kmol/m3)

950

950

950

TR0 (K)

300

300

300

Product 3

Table A2. Parameters of dynamic models describing reaction tasks in reactor RIII or RIV Product

PI

PII

PIII

k2 (1/h)

600

1000

3000

ER2 (K)

2000

2200

2500

CC0 (kmol/m3)

950

950

950

CD0 (kmol/m3)

0

0

0

CCF (kmol/m3)

47.5

47.5

47.5

CDF (kmol/m3)

902.5

902.5

902.5

TR0 (K)

300

300

300


Page 45 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table A3. Parameters of a scheduling model Symbol

Description

Value

Unit

C pBO

Unitary backorder cost of product p

1000

m.u./batch

C pIV

Unitary inventory cost of product p

200

m.u./batch

CP

Unitary utility cost

1

m.u./K

PCijkF

Fixed cost of mixing task

10

m.u.

PCijkF

Fixed cost of filtration task

30

m.u.

PCijkF

Fixed cost of separation task

100

m.u.

PTijkF

Fixed processing time of mixing task

1.0

hour

PTijkF

Fixed processing time of filtration task

1.0

hour

PTijkF

Fixed processing time of separation task

1.5

hour

UTi

Transition time in TI

0.2

hour

UTi

Transition time in RI and RII

0.5

hour

UTi

Transition time in FI

0.2

hour

UTi

Transition time in RIII and RIV

0.4

hour

UTi

Transition time in SI and SII

0.3

hour



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 48

Reference 1. Grossmann, I., Enterprise-wide optimization: A new frontier in process systems engineering. AIChE Journal 2005, 51, (7), 1846-1857. 2. Christofides, P. D.; Davis, J. F.; El-Farra, N. H.; Clark, D.; Harris, K. R. D.; Gipson, J. N., Smart plant operations: Vision, progress and challenges. Aiche Journal 2007, 53, (11), 2734-2741. 3. Quaglia, A.; Sarup, B.; Sin, G.; Gani, R., A systematic framework for enterprise-wide optimization: Synthesis and design of processing networks under uncertainty. Computers & Chemical Engineering 2013, 59, 47-62. 4. Wassick, J. M., Enterprise-wide optimization in an integrated chemical complex. Computers & Chemical Engineering 2009, 33, (12), 1950-1963. 5. Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F., Enterprise-wide modeling & optimization - An overview of emerging research challenges and opportunities. Computers & Chemical Engineering 2007, 31, (5-6), 692-711. 6. Grossmann, I. E., Advances in mathematical programming models for enterprise-wide optimization. Computers & Chemical Engineering 2012, 47, 2-18. 7. 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. Computers & Chemical Engineering 2012, 47, 157-169. 8. Méndez, C. A.; Cerdá, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M., State-of-the-art review of optimization methods for short-term scheduling of batch processes. Computers & Chemical Engineering 2006, 30, (6), 913-946. 9. Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A., Planning and Scheduling under Uncertainty: A Review Across Multiple Sectors. Industrial & Engineering Chemistry Research 2010, 49, (9), 3993-4017. 10. Tan, W.; Khoshnevis, B., Integration of process planning and scheduling - a review. Journal of Intelligent Manufacturing 2000, 11, (1), 51-63. 11. Chu, Y.; You, F.; Wassick, J. M.; Agarwal, A., Integrated planning and scheduling under production uncertainties: Bi-level model formulation and hybrid solution method. Computers & Chemical Engineering 2014, DOI: http://dx.doi.org/10.1016/j.compchemeng.2014.02.023. 12. Li, Z.; Ierapetritou, M. G., Integrated production planning and scheduling using a decomposition framework. Chemical Engineering Science 2009, 64, (16), 3585-3597. 13. Yue, D. J.; You, F. Q., Planning and Scheduling of Flexible Process Networks Under Uncertainty with Stochastic Inventory: MINLP Models and Algorithm. AIChE Journal 2013, 59, (5), 1511-1532. 14. Erdirik-Dogan, M.; Grossmann, I. E., Simultaneous planning and scheduling of single-stage multiproduct continuous plants with parallel lines. Computers & Chemical Engineering 2008, 32, (11), 2664-2683. 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. Computers & Chemical Engineering 2014, DOI: 10.1016/j.compchemeng.2014.02.023. 16. Birewar, D. B.; Grossmann, I. E., Simultaneous Production Planning and Scheduling in Multiproduct Batch Plants. Industrial & Engineering Chemistry Research 1990, 29, (4), 570-580. 17. Erdirik-Dogan, M.; Grossmann, I. E., A decomposition method for the simultaneous planning and scheduling of single-stage continuous multiproduct plants. Industrial & Engineering Chemistry Research 2006, 45, (1), 299-315. 18. Chu, Y.; You, F., Integrated Scheduling and Dynamic Optimization of Sequential Batch Processes with Online Implementation. AIChE Journal 2013, 59, (7), 2379-2406. 19. 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. Industrial & Engineering Chemistry Research 2013, 52, (23), 7867-7885. 20. Mishra, B. V.; Mayer, E.; Raisch, J.; Kienle, A., Short-term scheduling of batch processes. A comparative study of different approaches. Industrial & Engineering Chemistry Research 2005, 44, (11), 4022-4034. 46 ACS Paragon Plus Environment

Page 47 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


21. Nie, Y. S.; Biegler, L. T.; Wassick, J. M., Integrated scheduling and dynamic optimization of batch processes using state equipment networks. Aiche Journal 2012, 58, (11), 3416-3432. 22. Capon-Garcia, E.; Guillen-Gosalbez, G.; Espuna, A., Integrating process dynamics within batch process scheduling via mixed-integer dynamic optimization. Chemical Engineering Science 2013, 102, 139-150. 23. Flores-Tlacuahuac, A.; Grossmann, I. E., Simultaneous cyclic scheduling and control of a multiproduct CSTR. Industrial & Engineering Chemistry Research 2006, 45, (20), 6698-6712. 24. Chu, Y.; You, F., Integration of Scheduling and Dynamic Optimization of Batch Processes under Uncertainty: Two-Stage Stochastic Programming Approach and Enhanced Generalized Benders Decomposition Algorithm. Industrial & Engineering Chemistry Research 2013, 52, (47), 16851-16869. 25. 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. Computers & Chemical Engineering 2013, 58, (0), 315-333. 26. Zhuge, J.; Ierapetritou, M. G., Integration of scheduling and control with closed loop implementation. Industrial & Engineering Chemistry Research 2012, 51, (25), 8550-8565. 27. Chu, Y.; You, F., Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and large-scale global optimization algorithm. Computers & Chemical Engineering 2012, 47, 248-268. 28. Chu, Y.; You, F., Moving horizon approach of integrating scheduling and control for sequential batch processes. AIChE Journal 2014, 60, (5), 1654-1671. 29. Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Ross, R.; van Schijndel, J. M. G., Simultaneous design and control optimisation under uncertainty. Computers & Chemical Engineering 2000, 24, (2-7), 261-266. 30. Yuan, Z.; Chen, B.; Sin, G.; Gani, R., State ‐ of ‐the ‐ art and progress in the optimization ‐ based simultaneous design and control for chemical processes. AIChE Journal 2012, 58, (6), 1640-1659. 31. Bhatia, T.; Biegler, L. T., Dynamic optimization in the design and scheduling of multiproduct batch plants. Industrial & Engineering Chemistry Research 1996, 35, (7), 2234-2246. 32. Birewar, D. B.; Grossmann, I. E., Incorporating scheduling in the optimal design of multiproduct batch plants. Computers & Chemical Engineering 1989, 13, (1–2), 141-161. 33. Engell, S.; Harjunkoski, I., Optimal operation: Scheduling, advanced control and their integration. Computers & Chemical Engineering 2012, 47, 121-133. 34. Harjunkoski, I.; Nystrom, R.; Horch, A., Integration of scheduling and control-Theory or practice? Computers & Chemical Engineering 2009, 33, (12), 1909-1918. 35. Shobrys, D. E.; White, D. C., Planning, scheduling and control systems: why cannot they work together. Computers & Chemical Engineering 2002, 26, (2), 149-160. 36. Jans, R.; Degraeve, Z., Modeling industrial lot sizing problems: a review. International Journal of Production Research 2008, 46, (6), 1619-1643. 37. Karimi, B.; Ghomi, S.; Wilson, J. M., The capacitated lot sizing problem: a review of models and algorithms. Omega-International Journal of Management Science 2003, 31, (5), 365-378. 38. Chachuat, B.; Singer, A. B.; Barton, P. I., Global mixed-integer dynamic optimization. AIChE Journal 2005, 51, (8), 2235-2253. 39. Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S., Dynamic optimization in a discontinuous world. Industrial & Engineering Chemistry Research 1998, 37, (3), 966-981. 40. Biegler, L. T.; Grossmann, I. E., Retrospective on optimization. Computers & Chemical Engineering 2004, 28, (8), 1169-1192. 41. Cuthrell, J. E.; Biegler, L. T., On the optimization of differential‐algebraic process systems. AIChE Journal 1987, 33, (8), 1257-1270. 42. Gorissen, D.; Couckuyt, I.; Demeester, P.; Dhaene, T.; Crombecq, K., A Surrogate Modeling and Adaptive Sampling Toolbox for Computer Based Design. Journal of Machine Learning Research 2010, 11, 20512055. 43. Ong, Y. S.; Nair, P. B.; Keane, A. J., Evolutionary optimization of computationally expensive problems via surrogate modeling. Aiaa Journal 2003, 41, (4), 687-696. 47 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 48

44. Lange, K.; Hunter, D. R.; Yang, I., Optimization transfer using surrogate objective functions. Journal of computational and graphical statistics 2000, 9, (1), 1-20. 45. Demir, Y.; Isleyen, S. K., Evaluation of mathematical models for flexible job-shop scheduling problems. Applied Mathematical Modelling 2013, 37, (3), 977-988. 46. Pinto, J. M.; Grossmann, I. E., A continuous time mixed integer linear programming model for short term scheduling of multistage batch plants. Industrial & Engineering Chemistry Research 1995, 34, (9), 3037-3051. 47. Chu, Y.; You, F., Integrated scheduling and dynamic optimization by stackelberg game: Bilevel model formulation and efficient solution algorithm. Industrial & Engineering Chemistry Research 2014, 53, (13), 55645581. 48. Colson, B.; Marcotte, P.; Savard, G., An overview of bilevel optimization. Annals of Operations Research 2007, 153, (1), 235-256. 49. Chu, Y.; You, F.; Wassick, J. M., Hybrid method integrating agent-based modeling and heuristic tree search for scheduling of complex batch processes. Computers & Chemical Engineering 2014, 60, (0), 277-296. 50. Chu, Y. F.; Wassick, J. M.; You, F. Q., Efficient scheduling method of complex batch processes with general network structure via agent-based modeling. AIChE Journal 2013, 59, (8), 2884-2906. 51. Rosenthal, R. E., GAMS - A User's Guide. GAMS Development Corporation: Washington, DC, USA, 2013. 52. Caccavale, F.; Iamarino, M.; Pierri, F.; Tufano, V., Control and Monitoring of Chemical Batch Reactors. In Control and Monitoring of Chemical Batch Reactors, Springer-Verlag London Ltd: Godalming, 2011; pp 1184. 53. Luyben, W. L., Chemical reactor design and control. Wiley-Interscience: 2007. 54. Myers, R. H.; Anderson-Cook, C. M., Response surface methodology: process and product optimization using designed experiments. Wiley. com: 2009; Vol. 705. 55. Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E., Lagrangean heuristic for the scheduling and control of polymerization reactors. AIChE Journal 2008, 54, (1), 163-182.


Integrated Planning, Scheduling, and Dynamic Optimization for Batch

Recommend Documents