Optimal Short-Term Scheduling of Large-Scale ... - ACS Publications

11002

Ind. Eng. Chem. Res. 2009, 48, 11002–11016

Optimal Short-Term Scheduling of Large-Scale Multistage Batch Plants Pedro M. Castro,*,†,§ Iiro Harjunkoski,‡ and Ignacio E. Grossmann§ Departamento de Modelaçaõ e Simulaçaõ de Processos, Instituto Nacional de Engenharia, Tecnologia e InoVaçaõ, 1649-038 Lisboa, Portugal, ABB Corporate Research Center, Wallstadter Strasse 59, 68526 Ladenburg, Germany, and Department of Chemical Engineering, Carnegie Mellon UniVersity, Pittsburgh, PennsylVania 15213

This paper presents a new decomposition algorithm for the optimal scheduling of large-scale multiproduct plants containing a large number of orders. Rather than tackling highly complex, full-space problems that cannot be solved in reasonable time, the complete set of orders is scheduled sequentially by considering one, or a couple of them, at a time. As we proceed through the iterations, previously scheduled orders can be partly rescheduled to allow for some flexibility while keeping the combinatorial complexity at a manageable level. Once a complete schedule is obtained, the same concept is applied to improve the schedule locally. The user can choose to rely on either a unit-specific or a sequencing variable based continuous-time mixedinteger linear programming model. In addition, there are other parameters that affect how the decomposition is carried out, so the algorithm is highly versatile and adaptable to problems of varying sizes. The largest problem solved is a real-life, 50-order, 17-unit, 6-stage problem, for which a very good solution can be found in less than 1 min of computational time. 1. Introduction The increasingly large literature in the scheduling area highlights the successful application of optimization approaches to an extensive variety of challenging problems.1 This important achievement comes mainly from the remarkable advances in modeling techniques, algorithmic methods, and computational technologies that have been made in the past decade or so. However, there is still a significant gap between theory and practice. New academic developments are mostly tested on relatively small problems, whereas real-world applications consist of hundreds of batches, dozens of pieces of equipment, and long scheduling horizons. In order to make exact methods more attractive for real-world applications, efforts have been increasingly oriented toward the development of systematic techniques that allow maintaining the number of decisions at a reasonable level, even for large-scale problems. Manageable model sizes may be obtained by applying heuristic model reduction methods, decomposition methods, or aggregation techniques. Once an initial solution is generated within reasonable computational time, gradual improvement through optimization based techniques can be obtained with modest computational effort. Although these techniques can no longer guarantee optimality, this may not be critical in practice due to the following: (i) very short time is available to generate a solution; (ii) a theoretical optimality easily gets lost due to the dynamic nature of industrial environments; (iii) implementing the schedule as such is often limited by the real process; (iv) only a subset of the actual scheduling goals is taken into account. Roslo¨f et al.2 proposed a systematic method for complex industrial scheduling and rescheduling. It consists of an easyto-implement reordering procedure that improves schedules toward optimality. Starting with a feasible schedule, a subset of the jobs is released and then inserted back into the schedule between the other fixed jobs. Different criteria can be used for * To whom correspondence should be addressed. Tel.: +351210924643. Fax: +351-217167016. E-mail: [email protected].. † Instituto Nacional de Engenharia, Tecnologia e Inovaçaõ. ‡ ABB Corporate Research Center. § Carnegie Mellon University.

job selection, including random picks, and the computational complexity can be set to a desired level in relation to the total number of jobs present in the system and the available computer capacity, by adjusting the number of simultaneously released jobs. The authors used a continuous-time model with sequencing variables and illustrated their algorithm with an industrial case study from a Finnish paper-converting mill, involving a single processing unit and sequence-dependent changeovers. They also gave guidelines on how to adapt the algorithm to account for parallel units. With parallel units, rescheduling may involve reordering as well as reallocation, as in the algorithm of Meńdez and Cerda´.3 The authors also relied on a continuous-time, mixed-integer linear programming model (MILP) with sequencing variables, suitable for single stage batch plants with no resource constraints other than equipment. Following the strategy used by Pinto and Grossmann4 in their slot based continuous-time model for multistage plants, they noted that compact models can be derived after embedding preordering rules in the full-space mathematical formulation, in what may represent a suitable alternative to quickly generate good production schedules for large-scale facilities. Indeed, by applying an increasing slack time rule, they were able to generate solutions with lower total earliness than those previously reported. The concept behind previous2,3 rescheduling approaches of restricting the degrees of freedom of all but a subset of orders can be also used to derive a good initial schedule. Such a procedure has been called constructive scheduling by Castro et al.,5 who proposed two alternatives for scheduling multistage plants with sequence-dependent changeovers. In the first, the full set of orders (e.g., 30) was considered with the best five according to makespan being scheduled. The solution was then used to fix the unit release times for the next iteration. In other words, the approach completely freezes the schedule for subsequent iterations. The second approach was identical to that of Meńdez and Cerda´.3 Orders were allocated to iterations by following the lexicographic sequence. Those being considered for the first time could be assigned to every unit and take any order (position) in the sequence, while previously scheduled

10.1021/ie900734x CCC: $40.75  2009 American Chemical Society Published on Web 11/06/2009

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009

orders could change their starting times but not their unit assignments nor their relative positions with respect to other previously scheduled orders. Overall, the latter method had more flexibility and was found to be more robust and better suited for the problem at hand. Nevertheless, there is always a limit where it is no longer possible to account for all orders, where schedule freezing cannot be avoided. It can be viewed as the boundary between short- and medium-term scheduling. Medium-term scheduling often relies on rolling-horizon algorithms, which were first proposed by Dimitriadis et al.6 They separate the overall scheduling horizon into detailed and aggregate time blocks. The problem is then solved in a sequence of iterations, where the previous solution of the detailed time block is fixed and part of the remaining aggregate time block gets to be scheduled in detail. The most critical issue, besides specifying the length of the detailed time block, is finding a simplified (aggregate) model that leads to decisions that can actually be implemented in practice when considering the detailed model. In a multiproduct plant such decisions involve selecting the products to be considered on a particular subhorizon. In the approach of Lin et al.7 this is done accounting for demand distribution, unit utilization, and the limits on the complexity of the resulting scheduling problem. Their work was extended by Janak et al.,8 which allows violation of the complexity limit to ensure that each subhorizon contains a minimum number of days. Products are selected if (1) there is an order for the product in that subhorizon, or within a given amount of time into the future; (2) the product is used as a raw material for another product that is also included; (3) the product was still processing in the previous subhorizon; and (4) it is a campaign product included in a campaign for the current horizon. The underlying scheduling model is an enhanced version of the well-known, state-task network, unit-specific continuous-time formulation of Ierapetritou and Floudas.9 Overall, a large-scale industrial plant consisting of over 80 pieces of equipment and including processing recipes of hundreds of different products was effectively scheduled. Rescheduling of the same large-scale plant was discussed in ref 10. The goal was to provide an immediate response to an unexpected event such as equipment breakdown or the addition/ modification of orders, while taking into account the schedule currently in progress, particularly those parts that were not affected. Starting with a nominal schedule from the decomposition framework,8 a set of rules is proposed to fix tasks in the subhorizon that is in production at the time of the event. The simplified MILP model is then used to perform rescheduling, while the remaining time horizon is rescheduled from scratch in the same manner as when generating the nominal schedule. This is required, since the changes in inventory, demand satisfaction, and unit availability can make the previous solution infeasible or highly suboptimal. Overall, an updated production schedule could be determined in a reasonable amount of CPU time. Metaheuristic methods such as genetic algorithms are an alternative to mathematical programming for scheduling and rescheduling. He and Hui11 presented an intelligent evolutionary algorithm for single stage batch plants that combines metaheuristic methods with heuristic rules to solve large-scale problems with modest computational effort. At a given iteration, the algorithm selects the best heuristic rule from a set of prespecified ones (e.g., earliest completion time, earliest start time, shortest process time), instead of relying on a single rule throughout the search for the optimal solution. Their automatic rule selection and rule sequence evaluation methods were later improved,12

11003

