A Three-Stage Solution Algorithm for Chemical Production Scheduling

12th IFAC Symposium on Dynamics and Control of 12th IFAC Symposium on Dynamics and Control of Process including Biosystems 12th IFACSystems, Symposium on Dynamics Dynamics and Control Control of of 12th IFAC Symposium on and Process Systems, including Biosystems Available online at www.sciencedirect.com Florianópolis - SC,including Brazil, April 23-26,and 2019 Process Biosystems 12th IFACSystems, Symposium on Dynamics Control of Process Systems, Biosystems Florianópolis - SC,including Brazil, April 23-26, 2019 Florianópolis SC, Brazil, April 23-26, 2019 Process Systems, including Biosystems Florianópolis - SC, Brazil, April 23-26, 2019 Florianópolis - SC, Brazil, April 23-26, 2019

ScienceDirect

IFAC PapersOnLine 52-1 (2019) 838–843

A A A A

Three-Stage Solution Algorithm for Three-Stage Three-Stage Solution Solution Algorithm Algorithm for for Chemical Production Scheduling Three-Stage Solution Algorithm for Chemical Production Scheduling Chemical Production Scheduling Chemical Production Scheduling Hojae Lee ∗∗ Christos T. Maravelias ∗∗ ∗∗

Hojae Lee ∗∗ Christos T. Maravelias ∗∗ ∗∗ Hojae Hojae Lee Lee ∗ Christos Christos T. T. Maravelias Maravelias ∗∗ ∗ Hojae Lee Christos T. Maravelias ∗ University of Wisconsin-Madison, Madison, WI 53706 USA (Tel: ∗ University of Wisconsin-Madison, Madison, WI 53706 USA (Tel: ∗ University of Wisconsin-Madison, Madison, 608-772-9125; e-mail: Wisconsin-Madison, Madison, WI WI 53706 53706 USA USA (Tel: (Tel: 608-772-9125; e-mail: [email protected]) [email protected]) ∗ University of ∗∗ of Wisconsin-Madison, Madison, WI 53706 USA (Tel: 608-772-9125; e-mail: [email protected]) University of Wisconsin-Madison, Madison, WI 53706 USA (Tel: ∗∗University 608-772-9125; e-mail: [email protected]) University of Wisconsin-Madison, Madison, WI 53706 USA (Tel: ∗∗ 608-772-9125; e-mail: [email protected]) ∗∗ University of Wisconsin-Madison, Wisconsin-Madison, Madison, WI WI 53706 53706 USA USA (Tel: (Tel: 608-265-9026; e-mail: [email protected]) University of Madison, 608-265-9026; e-mail: [email protected]) ∗∗ University of Wisconsin-Madison, Madison, WI 53706 USA (Tel: 608-265-9026; e-mail: 608-265-9026; e-mail: [email protected]) [email protected]) 608-265-9026; e-mail: [email protected]) Abstract: A A novel novel algorithm algorithm that that exploits exploits the the strengths strengths of of discretediscrete- and and continuous-time continuous-time Abstract: Abstract: A novel algorithm that exploits the strengths of discreteand continuous-time scheduling mathematical programming formulations is proposed. It consists of three stages, stages, in in Abstract: A novel algorithm that exploits the strengths of discreteand continuous-time scheduling mathematical programming formulations is proposed. It consists three Abstract: A novel algorithm that exploits the strengths of discreteandof scheduling programming formulations proposed. It consists of three stages, which (i) (i) an anmathematical approximate solution is obtained obtained using a is discrete-time mixed-integer programming scheduling mathematical programming formulations isdiscrete-time proposed. Itmixed-integer consists of continuous-time three stages, in in which approximate solution is using a programming scheduling programming proposed. Itmixed-integer consists three stages, in which (i) approximate solution is using aa isdiscrete-time programming scheduling model, (ii) the the solution is formulations mapped onto continuous-time gridsof via via mapping which (i) an anmathematical approximate solution is obtained obtained usingonto discrete-time mixed-integer programming scheduling model, (ii) solution is mapped continuous-time grids aa mapping which (i) anand approximate solution isthe obtained using a discrete-time mixed-integer scheduling model, (ii) accuracy the solution issolution mapped onto continuous-time grids via viaprogramming a mapping mapping algorithm, (iii) the of is improved by solving a continuous-time linear scheduling model, (ii) the solution is mapped onto continuous-time grids a algorithm, and (iii) the accuracy of theissolution improved by solving a continuous-time linear scheduling model, (ii)Two the types solution mappedis continuous-time gridsnamely via a unitmapping algorithm, and (iii) accuracy of solution is improved by aa continuous-time linear programming model. ofthe continuous-time grids aresolving introduced, and algorithm, and (iii) the the accuracy ofof the solution is onto improved by solving continuous-time linear programming model. Two types continuous-time grids are introduced, namely unitand algorithm, and (iii) the accuracy of the solution is improved by solving a continuous-time programming model. Two types of continuous-time grids are introduced, namely unitand material-specific grids, to ensure feasible unit utilization and material balance in the model. An programming model. Two typesfeasible of continuous-time grids are introduced, namely unit-linear and material-specific grids, to ensure unit utilization and material balance in the model. An programming model. types of continuous-time grids are introduced, namely unitand material-specific grids,Two to study ensure feasible unit utilization utilization andthe material balance in the the model. model. An extensive computational is performed, showing that algorithm is capable of finding material-specific grids, to ensure feasible unit and material balance in An extensive computational isfeasible performed, showing that algorithm is capable of finding material-specific grids,orders to study ensure unitfaster utilization andthe material balance in the model. An extensive computational study is performed, showing that the algorithm is capable capable of finding finding high quality quality solutions of is magnitudes thanthat traditional methods. extensive computational study performed, showing the algorithm is of high solutions orders of magnitudes faster than traditional methods. extensive computational study is performed, showing that the algorithm is capable of finding high quality solutions orders of magnitudes faster than traditional methods. high quality solutions orders of magnitudes faster than traditional methods. © 2019, IFACsolutions (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. high quality orders of magnitudes than traditional methods. Keywords: Chemical production production scheduling, faster solution method, scheduling algorithm, Keywords: Chemical scheduling, solution method, scheduling algorithm, network network Keywords: Chemical production production scheduling, scheduling, solution solution method, method, scheduling scheduling algorithm, environment Keywords: Chemical algorithm, network network environment Keywords: Chemical production scheduling, solution method, scheduling algorithm, network environment environment environment 1. INTRODUCTION INTRODUCTION One of of the the main main strengths strengths of of discrete-time discrete-time models models is is due due 1. One 1. INTRODUCTION One the main strengths of discrete-time models is due to the theof fact that the locations the time points are 1. INTRODUCTION One of fact the that mainthe strengths of of discrete-time models isfixed due to locations of the time points are fixed 1. INTRODUCTION One of fact the that main strengths of of discrete-time due to the locations of the time time points points areisfixed fixed andthe known, which allows direct modeling ofmodels intermediate to the fact that the locations the are and known, which allows direct modeling of intermediate to the fact that the locations of the time points are fixed and known, which allows direct modeling of intermediate Chemical production production scheduling scheduling is is one one of of the the important important events, events, suchwhich as intermediate intermediate shipments and deliveries, and known, allows directshipments modeling of intermediate Chemical such as and deliveries, and known, which allows direct modeling of intermediate Chemical production scheduling is one one of the important events, such as intermediate shipments and deliveries, layers of decision making involved in the operation of as well as linear modeling of inventory/backlog profiles Chemical production scheduling is of the important events, such as intermediate shipments and deliveries, layers of decision making involved in of thethe operation of as well such as linear modeling of shipments inventory/backlog profiles Chemical production scheduling is one important events, as intermediate and deliveries, layers of decision making involved in the operation of as well as linear modeling of inventory/backlog profiles chemical plants. and availability/cost of shared utilities, both constant layers of plants. decision making involved in the operation of and as well as linear modeling of inventory/backlog profiles chemical availability/cost of shared utilities, both constant layers of plants. decision making involved in the operation of and as as linear modeling of inventory/backlog profiles chemical availability/cost of utilities, constant time-varying. Importantly, extensive computational chemical plants. andwell availability/cost of shared shared utilities, both both constant and time-varying. Importantly, extensive computational Scheduling involves the allocation of shared resources to chemical plants. availability/cost of shared utilities, both constant time-varying. Importantly, extensive computational Scheduling involves the allocation of shared resources to and studies showed that thatImportantly, discrete-timeextensive models are are faster than than and time-varying. computational showed discrete-time models faster Scheduling involves the of resources to tasks over over time time (Pinedo, 2012). Plant Plant operation based on on Scheduling involves the allocation allocation of shared shared resources to studies and time-varying. Importantly, extensive computational studies showed that discrete-time models are faster than tasks (Pinedo, 2012). operation based their continuous counterparts when solving large-scale studies showed that discrete-time models are faster than Scheduling involves the allocation of shared resources to continuous counterparts when solving large-scale tasks over 2012). Plant on optimization-based scheduling brings significantbased benefits tasks over time time (Pinedo, (Pinedo, 2012). brings Plant operation operation based on their studies showed that discrete-time models are2011). faster than continuous counterparts solving large-scale optimization-based scheduling significant benefits instances (Sundaramoorthy and when Maravelias, their continuous counterparts when solving large-scale tasks overoftime (Pinedo, 2012). Plant operation based on their instances (Sundaramoorthy and Maravelias, 2011). optimization-based scheduling brings significant benefits in terms increased plant throughput, customer demand optimization-based scheduling brings significant benefits their continuous counterparts when solving large-scale instances (Sundaramoorthy and Maravelias, 2011). in terms of increased plant throughput, customer demand instances (Sundaramoorthy and Maravelias, 2011). optimization-based scheduling brings the significant benefits in terms plantetc. throughput, customer demand satisfaction, cost savings, savings, Among different meth- On On the the other other hand, continuous-time continuous-time models require require signifiin terms of of increased increased plant throughput, customer demand instances (Sundaramoorthy and Maravelias, 2011). signifisatisfaction, cost etc. Among the different methhand, models in terms of increased plant throughput, customer demand satisfaction, cost savings, etc. Among the different methOn the other hand, continuous-time models require signifiods used for scheduling, methods based on mathematicantly fewer time points to represent any solution, leading satisfaction, savings, methods etc. Among the on different meth- cantly On thefewer othertime hand,points continuous-time require leading signifiods used forcost scheduling, based mathematito representmodels any solution, satisfaction, savings, etc. the on different meth- cantly On thefewer other hand,points continuous-time require signifiods used scheduling, methods based mathematitime to any solution, leading cal programming programming have been been theAmong most prominent (Mendez to smaller smaller model sizes. Furthermore, they can represent represent ods used for forcost scheduling, methods based on mathematicantly fewer time points to represent representmodels any solution, leading cal have the most prominent (Mendez to model sizes. Furthermore, they can ods used for scheduling, methods based on mathematicantly fewer time points to represent any solution, leading cal programming have been the most prominent (Mendez to smaller model sizes. Furthermore, they can represent et al., al., 2006; Harjunkoski Harjunkoski et al., al., 2014). thesmaller solutions moresizes. accurately, as the the timing timing of events events are cal programming have been the 2014). most prominent (Mendez the to model Furthermore, they can represent et 2006; et solutions more accurately, as of are cal programming have been the 2014). most prominent (Mendez the to sizes. Furthermore, they can represent et al., 2006; Harjunkoski et al., solutions more accurately, as the of events are partsmaller of the themodel optimization decision andtiming are not not bound by et al., 2006; Harjunkoski et al., 2014). the solutions more accurately, as the timing of events are part of optimization decision and are bound by One of2006; the main main challengesetfor for mathematical programming part et al.,of Harjunkoski al.,mathematical 2014). the solutions more accurately, as the timing of events are of the optimization decision and are not bound by One the challenges programming predetermined location of the time points. part of the optimization are not bound by location of decision the timeand points. One of challenges for programming models formain scheduling is that that they are are difficult difficult to solve solve predetermined One of the the main challenges for mathematical mathematical programming part of the optimization are not bound by predetermined location of decision the time timeand points. models for scheduling is they to predetermined location of the points. One of the main challenges for mathematical programming models for scheduling is that they are difficult to solve due to the combinatorial nature of the problem, especially In this work, we propose an algorithm that harnesses the the models forcombinatorial scheduling is nature that they are difficultespecially to solve In predetermined location ofan thealgorithm time points. due to the of the problem, this work, we propose that harnesses models for scheduling is that they are difficult to solve due to the combinatorial nature of the problem, especially In this work, we propose an algorithm that harnesses the for the industrial-scale instances. Fast solution is essenadvantages of both discreteand continuous-time formuladue the to the combinatorialinstances. nature of the In this work, propose an and algorithm that harnesses the for industrial-scale Fastproblem, solutionespecially is essen- advantages of we both discretecontinuous-time formuladue to the combinatorial nature the especially In thiswhile work, we propose an and algorithm that harnesses the for industrial-scale instances. Fast solution is of both discretecontinuous-time formulatial,the however, since it it allows allows theofoperators operators to frequently frequently tions, overcoming their limitations, namely Discretefor the industrial-scale instances. Fastproblem, solution is essenessen- advantages advantages of both discreteand continuous-time formulatial, however, since the to tions, while overcoming their limitations, namely Discretefor the industrial-scale instances. Fast solution is essenadvantages of both discreteand continuous-time formulatial, however, since it allows the operators to frequently tions, while overcoming their limitations, namely Discretereschedule (e.g.since online scheduling), giving them them the ability Continuous Continuous Algorithm their (DCA). The proposed proposed algorithm tial, however, it scheduling), allows the operators to the frequently tions, while overcoming limitations, namelyalgorithm Discretereschedule (e.g. online giving ability Algorithm (DCA). The tial, however, since it scheduling), allows thedemand operators to the frequently tions, while overcoming their limitations, namelyalgorithm Discretereschedule online giving them ability Algorithm (DCA). The proposed to react react to (e.g. disturbances such as as fluctuations. consists of three three stages: (i) (i) obtain an approximate approximate solution reschedule (e.g. online scheduling), giving fluctuations. them the ability Continuous Continuous Algorithm (DCA). The proposed algorithm to to disturbances such demand consists of stages: obtain an solution reschedule (e.g. online scheduling), giving fluctuations. them the ability consists Continuous Algorithm (DCA). The proposed algorithm to react to disturbances such as demand of three stages: (i) obtain an approximate solution using a discrete-time mixed-integer programming (MIP) to react to disturbances such as demand fluctuations. consists three stages:mixed-integer (i) obtain an approximate a of discrete-time programmingsolution (MIP) To react address the aforementioned aforementioned challenges, the process process using to to disturbances such as demand fluctuations. consists three stages: (i) obtain an approximate solution using aa of discrete-time mixed-integer programming (MIP) To address the challenges, the formulation with a large discretization step length, (ii) using discrete-time mixed-integer programming (MIP) with a large discretization step length, (ii) To address the challenges, the systems engineering (PSE) community community has studied studied differ- formulation To address the aforementioned aforementioned challenges, the process process using a discrete-time mixed-integer programming (MIP) formulation with aa solution large discretization step length, (ii) systems engineering (PSE) has differmap the obtained onto the two novel type of formulation with large discretization step length, (ii) To address the aforementioned challenges, the process the obtained onto the two novel type(ii) of systems engineering (PSE) community has different modeling modeling approaches. How to represent represent and model model systems engineering (PSE) How community has studied studied differ- map formulation withgrids a solution large step length, map the obtained solution the two novel type of ent approaches. to and continuous-time usingdiscretization aonto mapping algorithm, and (iii) map the obtained solution onto the two novel type of systems engineering (PSE) community has studied differcontinuous-time grids using a mapping algorithm, and (iii) ent modeling approaches. How to represent and model timemodeling is one one of ofapproaches. them; whether whether totobase base it on on and precedence ent Howto represent model continuous-time map the using two novel type(iii) of grids using aaonto mapping algorithm, and time is them; it precedence refine the the obtained solution tosolution higher quality continuouscontinuous-time grids mapping algorithm, and (iii) ent modeling approaches. How tobase represent and model refine the solution to aausing higher quality using aa continuoustime is one of them; whether to it on precedence relationships or time-grids, whether to adopt global/local time is one of them; whether to base it on precedence continuous-time grids using a mapping algorithm, and (iii) refine the solution to a higher quality using a continuousrelationships or them; time-grids, whether to adopt time linear linear programming (LP) quality formulation. refine the solution to a higher using a continuoustime is one relations of to base it on global/local precedence programming (LP) formulation. relationships or whether to global/local precedence orwhether common/unit-specific time-grids, time relationships or time-grids, time-grids, whether to adopt adopt global/local refine the solution to a higher quality using a continuoustime linear programming (LP) formulation. precedence relations or common/unit-specific time-grids, time linear programming (LP) formulation. relationships or time-grids, whether to adopt global/local precedence common/unit-specific whether to to relations representor time in aa discrete discrete or or aa time-grids, continuous time linear programming (LP) formulation. precedence relations ortime common/unit-specific time-grids, whether represent in continuous precedence ortime common/unit-specific 2. PROBLEM PROBLEM STATEMENT STATEMENT whether to represent in or continuous fashion, etc. etc. (Maravelias, 2012). In general, general, if time-grids, model is is whether to relations represent time in aa discrete discrete or aaif continuous 2. fashion, (Maravelias, 2012). In aa model 2. PROBLEM STATEMENT whether to represent time in a discrete or a continuous fashion, etc. (Maravelias, 2012). In general, if a model is based on discrete representation of time, it is referred as 2. PROBLEM STATEMENT fashion, (Maravelias, 2012). In general, if referred a model as is based onetc. discrete representation of time, it is PROBLEM STATEMENT fashion, (Maravelias, 2012). general, if continuous a model as is We study the 2.short-term based on discrete representation of is referred a discrete-time discrete-time model, and if it it In is time, basedit on chemical production scheduling scheduling onetc. discrete representation of time, iton is referred as abased model, and if is based continuous We study the short-term chemical production based on discrete representation of time, it is referred as aarepresentation discrete-time model, and if it is based on continuous We study the short-term chemical production scheduling of time, it is referred as a continuous-time problem in network environments. Specifically, given (i) fafadiscrete-time model, and if it is based on continuous We study the short-term chemical production scheduling representation ofmodel, time, itand is referred as a continuous-time inthe network environments. Specifically, (i) arepresentation discrete-time if ithave is based on advantages continuous problem We study short-term chemical productiongiven scheduling representation of time, it is referred as aa continuous-time problem in network environments. Specifically, given (i) famodel. The two types of models different cility data (e.g plant layout, unit connectivity), (ii) producof time, it is referred as continuous-time problem in network environments. Specifically, given (i) famodel. The twooftypes ofitmodels haveas different advantages cility data plantenvironments. layout, unit connectivity), (ii) producrepresentation time, is referred a continuous-time problem in(e.g network Specifically,routes), given (i)(iii) famodel. The of different advantages data (e.g plant layout, (ii) and disadvantages disadvantages (Floudas and have Lin, 2004). 2004). tion recipe recipe (e.g. processing time,connectivity), production model. The two two types types of models models have different advantages cility cility data (e.g plant layout, unit unit connectivity), (ii) producproducand (Floudas and Lin, tion (e.g. processing time, production routes), (iii) model. The two types of models have different advantages cility data (e.g plant layout, unit connectivity), (ii) producand disadvantages (Floudas and Lin, 2004). tion recipe (e.g. processing time, production routes), (iii) and disadvantages (Floudas and Lin, 2004). tion recipe (e.g. processing time, production routes), (iii) and disadvantages (Floudas and Federation Lin, 2004). tion recipe (e.g. Ltd. processing production routes), (iii) 2405-8963 © 2019, IFAC (International of Automatic Control) Hosting by Elsevier All rightstime, reserved.

Copyright © 2019 IFAC 838 Copyright 2019 IFAC 838 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 838 Copyright © 2019 2019 IFAC IFAC 838 10.1016/j.ifacol.2019.06.166 Copyright © 2019 IFAC 838

2019 IFAC DYCOPS Florianópolis - SC, Brazil, April 23-26, 2019

Hojae Lee et al. / IFAC PapersOnLine 52-1 (2019) 838–843

839

sequence between batches, and maps the discrete-time solution onto the two types of continuous-time grids that we newly introduce: material- and unit-specific grids. Here, all the necessary information including assignments of batches to units, sequencing of tasks assigned to a unit and tasks that produce and consume the same material is extracted in the form of algorithmic sets and parameters (see Section 3.3 for details). Finally, in the third stage model, an LP is solved to find a more accurate solution. Specifically, the information from the second stage is utilized to enforce feasible material balance and unit utilization. See Figure 1 for an overview. 3.2 First stage discrete-time MIP scheduling model We use State-Task Network (STN) model proposed by Shah et al. (1993). The model adopts a common time grid that is uniformly discretized into intervals of width, δ. We introduce set N  n to denote the resulting set of time points. We use the over bar notation to denote the converted problem data in terms of the discrete-time grid (e.g. τ¯ij = (τijF + βjmax τijmin )/δ). The model consists of unit utilization (1), unit capacity (2), and material balance constraints (3). n 



i∈Ij n ≥n−¯ τij +1

Xijn ≤ 1 ∀ j, n

D ≤ βjmax Xijn ∀ i, j ∈ Ji , n βjmin Xijn ≤ Bijn D D Sk(n+1) = Skn +

Fig. 1. Overview of Discrete-Continuous Algorithm

+

resource availability (e.g. initial inventory, utility availability), and (iv) demand information (e.g. demand amount and dates), we aim to solve for the optimal schedule in terms of the three objective functions: makespan minimization, cost minimization and profit maximization. The following sets and parameters are used to describe the problem data. We introduce a set of tasks i ∈ I, a set of units j ∈ J, and a set of materials k ∈ K. We also define subsets: Ij denotes tasks that can be performed in − unit j, I+ k /Ik denote tasks that produce/consume material k, and Ji denotes units that can perform task i. The fixed/variable processing cost and time of task i performed F V /αij and τijF /τijV , respectively. The in unit j are given by αij fraction of material k produced (>0) or consumed (<0) by task i is denoted by ρik . The minimum/maximum capacity of unit j is given by βjmin /βjmax . Finally, the demand amount of material k at the end of horizon, H, is given by ξk . 3. DISCRETE-CONTINUOUS ALGORITHM





ρik



j∈Ji

839

(2)

D Bij(n−¯ τij )

D Bijn − ξ¯kn ∀ k, n

(3)

where Xijn ∈ {0, 1} denotes the start of task i in unit j at D time point n, Bijn ∈ R+ denotes the batch size of task i D ∈ R+ processed in unit j starting at time point n, and Skn denotes the inventory of material k during time period n. Parameter ξ¯kn denotes the delivery (<0) or demand (>0) amount of material k at time point n. The objective functions considered are makespan minimization (4), cost minimization (5), and profit maximization (6). min M S : M S ≥ min



i∈I j∈Ji

 

(n + τ¯ij )Xijn ∀ n

F V D (αij Xijn + αij Bijn )

(4) (5)

i∈I j∈Ji n∈N



D πk (Sk( + ¯ H+1)

k∈K

The first stage model aims to obtain an approximate solution using a relatively large discretization step length. The second stage mapping algorithm identifies the relative



j∈Ji

i∈I+ k

i∈I− k

max

3.1 Overview

ρik

(1)

−

 



ξ¯kn )

n∈N F V D (αij Xijn + αij Bijn )

(6)

i∈I j∈Ji n∈N

where M S is a nonnegative continuous variable denoting makespan.

2019 IFAC DYCOPS 840 Florianópolis - SC, Brazil, April 23-26, 2019

Hojae Lee et al. / IFAC PapersOnLine 52-1 (2019) 838–843

Fig. 2. Illustration of mapping batches of tasks to (a) unitspecific grid, and (b) material-specific grid It is important to note that any discrete-time model, such as models based on resource task network representation (Pantelides, 1994), can be used for the first stage model. In addition, any solution methods for discrete-time models, such as reformation, tightening constraints, and decomposition methods can be used (Velez and Maravelias, 2014). 3.3 Second stage mapping algorithm The solution obtained in the first stage is mapped onto two types continuous-time grid, namely unit- and materialspecific grids, using the mapping algorithm presented in this subsection. We introduce set M to denote the set of grid points on the aforementioned grids. Index pair (j, m) represents the grid point m mapped onto the unit j grid. As shown in Figure 2a, the start of the batches is mapped to a grid point. Index pair (k, m) represent the grid point m mapped onto the material k grid. Any batches of tasks that produce or consume material k are mapped on to the k grid. We introduce subsets C MP k /Mk to identify grid points that correspond to batches that produce/consume material k. As shown in Figure 2b, the start/end of the batches that consume/produce material k are mapped onto the grid. The necessary information extracted from the first stage solution (e.g. batch-unit assignment, sequencing of batches of tasks within each material and unit grid) are recorded using algorithmic parameters. Specifically, the information regarding assignments of batches of tasks to unit- and 1 2 , fjkm and material-specific grids are represented by fikm 4 fijkmm . Sequencing relations (e.g. preceding and succeeding batches) between batches that consume and produce 3 material k are represented by fkmm  . An illustrative example is shown in Figure 3. The full mapping algorithm is given in Appendix A.

Fig. 3. Illustration of the mapping algorithm: (a) a simple network, and (b) first stage solution being mapped and the resulting algorithmic sets and parameters It is important to note that the other characteristics of the solution, such as batch sequencing within units and assignments of batches to units remain unchanged, which may limit the improvements in solution accuracy in some cases. We introduce the following continuous nonnegative C to denote the batch size of the task mapped variables: Bkm C to material point (k, m), Skm to denote the inventory level U M of material k at material point (k, m), and Tjm and Tkm to denote the timing of unit point (j, m) and material point (k, m), respectively. We ensure that the finish time of the preceding batches is less or equal to the timing of the succeeding batches. 

M M 3 Tkm ∀k, m, m : fkmm ≤ Tkm   = 1

(7)

We enforce the sequencing between the batches assigned to the same unit.

3.4 Third stage continuous-time LP model Based on the first stage solution and the extracted information, we solve a continuous-time LP model which re-solves for the timing of events, inventory levels and the batch sizes, ultimately obtaining a more accurate solution. 840

U C U Tjm + (τijF + τijV Bkm  ) ≤ Tj(m+1)



4 ∀i, j, k, m, m : fijkmm  = 1

(8)

2019 IFAC DYCOPS Florianópolis - SC, Brazil, April 23-26, 2019

Hojae Lee et al. / IFAC PapersOnLine 52-1 (2019) 838–843

841

Since a single batch is mapped onto multiple points across different unit- and material-specific grids, we need to ensure that the timing of all the relevant grid points are in coordination. 

U M 4 Tjm = Tkm ∀i, j, k, m, m ∈ MC  k : fijkmm = 1

U Tjm

=

M F Tkm  −(τij

+

C τijV Bkm )

(9) (10)



4 ∀i, j, k, m, m ∈ MP k : fijkmm = 1

The batch size has to be bounded by the capacity of the unit it is assigned to. Also, since a single batch can be mapped to multiple material-specific grids, the corresponding batch sizes should be identical to one another.

C Bkm 

C 2 ≤ βjmax ∀j, k, m : fjkm =1 βjmin ≤ Bkm

=



BkC m



(11)



∀i, j, k, k , m, m , m :

4 4 = 1, k < k fijkmm  = f ijk mm



(12)

3.5 Selection of key model parameters

We ensure the inventory balance constraint by enforcing that the inventory level of material k at grid point m is equal to the summation of the initial inventory and the amount produced, minus the amount consumed up to that point. 



C Skm =γk +

+ 1 m>m ∈MP k i∈Ik :f

+

C ρik Bkm  

 ikm



− 1 m≥m ∈MC k i∈Ik :f

 ikm

=1

=1

C ρik Bkm ∀k, m ∈ MC  k

(13) where γk is the initial inventory of material k. For makespan minimization, we make sure that the makespan is greater or equal to each time a material is produced. M ∀k, m ∈ MP minM S : M S ≥ Tkm k

(14)

For cost minimization, we minimize the total processing cost: min

 



i∈I j∈Ji k∈K m∈MP :f 1 =f 2 =1 k ikm jkm

1 F V C + αij Bkm ) (αij |K+ | i

For profit maximization, we maximize the total revenue minus the cost: 

k∈K

−

C πk Sk,M M k

 



i∈I j∈Ji k∈K m∈MP :f 1 k

ikm

2 =fjkm =1

DCA is flexible in that some of the key model parameters can be adjusted according to the needs of the user. Specifically, the discretization step length used in the first stage model can be adjusted if, for example, a more accurate approximate solution is preferred (i.e. shorter step length) or faster solution time is preferred (i.e. longer step length). Furthermore, for cost minimization or profit maximization instances, the horizon length can be relaxed (i.e. made longer), to account for the discretization errors present in the first stage model. Although not discussed in detail, we note that appropriately chosen key model parameters may lead to significantly improved computational speed without compromising the solution quality. 4. COMPUTATIONAL STUDY All the runs in this study are performed using CPLEX 12.6.2 modeled using GAMS 24.9.2 on a cluster with Intel Xeon (E5520) processors at 2.27 GHz and 16 GB of RAM running on CentosOS Linux 7. Resource limit of 1 hr is used, and relative gap of 1% and 3% are allowed for cost minimization and profit maximization instances, respectively. Otherwise, default options are used. 4.1 Illustrative example

(15) | is necessary to avoid counting the The term 1/|K+ i processing cost of a single batch more than once; recall that a single batch can be mapped to more than one grid.

min

Fig. 4. STN representation of the example network

1 F V C + αij Bkm ) (αij |K+ i | (16) 841

The STN representation of a modified network from Maravelias and Papalamprou (2009) is shown in Figure 4, It consists of 37 tasks, 12 units, and 28 states. The processing times of the tasks are irregular (e.g. 3.25, 5.125 h) and are dependent on the batch size. Given no demand, the objective is to maximize the profit throughout 96 h of horizon. The example is solved using DCA. The resulting model has 4583 constraints and 1480 discrete variables. A solution with objective function value of $ 9,205 is obtained in 25 seconds (see Figure 5 for the Gantt chart). Note that when solved with a continuous-time model (Sundaramoorthy and Karimi, 2005), a solution of $ 5,690 is obtained in 1000 seconds with 85% gap.

2019 IFAC DYCOPS 842 Florianópolis - SC, Brazil, April 23-26, 2019

Hojae Lee et al. / IFAC PapersOnLine 52-1 (2019) 838–843

Fig. 5. Gantt chart for the solution obtained by DCA; numbers in the box represent task/batch size 4.2 Computational performance

We compare the performance of DCA against two traditional scheduling models: (i) a discrete-time model (Shah et al., 1993) denoted as DM, and (ii) a continuous-time model (Sundaramoorthy and Karimi, 2005) denoted as CM. A total of 630 runs are included in this study, considering a combination of 7 different networks from the literature, 10 instances for each network, 3 objective functions, and 3 different approaches. The discretization step length of DM is chosen such that the resulting MIP model is tractable (e.g. 1 hr or 2 hr for the instances with 120h horizon length), while the number of time slots for CM is chosen such that additional time slots do not lead to better solutions. In Figure 6a, we show the performance of DCA (black dots), CM (red rectangles), and DM (blue triangles) in terms of solution quality (i.e. q := |BestF oundObj − Obj|/BestF oundObj) and CPU time normalized with respect to the CPU time of the fastest approach. The results show that in most instances, DCA provides the highest quality solution the fastest, showing up to 4 orders of magnitude speedup. In Figure 6b, where the cumulative distribution of instances as a function of q is shown, it can be observed that DCA obtains the best solution in around 85% of the instances, and obtains high quality solution (∼4% inferior than the best found) in the rest of the instances. Figure 6c shows the distribution of instances in which DCA finds solutions with better, same, or worse objective function values compared to the other two methods, shown for each objective function. Numbers in the blocks represent either the average difference in the quality of solution (green/red box), or the amount of speedup observed (gray box). DCA finds the best quality solution in the majority of cost minimization instances, while it finds significantly higher quality solutions (∼6%) in a large portion of makespan minimization and profit maximization instances. In cases of small-scale instances in which traditional methods can find the optimal solutions very fast, DCA finds the optimal solution in most cases, but the advantage in CPU time is not as significant (results not shown). 842

Fig. 6. Computational results: (a) solution quality vs CPU, (b) cumulative distribution based on solution quality, and (c) distribution of better/same/worse quality solutions by objective function

2019 IFAC DYCOPS Florianópolis - SC, Brazil, April 23-26, 2019

Hojae Lee et al. / IFAC PapersOnLine 52-1 (2019) 838–843

843

5. CONCLUSION

Sundaramoorthy, A. and Karimi, I.A. (2005). A simpler better slot-based continuous-time formulation for short-term scheduling in multipurpose batch plants. A novel solution algorithm for chemical production scheduling, which combines the strengths of discrete- and continuous- Chemical Engineering Science, 60(10), 2679–2702. doi: 10.1016/j.ces.2004.12.023. time scheduling formulations, is proposed. It consists of three stages in which (i) an approximate solution is ob- Sundaramoorthy, A. and Maravelias, C.T. (2011). Computational Study of Network-Based Mixed-Integer tained using a discrete-time MIP scheduling model, (ii) Programming Approaches for Chemical Production the solution is mapped onto the unit- and material-specific Scheduling. Industrial & Engineering Chemistry Regrids, and (iii) the accuracy of the solution is improved by search, 50(9), 5023–5040. doi:10.1021/ie101419z. solving the continuous-time LP model. The algorithm is flexible in that some of the key model parameters (e.g. Velez, S. and Maravelias, C.T. (2014). Advances in MixedInteger Programming Methods for Chemical Production discretization step length) can be adjusted to further imScheduling. Annual Review of Chemical and Biomolecprove the performance. Finally, through an extensive comular Engineering, 5, 97–121. doi:10.1146/annurevputational study, it was shown that the proposed approach chembioeng-060713-035859. is capable of obtaining high quality solutions, up to four orders of magnitude faster than traditional approaches. Appendix A. MAPPING ALGORITHM ACKNOWLEDGEMENTS

Table A.1. Algorithmic parameters

H. Lee would like to acknowledge support from University of Wisconsin - Wisconsin Distinguished Graduate Fellowship, as well as the Kwanjeong Educational Foundation, South Korea. REFERENCES Floudas, C.A. and Lin, X.X. (2004). Continuoustime versus discrete-time approaches for scheduling of chemical processes: a review. Computers & Chemical Engineering, 28(11), 2109–2129. doi: 10.1016/j.compchemeng.2004.05.002. Harjunkoski, I., Maravelias, C.T., Bongers, P., Castro, P.M., Engell, S., Grossmann, I.E., Hooker, J., Mendez, C., Sand, G., and Wassick, J. (2014). Scope for industrial applications of production scheduling models and solution methods. Computers & Chemical Engineering, 62, 161–193. doi:10.1016/j.compchemeng.2013.12.001. Maravelias, C.T. and Papalamprou, K. (2009). Polyhedral results for discrete-time production planning mip formulations for continuous processes. Computers & Chemical Engineering, 33(11), 1890–1904. doi:10.1016/j.compchemeng.2009.05.015. URL ://WOS:000271370500010. Maravelias, C.T. (2012). General framework and modeling approach classification for chemical production scheduling. Aiche Journal, 58(6), 1812–1828. doi: 10.1002/aic.13801. Mendez, C.A., Cerda, J., Grossmann, I.E., Harjunkoski, I., and Fahl, M. (2006). State-of-the-art review of optimization methods for short-term scheduling of batch processes. Computers & Chemical Engineering, 30(6-7), 913–946. doi:10.1016/j.compchemeng.2006.02.008. Pantelides, C.C. (1994). Unified frameworks for optimal process planning and scheduling. In Proceedings on the second conference on foundations of computer aided operations, 253–274. Cache Publications New York. Pinedo, M.L. (2012). Scheduling: theory, algorithms, and systems. Springer Science & Business Media. Shah, N., Pantelides, C.C., and Sargent, R.W.H. (1993). A General Algorithm for Short-Term Scheduling of Batch-Operations – II. Computational Issues. Computers & Chemical Engineering, 17(2), 229–244. doi: 10.1016/0098-1354(93)80016-g. 843

MjU MkM

Total number of tasks performed in unit j Total number of tasks that produced/consumed material k 1 One if task i is mapped to material grid point (k, m) fikm 2 One if a task mapped to material grid point (k, m) is fjkm performed in unit j One if point where material k is produced (i.e. m) is f3  kmm  before point where it is consumed (i.e. m ) 4 f  One if task i is mapped to unit grid point (j, m) and ijkmm  material grid point (k, m )

Table A.2. Mapping algorithm 0: 1: 2: 3: 4: 5: 6: 7: 8: 9:

Initialize: MjU = 0 ∀j, MkM = 0 ∀k, f 1∼4 = 0, MP k = MC = ∅ ∀k k Loop n ∈ N Loop k ∈ K Loop i ∈ I+ k Loop j ∈ Ji If Xij(n−¯ τij ) = 1 MkM ←  MkM + 1 U Mj ← Xi jn   n
10:

f2

11:

f4

12: 13: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24:

k

j

)

M) jk,m(=Mk

←1 

ijk,m(=MjU ),m (=M M ) k

←1

EndIf EndLoop x2 Loop i ∈ I− k Loop j ∈ Ji If Xijn = 1 MkM + 1 MkM ←  MjU ← Xi jn   n ≤n i ∈I C M MC k ← Mk ∪ {Mk } f1 M ← 1 ik,m(=Mk )

f2

←1

j

M) jk,m(=Mk f3 ← 1 ∀m < MkM , m ∈  kmm 4 f ←1  ijk,m(=MjU ),m (=M M ) k

25: EndIf 26: EndLoop x4



M MP k , m = Mk