- Email: [email protected]
A heuristic approach for solution of MINLP problem for production scheduling in preemptive mode
Chemical Engineering and Processing 44 (2005) 933–940
A heuristic approach for solution of MINLP problem for production scheduling in preemptive mode Mehmet Yuceer, Ridvan Berber∗ Department of Chemical Engineering, Faculty of Engineering, Ankara University, Tandogan, 06100 Ankara, Turkey Received 8 March 2004; received in revised form 15 October 2004; accepted 3 November 2004 Available online 12 February 2005
Abstract Chemical engineering problems formulated as mixed integer nonlinear programming (MINLP) model are difficult to solve when the model is non-convex. In order to overcome this difficulty, a semi-heuristic algorithm for production scheduling was developed in this work. Using this approach, the non-convex MINLP problem is first considered as an MILP problem without dividing the orders into units. The order causing prolonged delivery time is thus identified, constraints for this order are then relaxed and MILP problem is re-solved using the new constraints. Having reached the new schedule, the quantitative distribution of the specific order to different units is determined by solving the LP problem that does not contain integer variables since allocation of orders to the units, and processing order are known. The results obtained with three example problems indicate improvements over previously reported schedules and therefore, give promise that the suggested strategy may be used in moderately sized industrial applications. © 2004 Elsevier B.V. All rights reserved. Keywords: Production scheduling; MILP; MINLP
1. Introduction The single-stage multi-product plants with non-identical parallel production units are common to processing pharmaceuticals, polymer, food industries and also machinery. Hard industrial competition requires that customer orders not only fulfilled on time, but also manufactured with utmost economic benefit in such plants. This, in turn dictates a strict scheduling of orders in available equipment. If one adds many constraints with regard to the scheduling, such as predecessors/successors, cycles to be avoided and the scheduling horizon; the problem becomes a huge one, right at the ‘formulation’ phase, even before attempting to solve. The solutions provided until the late 1980s had not been considered satisfactory to prevent production tardiness. After, ∗
Corresponding author. Fax: +90 312 2127464. E-mail address: [email protected] (R. Berber).
0255-2701/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cep.2004.11.001
probably, the first attempt towards the solution of this NPcomplete problem by Pekny and Miller [1], some approaches using heuristics [2], MILP formulation in discrete time [3,4] and continuous time with heuristic [5] have been proposed. However, one major difficulty with MILP formulations is that the problem size increases exponentially with the number of units/orders. In a most recent work, Yuceer et al. [6] identified one deficiency and a misconception about the formulation of such problems. They claimed that, contrary to the foregoing statement, the fact that every order can be processed in every unit or after every order in one unit is more contributive to the size of such problems. An overview of developments in the scheduling of multiproduct, multipurpose batch and continuous processes has been recently presented by Floudas and Lin [7]. Besides this review, the two issues can be cited to indicate the direction that the recent research work is moving towards; namely, robustness and cyclic scheduling. Robustness, however currently limited to linear systems, accommodates un-
934
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
certainty in the scheduling algorithm. This has been tackled, for example, by Lin et al. [8], who used a robust optimization methodology previously developed for linear programming problems and extended it to MILP problems. A two-stage stochastic integer program for scheduling under demand uncertainty is also presented by Engell et al. [9] by using scenario decomposition method based on Lagrangian relaxation and some heuristics. Cyclic scheduling concerns with the periodic execution of sub-schedules with a smaller time horizon, particularly for the case when the time horizon is long compared to the duration of individual tasks. Two MINLP formulations and case studies for such problems have been presented by Alle et al. [10] and Wu and Ierapetritou [11]. In our previous work [6], an automatically generated MILP formulation was used to schedule a set of orders in non-identical parallel processing units without dividing any order. However, if the process time of any of the orders in one unit is very long, then it becomes necessary to divide this order into two or more pieces so that other units may also be used to process it. This strategy, which can also be defined as the preemptive mode, will undoubtedly help meet the customer expectations and prevent tardiness. The special difficulty associated with such a strategy, which results in an MINLP formulation, is that the objective function becomes nonlinear, and furthermore non-convex due the multiplicative appearance of integer and real variables. Since available MINLP algorithms are for convex functions only, the problem is seemingly very difficult to solve. This work tackles such problems by developing a semi-heuristic algorithm for the solution, and shows that it is effective for moderately sized cases.
2. Problem statement, model and solution The scheduling problem in single-stage, parallel unit, multi-product plants is considered. Customer orders, each composed of only one product, are to be processed in a subset of available processing units. Total production time is chosen as the objective to be minimized. The model accounts for a number of physical constraints, which are likely to occur in practice, such as follows: • orders may have predecessor and successors; • both the predecessor and the successor to a given order must be manufactured in the same processing unit; • each equipment has only one single order to be first processed in it; • each order can be first processed in only one processing unit; • scheduling horizon (i.e. the maximum period of time available to complete all production requirements) is considered as a constraint to ensure that a due date for the whole set of orders to be processed is met;
Fig. 1. Semi-heuristic algorithm for the solution of non-convex MINLP problem.
• cycle constraints were included such that a job can not be simultaneously both the predecessor and the successor to another one; • release time of units, i.e. the time needed for any unit to become available, is incorporated in the formulation of the objective function. The original problem is first handled as an MILP problem without any attention to dividing the orders. The solution to this problem identifies the orders that extend the processing time, therefore, indicating where the attention should be focused for dividing an order to be processed in multiple units. The orders that can be divided are indicated by the developed software. However, the choice from this list is made by the user. For the type of problems under consideration, dividing one or two orders appears to be satisfactory. In the next step, the constraints for the order to be divided are relaxed so that it is allowed to be processed in more than a unit, and resulting MILP problem is re-solved under these circumstances. Only one round of relaxation is done, and the details are explained in the following Section 2.3. The “relaxed” solution provides the schedule in which the particular order will be processed in multiple units in an optimal setting. As the allocation to the units and the processing order becomes known and the integer variables are eliminated; then the real amount to be processed in each unit remains to be determined, which is done by solving the corresponding linear problem (LP). The solution algorithm developed is presented in Fig. 1 as a flow diagram. 2.1. Objective function The objective of the problem we propose here is the minimization of a function comprising the total processing time of orders being considered, as represented by the following
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
935
mathematical statement: RTEl + J= (XFil + Xmil ) × Til × XMil i∈I
l∈L
+
m ∈ PRil
Xmil × TCLmil .
(1)
m ∈ PRil
2.2. Constraints The MILP formulation incorporates a number of constraints that might be encountered in practice. For some constraints we basically followed the same definition and formulation proposed by Cerda et al. [5], however added a new constraint and mathematically defined the scheduling horizon as a constraint. (a) Predecessor constraint. Each order has at most one unique predecessor being manufactured just before in the same unit. In other words, every order is to be manufactured either as the first order or after another one in only one unit. The order that precedes another one is called the predecessor. The binary variable XFil denotes that the order Oi is the first being processed in unit l. Xmil stands for the schedule of processing Om after order Oi in unit l. The predecessor constraint is represented mathematically by the following equation: Xmil + XFil = 1 i = 1, 2, . . . , n. l ∈ Li m ∈ PRil
l ∈ Li
(2)
(b) Successor constraint. Every order can be either last processed, or manufactured just before another one in the assigned equipment item, and therefore, features at most one single successor. Xikl ≤ 1 i = 1, 2, . . . , n (3) l ∈ Li k ∈ SUil
Predecessor and successor constraints assure the assignment of consecutive orders to the same unit. (c) Both the predecessor and the successor to a given order must be manufactured in the same processing unit. XFil + Xmil + Xiks ≤ 1 m ∈ PRil
s ∈ Li ,s=l k ∈ SUil
l ∈ Li , i = 1, 2, . . . , n
(4)
This inequality constraint dictates that if an order Oi is processed in one unit, there should not be any successor to this order in any other unit. (d) Each equipment has only one single order to be first processed in it. XFil = 1 l∈L (5) i ∈ Il
Fig. 2. Graphical representation of cyclic constraints [5].
(e) Each order can be first processed in only one processing unit. This constraint was added in this work, since the above constraints, which had also been employed [5], could lead to the results that one order might be assigned to more than one unit as a first order to be processed. We encountered such cases when moderately large problems were chosen. XFil ≤ 1 i = 1, 2, . . . , n (6) l ∈ Li
(f) The scheduling horizon is a criterion relating to the timely satisfaction of orders, and presents the maximum time period available to complete the production requirements of the current set of orders. Here it was mathematically defined and incorporated into the formulation as a constraint, in addition to other practical constraints. This constraint ensures that a due date for the whole set of orders under consideration is met. The definition is as follows: M = max RTEl + XFil Til + m ∈ PRil
×(Til + TCLmil )Xmil ,
m ∈ PRil
l = 1, 2, . . . , r)] (7)
To clarify this definition, we can consider an arbitrary unit number 1 and four products (orders) that can be produced in this unit, namely 1, 2, 3 and 5. We assume furthermore that product 1 has one predecessor (product 2) and one successor (product 5). Product 3 is the successor of product 5, and this set of orders is to be completed within 50 h. Under these circumstances, the scheduling horizon becomes: M = RTE1 + (TP11 × NB11 )XF11 + (TP21 × NB21 )XF21 + (TP31 × NB31 )XF31 + (TP51 × NB51 )XF51 + (TP51 × NB51 + TCL151 )X151
936
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
+ (TP21 × NB21 + TCL211 ) × X211 + (TP31
Table 1c Input data, predecessors and order quantities, Example 3
× NB31 + TCL531 )X531 ≤ 50 Depending on the optimum schedule to be found, some of the binary variables will be either 0 or 1, and the total production time in the unit will be less than the scheduling horizon, assuming 50 h, if the constraint is met. Had there been any other unit (say unit 2), we would have written a similar constraint and the maximum of those two numbers, (i.e. the total production times in units 1 and 2) would meet the criteria of being less than 50. (g) Cycle constraint (i). A job cannot be simultaneously both the predecessor and the successor to another one. (Xikl + Xkil ) ≤ 1 k ∈ PRil , l ∈ Lik
i = 1, 2, . . . , n
(8)
(h) Cycle constraint (ii). An order Oi cannot be processed just before the order Ok and immediately after the successor to Ok . Cerda et al. [5] represented this constraint as follows: (Xikl + Xkwl + Xwil ) ≤ 2 k ∈ SUil , l ∈ Likw
w ∈ [SUkl ∩ PRil ],
i = 1, 2, . . . , n
(9)
Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8 O9
O8, O9, O7, O2 O4, O1, O5, O7 O4, O2 O1, O9, O7, O2 O3, O9, O2 O3 O4, O1, O9, O2 O1 O4, O3, O7
770 720 1200 600 990 700 500 940 1150
2.3. Relaxation of the constraints According to the fourth step of the algorithm given in Fig. 1, the MILP problem was re-solved after allowing some relaxations in the constraints. In the new (i.e. relaxed) form of constraint equations, D represents the number of units in which divided orders are to be processed. The value of this number was taken as 2 in the specific examples considered here indicating that the orders were divided in two parts. These constraints, namely Eqs. (2)–(4) and (6), include the following relaxations: (a) Predecessor constraint. Xmil + XFil = D l ∈ Li m ∈ PRil
Likw is the subset of the processing units where orders Oi , Ok and Ow can all be processed. Ow represents the successor of Ok . This constraint prevents the occurrence of cyclic scheduling such as shown in Fig. 2.
l ∈ Li
for every ith order to be divided
(10)
Table 2a Production times of one batch of orders in units, Examples 1 and 2
Table 1a Input data, predecessors and order quantities, Example 1 Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8
O3, O6 O3, O7 O2, O6 O5, O6 O4, O7 O1, O5 O3, O5, O8 O5, O7
700 1050 900 1000 650 1350 950 850
Table 1b Input data, predecessors and order quantities, Example 2 Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
O3, O6, O9, O10 O3, O7 O2, O6, O9, O10 O5, O6, O10 O4, O7, O10 O1, O5, O9 O3, O5, O8 O5, O7 O1, O3, O6, O10 O1, O3, O4, O9
550 850 700 900 500 1050 950 850 450 650
Products
Unit 1
Unit 2
Unit 3
Unit 4
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
1.7 – 1.25 – – 2.4 – – 1.6 2.6
– – – 1.7 1.4 1.8 – – – 1.9
– 0.9 1.1 – – – 1.05 – – –
– – – – 0.85 – 1.65 2.1 – –
Table 2b Production times of one batch of orders in units (Ti day/batch), Example 3 Products
Unit 1
Unit 2
Unit 3
Unit 4
O1 O2 O3 O4 O5 O6 O7 O8 O9
1.25 2.6 – 1.7 – – 1.6 – 2.4
– 1.9 2.1 – 1.7 – – – 1.8
1.1 – 1.9 – – – – 0.9 –
– – 1.8 – – 2.1 – – –
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
(b) Successor constraint. Xikl ≤ D
Table 3a Input data, batch capacities, Examples 1 and 2 Products
Unit 1 (kg/batch)
Unit 2 (kg/batch)
Unit 3 (kg/batch)
Unit 4 (kg/batch)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
100 – 140 – – 280 – – 200 250
– – – 120 90 210 – – – 270
– 210 170 – – – 390 – – –
– – – – 130 – 290 120 – –
937
l ∈ Li k ∈ SUil
for every ith order to be divided
(11)
(c) Both the predecessor and the successor to a given order must be manufactured in the same processing unit. XFil + Xmil + Xiks ≤ D l ∈ Li m ∈ PRil
s ∈ Li ,s=l k ∈ SUil
for every ith order to be divided
(12)
(d) Each order can be first processed in only one processing unit. XFil ≤ D for every ith order to be divided
Table 3b Input data, batch capacities, Example 3
l ∈ Li
Products
Unit 1 (kg/batch)
Unit 2 (kg/batch)
Unit 3 (kg/batch)
Unit 4 (kg/batch)
O1 O2 O3 O4 O5 O6 O7 O8 O9
140 250 – 100 – – 200 – 280
– 270 110 – 120 – – – 210
170 – 115 – – – – 210 –
– – 120 – – 120 – – –
(13)
In real life situations, there is a possibility that some units be not available for job assignment because of either a previously assigned job being processed, or maintenance. Therefore, at any time when a set of orders is to be scheduled, such delays in the availability of some units need to be taken into account. The time needed for any unit to become available is called the release time of equipment (RTE). The model incorporates the release time of units in the formulation of the objective function, and the scheduling horizon. Although
Table 4a Cleaning times, Example 1 and 2 Predecessors (↓)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
Successors (→) O1
O2
O3
O4
O5
O6
O7
O8
O9
O10
– – 1.0 – – 1.4 – – 2.1 1.5
– – 0.15 – – – 1.8 – – –
– 1.1 – – – 0.3 – – 1.25 0.6
– – – – 0.3 0.7 – – – 0.75
– – – 0.05 – – 0.85 – – 0.5
0.65 – – – 0.7 – – – 0.8 –
– – 0.3 – 0.9 – – 1.65 – –
– – – – 0.6 – 0.45 – – –
0.85 – 1.6 – – 1.2 – – – 0.7
0.4 – 0.2 0.5 – – – – 0.65 –
Table 4b Cleaning times, Example 3 Predecessors (↓)
O1 O2 O3 O4 O5 O6 O7 O8 O9
Successors (→) O1
O2
O3
O4
O5
O6
O7
O8
O9
– 0.6 – – – – 1.25 1.1 0.3
0.2 – – 0.4 0.5 – 0.65 – –
0.4 0.5 – – 0.05 – – 0.9 –
1.0 1.5 – – – – 2.1 – 1.4
– 0.75 0.3 – – – – – 0.7
– – 0.6 – – – – – –
1.6 0.7 – 0.85 – – – – 1.2
0.15 – 0.8 – – – – – –
– – 0.7 0.65 – – 0.8 – –
938
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
for the current MILP models heuristic rules are generally resorted to for feasible solutions, the specific MILP model and solution proposed here avoid any such need. It is to be mentioned that only physical type of operations for production are considered here, and thus processing time depends on the batch size. When chemical reactions are involved, this assumption will not be valid and the case will get more complicated. In such instances, the final phase of the algorithm, where the quantitative sizes of divisions are determined, will no longer be a linear problem. A MATLAB program was developed to formulate the problem automatically so that it can be solved by a branch and bound technique in MILP form. The program takes the data related to the properties of the orders and production units, creates the model (i.e. objective function and all the constraint equations), produces the input files for MINOS package [12] to be run in Fortran environment for branch and bound iterations, and necessary other files for easy interpretation of the MINOS output. The user, then, easily extracts the optimum results, i.e. allocation of orders into the production units and the schedule in each unit, and draws the Gantt chart.
optimum scheduling out of this solution resulted in an objective function value of 77.15, 81.4 and 90.2 days for Examples 1, 2 and 3, respectively. The completion time for the set of orders was 24.55, 26.25 and 33 days. In all cases, the scheduling horizon was considered as a constraint. However, the results, as shown in Figs. 3a and 4a, indicated that dividing the order 6, coincidentally in the first two examples, would have resulted in even better optimal solution in terms of fulfilling the possible customer demand of earlier completion of the set of orders in Examples 1 and 2. Fig. 5a indicate that both dividing the two orders 3 and 9 would be benefical in terms of earlier completion of the set in Example 3. When the solution strategy proposed in this work was implemented for these examples, the optimum schedules shown in Gantt diagrams of Figs. 3b, 4b and 5b were found. The objective functions reached in this case were 78.45, 81.5 and 96.7 days for Examples 1, 2 and 3, and
3. Examples and results A set of example problems was solved in order to test the proposed method. Two of the examples presented here are some of the specific problems tackled by Cerda et al. [5]. These examples were previously reconsidered and solved by Yuceer et al. [6] with an objective function of minimizing the total production time. In addition to these two examples, a third one was added in order to demonstrate the possibility of dividing two orders by the suggested algorithm. Example 1 has 8 orders to be processed in 4 units whereas Example 2 considers a set of 10 orders for 4 units as indicated in Tables 1a and 1b, respectively. Example 3 has 4 units and a customer demand of 9 products as shown in Table 1c. All units are assumed to be available at the time of scheduling. Tables 2a and 2b show the production times of one batch for each product. Production capacities for orders in every unit are given in Tables 3a and 3b, while the cleaning times are indicated in Tables 4a and 4b. Data in Tables 2a, 3a and 4a hold for the first two examples, with the difference that in Example 1 only the first 8 of the orders are taken into account. There are limitations that the set of orders should be completed within 25 days in Example 1, and within 27 days in Example 2. Example 1 has 29 binary variables and 314 constraint equations, whereas Example 2 has 47 binary variables and 648 constraint equations. Example 3 involves 43 binary variables and 607 constraint equations. The objective for all is to minimize the total production time. These problems were first solved with an MILP formulation, i.e. without resorting to any division for orders. The
Fig. 3. (a) Solution of the Example 1 without dividing the orders [6]; (b) solution of the Example 1 when order 6 is divided between units (this work).
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
939
Fig. 4. (a) Solution of the Example 2 without dividing the orders [6]; (b) solution of the Example 2 when order 6 is divided between units (this work).
Fig. 5. (a) Solution of Example 3 without dividing the orders (this work); (b) solution of Example 3 when order 3 and 9 are divided between units (this work).
the orders were to be completed within 21.65, 23.85 and 26.25 days, respectively. These results are clearly indicative of the fact that, with relatively small sacrifices of 1.3 and 0.1 days in the objective functions, the set of orders would be completed 2.9 and 2.4 days earlier for example problems 1 and 2. For Example 3, the objective function sacrifices 6.5 days when the suggested algorithm is implemented for the division of two orders, but the set of orders are completed 6.75 days earlier. The CPU time for runs on a 2.4 GHz Pentium-4 processor, 256 MB memory and 40 GB hard disk capacity PC is shown in Table 5 for all of the example problems. The time indicated involves the total time spent for computing in three phase of the algorithm, i.e. solutions for MILP, relaxed MILP and LP problems.
Table 5 CPU time for example runs Example Number of number units orders
Number of Number of binary variables constraint equations
CPU time (s)
1 2 3
29 47 43
22.18 56.84 55.58
4–8 4–10 4–9
314 648 607
4. Conclusions This work considered optimal scheduling of a set of orders in multi-product batch plants with non-identical, single-stage, parallel processing units in preemptive operation mode. Particular emphasis was placed in dividing one or two orders into multiple units, and a semi-
940
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
heuristic algorithm was proposed to solve the MINLP problem with a non-convex objective function. The results obtained from some example problems indicate improvements over previous schedules where such emphasis was not considered, and therefore, show promise that the suggested strategy can be used for implementation in industrial environment.
Appendix A. Notation
D I Il J L Li Lik M n NBil PRi PRil
r RTEl SUi SUil
TCLmil Til TPil XFil
number of units in which divided orders are to be processed set of available orders set of orders that can be processed in unit l objective function set of available units set of available units for Oi set of available units for the pair of orders (Oi , Ok ) scheduling horizon (the finishing time for the whole set of orders) total number of orders (each order corresponding to a specific product) number of batches of Oi to be processed when the unit l is assigned to Oi whole set of predecessors of order Oi set of customer orders that can be processed just before Oi in unit l, i.e. feasible predecessors of Oi in unit l total number of processing units release time of unit l set of feasible successors to order O set of customer orders that can be processed immediately after Oi in unit l, i.e. feasible successors to Oi in unit l cleaning time for the ordered pair of jobs (Oi , Om ) in unit l total processing time of order Oi in unit l. processing time of a batch of Oi in unit l binary variable denoting that order Oi is the first being processed in unit l
Xmil
XMil
binary variable denoting that the processing of order Om takes place in unit l just before the campaign for Oi amount of the Oi order to be processed in unit l
References [1] J.F. Pekny, D.L. Miller, Exact solution of the no-wait flowshop scheduling problem with a comparison to heuristic methods, Comput. Chem. Eng. 15 (1991) 741–748. [2] R.F.H. Mussier, L.B. Evans, An approximate method for the production scheduling of industrial batch processes with parallel units, Comput. Chem. Eng. 13 (1989) 229–238. [3] E. Kondili, C.C. Pantelides, R.W.H. Sargent, A general algorithm for short-term scheduling of batch operations. I: MILP formulation, Comput. Chem. Eng. 17 (1993) 211–227. [4] N. Shah, C.C. Pantelides, R.W.H. Sargent, A general algorithm for short-term scheduling of batch operations. II: computational issues, Comput. Chem. Eng. 17 (1993) 229–244. [5] J. Cerda, G.P. Henning, I.E. Grossmann, A mixed integer linear programming model for short-term scheduling of single-stage multiproduct batch plants with parallel lines, Ind. Eng. Chem. Res. 36 (1997) 1695–1707. [6] M. Yuceer, R. Berber, Z. Ozdemir, Automatic generation of production scheduling problems in single-stage multi-product batch plants: some numerical examples, in: (Proceedings of International Conference of Computational Methods in Sciences and Engineering—ICCMSE 2004, Athens, Greece, 19–23 November 2004), Lecture Series on Computer and Computational Sciences, vol. 1, VSP International Science Publishers, 2004, pp. 657–660. [7] C.A. Floudas, X. Lin, Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review, Comput. Chem. Eng. 28 (2004) 2109–2129. [8] X. Lin, S.L. Janak, C.A. Floudas, A new robust optimization approach for scheduling under uncertainty: I. Bounded uncertainty, Comput. Chem. Eng. 28 (2004) 1069–1085. [9] S. Engell, A. Markert, G. Sand, R. Schultz, Aggregated scheduling of a multiproduct batch plant by two-stage stochastic integer programming, Optim. Eng. 5 (2004) 335–359. [10] A. Alle, L.G. Papageorgiou, J.M. Pinto, A mathematical programming approach for cyclic production and cleaning scheduling of multistage continuous plants, Comput. Chem. Eng. 28 (2004) 3–15. [11] W. Wu, M. Ierapetritou, Cyclic short term scheduling of multiproduct batch plants using continuous-time representation, Comput. Chem. Eng. 28 (2004) 3–15. [12] B.A. Murtagh, M.A. Saunders, MINOS version 5.4 User’s Guide, Systems Optimization Laboratory, Stanford University, Stanford, CA, 1975.
A heuristic approach for solution of MINLP problem for production scheduling in preemptive mode Mehmet Yuceer, Ridvan Berber∗ Department of Chemical Engineering, Faculty of Engineering, Ankara University, Tandogan, 06100 Ankara, Turkey Received 8 March 2004; received in revised form 15 October 2004; accepted 3 November 2004 Available online 12 February 2005
Abstract Chemical engineering problems formulated as mixed integer nonlinear programming (MINLP) model are difficult to solve when the model is non-convex. In order to overcome this difficulty, a semi-heuristic algorithm for production scheduling was developed in this work. Using this approach, the non-convex MINLP problem is first considered as an MILP problem without dividing the orders into units. The order causing prolonged delivery time is thus identified, constraints for this order are then relaxed and MILP problem is re-solved using the new constraints. Having reached the new schedule, the quantitative distribution of the specific order to different units is determined by solving the LP problem that does not contain integer variables since allocation of orders to the units, and processing order are known. The results obtained with three example problems indicate improvements over previously reported schedules and therefore, give promise that the suggested strategy may be used in moderately sized industrial applications. © 2004 Elsevier B.V. All rights reserved. Keywords: Production scheduling; MILP; MINLP
1. Introduction The single-stage multi-product plants with non-identical parallel production units are common to processing pharmaceuticals, polymer, food industries and also machinery. Hard industrial competition requires that customer orders not only fulfilled on time, but also manufactured with utmost economic benefit in such plants. This, in turn dictates a strict scheduling of orders in available equipment. If one adds many constraints with regard to the scheduling, such as predecessors/successors, cycles to be avoided and the scheduling horizon; the problem becomes a huge one, right at the ‘formulation’ phase, even before attempting to solve. The solutions provided until the late 1980s had not been considered satisfactory to prevent production tardiness. After, ∗
Corresponding author. Fax: +90 312 2127464. E-mail address: [email protected] (R. Berber).
0255-2701/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cep.2004.11.001
probably, the first attempt towards the solution of this NPcomplete problem by Pekny and Miller [1], some approaches using heuristics [2], MILP formulation in discrete time [3,4] and continuous time with heuristic [5] have been proposed. However, one major difficulty with MILP formulations is that the problem size increases exponentially with the number of units/orders. In a most recent work, Yuceer et al. [6] identified one deficiency and a misconception about the formulation of such problems. They claimed that, contrary to the foregoing statement, the fact that every order can be processed in every unit or after every order in one unit is more contributive to the size of such problems. An overview of developments in the scheduling of multiproduct, multipurpose batch and continuous processes has been recently presented by Floudas and Lin [7]. Besides this review, the two issues can be cited to indicate the direction that the recent research work is moving towards; namely, robustness and cyclic scheduling. Robustness, however currently limited to linear systems, accommodates un-
934
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
certainty in the scheduling algorithm. This has been tackled, for example, by Lin et al. [8], who used a robust optimization methodology previously developed for linear programming problems and extended it to MILP problems. A two-stage stochastic integer program for scheduling under demand uncertainty is also presented by Engell et al. [9] by using scenario decomposition method based on Lagrangian relaxation and some heuristics. Cyclic scheduling concerns with the periodic execution of sub-schedules with a smaller time horizon, particularly for the case when the time horizon is long compared to the duration of individual tasks. Two MINLP formulations and case studies for such problems have been presented by Alle et al. [10] and Wu and Ierapetritou [11]. In our previous work [6], an automatically generated MILP formulation was used to schedule a set of orders in non-identical parallel processing units without dividing any order. However, if the process time of any of the orders in one unit is very long, then it becomes necessary to divide this order into two or more pieces so that other units may also be used to process it. This strategy, which can also be defined as the preemptive mode, will undoubtedly help meet the customer expectations and prevent tardiness. The special difficulty associated with such a strategy, which results in an MINLP formulation, is that the objective function becomes nonlinear, and furthermore non-convex due the multiplicative appearance of integer and real variables. Since available MINLP algorithms are for convex functions only, the problem is seemingly very difficult to solve. This work tackles such problems by developing a semi-heuristic algorithm for the solution, and shows that it is effective for moderately sized cases.
2. Problem statement, model and solution The scheduling problem in single-stage, parallel unit, multi-product plants is considered. Customer orders, each composed of only one product, are to be processed in a subset of available processing units. Total production time is chosen as the objective to be minimized. The model accounts for a number of physical constraints, which are likely to occur in practice, such as follows: • orders may have predecessor and successors; • both the predecessor and the successor to a given order must be manufactured in the same processing unit; • each equipment has only one single order to be first processed in it; • each order can be first processed in only one processing unit; • scheduling horizon (i.e. the maximum period of time available to complete all production requirements) is considered as a constraint to ensure that a due date for the whole set of orders to be processed is met;
Fig. 1. Semi-heuristic algorithm for the solution of non-convex MINLP problem.
• cycle constraints were included such that a job can not be simultaneously both the predecessor and the successor to another one; • release time of units, i.e. the time needed for any unit to become available, is incorporated in the formulation of the objective function. The original problem is first handled as an MILP problem without any attention to dividing the orders. The solution to this problem identifies the orders that extend the processing time, therefore, indicating where the attention should be focused for dividing an order to be processed in multiple units. The orders that can be divided are indicated by the developed software. However, the choice from this list is made by the user. For the type of problems under consideration, dividing one or two orders appears to be satisfactory. In the next step, the constraints for the order to be divided are relaxed so that it is allowed to be processed in more than a unit, and resulting MILP problem is re-solved under these circumstances. Only one round of relaxation is done, and the details are explained in the following Section 2.3. The “relaxed” solution provides the schedule in which the particular order will be processed in multiple units in an optimal setting. As the allocation to the units and the processing order becomes known and the integer variables are eliminated; then the real amount to be processed in each unit remains to be determined, which is done by solving the corresponding linear problem (LP). The solution algorithm developed is presented in Fig. 1 as a flow diagram. 2.1. Objective function The objective of the problem we propose here is the minimization of a function comprising the total processing time of orders being considered, as represented by the following
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
935
mathematical statement: RTEl + J= (XFil + Xmil ) × Til × XMil i∈I
l∈L
+
m ∈ PRil
Xmil × TCLmil .
(1)
m ∈ PRil
2.2. Constraints The MILP formulation incorporates a number of constraints that might be encountered in practice. For some constraints we basically followed the same definition and formulation proposed by Cerda et al. [5], however added a new constraint and mathematically defined the scheduling horizon as a constraint. (a) Predecessor constraint. Each order has at most one unique predecessor being manufactured just before in the same unit. In other words, every order is to be manufactured either as the first order or after another one in only one unit. The order that precedes another one is called the predecessor. The binary variable XFil denotes that the order Oi is the first being processed in unit l. Xmil stands for the schedule of processing Om after order Oi in unit l. The predecessor constraint is represented mathematically by the following equation: Xmil + XFil = 1 i = 1, 2, . . . , n. l ∈ Li m ∈ PRil
l ∈ Li
(2)
(b) Successor constraint. Every order can be either last processed, or manufactured just before another one in the assigned equipment item, and therefore, features at most one single successor. Xikl ≤ 1 i = 1, 2, . . . , n (3) l ∈ Li k ∈ SUil
Predecessor and successor constraints assure the assignment of consecutive orders to the same unit. (c) Both the predecessor and the successor to a given order must be manufactured in the same processing unit. XFil + Xmil + Xiks ≤ 1 m ∈ PRil
s ∈ Li ,s=l k ∈ SUil
l ∈ Li , i = 1, 2, . . . , n
(4)
This inequality constraint dictates that if an order Oi is processed in one unit, there should not be any successor to this order in any other unit. (d) Each equipment has only one single order to be first processed in it. XFil = 1 l∈L (5) i ∈ Il
Fig. 2. Graphical representation of cyclic constraints [5].
(e) Each order can be first processed in only one processing unit. This constraint was added in this work, since the above constraints, which had also been employed [5], could lead to the results that one order might be assigned to more than one unit as a first order to be processed. We encountered such cases when moderately large problems were chosen. XFil ≤ 1 i = 1, 2, . . . , n (6) l ∈ Li
(f) The scheduling horizon is a criterion relating to the timely satisfaction of orders, and presents the maximum time period available to complete the production requirements of the current set of orders. Here it was mathematically defined and incorporated into the formulation as a constraint, in addition to other practical constraints. This constraint ensures that a due date for the whole set of orders under consideration is met. The definition is as follows: M = max RTEl + XFil Til + m ∈ PRil
×(Til + TCLmil )Xmil ,
m ∈ PRil
l = 1, 2, . . . , r)] (7)
To clarify this definition, we can consider an arbitrary unit number 1 and four products (orders) that can be produced in this unit, namely 1, 2, 3 and 5. We assume furthermore that product 1 has one predecessor (product 2) and one successor (product 5). Product 3 is the successor of product 5, and this set of orders is to be completed within 50 h. Under these circumstances, the scheduling horizon becomes: M = RTE1 + (TP11 × NB11 )XF11 + (TP21 × NB21 )XF21 + (TP31 × NB31 )XF31 + (TP51 × NB51 )XF51 + (TP51 × NB51 + TCL151 )X151
936
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
+ (TP21 × NB21 + TCL211 ) × X211 + (TP31
Table 1c Input data, predecessors and order quantities, Example 3
× NB31 + TCL531 )X531 ≤ 50 Depending on the optimum schedule to be found, some of the binary variables will be either 0 or 1, and the total production time in the unit will be less than the scheduling horizon, assuming 50 h, if the constraint is met. Had there been any other unit (say unit 2), we would have written a similar constraint and the maximum of those two numbers, (i.e. the total production times in units 1 and 2) would meet the criteria of being less than 50. (g) Cycle constraint (i). A job cannot be simultaneously both the predecessor and the successor to another one. (Xikl + Xkil ) ≤ 1 k ∈ PRil , l ∈ Lik
i = 1, 2, . . . , n
(8)
(h) Cycle constraint (ii). An order Oi cannot be processed just before the order Ok and immediately after the successor to Ok . Cerda et al. [5] represented this constraint as follows: (Xikl + Xkwl + Xwil ) ≤ 2 k ∈ SUil , l ∈ Likw
w ∈ [SUkl ∩ PRil ],
i = 1, 2, . . . , n
(9)
Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8 O9
O8, O9, O7, O2 O4, O1, O5, O7 O4, O2 O1, O9, O7, O2 O3, O9, O2 O3 O4, O1, O9, O2 O1 O4, O3, O7
770 720 1200 600 990 700 500 940 1150
2.3. Relaxation of the constraints According to the fourth step of the algorithm given in Fig. 1, the MILP problem was re-solved after allowing some relaxations in the constraints. In the new (i.e. relaxed) form of constraint equations, D represents the number of units in which divided orders are to be processed. The value of this number was taken as 2 in the specific examples considered here indicating that the orders were divided in two parts. These constraints, namely Eqs. (2)–(4) and (6), include the following relaxations: (a) Predecessor constraint. Xmil + XFil = D l ∈ Li m ∈ PRil
Likw is the subset of the processing units where orders Oi , Ok and Ow can all be processed. Ow represents the successor of Ok . This constraint prevents the occurrence of cyclic scheduling such as shown in Fig. 2.
l ∈ Li
for every ith order to be divided
(10)
Table 2a Production times of one batch of orders in units, Examples 1 and 2
Table 1a Input data, predecessors and order quantities, Example 1 Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8
O3, O6 O3, O7 O2, O6 O5, O6 O4, O7 O1, O5 O3, O5, O8 O5, O7
700 1050 900 1000 650 1350 950 850
Table 1b Input data, predecessors and order quantities, Example 2 Products
Predecessors
Quantity (kg)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
O3, O6, O9, O10 O3, O7 O2, O6, O9, O10 O5, O6, O10 O4, O7, O10 O1, O5, O9 O3, O5, O8 O5, O7 O1, O3, O6, O10 O1, O3, O4, O9
550 850 700 900 500 1050 950 850 450 650
Products
Unit 1
Unit 2
Unit 3
Unit 4
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
1.7 – 1.25 – – 2.4 – – 1.6 2.6
– – – 1.7 1.4 1.8 – – – 1.9
– 0.9 1.1 – – – 1.05 – – –
– – – – 0.85 – 1.65 2.1 – –
Table 2b Production times of one batch of orders in units (Ti day/batch), Example 3 Products
Unit 1
Unit 2
Unit 3
Unit 4
O1 O2 O3 O4 O5 O6 O7 O8 O9
1.25 2.6 – 1.7 – – 1.6 – 2.4
– 1.9 2.1 – 1.7 – – – 1.8
1.1 – 1.9 – – – – 0.9 –
– – 1.8 – – 2.1 – – –
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
(b) Successor constraint. Xikl ≤ D
Table 3a Input data, batch capacities, Examples 1 and 2 Products
Unit 1 (kg/batch)
Unit 2 (kg/batch)
Unit 3 (kg/batch)
Unit 4 (kg/batch)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
100 – 140 – – 280 – – 200 250
– – – 120 90 210 – – – 270
– 210 170 – – – 390 – – –
– – – – 130 – 290 120 – –
937
l ∈ Li k ∈ SUil
for every ith order to be divided
(11)
(c) Both the predecessor and the successor to a given order must be manufactured in the same processing unit. XFil + Xmil + Xiks ≤ D l ∈ Li m ∈ PRil
s ∈ Li ,s=l k ∈ SUil
for every ith order to be divided
(12)
(d) Each order can be first processed in only one processing unit. XFil ≤ D for every ith order to be divided
Table 3b Input data, batch capacities, Example 3
l ∈ Li
Products
Unit 1 (kg/batch)
Unit 2 (kg/batch)
Unit 3 (kg/batch)
Unit 4 (kg/batch)
O1 O2 O3 O4 O5 O6 O7 O8 O9
140 250 – 100 – – 200 – 280
– 270 110 – 120 – – – 210
170 – 115 – – – – 210 –
– – 120 – – 120 – – –
(13)
In real life situations, there is a possibility that some units be not available for job assignment because of either a previously assigned job being processed, or maintenance. Therefore, at any time when a set of orders is to be scheduled, such delays in the availability of some units need to be taken into account. The time needed for any unit to become available is called the release time of equipment (RTE). The model incorporates the release time of units in the formulation of the objective function, and the scheduling horizon. Although
Table 4a Cleaning times, Example 1 and 2 Predecessors (↓)
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
Successors (→) O1
O2
O3
O4
O5
O6
O7
O8
O9
O10
– – 1.0 – – 1.4 – – 2.1 1.5
– – 0.15 – – – 1.8 – – –
– 1.1 – – – 0.3 – – 1.25 0.6
– – – – 0.3 0.7 – – – 0.75
– – – 0.05 – – 0.85 – – 0.5
0.65 – – – 0.7 – – – 0.8 –
– – 0.3 – 0.9 – – 1.65 – –
– – – – 0.6 – 0.45 – – –
0.85 – 1.6 – – 1.2 – – – 0.7
0.4 – 0.2 0.5 – – – – 0.65 –
Table 4b Cleaning times, Example 3 Predecessors (↓)
O1 O2 O3 O4 O5 O6 O7 O8 O9
Successors (→) O1
O2
O3
O4
O5
O6
O7
O8
O9
– 0.6 – – – – 1.25 1.1 0.3
0.2 – – 0.4 0.5 – 0.65 – –
0.4 0.5 – – 0.05 – – 0.9 –
1.0 1.5 – – – – 2.1 – 1.4
– 0.75 0.3 – – – – – 0.7
– – 0.6 – – – – – –
1.6 0.7 – 0.85 – – – – 1.2
0.15 – 0.8 – – – – – –
– – 0.7 0.65 – – 0.8 – –
938
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
for the current MILP models heuristic rules are generally resorted to for feasible solutions, the specific MILP model and solution proposed here avoid any such need. It is to be mentioned that only physical type of operations for production are considered here, and thus processing time depends on the batch size. When chemical reactions are involved, this assumption will not be valid and the case will get more complicated. In such instances, the final phase of the algorithm, where the quantitative sizes of divisions are determined, will no longer be a linear problem. A MATLAB program was developed to formulate the problem automatically so that it can be solved by a branch and bound technique in MILP form. The program takes the data related to the properties of the orders and production units, creates the model (i.e. objective function and all the constraint equations), produces the input files for MINOS package [12] to be run in Fortran environment for branch and bound iterations, and necessary other files for easy interpretation of the MINOS output. The user, then, easily extracts the optimum results, i.e. allocation of orders into the production units and the schedule in each unit, and draws the Gantt chart.
optimum scheduling out of this solution resulted in an objective function value of 77.15, 81.4 and 90.2 days for Examples 1, 2 and 3, respectively. The completion time for the set of orders was 24.55, 26.25 and 33 days. In all cases, the scheduling horizon was considered as a constraint. However, the results, as shown in Figs. 3a and 4a, indicated that dividing the order 6, coincidentally in the first two examples, would have resulted in even better optimal solution in terms of fulfilling the possible customer demand of earlier completion of the set of orders in Examples 1 and 2. Fig. 5a indicate that both dividing the two orders 3 and 9 would be benefical in terms of earlier completion of the set in Example 3. When the solution strategy proposed in this work was implemented for these examples, the optimum schedules shown in Gantt diagrams of Figs. 3b, 4b and 5b were found. The objective functions reached in this case were 78.45, 81.5 and 96.7 days for Examples 1, 2 and 3, and
3. Examples and results A set of example problems was solved in order to test the proposed method. Two of the examples presented here are some of the specific problems tackled by Cerda et al. [5]. These examples were previously reconsidered and solved by Yuceer et al. [6] with an objective function of minimizing the total production time. In addition to these two examples, a third one was added in order to demonstrate the possibility of dividing two orders by the suggested algorithm. Example 1 has 8 orders to be processed in 4 units whereas Example 2 considers a set of 10 orders for 4 units as indicated in Tables 1a and 1b, respectively. Example 3 has 4 units and a customer demand of 9 products as shown in Table 1c. All units are assumed to be available at the time of scheduling. Tables 2a and 2b show the production times of one batch for each product. Production capacities for orders in every unit are given in Tables 3a and 3b, while the cleaning times are indicated in Tables 4a and 4b. Data in Tables 2a, 3a and 4a hold for the first two examples, with the difference that in Example 1 only the first 8 of the orders are taken into account. There are limitations that the set of orders should be completed within 25 days in Example 1, and within 27 days in Example 2. Example 1 has 29 binary variables and 314 constraint equations, whereas Example 2 has 47 binary variables and 648 constraint equations. Example 3 involves 43 binary variables and 607 constraint equations. The objective for all is to minimize the total production time. These problems were first solved with an MILP formulation, i.e. without resorting to any division for orders. The
Fig. 3. (a) Solution of the Example 1 without dividing the orders [6]; (b) solution of the Example 1 when order 6 is divided between units (this work).
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
939
Fig. 4. (a) Solution of the Example 2 without dividing the orders [6]; (b) solution of the Example 2 when order 6 is divided between units (this work).
Fig. 5. (a) Solution of Example 3 without dividing the orders (this work); (b) solution of Example 3 when order 3 and 9 are divided between units (this work).
the orders were to be completed within 21.65, 23.85 and 26.25 days, respectively. These results are clearly indicative of the fact that, with relatively small sacrifices of 1.3 and 0.1 days in the objective functions, the set of orders would be completed 2.9 and 2.4 days earlier for example problems 1 and 2. For Example 3, the objective function sacrifices 6.5 days when the suggested algorithm is implemented for the division of two orders, but the set of orders are completed 6.75 days earlier. The CPU time for runs on a 2.4 GHz Pentium-4 processor, 256 MB memory and 40 GB hard disk capacity PC is shown in Table 5 for all of the example problems. The time indicated involves the total time spent for computing in three phase of the algorithm, i.e. solutions for MILP, relaxed MILP and LP problems.
Table 5 CPU time for example runs Example Number of number units orders
Number of Number of binary variables constraint equations
CPU time (s)
1 2 3
29 47 43
22.18 56.84 55.58
4–8 4–10 4–9
314 648 607
4. Conclusions This work considered optimal scheduling of a set of orders in multi-product batch plants with non-identical, single-stage, parallel processing units in preemptive operation mode. Particular emphasis was placed in dividing one or two orders into multiple units, and a semi-
940
M. Yuceer, R. Berber / Chemical Engineering and Processing 44 (2005) 933–940
heuristic algorithm was proposed to solve the MINLP problem with a non-convex objective function. The results obtained from some example problems indicate improvements over previous schedules where such emphasis was not considered, and therefore, show promise that the suggested strategy can be used for implementation in industrial environment.
Appendix A. Notation
D I Il J L Li Lik M n NBil PRi PRil
r RTEl SUi SUil
TCLmil Til TPil XFil
number of units in which divided orders are to be processed set of available orders set of orders that can be processed in unit l objective function set of available units set of available units for Oi set of available units for the pair of orders (Oi , Ok ) scheduling horizon (the finishing time for the whole set of orders) total number of orders (each order corresponding to a specific product) number of batches of Oi to be processed when the unit l is assigned to Oi whole set of predecessors of order Oi set of customer orders that can be processed just before Oi in unit l, i.e. feasible predecessors of Oi in unit l total number of processing units release time of unit l set of feasible successors to order O set of customer orders that can be processed immediately after Oi in unit l, i.e. feasible successors to Oi in unit l cleaning time for the ordered pair of jobs (Oi , Om ) in unit l total processing time of order Oi in unit l. processing time of a batch of Oi in unit l binary variable denoting that order Oi is the first being processed in unit l
Xmil
XMil
binary variable denoting that the processing of order Om takes place in unit l just before the campaign for Oi amount of the Oi order to be processed in unit l
References [1] J.F. Pekny, D.L. Miller, Exact solution of the no-wait flowshop scheduling problem with a comparison to heuristic methods, Comput. Chem. Eng. 15 (1991) 741–748. [2] R.F.H. Mussier, L.B. Evans, An approximate method for the production scheduling of industrial batch processes with parallel units, Comput. Chem. Eng. 13 (1989) 229–238. [3] E. Kondili, C.C. Pantelides, R.W.H. Sargent, A general algorithm for short-term scheduling of batch operations. I: MILP formulation, Comput. Chem. Eng. 17 (1993) 211–227. [4] N. Shah, C.C. Pantelides, R.W.H. Sargent, A general algorithm for short-term scheduling of batch operations. II: computational issues, Comput. Chem. Eng. 17 (1993) 229–244. [5] J. Cerda, G.P. Henning, I.E. Grossmann, A mixed integer linear programming model for short-term scheduling of single-stage multiproduct batch plants with parallel lines, Ind. Eng. Chem. Res. 36 (1997) 1695–1707. [6] M. Yuceer, R. Berber, Z. Ozdemir, Automatic generation of production scheduling problems in single-stage multi-product batch plants: some numerical examples, in: (Proceedings of International Conference of Computational Methods in Sciences and Engineering—ICCMSE 2004, Athens, Greece, 19–23 November 2004), Lecture Series on Computer and Computational Sciences, vol. 1, VSP International Science Publishers, 2004, pp. 657–660. [7] C.A. Floudas, X. Lin, Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review, Comput. Chem. Eng. 28 (2004) 2109–2129. [8] X. Lin, S.L. Janak, C.A. Floudas, A new robust optimization approach for scheduling under uncertainty: I. Bounded uncertainty, Comput. Chem. Eng. 28 (2004) 1069–1085. [9] S. Engell, A. Markert, G. Sand, R. Schultz, Aggregated scheduling of a multiproduct batch plant by two-stage stochastic integer programming, Optim. Eng. 5 (2004) 335–359. [10] A. Alle, L.G. Papageorgiou, J.M. Pinto, A mathematical programming approach for cyclic production and cleaning scheduling of multistage continuous plants, Comput. Chem. Eng. 28 (2004) 3–15. [11] W. Wu, M. Ierapetritou, Cyclic short term scheduling of multiproduct batch plants using continuous-time representation, Comput. Chem. Eng. 28 (2004) 3–15. [12] B.A. Murtagh, M.A. Saunders, MINOS version 5.4 User’s Guide, Systems Optimization Laboratory, Stanford University, Stanford, CA, 1975.