and a new automatic rule combination method was proposed that exhibited an order of magnitude higher convergence speed with increased search ability. The purpose of this paper is to systematically address the short-term scheduling of multistage batch plants following the decision, for each particular product, of the number and size of batches to produce, resulting in a set of orders with fixed processing times and release/due dates. The details about the problem statement are given in section 2. A complex algorithm is proposed that can be parametrized for the fast and efficient solution of problems of varying size. The key idea is explained in section 3, which is essentially the one proposed by Ro¨slof et al.2 and further explored by Meńdez and Cerda´.3 The novel aspects are the following: (i) the use of a multiple time grid continuous-time model instead of one based on sequencing variables; (ii) the introduction of an intermediate level of complexity between the full-space problem and the option of fixing both the order-unit assignments and relative ordering of previously scheduled orders; (iii) validation of the algorithm through the solution of example problems of moderate size for which the optimal solution is known, in order to measure the optimality gap and difference in total computational effort. Nevertheless, and concerning the first novel aspect, the user can choose between the multiple time grid MILP model of Castro and Grossmann14 or the one with global precedence sequencing variables of Harjunkoski and Grossmann.15 The required model constraints are given in section 4, featuring dynamic domains so that they are suitable for the solution of the full-space as well as the subproblems that are generated by the decomposition algorithm. While the basic underlying idea is the same, the algorithm has been divided into a constructive scheduling part (section 5), which generates a good feasible schedule, and a rescheduling part (section 6), which normally improves the solution through local search. Full details are given that include the specification of the dynamic sets that define the constraints domain. Extensive computational studies are performed in section 7. Analysis of the results is performed from a comparative perspective in terms of (i) fullspace model vs new decomposition algorithm, (ii) fixed vs variable relative positioning strategy, (iii) alternative preordering heuristics for order-iteration allocation, (iv) time grid vs sequencing based model, (v) number of orders considered per iteration, (vi) performance of constructive vs rescheduling part. Finally, conclusions are given in section 8. 2. Problem Definition Consider a short-term scheduling problem of multistage, multiproduct batch plants. Given are a set I of orders that must follow a sequence of processing stages k ∈ K to reach the condition of final products. It is assumed that all orders go through all stages in the same sequence. A set m ∈ M of equipment units is available at the plant, each allocated to a given stage. Set Mk includes those belonging to stage k. The processing times pi,m, release ri, and due dates di are assumed to be known. It is also assumed that setup times are sequence independent (included in the processing time) and that there is unlimited intermediate storage (UIS). The objective is to minimize the makespan. The process representation is shown in Figure 1 in the form of a resource-task network16 (RTN). The plant resources (circles) are the equipment units and material states, while the processing tasks are represented as rectangles. The material state is directly associated with the stage where it is produced. Thus,

11004


Figure 1. RTN representation for order I1, featuring a total of |M| units and |K| stages.

Figure 2. Illustration of constructive scheduling algorithm for fixed relative positions (FRP) and one order at a time (NOS ) 1).

processing tasks belonging to stage k consume material at state k - 1 and produce material at state k. 3. Key Idea In large-scale scheduling problems the number of orders typically exceeds by far the number of equipment units. We acknowledge this fact by proposing a decomposition method that reduces the complexity by scheduling a subset of the orders at a time. Deriving a schedule for a multistage plant with parallel units can be viewed as involving two decision levels: (i) assigning orders to units; (ii) sequencing orders on every unit. We follow this hierarchy to set different degrees of freedom for the orders. Orders being considered for the first time have total freedom, i.e., they can be assigned to every possible unit and take any position in the sequence. On the other hand, the timing of events of previously scheduled orders is allowed to change, but not their assignments to units. Furthermore, depending on the method, some sequencing restrictions may apply. Overall, to construct the full schedule, several iterations (set J) are involved. Given the number of orders to schedule with total freedom, NOS, the number of iterations is determined by eq 1. |J| ) |I|/NOS

(1)

The other aspect concerns the order-iteration assignments, which can be done in a variety of ways ranging from preordering heuristics to a random choice. We will be assuming that the elements of set Ij are the orders being considered for the first time in iteration j. They will also be a part of subsequent iterations up to the last one, so the resulting mathematical

problems increase in complexity as we go from one iteration to the next one. 3.1. Fixed Relative Positions (FRP). The first alternative of the constructive scheduling procedure keeps complexity to a minimum by assuming that previously scheduled orders cannot change their relative position in the sequence. This is illustrated in Figure 2 for one order at a time (NOS ) 1) with a simple example involving four units (two units per stage). We rely on the concept of time intervals (slots), which are unit-specific; i.e., events occurring in different units are affiliated with different time grids. In particular, we assume that a single time slot per order and per stage is enough to account for all valid alternatives within the constrained scenario. To make the diagram compact, timing issues between slots of different units have been neglected (e.g., an order can only start to be processed at stage k + 1 after the completion of stage k). In the first iteration (Figure 2), there is a single order under consideration (blue) so generality is ensured by postulating one slot for every time grid. The solver will then decide to assign the order to either (a) M1 or M2 in stage 1 or (b) M3 or M4 in stage 2. Let us assume that M1 and M4 are chosen. For j ) 2, the time grids no longer look the same. While vacant units remain with a single time slot, two intervals are added to the grids of M1 and M4 to account for all possibilities: production of red before or after blue. Notice that, in both units, the blue order is preassigned to slot number 2. The outcome is now production of the red order in M2 and M4 (prior to the blue order). Moving on to the third iteration (j ) 3), there will be units with one, three, and five time intervals and, again, previously tackled orders are fixed to a certain time slot. In this respect, in unit M4, the red order is fixed to slot 2 and the blue


11005

Figure 3. Illustration of constructive scheduling algorithm for fixed relative positions (FRP) and two orders at a time (NOS ) 2).

Figure 4. Illustration of constructive scheduling algorithm for variable relative positions (VRP) and one order at a time (NOS ) 1).

to slot 4 so that their relative position remains unchanged and the new green order is given full flexibility (production prior/ after both orders or in between). An additional order (black) is included to end the illustration at the start of iteration five. The same concept can be applied to other values of parameter NOS. Minor changes are required as can be seen in Figure 3 for the same set of orders and NOS ) 2. Now we start to prepostulate two time slots for every unit to allow production of the red and blue orders on the same unit. Assuming that this is the optimal choice for unit M4, with red before blue, the red order is fixed to slot 3 and the blue to slot 6, in iteration 2. This allows for the green and black orders to be assigned to the two first/intermediate/last slots. Needless to say, the solution space increases with a larger value of NOS up to the full-space limit (NOS ) |I|, |J| ) 1). However, so does the complexity per iteration, as will be seen later, so a trade-off is involved. 3.2. Variable Relative Positions (VRP). The alternative to fixed relative positions is to keep unit assignments unchanged but let previously assigned orders to move freely around the grid. This strategy is illustrated in Figure 4 for one order at a time (NOS ) 1). The first iteration is identical to FRP and involves a single time slot per grid. Again we assign the blue order to M1 and M4, but now these units require two rather than three slots for j ) 2. This is because the blue order is no longer fixed to a particular time slot. While fewer slots are being used, the complexity in terms of binary variables of the type (order, unit, slot) has in fact doubled (6 and 0 binaries for the red and blue orders in Figure 2, respectively, versus 6 + 6 ) 12 in Figure 4). Nevertheless, not all combinations are possible in VRP since we force unoccupied slots to be the last ones in the grid, in order to reduce solution degeneracy. The iteration outcome is thus blue in slots 1 of M1 and 2 of M4, and red in slots 1 of M2 and M4. Iteration 3 further illustrates this point, while iteration 4 shows a change in the relative position of the

blue and red orders in M4. It now seems that the best option is to employ a black-blue-red sequence. Overall, using variable relative positions widens the solution space but is more demanding computationally than its fixed positions counterpart. 4. Underlying Mathematical Formulations The constructive scheduling algorithm can decompose the original problem into different levels of complexity depending on the number of orders that are considered at a time. One important property is that it can also handle the full-space problem as a special case, by making NOS ) |I|. In such a case, there will be a single iteration and we define for every time grid as many time slots as the number of orders. Even though this ensures optimality, it is not the best approach in problems featuring units in parallel since significantly fewer slots will generally suffice and because the size of the grid has a profound effect on computational effort. The typical approach used by the research community is to use an iterative procedure that starts with a small number of slots and keeps solving the problem following a single increase in the number of slots until there is no improvement in the objective function or the problem becomes intractable. Such a procedure has already been classified by some authors17 as being heuristic to highlight the fact that only in the limit of enough time slots can optimality be guaranteed. From another perspective, we can say that the main drawback of time grid based continuous-time formulations when compared to those relying on sequencing variables for event representation1 is the uncertainty in the required number of time slots. This is multiplied by the number of equipment units in the case of unit-specific1 approaches, which is overcome by specifying the same number of slots for every unit. However, this is not the most logical thing to do in situations where the number of parallel units per stage has a wide variability, for

11006


which a couple of slots may suffice for some equipment while others may need several slots. One major property of the proposed algorithm is that it eliminates the uncertainty as part of the solution strategy. Recall (section 3) that, at each iteration, we add enough time slots to account for all possibilities and then eliminate those that were not assigned to orders, so that the size of the grid is kept at a minimum in subsequent iterations. Furthermore, the grids can and will normally grow at different rates, to ensure maximum computational efficiency. While it can be argued that the uncertainty in the number of time slots is relevant only when we want to guarantee optimality, which is not pursued in this case anyway, we know from experience that, due to its strong influence on solution quality, good solutions can only be achieved with an adequate number of time slots on every iteration. In fact, the proposed algorithm is, to the best of our knowledge, the first relying on a time grid based continuoustime formulation that successfully overcomes the number of slots specification issue. While others8 have proposed closely related decomposition methods for large-scale scheduling, they have neglected this issue, thus leading to another possible cause for a suboptimal solution outcome. This is the main novelty of the paper in the context of methods that gradually construct a new or improve an existing schedule, which have been around for a long time. In a recent paper13 we have compared three conceptually different unit-specific approaches that are able to solve single batch problems. The one by Castro and Grossmann14 has been shown to require a single time slot per order and per stage to find the optimal solution, and for this reason its reformulated version13 was classified as a minimum event points approach. On the other hand, the ones by Shaik and Floudas18 and Castro and Novais13 can be applied to the more general multiple product batches case. However, they rely on timing constraints that directly relate event points belonging to different units that normally force the number of slots past the theoretical minimum. Since the decomposition strategies FRP and VRP assume that a single time slot per order and per stage is enough to account for all valid alternatives within the iteration constrained scenarios, the approach of Castro and Grossmann14 is the chosen one. Precedence based continuous-time models are an alternative to models relying on time grids. For comparative purposes, we also implement the general precedence model of Harjunkoski and Grossmann.15 The most relevant features are briefly described next, where the main differences with regard to the full-space models are stricter domains for variables and constraints and the introduction of slack variables to always ensure that a solution is returned. This is a critical issue in a decisionmaking tool that is to be used by industry. 4.1. Time Grid Based (TG). The continuous-time model of Castro and Grossmann14 employs multiple time grids, one for each equipment unit. Binary variables Ni,m,t identify the execution of order i in unit m at event point (slot) t. Variables Rm,t give the excess amount of unit m at point t and, although not required, are also defined as binary for performance reasons. The two sets are related through the excess resource balances, eq 3. The remaining variables are all continuous non-negative. The absolute time of event point t in grid m is given by Tt,m. Two consecutive event points must be spaced by at least as much time as the processing time of the order being executed, eq 4. Variable TDi,k gives the transfer time of order i in stage k. On the one hand, it must be greater than the finishing time of the

order in stage k, eq 5. On the other hand, it must be lower than the order’s starting time in stage k + 1, eq 6. The makespan is 1 and S2i , are slack variables. These three given by MS, while St,m variables are involved in the objective function, eq 2. It involves makespan minimization together with penalizing terms for the slack variables, which will be equal to zero in many cases, in particular those not involving due dates. The contribution of each term is controlled through parameters R and β. Equations 7 and 8 are the release and due date constraints, respectively. In cases where the due dates cannot be met, eq 8 is relaxed through its slack variable S1t,m. Due to the form of the constraint, the slack variables are needed for all pairs (event point, unit) and not just for those belonging to last stage units. Nevertheless, only the latter are penalized in the objective function (see eq 2). Finally, eq 9 ensures that each order is processed exactly once on each stage.

∑ ∑

min MS +

1 RSt,m +

∑N

Rm,t ) 1 -

i,m,t

∑

βS2i

(2)

∀m ∈ M, t ∈ T mactive

(3)

active m∈M|K| t∈T m

i∈Iactive

i∈Im,t

Tt+1,m | t∉Tmlast + MS| t∈Tmlast - Tt,m g

∑N

∀m ∈ M, t ∈ T mactive

i,m,tpi,m

(4)

i∈Im,t

TDi,k g Tt,m + Ni,m,tpi,m - H(1 -

∑

Ni,m,t')

active t'∈T m t'gt i∈Im,t'

∀k ∈ K, k * |K|, m ∈ Mk, t ∈ T mactive, i ∈ Im,t TDi,k-1 e Tt,m + H(1 -

∑

(5)

Ni,m,t') ∀k ∈ K, k *

active t'∈Tm t'et i∈Im,t'

1, m ∈ Mk, t ∈ Tmactive, i ∈ Im,t (6) Tt,m g

∑N

i,m,t(ri

∑

+

i∈Im,t

k'∈K k'k

∀k ∈ K, m ∈ Mk, t ∈ T mactive

(8)

i∈Im,t

∑ ∑

Ni,m,t ) 1 ∀i ∈ Iactive, k ∈ K

(9)

active m∈Mk t∈Tm i∈ Im,t

The following sets of constraints are redundant but were found to considerably improve the performance of the mathematical formulation.14 Equations 10 and 11 act as lower and upper bounds on the transfer times. Again, slack variables are used (S2i ) to allow for the orders due date to be violated. Nevertheless, it is important to highlight that variables TDi,k have more freedom than their Tt,m counterparts, so it is more 1 to be active. Thus, one should likely for slack variables St,m make R . β. Equation 12 avoids solutions where an order follows a forbidden path. Equations 13 and 14 are responsible for reducing the integrality gap. The number of symmetric


solutions is reduced by enforcing the allocation of tasks to event points with as low an index as possible (see eq 15, Figure 3, Figure 4, and section 3.2).

min MS +

∑

RS2i

11007

(16)

i∈Iactive

Tfi',k g Tfi,k + pi',m - max di''(3 - Xi,i',k - Yi,m i''∈Iactive

Yi',m) ∀i ∈ Iactive, i' ∈ Iactive, i' > i, k ∈ K, m ∈ Mk ∩ Mi ∩ Mi'

TDi,k g ri +

∑ ∑ ∑

Ni,m,tpi,m ∀i ∈ Iactive, k ∈ K, k * |K|

active m∈Mk' t∈Tm i∈Im,t

k'∈K k'ek

Tfi,k g Tfi',k + pi,m - max di''(2 + Xi,i',k - Yi,m i''∈Iactive active

(10)

, i' ∈ Iactive, i' > Yi',m) ∀i ∈ I i, k ∈ K, m ∈ Mk ∩ Mi ∩ Mi' Tfi,k e Tfi,k+1 -

TDi,k e di + S2i -

∑ ∑ ∑

(11)

∑

(18)

Yi,mpi,m ∀i ∈ Iactive, k ∈ K, k * |K|

m∈Mk+1 m∈Mi


(19)

active m∈Mk' t∈Tm i∈Im,t

k'∈K k'>k

(17)

Tfi,k g

∑

Yi,m(ri + pi,m) ∀i ∈ Iactive, k ∈ K, k ) 1

m∈Mk m∈Mi

(20)

∑

(Ni,m,t | i∈Im,t + Ni,m',t | i∈Im',t) e

active t∈Tm i∈Im,t∨i∈Im',t

1 ∀i ∈ Iactive, m ∈ Mi, m' ∈ Mi, (m, m') ∈ FP

∑N

Tt,m e MS -

∑

i,m,t(pi,m

(22)

∑

+

min pi,m') ∀k ∈ K, m ∈ Mk, t ∈

Tmactive

Yi,m ) 1 ∀i ∈ Iactive, k ∈ K

(23)

(13)

Yi,m + Yi,m' e 1 ∀i ∈ Iactive, m ∈ Mi, m' ∈ Mi, (m, m') ∈ FP

∑

i∈Iactive m∈Mi

TDi,k e MS m∈Mk'

MS g Tfi,|K| ∀i ∈ Iactive

m∈Mk m∈Mi

m'∈Mk' k'∈k m'∈Mi k'>k

k'∈k k'>k

(21)

(12)

i∈Im,t

∑ ∑ ∑

Tfi,k e di + S2i ∀i ∈ Iactive, k ∈ K, k ) |K|


(14)

active Tm i∈Im,t

Rm,t g Rm,t-1 ∀m ∈ M, t ∈ Tmdeg

Yi,mpi,m e MS - min (

∑

min pi,m') - min (ri +

i∈Iactive k'∈K m'∈Mk' m∈Mi k'>k

∑

k'∈K k'>k

(24)

min pi,m') ∀k ∈ K, m ∈ Mk,

m'∈Mk'

i∈Iactive m∈Mi

∪ Mi * L

i∈Iactive

(25)

5. Constructive Scheduling Algorithm

(15)

4.2. Sequencing Variables Based (SV). In the continuoustime model of Harjunkoski and Grossmann,15 order sequencing in stage k is explicitly defined through binary variables Xi,i′,k (with i′ > i). Binaries Yi,m assign the execution of order i in unit m, while positive continuous variables Tfi,k give the order’s ending time in stage k. When compared to the time grid based formulation, it is harder to relate orders being executed in the same stage (eqs 17 and 18) and easier to ensure proper material transfer between consecutive stages (eq 19). The release date constraint, eq 20, is very similar to eq 7 and the due date constraint does not require a big-M term (eq 21). Since there is a single set of timing variables, only the slack variables S2i are 1 needed. They have a similar importance to slack variables St,m in the previous model and are thus given the same weight R in the objective function with regard to the makespan, eq 16. The latter must be greater than the ending time of all orders in the last production stage, eq 22. Equations 23 and 24 have the same effect of eqs 9 and 12, respectively. Finally, eq 25 is not absolutely necessary, but it tightens the formulation.14

In the previous section, the mathematical formulations were presented in a general form, meaning that the domain of the variables and constraints was defined using dynamic sets. The proposed constructive scheduling algorithm will simply manipulate the dynamic sets.21 Nevertheless, it does need to access problem data and the value of some model variables from the preceding iteration. To distinguish the variable from its value (level), an attribute (.l) is employed for the latter. Four decisions need to be made prior to algorithm execution. The first is related to the extent of the decomposition process, where one selects how many orders are considered at a time through parameter NOS, and consequently sets the number of iterations. The second involves rank ordering to specify the way in which orders are distributed over the iterations (set Ij). The underlying mathematical formulation is the third choice. Last but not least, we have to decide on the constructive scheduling concept, where fixed and variable relative position options are available. While these were illustrated in section 3 solely for the time grid based formulation, the concept is the same for the sequence based model as one will see in section 5.2. 5.1. Details for Time Grid Based Model. A particular iteration j starts with the specification of the orders that are being considered for the first time Icurrent; see Figure 5. These are then added to the previous orders (Iprevious) to build the set of active orders, Iactive. The next step is to determine the number of orders already assigned to unit m, npm, by checking the current value

11008


Figure 5. Constructive scheduling algorithm using formulation TG.

of binary variables Ni,m,t. Afterward, we decide whether to go for a fixed relative position (FRP) strategy or its variable counterpart (VRP). With FRP, we need to determine the position in the sequence (posi,m) for every order i allocated to unit m. The required module M1 is given in Figure 6, where the tricky part is not to count the empty slots (see Figures 2 and 3). The dynamic sets updating block is shown on the middle left of Figure 5 for FRP. The first, Mi, gives the units that are processing or can process order i. The former applies to previous orders that have been assigned to the unit (posi,m g 1). Current orders can be allocated to the unit whenever their processing time, pi,m, is greater than zero. A critical aspect concerns the specification of the active time slots for unit m, Tmactive. We must leave one slot for each already assigned order, npm, plus NOS + npm(NOS) slots to allow new orders to be placed before and after the positions already occupied (refer again to Figures 2 and 3). The last of these slots will be the sole member of set Tlast m . The elements of set Im,t are the orders that can be processed in unit m at slot t. Previous orders have known positions, while current orders can be allocated to all other empty slots, provided that the unit is suitable. Finally, Tmdeg gives the vacant time slots in unit m whose predecessor is also free. It sets the time domain of the constraint that reduces solution degeneracy (eq 15), which will force current orders to be allocated to slots in the immediate vicinity of those already occupied. Naturally, it will be nonempty only for NOS g 2. For the VRP strategy, updating the elements of the dynamic sets is more straightforward. This is because the relative positions from previous iterations are no longer relevant since everything but the order-unit assignments is allowed to change. The definition of set Mi is essentially the same and it is now useful to define set Im, which gives the orders that can be processed in unit m. Each order will need exactly one time slot, Tmactive, and can be executed on each and every one, Im,t (refer to Figure 4 and section 3.2). As a consequence, all slots but the first (Tmdeg) will be subject to the fewer degenerate solutions constraint.

Following the dynamic sets update, the time grid based model can be solved to find the optimal schedule for the active set of orders. We then continue to the next iteration and repeat the process until the moment all the orders have been scheduled. After the last iteration, the algorithm reports the final solution. 5.2. Details for Sequence Based Model. The algorithm that uses the sequence based model is given in Figure 7. The module to determine the position of previously scheduled orders is included in the diagram, giving the impression that it is significantly more complex than its TG counterpart when in fact it is only slightly more complex. The added complexity is due to the concept of general precedence, making it harder to find which binary variables to fix for the fixed relative positions strategy (FRP). This is done explicitly by setting upper and lower bounds (attributes .up and .lo, respectively) instead of relying solely on dynamic sets. In contrast, the time grid based approach is closely related to the concept of immediate precedence, despite the fact that there can be vacant slots in between orders that complicate the algorithm. The part up to deciding the decomposition strategy is basically the same as before. For a particular iteration, one identifies the current, previous, and active set of orders, identifies the number of orders assigned to unit m, and restricts order-unit assignments through the set Mi. This is all that is needed for the variable relative positions strategy (VRP) prior to solving model SV. The solution yields the ending times of order i in stage k (Tfi,k) and the assignments (Yi,m), and with these values we calculate the finishing times on the units (TAUXi,m), which are useful to determine order positioning. We then proceed to the next iteration or terminate the algorithm by reporting the final schedule. Before solving the mathematical programming problem corresponding to the FRP option, we need to determine the position of every order in the sequence for each unit (module M2) and, based on this, fix the proper sequencing variables (M3). These modules are now described for an intermediate iteration (j > 1) and for a certain unit m. Starting with position


11009

M2, index j′ is not necessarily equal to the position under consideration. Nevertheless, it is still required to identify when all positions have been searched. Afterward, we proceed to the next equipment unit. Once all units have been considered, the model is solved to find the optimal schedule. 6. Rescheduling Algorithm

Figure 6. M1 module to determine the position in the sequence of previously scheduled orders (FRP option).

1 in the sequence (j′ ) 1), point 1 in Figure 7, and assuming that order i′ ) 1 has been previously assigned to m, the expression in point 2 holds true. If the order’s ending time in the stage to which m belongs to is not the lowest among all orders i′′ allocated to m, point 3, we proceed to the next order, point 4. We will eventually return to the same decision point 5 with the actual first order. When this occurs, we set its position posi′,m ) j′, and update the value of the auxiliary parameters, point 6, so that we can move to the next position in the sequence, point 7. After all positions have been found for unit m, we can go on and fix the sequencing variables. Starting with j′ ) 1, point A, we are again looking for the order i′ that has the first position in the sequence, location B (note that AUX ) 0). The binary variables will be fixed for all orders i′′ processed after i′ in unit m, point C. In point D we evaluate if i′′ has a higher index than i′. If the answer is true, the binary Xi′,i′′,k is in the model and it must be set to 1 since i′ globally precedes i′′, block E. Otherwise, it is its counterpart Xi′′,i′,k that is part of the model. It is set to zero in block F since i′′ does not precede i′. After all sequencing variables relating i′ with the other orders already assigned to the unit have been set, we increment the auxiliary parameter in block G. It may happen that the order allocated to the next position in the sequence has a higher index than i′. If this is true, point H, we can set more binary variables without further increasing index j′, location I. In other words, contrary to what happens in module

The constructive scheduling algorithm is often capable of providing a good solution to the problem, which might be considered suitable for practical purposes. However, it may happen that the returned solution does not meet all due date constraints, which is a serious flaw in cases where such constraints can actually be met. In addition, the solution is somewhat dependent on the order-iteration selection (subsets Ij), as will be seen in the next section, so performing a local search around the best schedule seems perfectly logical. To overcome these weaknesses, we now present a rescheduling step that gradually improves the previous solution at low computational cost. It is based on the ideas first proposed by Roslo¨f et al.2 and employed by Meńdez and Cerda´3 that are the ones behind the fixed relative positions strategy. The algorithm uses the sequence based model because it typically finds the optimal solution faster than the time grid based formulation. Speed is critical whenever multiple trials and limited computational resources are involved. The most complex part of the algorithm is finding order positioning and fixing the proper sequencing variables for the mathematical problem. Fortunately, the hard work has already been done for the constructive algorithm, so now we can simply call modules M2 and M3 (detailed in Figure 7). It is, however, important to highlight that in the constructive process the number of orders to consider increases from one iteration to the next, while now modules M2 and M3 will consider the full set of orders with the exception of those being rescheduled, regardless of the iteration number. The algorithm is shown in Figure 8. The basic idea is to release one or a couple of orders from the schedule, and try to reassign or resequence them to see if a better solution can be found. For simplicity, we assume the same order-iteration ranking as in the constructive step, given by Ij. This is because earlier orders took into account fewer orders than those that were scheduled later, when deciding unit allocation, so they have a higher potential for improvement and should be rescheduled first. Given are also the values of some SV model variables corresponding to the final schedule from the constructive algorithm. If model TG was used instead, the correspondence can be made using eqs 26 and 27. Yi,m.l )

∑N

i,m,t.l

∀i ∈ I, m ∈ M

(26)

t∈T

Tfi,k.1 )

∑ ∑ [N

i,m,t.l(Tt,m.1

+ pi,m)] ∀i ∈ I, k ∈ K

m∈Mk t∈T

(27) The initialization step basically consists of setting the current schedule as the best one. This involves storing the value of the objective function (OBJ), the orders assigned to unit m (Imbest), best , and (b) on and their finishing times (a) on every stage, Tfi,k every unit, TAUXi,m, which will be required in module M2. Next, we define the sets of (i) current orders to be rescheduled with total freedom and (ii) remaining orders (named again previous to allow modules M2 and M3 to be reused), which will be allowed to change their starting times but not their unit assignments nor their relative positions with respect to orders

11010


Figure 7. Constructive scheduling algorithm using formulation SV.

of the same set. The auxiliary parameters, TAUXi,m and posi,m, are then reset to zero for current and all orders, respectively, while npm gives the number of orders assigned to a particular unit without counting those to be rescheduled. The following steps involve resetting the domain of all the binary variables and determining for every order the elements of set Mi. Once the order positioning has been found (M2), and the proper sequencing variables have been fixed (M3) for every

equipment unit, the mathematical problem can be solved. Hopefully, we will be able to find a better solution and will need to update the best solution parameters. Otherwise, we will need to update the value of the model variables with those corresponding to the best solution. Either way, we need to update parameter TAUXi,m. We then proceed to the next iteration or terminate the algorithm after reporting the best found solution.


11011

Figure 8. Rescheduling algorithm using formulation SV.

7. Computational Results The performance of the constructive and rescheduling algorithms is evaluated through the solution of 13 example problems. The first eight (P3-P10) have been thoroughly studied13,14 and can be tackled by the full-space versions of time grid (TG) and sequencing variables (SV) models. The data for the remaining five (P11-P15) are given as Supporting Information. Problem P13 is adapted from a real industrial case involving 50 orders, 17 units, and 6 stages, while P12 considers just the first 30 orders. Neither release nor due dates are specified in P12 and P13. Data for P11, P14, and P15 were generated randomly. In order to ensure that priority is given to finding schedules that meet the due date constraints for all orders, we have set R ) 10 and β ) 1 (see eqs 2 and 16). While it is important to emphasize that such values affect the performance of the fullspace MILP models and decomposition algorithm, it is beyond the scope of the paper to evaluate their influence in computational times. The algorithms and MILP models were implemented and solved in GAMS 22.8, using CPLEX 11.1 with default options. The termination criterion was either a relative optimality tolerance ) 10-6 or a maximum computational time equal to (i) 60 000 CPU s for the full-space models, (ii) 3600 CPU s per iteration for the constructive scheduling algorithm, and (iii) 60 CPU s per iteration for the rescheduling algorithm. The hardware consisted of a laptop with an Intel Core2 Duo T9300 2.5 GHz processor, with 4 GB of RAM and running Windows Vista Enterprise. 7.1. Full-Space Models (FS) vs New Decomposition Algorithm (DA). Decomposition algorithms are useful for largescale problems, those for which current hardware capabilities prevent full-space models from finding a good or even a feasible

solution in a reasonable time. As can be seen from Table 1, the full-space sequencing variable model found no solution for P14 and P15 before an out-of-memory termination by the solver. Although it performs better for P13, the solution found up to 16.7 h of computational time is still distant (16.8% higher makespan) from the best solution found by the proposed decomposition method. On the other hand, the sequencing variable model (SV) is still capable of finding a good solution for P12 (8% longer makespan), while the time grid based model could not even find a feasible solution. This behavior is not new13,14 and was the main motivation for adopting SV as the underlying model for the rescheduling algorithm. However, if the time grid based model (TG) is able to find a feasible solution, it can be orders of magnitude faster at proving optimality, as can be seen for P7-P11. For P11, it could even find the optimal solution, contrary to SV. Note that a makespan of 207 was found using |T| ) 4 (see the schedule in Figure 9) and confirmed for |T| ) 5 in 15544 CPU s. It is also worthwhile to mention that CPLEX 10.2 was unable13 to find the optimal solution for P9 and P10 with TG up to the same maximum resource limit, whereas now (CPLEX 11.1) 10 min are enough to solve both problems to optimality. While these results illustrate the influence of solver options on computational performance, it is beyond the scope of this paper to test other values since there will always be a point where problems arising from full-space models become intractable. The new decomposition strategy has proven useful for the larger problems (P12-P15). As will be seen in detail in the next subsections, it can find close-to-best known solutions in a few minutes of computational time. Nevertheless, the best solutions for P12, P13 (shown in Figure 10), and P15 were found in respectively 2.5, 1.6, and 4.1 CPU h. The full-space model

11012


Table 1. Comparative Results for Full-Space Models and New Decomposition Algorithm (Optimal Solution in Bold) problem

features (|I|,|M|,|K|)

P3

(8,6,2)

P4

(8,6,2)

P5

(10,6,3)

P6

(10,6,3)

P7

(12,8,3)

P8

(12,8,3)

P9

(15,6,2)

P10

(10,8,4)

P11

(16,8,2)

P12

(30,17,6)

P13

(50,17,6)

P14

(50,12,6)

P15

(50,20,5)

a

method

model

FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS FS DA FS DA FS DA FS DA

TG SV TG,SV TG SV SV TG SV TG,SV TG SV TG,SV TG SV TG,SV TG SV TG,SV TG SV TG TG SV TG TG SV SV TG SV TG SV TG SV SV SV SV

best solution 793 796 793

1435

1394

1456

1655

235 251 252 254 207 219 230 no solution 21.791 20.168 36.136 30.053 no solution 984 no solution 535

CPU s 1.40 0.35 1.29 0.45 0.38 1.11 5.00 32.6 1.83 1.59 0.87 1.89 602 57 802 2.40 6.97 2 547 2.39 448 17 650 7.83 117 3 748 162 679 7 277a 3.29 60 000 60 000 8 855 60 000 5 833 33 061a 350 26 949a 14 776

Out-of-memory termination.

TG should be preferred in the remaining cases. Despite the 13 trials resulting from different combinations of (i) ordering heuristics, (ii) relative position strategy, (iii) values of parameters NOS, and (iv) underlying model, the decomposition algorithm (DA) was unable to find the optimal solution for P3 and P9-11. The one for P11 is particularly poor, 11.1% above the optimum. Note that in the rows belonging to DA (Table 1) we list the models that were able to find the best solution reported. In examples where both TG and SV were successful, the fastest run of SV (CPU time listed in the last column) was always faster than the fastest of TG. More detailed results can be found in the following tables.

Figure 10. Best found solution for P13.

7.2. Fixed (FRP) vs Variable Relative Positions (VRP). We start evaluating the performance of the proposed algorithm by looking into the two alternative relative position strategies. Recall that the variable relative position strategy (VRP) is not as restrictive as its fixed relative position counterpart, since only order-unit assignments are fixed, thus keeping most of the degrees of freedom. The advantage is that better solutions do result in the constructive step but only if each individual problem does not become too complex. Clearly, the results in Table 2 show that the VRP strategy is not appropriate for large-scale problems since it takes roughly 21 and 36 h for P12 and P13 to get to really poor solutions (58 and 171% above the best ones found). Note that the number of discrete variables (DV) is only about 40% of the full-space model (compare with Table 7), when for FRP with NOS ) 1 it was 1 order of magnitude lower. Interestingly, such solutions were significantly improved (20 and 58%) by the rescheduling part of algorithm, in less than 1 min. Note that the FRP strategy was used because the significantly higher computational times make VRP a bad option for rescheduling. In view of the above results, VRP is abandoned for the remainder of the paper. 7.3. Earliest Due Date (EDD) vs Minimum Slack Times (MST). The first decision in the constructive and rescheduling algorithms involves selecting which orders of the portfolio are considered in a particular iteration. Preordering heuristics have been used by some authors3,4,19 in continuoustime models to reduce the complexity, and were shown20 to compromise optimality by at least 4% on a single stage problem with parallel units. In the context of this paper, preordering is not as restrictive since it is still possible for order i′ to be processed before order i despite being scheduled in a later iteration. We have used two common heuristics, earliest due date (EDD), and increasing slack times (MST). The concept in both is to schedule the orders with a smaller time window (spani) first, with the latter using information concerning not only due dates but also release dates as well as the minimum processing times on every stage, see eq 28. While other heuristics could have been studied (e.g., most work remaining, shortest processing time, modified due date), we feel that those chosen clearly illustrate the importance of preordering heuristics in both solution quality and computational effort. spanMST ) d i - ri i

Figure 9. Optimal solution for P11.

∑

min pi,m ∀i ∈ I

m∈M k∈K m∈M i k

(28)

The solutions obtained after the constructive and rescheduling steps are shown in Table 3 for the time grid model and in Table


11013

Table 2. Computational Statistics for VRP Strategy with MST Heuristic, TG Model, and NOS ) 1 constructive step

rescheduling

problem

discrete variables

continuous variables

equations

CPU s

solution

solution

CPU s

P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13

92 88 220 220 330 328 232 300 222 2940 7058

150 146 318 318 450 446 332 430 232 3498 7978

229 225 532 532 725 722 472 777 485 6 343 13 887

1.57 1.45 5.23 3.62 4.28 6.28 8.15 210 5.19 74 926 130 252

844 848 1476 1433 1456 1655 261 286 244 31.893 81.386

844 826 1476 1433 1456 1655 261 270 244 25.050 51.933

0.92 1.09 1.69 1.80 2.42 3.05 2.91 2.56 2.74 20.0 58.2

Table 3. Solution after Constructive and Rescheduling Steps Using FRP Strategy and Model TG preordering heuristic EDD

preordering heuristic MST

problem

optimum/best

NOS ) 1

NOS ) 2

NOS ) 3

NOS ) 1

NOS ) 2

NOS ) 3

P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

793 793 1435 1394 1456 1655 235 252 207 20.168 30.053 984 535

844 826 1476 1453 1456 1655 258 263 272 21.333 30.640 1044 550

803 808 1435 1394 1456 1655 264 265 238 20.433 30.053 996 570

808 803 1502 1410 1456 1655 258 275 239 20.168 30.786 1010 568

844 826 1476 1433 1456 1655 272 267 238 21.699 30.733 1025 561

796 808 1435 1394 1456 1655 251 260 221 21.333 30.219 994 566

803 803 1502 1394 1456 1655 252 264 233 21.333 31.794 1011 566

Table 4. Solution after Constructive and Rescheduling Steps Using FRP Strategy Model SV preordering heuristic EDD

preordering heuristic MST

problem

optimum/best

NOS ) 1

NOS ) 2

NOS ) 3

NOS ) 1

NOS ) 2

NOS ) 3

P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

793 793 1435 1394 1456 1655 235 252 207 20.168 30.053 984 535

849 849 1476 1449 1456 1655 271 326 225 21.022 30.405 1026 563

803 848 1435 1394 1456 1655 250 262 239 20.433 30.218 995 547

808 793 1502 1410 1456 1655 255 263 246 20.433 30.133 996 535

849 849 1476 1433 1475 1655 260 267 226 21.983 30.133 1010 557

796 848 1435 1394 1456 1655 260 260 235 21.613 30.101 984 556

803 793 1502 1394 1456 1655 255 266 230 20.523 30.879 998 547

4 for its sequencing variable counterpart. In terms of success rate in finding the optimal/best known solution (identified in bold), the MST heuristic was successful in 19 trials as opposed to 20 by EDD, out of 78. However, if one calculates the quadratic sum of the relative errors, one finds a value of 0.231 for MST and 0.441 for EDD. In terms of computational time, the results for the largest problems (Tables 5 and 6) favor the EDD heuristic for TG and MST for SV. Thus, one can conclude that it is on average better to use the minimum slack times heuristic. 7.4. Time Grid (TG) vs Sequencing Variable (SV) Model. The sequencing variable model has a slightly better solution-quality performance than the time grid model with 20 vs 19 successful runs and total quadratic errors equal to 0.333 (SV) vs 0.339 (TG). Often, it also has a lower computational effort than TG even though the latter was faster for the larger problems with the EDD heuristic (see Tables 5 and 6). Overall, the differences are minor, so it is fair to say that the models have identical performances and so both can be employed. Finding out that the heuristic has a significant impact not only

Table 5. Total Computation Effort (CPU s) in Constructive Step Using FRP Strategy and MST Heuristic model TG

model SV

problem NOS ) 1 NOS ) 2 NOS ) 3 NOS ) 1 NOS ) 2 NOS ) 3 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

2.47 2.50 3.25 3.12 3.87 3.80 4.73 3.37 5.13 15.5 40.0 43.0 28.8

1.66 1.37 2.50 2.26 2.38 2.31 3.11 2.69 3.20 350 19 017 2 024 108

1.35 1.32 2.80 2.51 1.86 1.76 3.16 16.0 3.13 24 956 42 056 41 686 3 818

1.45 1.15 1.59 1.62 1.62 2.43 2.03 1.59 2.37 11.1 58.1 47.8 36.3

0.69 0.49 0.76 1.13 1.10 1.15 1.41 1.06 1.75 412 8141 76.8 135

0.58 0.52 0.94 0.75 0.83 0.96 1.93 1.29 1.44 15 085 39 575 27 393 708

on solution quality but also on total computational effort was very interesting. While the decomposition algorithm generates subproblems of identical size at certain iterations for the same

11014


Table 6. Total Computation Effort (CPU s) in Constructive Step Using FRP Strategy and EDD Heuristic model TG

model SV

problem NOS ) 1 NOS ) 2 NOS ) 3 NOS ) 1 NOS ) 2 NOS ) 3 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

2.55 2.51 3.33 3.21 3.76 3.84 4.90 3.27 5.00 15.5 40.0 37.4 27.0

1.69 1.31 2.52 2.44 2.15 2.22 3.36 2.37 3.37 59.9 4973 1326 103

1.32 1.34 3.04 2.68 1.95 1.87 3.13 3.79 3.24 8 397 27 621 40 558 3 731

1.18 0.98 1.39 1.78 1.92 2.01 2.33 1.48 2.69 10.3 46.9 59.8 26.8

0.61 0.76 0.88 0.89 0.96 1.43 1.47 0.97 1.47 247 18257 64.3 4014

0.37 0.48 0.81 0.74 0.69 1.11 1.57 1.17 1.44 14 862 35 539 26 330 14 444

value of NOS, it seems that some are just more difficult to solve to optimality than others. Another interesting result comes from the fact that the solutions from TG and SV are normally different for the same ordering heuristic and NOS value. It highlights another disadvantage of decomposition strategies. Since several problems are solved in sequence up to the generation of the final schedule, it is highly likely that degenerate solutions will occur. In other words, there will typically be more than one partial schedule with the same value of the objective function. The one selected by the solver will be used to fix certain assignment/sequencing variables, so it will affect the final outcome. Thus, while the values of the objective function may be the same up to some iteration, they can diverge so that alternative schedules are generated, one of them featuring a minimum makespan. It is beyond the scope of the paper to study other objective functions that are able to minimize this effect, but it may be an interesting research subject. For now, we just acknowledge that it is good to have alternative models to derive a few good schedules. 7.5. Influence of NOS Parameter. The decision on the number of orders to first consider (the value of parameter NOS) has a major influence on problem size and consequently on computational effort. All problems generated in the constructive step can be solved to optimality, with the exceptions being the last iterations of P12 and P14 (NOS ) 3) and P13 (NOS g 2). While an increase in NOS widens the feasible region associated with each subproblem, augmenting the probability of finding the optimal solution, it may also make it more difficult to solve the problems on each iteration to optimality due to the consideration of larger, more complex problems. Usually, increasing the value from 1 to 2 led to better solutions, but not from 2 to 3. In fact, the total quadratic error for TG, considering

both heuristics, was equal to 0.193, 0.059, and 0.086 (for SV we got 0.182, 0.077, and 0.074), respectively, for NOS ) 1, 2, and 3. Thus, in terms of the solution quality/computational effort trade-off, it is better to use NOS ) 2. Note, however, that this decision ultimately depends on the characteristics of the problem at hand and on hardware/solver efficiency. The whole purpose of the decomposition algorithm is to generate smaller mathematical problems than the full-space models. The actual values are given in Table 7 for the last problem solved (TG), which considers all problem orders. Model SV gives rise to smaller mathematical problems. Note that the last problem does not correspond to the largest one, in cases where |I| is not a multiple of NOS, which is the reason why the increase in the number of model entities is not roughly proportional to the increase in NOS. It is particularly interesting to find out that, for P3-P6 and P9-P11, the numbers of variables and equations of the last generated problems for NOS ) 3 are even higher than those resulting from the full-space model (FS). It is thus not totally surprising that the former set of problems is solved faster by the full-space models (recall the results in Table 1). Naturally, the size reduction becomes more evident as the number of orders, units, and stages increases. For the larger problems, P12-P15, one gets roughly 1 order of magnitude reduction in the number of binary variables for NOS ) 1, thus keeping the complexity at a manageable level. 7.6. Rescheduling Step. The purpose of the rescheduling step is to perform local searches to see if the solution from the constructive step can be improved. The results in Table 8 show that this was accomplished in most of the cases (67%). It is clear that the improvement rate increases with an increase in the value of NOS, which is not surprising since a wider feasible region is being considered. Given the observed behavior of the constructive step where the gains in solution quality with an increase in NOS do not compensate the much higher computational effort, for large-scale problems it is perhaps a better option to use the lowest possible value to generate an initial schedule, and a higher value in the rescheduling step. The drawback is that fewer iterations are involved, but this could be overcome by repeating the outermost loop in Figure 8, eventually using different ordering heuristics, possibly including a random order/iteration selection. Nevertheless, it was observed that the better the initial point from the constructive step, the better the solution after rescheduling, so the former should not be totally neglected. One relevant aspect of the rescheduling algorithm is that it successfully overcame all due date violations in some of the initial schedules from P7. The large penalties in the objective function led to major improvements during rescheduling (shown

Table 7. Influence of Parameter NOS on Computational Size of Last Problem Tackled in Constructive Step and Comparison to Full-Space Model TG discrete variables problem P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

continuous variables

equations

FS

NOS ) 1

NOS ) 2

NOS ) 3

FS

NOS ) 1

NOS ) 2

NOS ) 3

FS

NOS ) 1

NOS ) 2

NOS ) 3

159 159 320 320 1224 1224 576 440 544 6930 19 275 15 300 13 260

68 68 120 120 148 148 124 160 136 700 1173 2000 1020

132 132 228 228 288 280 192 304 272 1330 2290 2376 2040

186 186 252 252 474 465 408 336 288 2250 3150 3276 2820

186 186 371 371 1345 1345 628 511 593 7336 19 951 15 851 13 721

154 154 272 272 334 334 280 362 306 1612 2685 2702 2292

246 246 428 428 538 530 416 570 506 2588 4388 4454 3812

336 336 536 536 776 767 668 714 610 3830 5858 5954 5112

298 298 664 664 2138 2138 825 962 776 13 243 35 033 27 838 23 190

247 247 472 472 581 581 448 651 489 3005 5011 5051 4217

387 387 735 735 925 917 686 1019 789 4790 8130 8271 6959

531 531 956 956 1301 1292 1000 1311 1021 7001 10 868 11 071 9351

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009 Table 8. Solution Improvement (%) after Rescheduling Step for FRP Strategy and MST Heuristic NOS ) 1

NOS ) 2

NOS ) 3

problem

TG

SV

TG

SV

TG

SV

P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

0 3.28 0 0 0 0 0 7.61 4.80 5.17 1.36 0 1.06

0 0.59 0 0 33.7 0 4.06 6.97 4.24 3.51 0.74 0.69 1.76

1.49 0 0 0 0 0 4.20 0 2.64 1.16 2.46 0 0.35

0 0 0 0 27.8 0 1.89 0 0.84 5.14 1.62 0.10 0.71

0.62 0.62 0.07 3.13 43.5 0 2.70 1.86 7.54 1.61 1.34 0.69 0.88

0.62 7.03 0.07 3.13 43.5 0 3.77 1.12 4.56 5.12 0.34 0.80 1.26

Table 9. Total Computational Effort in Rescheduling Step for FRP Strategy, MST Heuristic, and TG Model in Constructive Step problem

NOS ) 1

NOS ) 2

NOS ) 3

P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

0.87 0.88 1.22 1.14 1.77 1.70 2.20 1.77 1.92 18.9 216 153 69.0

0.69 0.50 1.00 0.87 1.46 1.47 1.83 1.24 1.62 263 752 396 292

0.62 0.62 1.04 1.11 1.42 1.20 1.70 1.06 3.23 528 978 1029 657

in bold in Table 8). Overall, rescheduling does not involve too much time if a small value is set for the maximum resource limit (e.g., 60 CPU s); see Table 9. Similar results can be found when using model SV in the constructive step. Worst case, it is equal to such value times the number of iterations. Since the latter is strongly dependent on NOS, one can actually reduce the total time while achieving better solutions. It is also important to emphasize that when releasing a single order (i) all generated MILPs could be solved to optimality and (ii) the total computational time was under 4 min, regardless of the model and heuristic used. 8. Conclusions This paper has presented a flexible decomposition algorithm for the optimal scheduling of large-scale multistage batch plants with unlimited intermediate storage. Given several orders, the key idea is to first schedule one, two, or three at a time rather than tackling the full set simultaneously. Such orders will keep their unit assignments fixed in subsequent iterations, but their timing is allowed to change somewhat depending on the relative position strategy employed. Two alternatives have been proposed: fixed and variable. The former does not allow previously assigned orders to change their relative positions in the sequence, while the latter keeps most of the degrees of freedom for scheduling, potentially leading to better solutions. The drawback of the variable position strategy is that previous orders are fully resequenced over and over again, which has been shown to have a high computational cost. Fortunately, the fixed position strategy is able to keep the complexity at a manageable level, allowing generation of good solutions to the problem very rapidly. More specifically, we were able to tackle a real-life 50-order, 17-unit, 6-stage problem in less than 1 min.

11015

The decomposition strategy has been illustrated with two conceptually different continuous-time models, one relying on multiple time grids for event representation, and the other relying on global precedence sequencing variables. The former has a well-known disadvantage when compared to the latter: the need to specify the number of time slots in the grid, which strongly affects the quality of the solution and cannot be predicted accurately. A major breakthrough has been to remove this limitation, by specifying per iteration the minimum number required to consider all possibilities without compromising tractability of the mathematical problem. The algorithm has been shown to be very efficient and identical in performance to its sequencing variables counterpart. Ensuring feasibility is key for decomposition algorithms, which is typically achieved through the use of slack variables that relax some of the constraints. The dilemma is that one does not know if such variables are active due to the decomposition process or the problem itself. A local search strategy has been developed to generally improve the solution from the so-called constructive step, which succeeded in removing all pseudoinfeasibilities found. Other parameters affecting the algorithm’s performance that can be controlled by the user are the preordering heuristic and the number of orders first considered in the constructive step. Extensive computational studies have been performed, leading to the conclusion that it is better to use the minimum slack time (MST) rather than the earliest due date (EDD) heuristic for order-iteration assignment and that the best trade-off between solution quality and total computational effort is achieved for two orders per iteration. Work is underway to widen the scope of the algorithm to deal with sequence-dependent changeovers, which is often a very important issue in multistage multiproduct plants. While it is straightforward to adapt the sequencing variables model to such cases, the same is not true for the time grid based model. Two alternatives exist depending on the consideration of explicit or implicit changeovers, requiring respectively four- or threeindex binary variable models.22 Either way, the constructive step needs to be changed significantly. Acknowledgment The authors gratefully acknowledge financial support from Fundaçaõ Luso-Americana and the Center for Advanced Process Decision-making at Carnegie Mellon University. Supporting Information Available: Tables of data for problems P11-P15. This information is available free of charge via the Internet at http://pubs.acs.org. Nomenclature Sets FP ) forbidden paths, pairs of units (m,m′) where orders cannot go from unit m to m′ I/i, i′ ) orders Iactive ) active orders for scheduling model Imbest ) orders assigned to unit m in best found solution Icurrent ) orders being currently considered Iprevious ) previously scheduled orders Ij ) orders being considered for the first time in iteration j Im ) orders that can be executed in unit m Im,t ) orders that can be processed in unit m at interval t J/j, j′ ) iterations of constructive scheduling algorithm K/k, k′ ) stages

11016


M/m, m′ ) units Mi ) units that can process order i Mk ) units belonging to stage k T/t, t′ ) time slots (event points) Tmactive ) active event points for unit m Tmdeg ) vacant time slots in unit m whose predecessor is also free Tfix ) event points with fixed orders Tmlast ) last active event point in unit m Parameters AUX ) auxiliary parameter BEST ) objective of best found solution di ) due date of order i H ) time horizon NOS ) number of orders to schedule with total freedom npm ) orders already assigned to unit m pi,m ) processing time of order i in unit m posi,m ) position of order i in the sequence of unit m ri ) release date of order i TAUXi,m ) ending time of order i in unit m best Tfi,k ) ending time of order i in stage k in best found solution R ) contribution of first type of slack variables to the objective function β ) contribution of second type of slack variables to the objective function Variables MS ) makespan Ni,m,t ) binary variable that identifies the execution of order i in unit m at interval t OBJ ) objective function Tt,m ) absolute time of event point t on unit m TDi,k ) transfer time of order i in stage k Rm,t ) excess amount of unit m at event point t 1 St,m ) slack variable associated with event point t on unit m 2 Si ) slack variable associated with order i Tfi,k ) ending time of order i in stage k Xi,i′,k ) binary variable indicating that order i globally precedes order i′ in stage k Yi,m ) binary variable that identifies the execution of order i in unit m Variable Attributes .l ) current value (level) .lo ) lower bound .up ) upper bound

Literature Cited (1) Meńdez, C. A.; Cerda´, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art Review of Optimization Methods for Short-Term Scheduling of Batch Processes. Comput. Chem. Eng. 2006, 30, 913. (2) Roslo¨f, J.; Harjunkoski, I.; Bjo¨rkqvist, J.; Karlsson, S.; Westerlund, T. An MILP-based reordering algorithm for complex industrial scheduling and rescheduling. Comput. Chem. Eng. 2001, 25, 821.

(3) Meńdez, C.; Cerda´, J. Dynamic scheduling in multiproduct batch plants. Comput. Chem. Eng. 2003, 27, 1247. (4) Pinto, J. M.; Grossmann, I. E. A continuous time mixed integer linear programming model for short term scheduling of multistage batch plants. Ind. Eng. Chem. Res. 1995, 34, 3037. (5) Castro, P.; Meńdez, C.; Grossmann, I.; Harjunkoski, I.; Fahl, M. Efficient MILP-based solution strategies for large-scale industrial batch scheduling problems. In Computer Aided Chemical Engineering; Marquardt, W., Pantelides, C., Eds.; Elsevier: New York, 2006; Vol. 21 (part 2), p 2231. (6) Dimitriadis, A. D.; Shah, N.; Pantelides, C. C. RTN-Based Rolling Horizon Algorithms for Medium Term Scheduling of Multipurpose Plants. Comput. Chem. Eng. 1997, 21, S1061. (7) Lin, X.; Floudas, C. A.; Modi, S.; Juhasz, N. M. Continuous-Time Optimization Approach for Medium-Range Production Scheduling of a Multiproduct Batch Plant. Ind. Eng. Chem. Res. 2002, 41, 3884. (8) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production Scheduling of a Large-Scale Industrial Batch Plant I. Short-Term and Medium-Term Scheduling. Ind. Eng. Chem. Res. 2006, 45, 8234. (9) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 1. Multipurpose Batch Processes. Ind. Eng. Chem. Res. 1998, 37, 4341. (10) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production Scheduling of a Large-Scale Industrial Batch Plant. II. Reactive Scheduling. Ind. Eng. Chem. Res. 2006, 45, 8253. (11) He, Y.; Hui, C. W. Rule-Evolutionary Approach for Single-Stage Multiproduct Scheduling with Parallel Units. Ind. Eng. Chem. Res. 2006, 45, 4679. (12) He, Y.; Hui, C. W. Automatic Rule Combination Approach for Single Stage Process Scheduling Problems. AIChE J. 2007, 53, 2026. (13) Castro, P. M.; Novais, A. Q. Short-term scheduling of multistage batch plants with unlimited intermediate storage. Ind. Eng. Chem. Res. 2008, 47, 6126. (14) Castro, P. M.; Grossmann, I. E. New Continuous-time MILP model for the short-term scheduling of multistage batch plants. Ind. Eng. Chem. Res. 2005, 44, 9175. (15) Harjunkoski, I.; Grossmann, I. Decomposition Techniques for Multistage Scheduling Problems using Mixed-integer and Constraint Programming Methods. Comput. Chem. Eng. 2002, 26, 1533. (16) Pantelides, C. C. Unified Frameworks for the Optimal Process Planning and Scheduling. In Proceedings of the Second Conference on Foundations of Computer Aided Operations; Cache Publications: New York, 1994; p 253. (17) Liu, Y.; Karimi, I. A. Scheduling multistage, multiproduct batch plants with nonidentical parallel units and unlimited intermediate storage. Chem. Eng. Sci. 2007, 62, 1549. (18) Shaik, M. A.; Floudas, C. A. Unit-specific event-based continuoustime approach for short-term scheduling of batch plants using RTN framework. Comput. Chem. Eng. 2008, 32, 260. (19) Ierapetritou, M. G.; Hene´, T. S.; Floudas, C. A. Effective continuous-time formulation for short-term scheduling, 3 multiple intermediate due dates. Ind. Eng. Chem. Res. 1999, 38, 3446. (20) Castro, P. M.; Grossmann, I. E. An efficient MILP model for the short-term scheduling of single stage batch plants. Comput. Chem. Eng. 2006, 30, 1003. (21) Rosenthal, R. E. GAMSsa user’s guide. GAMS DeVelopment Corp.: Washington, DC, 2008. (22) Castro, P. M.; Grossmann, I. E.; Novais, A. Q. Two New Continuous-time models for the scheduling of multistage batch plants with sequence dependent changeovers. Ind. Eng. Chem. Res. 2006, 45, 6210.

ReceiVed for reView May 6, 2009 ReVised manuscript receiVed September 4, 2009 Accepted October 21, 2009 IE900734X

Optimal Short-Term Scheduling of Large-Scale ... - ACS Publications

Recommend Documents