Ind. Eng. Chem. Res. 2006, 45, 7603-7617
7603
Optimization of Aluminum Smelter Casthouse Operations Pradeep Prasad and Christos T. Maravelias* Department of Chemical and Biological Engineering, UniVersity of Wisconsin, Madison, Wisconsin 53705
Jeffrey Kelly Honeywell Process Solutions, 300 Yorkland BouleVard, Toronto, Ontario, Canada M2J 1S1
This paper presents a mixed-integer linear programming (MILP) model for the scheduling of a multistage process for the production of aluminum casts of different alloys, using parallel furnaces and casters. In contrast to the common approach in multistage models of considering a fixed set of orders being processed sequentially in the stages, the modeling approach in this paper accounts for actual material flows and, thus, provides flexibility with respect to the actual number of batches to be processed for meeting a given demand. This enables the model to handle parallel nonuniform units with variable capacities, where it is difficult to a priori decide on the number of batches that are required to satisfy the orders. The model also features a material balance over the furnace section, to capture processing details. A decomposition scheme that consists of a master problem and a sub-problem is developed and is used in an iterative algorithm to solve medium- to large-sized problems in reasonable computational time. 1. Introduction The casting of aluminum into bulk geometries of ingots, billets, and slabs is the primary source of refined aluminum alloys to secondary manufacturers. Traditionally, the planning of aluminum smelter casthouses is performed using material requirement planning (MRP) systems, which balance the supply of hot-metal feedstock to the casthouses from the pot lines and the demand for alloyed products. These tools are valuable for the long-term planning of orders with resource requirements (such as raw materials, utilities, etc). However, they are not suitable for the shorter-term scheduling of production in casthouses, given their inability to model production details such as equipment availabilities and capacities. This paper presents a mixed-integer linear programming (MILP) model for the shortterm scheduling of such aluminum ingot production facilities. The goal of the scheduling model is to assign tasks to units and sequence them, while respecting the process requirements and availability of appropriate resources. Although the optimization of aluminum and other similar metal-casting facilities has been studied earlier by other researchers,1-10 the problem considered in this paper has not been addressed. Gravel et al.1,2 obtained feasible schedules for an aluminum casting plant, using genetic algorithm and ant colony metaheuristic approaches. Pacciarelli and Pranzo3 studied a steel casting example, using an alternative graph technique. Nicholls4,5 used an MILP approach to calculate equipment throughputs over a given shift after considering equipment availability and resource limitations, but he did not study the detailed sequencing of tasks and changeovers. Moon and Hrymak6 developed an MILP formulation for the various stages during the steel annealing process and also modeled crane requirements. Harjunkoski and Grossmann7 studied scheduling in a steel production plant with multiple stages such as the problem described in this paper, but only one of the stages had two parallel units: the remainder had a single unit. They utilized a decomposition-based strategy that involved four optimization steps wherein the orders were first grouped based on similar * To whom correspondence should be addressed. Tel.: +1-608-2659026. Fax: +1-608-262-5434. E-mail:
[email protected].
processing requirements. Scheduling was performed for each individual group, followed by sequencing of the groups. Naphade et al.8 used an MILP formulation and a heuristic solution procedure to schedule the casting of steel melts into ingots. They used a heuristic procedure to plan the production of ingots on a weekly basis. Tang et al.9 treated the scheduling problem in the steel-making process as a hybrid flowshop problem with associated process constraints, and used a Lagrangian relaxation-based method to solve the integer programming formulation. Schwindt and Trautmann10 studied the scheduling of the aluminum casting process, using parallel casting systems and sequence-dependent changeovers. Their process mainly involved the casting section, and their model did not include holding furnaces and restrictions on hot metal availability. Also, the casters were structurally somewhat different from those considered in the present study. The short-term scheduling of batch processes using mixedinteger programming methods has received considerable attention during the past decade. (General reviews can be obtained in Reklaitis,11 Pinto and Grossmann,12 Shah,13 and Mendez et al.14) Multiproduct batch processes can be classified based on the structure of the production network as multistage or multipurpose. Multistage processes feature single or parallel units in multiple stages, with each order passing through a specific sequence of stages. Multipurpose plants are more general, allowing the splitting, mixing, and recycle of streams. The current process is a multistage process but also exhibits some characteristics of multipurpose processes, such as the need to monitor batch sizes and inventory levels. Modeling approaches for multistage plants usually involve the assigning and sequencing of a fixed set of orders in units and stages. There is no splitting or mixing of batches, and the relative sequencing between different orders is used to enforce timing constraints. The availability of unlimited intermediate storage is commonly assumed. Overall equipment capacities are considered, but it is normally not required to rigorously monitor inventory levels of the states (chemicals). Pinto and Grossmann15 developed a slot-based continuous-time formulation for multistage batch plants, and they used a decomposition approach to solve large-scale problems. Me´ndez and Cerda´16 developed an
10.1021/ie060652h CCC: $33.50 © 2006 American Chemical Society Published on Web 10/03/2006
7604
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
Figure 1. Schematic of the casting process.
MILP formulation for a candy manufacturing plant that has several features, such as assorted products, order due dates, and sequence-dependent changeovers. Gupta and Karimi17 developed an MILP formulation for a multistage batch plant with parallel units. In a recent study, Castro and Grossmann18 evaluated discrete-time, continuous-time, and constraint-programming (CP)-based modeling approaches for the scheduling of multiproduct, multistage batch plants. Multipurpose plants require the monitoring of resources and chemicals (states) over the scheduling horizon and are commonly modeled using the state-task network (STN)19,20 or resource-task network (RTN)21 approaches, which utilize a time grid and decision variables that are indexed by time. Based on the time grid that is used, these formulations can be classified as discrete-, continuous-, or mixed-time approaches. The discrete-time approach19-21 divides the time horizon into uniform time intervals. As an alternative, several researchers have studied continuous-time formulations.22-29 Continuoustime formulations utilize fewer discrete variables and are more flexible, but are also more complex and use big-M constraints for matching time points with tasks. In an attempt to utilize the strengths of the two time representations, Maravelias30 presented a mixed-time representation for multipurpose plants. Floudas and Lin31 and Burkard and Hatzl32 have given reviews that contrast the formulations based on different time representations. Decomposition and heuristic algorithms are often used to address problems of industrial size.33-39 The casting process can produce several products (ingots of different sizes and alloy grades), and has multiple stages with parallel units at each stage. Other process requirements include variable processing times, sequence-independent changeovers, no intermediate storage between the stages, and release and due times for feed availability and demand satisfaction. A description of the casting process is given in Section 2. The modeling challenges posed by the process, and our proposed approach, are discussed in Section 3. Section 4 presents the detailed mathematical model, and Section 5 discusses an example based on the model. In Section 6, we discuss a decomposition-based iterative approach for solving larger systems. Computational results from the iterative algorithm are presented in Section 7, and Section 8 presents the conclusions. 2. Process Description The casting of aluminum into molds is a multistage process with several parallel holding furnaces that are fed with hot metal from a collection of pots or cells. The hot metal from the holding furnaces is fed to casters to be cast into molds of different shapes and sizes. A schematic of the production process inside an aluminum smelter casthouse is shown in Figure 1.
Figure 2. Operations within furnace.
The feed sources are electrolytic cells wherein alumina is reduced to molten aluminum. The hot metal is collected in crucibles, transported, and poured into a holding furnace. In practice, there can be several different pot lines feeding one or more casthouses within the same plant, so the requirement for hot metal can be satisfied from diverse sources, and this must be appropriately managed and coordinated. The tapping of the pots is described as a set-partitioning problem by Ryan.40 Scrap metal can be generated by repetitive cleaning, preventive, and corrective maintenance activities, and it can constitute another source of metal feed to the process. The feed sources often have release and due times. The release time specifies the earliest time that the production can consume the hot metal and the due time specifies the latest time when it is available. The furnace has three sub-operations, which are called filling, holding, and drawing, as illustrated in Figure 2. Filling can occur only when hot metal is available from the source. When filling the furnace, a draw to a downstream caster cannot occur and vice versa. The holding operation provides enough delay between the filling and drawing of a furnace so that alloying materials (such as magnesium) can be added within required proportions, depending on the desired alloy qualities. It also allows for sufficient time for the stock to mix and stabilize inside the furnace, as well as for remelting any scrap additions. Thus, there exists a requirement for some minimum holding time (τhold) after filling the furnace exists. There can also be a need for an upper bound (τhold,max) on the holding time, given the energy requirements for keeping the aluminum molten. In the general case, there can be multiple draws from a single charge to a furnace with a corresponding time delay between the draws. There is a maximum inventory (FCmax) that a furnace can hold. The fill-to-full restriction requires the filling quantity (Bfin) to take the furnace inventory at least above a certain quantity that is called fill-to-full (FTF). When drawing hot metal from a furnace, the quantity drawn (Bfout) cannot be such that the furnace inventory falls below the minimum allowable inventory
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7605
(FCmin). Also, there is a draw-to-empty (DTE) quantity, which does not allow hot-metal fills until the inventory in the furnace is below the DTE amount. A cleaning or repetitive maintenance operation (washing) of the furnace is required whenever there is a changeover of alloys. It dilutes the remaining hot metal in the furnace from the previous alloy and transitions the furnace to the next alloy. The washing operation is performed by filling the furnace with hot metal, in addition to the metal left in the furnace from the previous batch, and then disposing the hot metal via a casting operation in a caster (“washcast”). The casting operation receives hot metal from the furnaces into molds of different shapes and sizes. Only one furnace can feed a caster at a time. Most casters are of the direct chilled type, which implies a certain controlled procedure for cooling and forming the hot metal alloy into logs of various geometries and sizes. The typical cycle time for a caster is 2-4 h. Whenever there is a change in the mold, a down time needs to be scheduled for proper setup and retooling of the unit. Sequence-independent changeover times for casters are typically of the order of 2-8 h. After the caster, there can be a homogenizing or treating operation that is similar to annealing and then a sawing operation, which takes the 4000-7000 mm log lengths and cuts them into more manageable lengths, according to the customer order, thus leading to a cutting stock problem. It is assumed here that there are no capacity bottlenecks with the homogenizing and sawing operations and these can be included within the cycle time of the caster. Scrap or trim loss after sawing is recycled back to the remelt furnaces. The finished ingots are used to satisfy customer demand. The orders (demand) have a specified amount of alloy type and mold type (cross section). They also have release and due times. The release time is the time at which the storage facility becomes available for a certain demand order, and it constitutes the earliest time at which a caster can finish a batch toward that order. The due time is the latest time at which the demand must be satisfied. Often, the washcast product also has an external demand. 3. Proposed Modeling Approach The aluminum casting process has several complicating features that motivate the need for a new and more-flexible modeling approach. These features include the following: (1) A need for material balance. The fact that the furnace can carry over molten metal between batches implies that a material balance is required to account for the hot metal alloy in the furnace. This carryover of molten metal can have implications on the feasibility of schedules when the supply of hot metal is limited by release and due times. Also, during the wash operation, any remaining hot metal from a previous furnace operation becomes incorporated into the “wash” state in the same unit. (2) Assignment of batches to orders. Order quantities can be such that they require multiple caster batches to be satisfied. It could be argued that large orders can be split into several smaller orders; however, this may not be possible to accomplish in any unique way with parallel furnaces and casters of different sizes. (3) Single connectiVity. At any time, each unit can connect to only one unit from a preceding stage. For any operation, the entire quantity in the unit must be transferred from a single preceding unit and also a no intermediate storage (NIS) policy must be enforced between furnaces and casters. (4) Release and due times. Orders, as well as feed supply, have release and due time restrictions, in addition to quantity specifications.
(5) Washcast operation. The implication of washcast operation is that a changeover in the furnace stage impacts the caster stage, because a caster must accept the wash alloy. This is another feature that makes the proposed model different from most scheduling models in the literature. The STN/RTN approach19-21 easily allows the inventory of states to be monitored. However, in the discrete-time approach, the number of binary variables required becomes large as the problem size increases. The continuous-time approach is also inefficient in this case, because of the large number of event points expected for modeling release and due times. Also, the requirement of unique connectivity between units during the transfer of states causes the need for additional binary variables that represent the active connections between units and must be indexed by the units as well as time/event points. Thus, both discrete- and continuous-time representations result in a large number of binary variables. On the other hand, multistage formulations with time slots can be expected to use fewer binary variables, but they do not typically provide flexibility in the number of batches for meeting orders and do not include material balances. Therefore, we have combined the aforementioned modeling approaches and developed a slot-based MILP formulation with a material balance over hot metal quantities in the furnaces. Traditionally, slot-based models monitor orders as they move through the production sequence. However, we have modeled quantities of the states (feed, hot metal in furnaces, and finished casts) instead of orders. The material balance over the hot metal quantities in the furnace allows the carryover of some amount of hot metal from batch to batch. The possibility of multiple drops from a single charge to a furnace is not considered here. The casting operation is modeled as a batch process that receives hot metal from only one furnace at a time, and the washcast operation is modeled by defining a “wash” alloy. Because of the single connectivity requirement previously mentioned, we use binary variables to define the connection between the individual furnace and caster slots. Although casthouses can produce hundreds of different alloy grades from each furnace, the production model presented here assumes that individual grades can be strategically set into groups such that the alloys within the same group have similar processing requirements. In this paper, we refer to the alloy groups simply as alloys. Figure 3 depicts the framework of slots used to model the process. Independent slots are defined for the furnaces and casters, and this is possible because of the single connectivity requirement. The furnace slots are indexed by k, whereas caster slots are indexed by l. The slots can float within the available production time but must follow the numbering sequence in each individual unit. The duration of slots is variable in the furnaces but is fixed in the casters. The wash alloy is defined as a unique state and must be processed once for an alloy changeover to occur. Thus, furnace slots with different alloys should always have a slot between them that is assigned to the wash alloy. A separate wash state and slot is required for modeling changeover in the furnace, because the product of the changeover operation gets cast in the caster as the washcast and can even have external demand. There is no intermediate storage between the furnace and the caster; therefore, if there is a flow of molten metal between a furnace and caster slot, the finish time for the furnace slot is equal to the start time of the caster slot. Changeovers between different mold types in the caster are simply a time lag between the corresponding slots and are, therefore, modeled without a separate slot. Not using a separate slot for changeover is advantageous, because it
7606
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
Figure 3. Schematic depicting the features of the slot formulation used in the model.
potentially makes it easier to estimate the total caster slots required based on the batches of finished products instead of also having to anticipate the number of mold changes. 4. Mathematical Formulation We assume that the following data is known: (1) A set of furnaces ( j ∈ JF) and casters ( j ∈ JC), and their specifications (FC max j , FTFj and DTEj for the furnaces, and max minimum CC min j,c and maximum CC j,c sizes for the casters). (2) Minimum holding time (τhold) in the furnaces, cycle times (τcast) for the casters, and changeover durations in the furnace (τwash) and caster (τch). (3) A set of feed pots (j ∈ JP) with minimum (Φ min j ) and maximum (Φ max j ) quantity specifications, release times (rj), and due times (dj). (4) A set of alloys (s ∈ S) and mold types (c ∈ C). (5) A set of orders (m ∈ M) with the corresponding alloy S(m) ∈ S, mold type C(m) ∈ C, minimum due amount (Demm), release times (Fm), and due times (δm). The goal is to determine the assignment of alloys to furnace slots, mold types to caster slots, the batch sizes, timing of slots, the connectivity between furnace and caster slots, and the connectivity between caster slots and orders. 4.1. Feed Source to Furnace Flow. Each furnace j ∈ JF is allowed to have a different number of available slots K( j). For each furnace j ∈ JF and furnace slot k e K( j), we define a binary assignment variable xj,k,s that denotes the active alloy s ∈ S within the slot. The set of alloys S includes the wash (φ) alloy. Equation 1 enforces that, at most, one alloy can be assigned to a particular slot.
xj,k,s e 1 ∑ s∈S
∀ j ∈ JF, k e K( j)
(1)
Note that it does not force all slots to be assigned an alloy. There can exist a number of sources (feed pots) for hot metal, but only one of them feeds any particular slot in a furnace. This is enforced through eqs 2-4,
xj,k,s ∑ qj′, j,k ) ∑ s∈S
∀ j ∈ JF, k e K( j)
(2)
j′ ∈ JP
Bfinj,k,s ∑ Fpfj′, j,k ) ∑ s∈S
rj′qj′, j,k e Tfsj,k e dj′qj′, j,k + H(1 - qj′, j,k) ∀ j′ ∈ JP, j ∈ JF, k e K( j) (5) Equation 6 ensures that the minimum and maximum amount of hot metal available from a feed source is respected.
Φ j′min e
Fpfj′, j,k e Φ j′max ∑ ∑ j∈JF
∀j′ ∈ JP
(6)
keK(j)
4.2. Furnace. The assignment of alloys to furnace slots is achieved through eq 1. Equation 7 forbids the assignment of different alloys on consecutive slots.
xj,k,s +
∑
xj,k-1,s′ e 1
∀ j ∈ JF, s ∈ S*, 1 < k e K( j)
s′∈S*\{s}
(7)
In eq 7, S* is the set of alloys other than wash (φ). Equation 8 prevents an empty slot from occurring between different alloys.
xj,k,s e ∑xj,k-1,s ∑ s∈S s∈S
∀ j ∈ JF, 1 < k e K( j)
(8)
Therefore, whenever a change in the alloy is required, eqs 7 and 8 together enforce a wash slot (xj,k,φ) between the different alloys. Equation 8 also stacks any superfluous empty slots at the end and after the active slots. In addition, wash operations are not allowed to occupy the first or last slot in every furnace. Equations 9-11 monitor the furnace inventory at the beginning of the slots.
Sfj,k,s ) Sf0j,s + Bfinj,k,s
∀ j ∈ JF, s ∈ S, k ) 1 (9)
Sfj,k,s ) Sfj,k-1,s + Bfinj,k,s - Bfoutj,k-1,s - Bwashj,k,s ∀ j ∈ JF, s ∈ S*, 1 < k e K( j) (10) Sfj,k,φ ) Sfj,k-1,φ + Bfinj,k,φ - Bfoutj,k-1,φ +
∑ Bwashj,k,s
s∈S*
∀ j ∈ JF, k e K( j)
(3)
j′ ∈ JP
Fpfj′, j,k e FC max j qj′, j,k
the feed source j′ to the furnace slot (j,k) and is equal to the input flow (Bfinj,k,s) to the slot. The assignment variable xj,k,s is used to enforce the correct alloy s in variable Bfinj,k,s. The binary variable qj′,j,k is used in eq 5 to ensure that the start time of each furnace slot (Tfsj,k) is within the corresponding release and due times for the availability of feed source j′.
∀ j′ ∈ JP, j ∈ JF, k e K( j) (4)
where qj′,j,k denotes a connection between feed source j′ and furnace slot (j,k). Variable Fpfj′,j,k is the corresponding flow from
∀ j ∈ JF, 1 < k e K( j) (11) Equation 9 is the material balance for the first furnace slot only, where Sf0j,s is the opening inventory. Equation 10 monitors the starting inventory (Sfj,k,s) in all other furnace slots, using the starting inventory during the previous slot, minus the amount withdrawn from the furnace during the previous slot (Bfoutj,k-1,s),
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7607
plus the amount fed to the furnace in the current slot (Bfinj,k,s), minus the amount (Bwashj,k,s) that is to be incorporated into the wash state (if any) in the next slot. The wash operation dilutes the hot metal remaining from the earlier batch (Bwashj,k,s) with more hot metal, so that it can be disposed as a cast through the washcast operation. Equation 11 is the corresponding material balance equation for the wash operation and uses a summation over alloys for Bwashj,k,s, because the wash operation incorporates the metal remaining from any previous alloy as part of the wash state. An additional variable (SfTj,s) is defined through eq 12 to obtain the closing inventory.
SfTj,s ) Sfj,k,s - Bfoutj,k,s
∀ j ∈ JF, s ∈ S, k ) K( j) (12)
The maximum furnace capacity and fill-to-full requirement are expressed by eq 13, which, along with eq 1, also ensures that a furnace cannot hold an inventory of different alloys at the same time.
FTFjxj,k,s e Sfj,k,s e FC max j xj,k,s
∀ j ∈ JF, s ∈ S, k e K( j) (13)
Equations 14 and 15 ensure that the Bwashj,k,s variable is nonzero only when both the binary for the wash operation in the current slot and the binary for the corresponding alloy in the previous slot are activated.
∑ Bwashj,k,s e DTEjxj,k,φ
∀ j ∈ JF, k e K( j) (14)
s∈S*
Bwashj,k,s e DTEjxj,k-1,s
∀j ∈ JF, s ∈ S*, 1 < k e K( j) (15)
Equation 16 ensures that we do not draw more than the amount present in the furnace.
0 e (Sfj,k,s - Bfoutj,k,s) e DTEjxj,k,s ∀ j ∈ JF, s ∈ S, k e K( j) (16) Equation 17 allows the flow to a furnace slot be nonzero only if the corresponding alloy is active in the slot.
Bfinj,k,s e FC max j xj,k,s
∀ j ∈ JF, s ∈ S, k e K( j)
(17)
The start (Tfsj,k) and finish (Tffj,k) times for a furnace slot (j,k) are constrained by eqs 18 and 19. Equation 18 enforces that a slot can start only after the previous slot finishes and eq 19 enforces the minimum duration of a furnace slot, depending on the alloy that is active.
Tfsj,k g Tffj,k-1 Tffj,k g Tfsj,k + τ hold j
∀ j ∈ JF, 1 < k e K( j)
(18)
xj,k,s + τwashxj,k,φ ∑ s∈S*
∀ j ∈ JF, k e K( j) (19)
4.3. Furnace to Caster Flow. The assignment in the caster is made by the binary variable wj,l,c, which takes a value of 1 when the mold type c ∈ C is assigned to the slot l e L(j) of caster j ∈ JC, where L(j) is the number of available slots for caster j. The assignment of mold types to caster slots is enforced by eq 20.
∑wj,l,c e 1
c∈C
∀ j ∈ JC, l e L( j)
(20)
The process description requires that each batch in the caster is derived completely from one furnace; i.e., whenever a connection is made between a furnace and a caster slot, it should be exclusive. This is modeled using a binary variable uj,k,j′,l that takes a value of 1 if there is a flow of hot metal from the furnace slot (j,k) to the caster slot (j′,l). As seen in eqs 21 and 22, a connection is feasible only when corresponding assignment binaries in the furnace and caster slots are 1.
∑ ∑
uj,k, j′,l )
j∈JF k e K(j)
∑wj′, l,c
∀ j′ ∈ JC, l e L(j′)
(21)
c∈C
xj,k,s ∑ ∑ uj,k,j′,l ) ∑ s∈S
∀ j ∈ JF, k e K( j) (22)
j′∈JC leL( j′)
Equation 23 then ensures flow between a furnace slot and caster slot only if the connecting variable (uj,k,j′,l) is 1. max Ffcs,j,k,c, j′,l e max{CC c,j′ }uj,k, j′,l ∑ ∑ c∈C s∈S c∈C
∀ j ∈ JF, j′ ∈ JC, k e K( j), l e L( j′) (23)
To represent the problem with fewer discrete variables, we have chosen to represent the attributes of a caster slot only with the mold type. To monitor the specific alloy that is active within the selected mold type, we use the continuous variable Ffcs,j,k,c,j′,l that represents the flow of alloy s from furnace slot (j,k) to the caster slot (j′,l) with mold type c. Equations 24 and 25 ensure that the aforementioned flow equals the size of the active mold type (Bcj′,l,c) and that the correct amount of the corresponding alloy is drawn from the furnace.
∑ ∑ ∑ s∈S j∈JF
Ffcs, j,k,c, j′,l ) Bcj′,l,c
k e K(j)
∑∑ ∑
∀ j′ ∈ JC, c ∈ C, l e L( j′) (24) Ffcs, j,k,c, j′,l ) Bfoutj,k,s
c∈C j′∈JC l e L(j′)
∀ j ∈ JF, s ∈ S, k e K( j) (25)
We have assumed that the time required to transfer metal from a furnace to a caster is included in the cycle time. Thus, we ensure through the variable uj,k,j′,l in the big-M relations (eq 26),
Tcsj′,l - H(1 - uj,k, j′,l) e Tffj,k e Tcsj′,l + H(1 - uj,k, j′,l) ∀ j ∈ JF, j′ ∈ JC, k e K( j), l e L( j′) (26) that the finish time of the furnace slot (Tffj,k) coincides with the start time of a caster slot (Tcsj′,l) if they are “connected”. 4.4. Caster. The assignment constraint in eq 20 is supported by eqs 27 and 28. Equation 27 stacks all the empty slots at the end, whereas eq 28 activates the changeover operation in the caster through the binary variable wˆ j,l if the mold types in consecutive slots are different.
∑wj,l,c e c∈C ∑wj,l-1,c
∀ j ∈ JC, 1 < l e L( j)
(27)
c∈C
wˆ j,l g wj,l,c +
∑
wj,l-1,c′ - 1
c′∈C\{c}
∀ j ∈ JC, c ∈ C, 1 < l e L( j) (28) Unlike the furnace, separate slots are not used for handling changeovers in casters, and the changeover binary for a caster slot (j,l) can be equal to 1, in addition to the assigned mold type binary. Changeovers are not allowed in the first slot of
7608
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
each caster. The batch size in a caster slot (Bcj,l,c) is dependent on the type of the mold that is active. For each mold type and caster, the batch size is between a minimum and maximum amount as given by eq 29. The changeover duration is added before the start of the slot with the active wˆ j,l variable through eq 30. Equation 31 enforces the cycle time (τcast) in each active caster slot. In addition to these constraints, the start time for caster slots (Tcsj,l) is bounded below by the smallest holding time among the furnaces. max CC min c,j wj,l,c e Bcj,l,c e CC c,j wj,l,c ∀j ∈ JC, c ∈ C, l e L( j) (29)
∑wj,l,c c∈C
∀ j ∈ JC, l e L( j)
(31)
4.5. Demand Satisfaction. Each order m ∈ M is characterized by the alloy S(m), mold type C(m), amount (Demm), the release time (rm) and due time (dm). Production from each caster slot is used to satisfy only one order. A large order can be satisfied using multiple caster slots. The binary variable zj,l,m has a value of 1 if a batch in caster slot (j,l) is used to satisfy order m. As given by eq 32, satisfaction of order m can occur only from a caster slot (j,l) processing the correct mold type.
∀ j ∈ JC, m ∈ M, c ) C(m), l e L( j) (32)
zj,l,m e wj,l,c
The variable Ams,c,j,l,m gives the amount of the order m for a product of alloy s and mold type c, met from the production in caster slot (j,l). It can take a nonzero value only when zj,l,m is 1 in eq 33, where the yield factor ηj,c adjusts for metal losses during casting and finishing. max )zj,l,m Ams,c, j,l,m e max(ηj,c′CC j,c′ c′∈C
∀ j ∈ JC, m ∈ M, s ) S(m), c ) C(m), l e L( j) (33) The correct matching of alloy and mold type combination between a caster slot and an order is enforced in eq 34 through the flow variable for alloys between the furnace and the caster (Ffcs,j,k,c,j′,l),
∑
∑ ∑ zj,l,m g 1
Ams,c,j′,l,m ) ηj,c
∑ ∑ Ffcs,j,k,c,j′,l
∀ m ∈ M: Demm > 0
(37)
j∈JC leL(j)
∀ j ∈ JC, 1 < l e L( j) (30)
Tcsj,l g Tcfj,l-1 + τ w ˆ j,l ch
Tcfj,l ) Tcsj,l + τcast
The positive slack variables s1j,l,m and s2j,l,m allow the caster slots to violate the release and due times for orders, preventing the model from becoming infeasible, because of tight deadlines. The slack variables are minimized in the objective function. 4.6. Additional Constraints. Some additional constraints were determined to be useful for finding good solutions quickly. Equation 37 states that each order that has a nonzero minimum demand amount must have at least one connection to a caster slot.
Equation 38 states that each active caster slot processing mold type c should connect to an order that can accept the same mold type.
wj,l,c )
∑
∀ j ∈ JC, c ∈ C, l e L( j) (38)
zj,l,m
m∈MC(c)
Equation 39 activates the wash binary if the alloy changeover occurs in a furnace, and it is obtained via the addition of two changeover constraints of the form xj,k+1,φ g xj,k,s + xj,k+2,s′ - 1.
xj,k+1,s + xj,k+2,φ g xj,k,s + xj,k+1,s +
∑
(xj,k+2,s′ + xj,k+3,s′) - 2
s′∈S*\{s}
∀ j ∈ JF, s ∈ S*, 1 < k < K( j) - 2 (39) Equation 40 provides a lower bound to the flow from a furnace slot to a caster slot, if there is a connection between them.
Ffcs,j,k,c,j′,l g min(CC min ∑ ∑ j,c )uj,k,j′,l c∈C s∈S c∈C
∀ j ∈ JF, j′ ∈ JC, k e K(j), l e L( j′) (40)
4.7. Objective Function. The objective function (eq 41) used in our investigations minimizes the number of changeovers required in the furnaces and casters, and the sum of positive slack variables (earliness/lateness).
j∈JF keK(j)
m∈MS(s)∩MC(c)
∀ j′ ∈ JC, s ∈ S, c ∈ C, l e L( j′) (34) where MS(s) and MC(c) are sets of orders that can accept alloy s and mold type c respectively. Note that often washcast itself is a product with demand. Equation 35 ensures that we meet the orders, by specifying that the sum of Ams,c,j,l,m over all caster slots should be greater than the minimum demand amount for order m.
∑ ∑ Ams,c,j,l,m g Demm j∈JC leL(j)
∀ m ∈ M, s ) S(m), c ) C(m): Demm > 0 (35) The variable Ams,c,j,l,m is fixed to zero for all s ∉ S(m) and c ∉ C(m). Demand satisfaction for a product occurs when the caster slot finishes and must respect the release and due times for the orders through eq 36.
Fmzj,l,m - s1j,l,m e Tcfj,l e δmzj,l,m + H(1 - zj,l,m) + s2j,l,m ∀ j ∈ JC, m ∈ M, l e L( j) (36)
k1
j∈JF k>1
ω 2(
∑ ∑ ∑ (s1j,l,m + s2j,l,m))}
(41)
j∈JC leL(j) m∈M
qj′′,j,k, xj,k,s, uj,k,j′,l, wj′,l,c, wˆ j′,l, zj′,l,m ∈ { 0,1} Ams,c,j′,l,m, Bcj′,l,c, Bfinj,k,s, Bfoutj,k,s, Bwashj,k,s, Fpfj′′,j,k, Ffcs,j,k,c,j′,l, Sfj,k,s, Tcsj′,l, Tfsj,k, Tcfj′,l, Tffj,k g 0 (42) The coefficients ω1 and ω2 are dependent on the specific problem. The value of the objective function can only change in steps equal to the greatest common factor (GCF) of the individual terms. In this case, the objective function consists of both integer and continuous variables, of which the continuous slack variables are dependent on the time-related parameter values. If the GCF of the time-related parameters (Fj, δj, rj, dj, τcast, τch, τwash, τjhold) is γ, the smallest possible change in the
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7609
objective function will be the GCF of ω1 and γω2. This fact will be used later to help accelerate the search for the optimal solution. Equations 1-42 constitute the MILP model, which, henceforth, is referenced as the full model (FS). 4.8. Remarks. Additional features can easily be included to extend the model. Sometimes, a maximum holding time is required in the furnace to limit energy consumption when the furnaces are heated to keep the metals molten. Forbidden sequences of tasks in a unitsor forbidden connections between unitsscan be incorporated. It is also possible to model sequence-dependent changeovers instead of sequence-independent changeovers. The selection of the number of slots is a tricky but important issue in all slot-based models, because it is imperative to select the correct number of slots to obtain a good solution; however, at the same time, defining too many slots has an adverse effect on computational performance. The common approach to determining the number of slots starts with an estimation of the number of batches required to meet the demand. In this paper, the estimated slots are divided equally between the units, because the problem data does not include unit-specific costs or processing times, and, hence, there is no anticipated reason for assignments to favor any particular unit. The heuristic that is followed in this paper is to divide the estimated number of batches by the number of parallel units, round upward, and then add more slots to account for changeovers, depending on the difference between the number of units in a stage and the number of possible modes (e.g., furnaces and alloys). 5. Example An example (P1) is used here to illustrate the features of the model (data are presented in Appendix A). The example consists of a set of seven orders to be satisfied in a 30-h production horizon, using two furnaces, two casters, and two sources for hot metal. The furnaces have different sizes, but they have the same holding time requirements and can process three different alloys. The casters also have different sizes and can process three different mold types. As seen from the data, there are some incompatible furnace-to-caster connections, because of the size requirements. For example, Furnace 1 cannot feed Caster 1 if it is processing mold type “C3”. Note that there are two big orders that cannot be satisfied with just one caster batch: e.g., order 5 can be satisfied with either two large batches or three small batches, depending on the specific furnaces and casters used. The weight coefficients in the objective function were considered as follows: ω1 ) 10 and ω2 ) 100. The MILP model (FS) was solved for these data, using GAMS 22.0/CPLEX 10.0 on a personal computer (PC) that was using the Windows XP operating system with a 2.8 GHz processor and 512 MB RAM. The CPLEX options that were applied include integer solution feasibility emphasis and objdif. The nature of the objective function allows us to utilize the “objdif” feature in CPLEX effectively, where the solver searches for nodes with solutions that are better by the objdif value than the current solution. Because the problem data related to time have a greatest common factor of 1, the slack term in any integer solution will also take integer values. Therefore, the smallest possible change in the objective value between integer solutions (see section 4.7) is simply the GCF of ω1 and ω2, which is equal to 10. Because the objdif value used is 9.9999, the use of objdif will not cause nodes with potentially optimal solutions to be pruned in the branch-and-bound tree. Six slots were defined for each
Figure 4. Gantt chart for example problem P1.
furnace and caster in this example. The minimum number of batches required is nine, so there may be some superfluous slots in each unit. The model has 378 binary variables, 2385 continuous variables, and 2029 equations and is solved to optimality within ∼3 h. Figure 4 presents the resulting Gantt chart and shows many of the features that have been discussed previously, such as changeovers, washcast, and batches merging toward an order. The different shades of gray represent alloy grades, and the patterns represent different mold types. Furnace 1 processes six batches of the same alloy (A1). Furnace 2 starts with two batches of alloy A2 and then changes over to alloy A3. Slot 3 in Furnace 2 performs the wash operation, to allow the furnace to change over to a new alloy. The second and fifth slots in Furnace 1, and the fourth slot in Furnace 2, exhibit cases where the processing time is greater than the minimum of 3 h. The finish time of each furnace slot coincides with the start time of a caster slot. Caster 1 only processes mold type C2 for different alloys. Slot 3 in caster 2 is used for the washcast operation. Caster 2 starts with four batches of mold type C1 and then changes over to mold type C3. Accordingly, there is a changeover time imposed before the fifth slot in Caster 2. All the caster slots are of same duration (3 h). The time window within which demand satisfaction must occur for each order is shown in the figure. The dark triangles represent the actual time at which demand satisfaction occurs, and it is equal to the time at which a caster slot of the corresponding type finishes operation. All the orders are satisfied within the release and due times in this case. As shown in Figure 4, some of the orders were satisfied at two or more points (i.e., with more than one caster batch). Order 5 is satisfied here using three batches in Caster 1. All the dummy orders for washcast (orders 8-10) have a minimum demand of zero. Hence, there is only one washcast product here and it goes toward order 8. Figure 5 shows the profile of the inventory of molten metal in Furnace 2. It shows that the alloy is filled to the FTF level for each batch. It also shows molten metal being carried over between the second and third batches. The furnace is forced to be emptied after the third batch, because it is used for washing.
7610
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
xj,k,s )
∑ ∑
∀ j ∈ JF, k e KM( j), s ∈ S
Vj,k, j′,m
j′∈JC m∈MS(s)
(44)
∑ ∑ ∑ Vj,k, j′,m e LM( j′)
∀ j′ ∈ JC
(45)
j∈JF keKM(j) m∈M
Equation 46 ensures that the material flow (Am′j,k,j′,m) is zero if the corresponding variable Vj,k,j′,m is zero.
Am′j,k, j′,m e ηj′,cCC c,max j′ Vj,k, j′,m ∀ j ∈ JF, j′ ∈ JC, m ∈ M, c ) C(m), k e KM( j) (46) Figure 5. Furnace inventory profile for Furnace 2.
6. Decomposition Approach The FS model is computationally limited, and, therefore, in this section, we develop a decomposition approach that enables us to solve problem sizes of industrial importance effectively. The problem is decomposed into two smaller problems: (1) A master problem (MP) is solved to obtain the assignment, sequencing, and sizing of batches in the furnaces, the assignment of furnace batches to orders, and the routing of these batches via casters and toward orders. (2) Given a master solution, the sub-problem (SP) is solved to obtain a feasible solution. The SP is obtained using the FS model by incorporating the information available from the current solution of the MP. The two models are combined in an iterative scheme, where the master problem provides candidate (partial) solutions and an upper bound and the subproblem provides (if feasible) a detailed schedule and a lower bound. The algorithm terminates when the two bounds converge. The MP model, the SP model, and the iterative algorithm are described in sections 6.1, 6.2, and 6.3, respectively. 6.1. Master Problem (MP). The master problem provides a partial solution for the full problem by relaxing the restrictions on feed availability and approximating the furnace-to-caster connections as a simple time delay, thus implying ready availability of a suitable caster as soon as each furnace batch finishes. Variables xj,k,s, Tfsj,k, and Tffj,k and the constraints described by eqs 1, 7, 8, 18, 19, and 39 are common between the FS and MP formulations. The binary variable Vj,k,j′,m is introduced to model furnace f caster f order routings, and it takes a value of 1 if order m is satisfied via furnace slot (j,k) and caster j′. The variable Am′j,k,j′,m expresses the material flow for the corresponding Vj,k,j′,m. Binary variable yj′,c tracks the mold types that are active in each caster, and it takes a value of 1 if mold type c is processed at least once in the caster j′. The variable moldsj provides the number of mold types active in caster j. A single batch size (Bfj,k,s) is used for the furnace slots, instead of the incoming flow (Bfinj,k,s) and outgoing flow (Bfoutj,k,s) variables. The data for the master problem are the same as those of the FS model. Equations 1, 7, 8, and 43 are used for the assignment, sequencing, and sizing of batches in the furnaces, where KM(j) is the number of allowed slots in furnace j in the MP.
Bfj,k,s e FC max j xj,k,s
∀ j ∈ JF, k e KM( j), s ∈ S (43)
Equation 44 assigns each active furnace slot to one order via one caster, and eq 45 ensures that the number of orders passing through a caster does not exceed the number of slots allowed (LM(j)).
Equation 47 matches the batch size of each furnace slot (Bfj,k,s) to the corresponding material flow toward demand satisfaction, while accounting for the correct yield factor.
∑ ∑ ∑ (Am′j,k, j′,m/ηj′,c) ) Bfj,k,s
m∈MS(s) j′∈JC c∈C(m)
∀ j ∈ JF, k e KM( j), s ∈ S (47) Equation 48 is a demand satisfaction constraint, and eq 49 ensures that each order is met between the release and due times, with violations allowed through the variables s1′j,k,m and s2′j,k,m.
∑ ∑ ∑ Am′j,k, j′,m g Demm
∀m∈M
(48)
j∈JF keKM( j) j′∈JC
∑ Vj,k, j′,m - s2′j,k,m e Tffj,k e (δm - τ ) ∑ Vj,k, j′,m + H(1 - ∑ Vj,k, j′,m) + s1′j,k,m j′∈JC j′∈JC
(Fm - τcast)
j′∈JC cast
∀ j ∈ JF, k e KM( j), m ∈ M (49)
Equation 50 “activates” a certain minimum number of connections (Vj,k,j′,m) toward each order, based on the minimum demand amount:
∑ ∑ ∑ ηj′,cCC c,maxj′ Vj,k, j′,m g Demm
j∈JF keKM(j) j′∈JC
∀m ∈ M, c ) C(m): Demm > 0 (50) although it is not necessary, it enhances the solution of the MP. The furnace slot start and finish times are given by eqs 18 and 19, as described previously in the FS model. In addition, eq 39 is used to support the assignment constraints in activation of the wash binary. A conservative estimate of the number of changeovers required in caster j′ can be obtained as (moldsj′ - 1), where the number of mold types active in caster j′ (moldsj′) is obtained from eqs 51-53.
∑ ∑
∑
∀ j′ ∈ JC, c ∈ C
Vj,k, j′,m e LM( j′)yj′,c
j∈JF keKM( j) m∈MC(c)
(51)
moldsj g
∑yj,c
∀ j ∈ JC
(52)
c∈C
moldsj g 1
∀ j ∈ JC
(53)
An estimate of the makespan of the furnaces in terms of the slots used is calculated in eq 54 via variable MS.
MSF g
∑ ∑xj,k,s
keKM( j) s∈S
∀ j ∈ JF
(54)
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7611
The objective function (eq 55) minimizes the violation of release and due times, the number of changeovers in the furnaces, the estimated number of changeovers in the casters, and the estimated makespan.
obj′ ) min(ω1
∑ ∑
xj,k,φ +
j∈JF keKM(j)
ω2
∑ ∑ ∑
(s1′j,k,m + s2′j,k,m) + ω3MSF + ω4
j∈JF keKM(j) m∈M
∑
j∈JC
K(j) ) moldsj)
∀i
(56)
keKM(j)
where i i-1 R j,k, j′,m ) V j,k, j′,m
∀i
and
B ) i
∑ ∑ ∑∑
i-1 V j,k, j′,m
-1
∀i
MP x j,k,s ∑ ∑ s∈S
L(j′) )
MP V j,k, ∑ ∑ ∑ j′,m j∈JF m∈M
x MP j,k,s
) 0 w Sfj,k,s ) Bfinj,k,s ) Bfoutj,k,s )
∑ ∑ ∑ Ffcs, j,k,c, j,l ) 0 c∈C j′∈JC leL(j′)
∀ j ∈ JF, k ∈ K( j), s ∈ S
(57)
∀ j′ ∈ JC
(58)
keKM( j)
All discrete and continuous variables are fixed to zero for indices that are greater than K(j) and L(j). All constraints and all summations in the FS model are now written with the updated values for K(j) or L(j). Knowing the exact number of active caster slots in the SP, eq 20 is written as an equality and eq 27 is dropped. 6.2.3. Allowed Furnace-Caster Connections. Based on the MP variable V j,k, j′,m, a new set is introduced to define allowed connections between individual furnace slots and casters. The set FC(j,k,j′) is the set of all active connections between furnace slots (j,k) and caster j′. Note that this set defines connections to each caster unit and not to individual caster slots. MP V j,k, ∑ j′,m ) 1 w ( j,k, j′) ∈ FC( j,k, j′) m∈M
∀ j ∈ JF, j′ ∈ JC, k e K( j) (59)
j∈JF keKM( j) j′∈JC m∈M i-1 Here, V j,k, j′m is the solution of the model MP in the (i - 1)th iteration. Any iteration i of the MP contains (i - 1) integer cuts, based on the solution from previous iterations. Note that the integer cuts force the MP to give new integer solutions for Vj,k,j′,m, but the furnace assignments (xj,k,s) are allowed to be the same between iterations. This is because each furnace assignment can lead to several solutions for Vj,k,j′,m, some of which are better than others. Also note that the integer cut does not include the variables that had the value zero in the previous iteration. This prevents a “different” integer solution from being constituted simply by a similar solution as the old one but with the orders being met using a larger number of batches. Any particular iteration of the master problem consists of eqs 1, 7, 8, 18, 19, 39, and 43-56. Additional constraints are included in the MP as part of the iterative algorithm and are explained in section 6.3. 6.2. Sub-problem (SP). The MP is used to obtain the MP assignment of alloys to furnaces (x j,k,s ), the routing of alloy MP batches through the casters toward order satisfaction (V j,k, j′m), MP and the mold types active in each caster (y j′,c ). This information is used to simplify the FS model and solve it as the subproblem (SP). The solution of the SP is thus constrained by information from the MP. An outline of the simplifications performed on the FS in each iteration is provided below. Further details including the modified constraints are provided in Appendix B. 6.2.1. Fix Furnace Slot Assignments. The furnace slot assignments in the sub-problem (SP) are fixed to the MP values. The fixing of the furnace assignments also allows eqs 1, 7, and 8 to be dropped in the SP. Furthermore, all continuous variables that refer to the alloys that are not assigned to a particular slot are fixed to zero for that slot, e.g.,
∀ j ∈ JF
keKM( j)
(55)
Note that we can apply the same reasoning as in section 4.7 to prove that the smallest possible change in the objective function is the GCF of ω1, γω2, ω3, and ω4. Equation 56 defines an integer cut that prevents the MP from repeating integer solutions during the iterative process. i i R j,k, ∑ ∑ ∑ ∑ j′,mVj,k, j′,m e B j∈JF j′∈JC m∈M
6.2.2. Update Number of Active Slots. Because the MP solution is known, it is possible to use the exact number of active furnace slots K(j) in the SP. Similarly, the number of active caster slots L(j) can be updated based on the MP variable MP V j,k, j′m:
A connection between a furnace slot and a caster is allowed in the SP only if such a connection exists in the MP. Therefore, the variables for furnace-caster connections are simplified by setting uj,k, j′,l ) 0 and Ffcs, j,k,c, j′,l ) 0 if (j,k, j′) ∉ FC(j,k, j′). Equations 21, 22, 24, and 25 are tightened by expressing the summations on the left-hand side over the set FC(j,k, j′). Also, eqs 23, 26, and 40 are expressed only for the indices (j,k, j′) ∈ FC(j,k, j′). 6.2.4. Allowed Caster-Order Connections. The MP variMP able V j,k, j′m is also used to define CO(j′,m), which is the set of allowed caster (j′) to order (m) connections. Note again that the set refers to casters and not individual slots.
∑ ∑
MP V j,k, j′,m ) 1 w ( j′,m) ∈ CO( j′,m)
j∈JF keKM( j)
∀ j′ ∈ JC, m ∈ M (60)
Variables that denote inactive caster-order connections are fixed by setting zj′,l,m ) 0 and Am′j,k, j′,m ) 0 if (j′,m) ∉ CO(j′,m). The summations on the left-hand side of eqs 34 and 35 are expressed only over the set of allowed caster-order connections CO(j′,m). Furthermore, the order satisfaction constraints (eqs 33 and 36) are only expressed for (j,m) ∈ CO(j,m). 6.2.5. Additional Constraints. New constraints are added to exploit information available from the MP. Equation 61 states that a connection between furnace slot (j,k) and caster j′ in the MP implies such a connection to the caster (one of the slots) in the SP.
∑ uj,k, j′,l ) 1
leL( j′)
∀ j ∈ JF, j′ ∈ JC, k e K( j): (j,k,j′) ∈ FC( j,k,j′) (61) Equation 62 states that the number of caster-to-order connections in the MP is preserved in the SP. MP ∑ zj′,l,m ) j∈JF ∑ ∑ V j,k, j′,m
leL( j′)
keK( j)
∀ j ∈ JC, m ∈ M (62)
7612
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
Figure 6. Iterative algorithm flowchart.
Equation 63 states that the number of caster-order connections in the MP, for each individual mold type, should be also preserved in the SP.
∑ wj′,l,c ) j∈JF ∑ ∑ ∑
leL( j′)
MP V j,k, j′,m
∀ j ∈ JC, c ∈ C
keK( j) m∈MC(c)
∀ j′ ∈ JC, c ∈ C, l e L( j′)
(64)
6.3. Iterative Algorithm. We use an iterative approach, where the MP is used to generate a succession of candidate partial solutions for the SP (see Figure 6). Note that K(j), L(j), FC(j,k,j′), and CO(j,m), as defined in section 6.2, change between iterations, whereas KM(j) and LM(j) remain unchanged. The algorithm stops either when an MP becomes integerinfeasible or when the termination criterion (such as the number of iterations, the resource limit, the optimality gap, etc.) is satisfied. An infeasible MP indicates that no new integer solution exists that is different from the previous iterations while also meeting the bounds imposed; in other words, the last feasible SP solution cannot be improved any further. The last feasible SP solution is also optimal if all SPs are solved to optimality. Additional constraints are added to the MP and SP to improve computational efficiency by reducing the search space. 6.3.1. Additions to the MP. Equation 65 defines an upper bound on the number of estimated changeovers allowed in the master problem.
∑ ∑
xj,k,φ +
j∈JF keKM( j)
∑ (moldsj′ - 1) e chngs - 1
j′∈JC
i′ i′ x j,k,φ + ∑ ∑ w ˆ j′,l ∑ ∑ j∈JF j′∈JC keK( j)
(65)
(66)
The value of LB is updated every time a master problem is solved to optimality. 6.3.2. Additions to the SP. Equation 67 provides an upper bound on the slack variables, based on the best available SP solution.
∑ ∑ ∑ (s1j,l,m + s2j,l,m) e slack
(67)
j∈JC leL( j) m∈M
where
slack )
i′ i′ (s1 j,l,m + s2 j,l,m ) ∑ ∑ ∑ j∈JC m∈M leL( j)
The parameter “slack” during any iteration is obtained from the sum of the slack variables in the previous feasible SP (in iteration i′). Equation 68 provides an upper bound on the number of changeovers in the SP solution.
∑ ∑ xj,k,φ + j′∈JC ∑ ∑ wˆ j′,l e chngs - 1
j∈JF keK( j)
(68)
leL( j′)
It is activated only after an SP solution with zero slack variables is determined. The parameter chngs during any iteration is obtained from the total number of changeovers in the previous feasible SP. 7. Computational Results
where
chngs )
obj′ g LB
(63)
Equation 64 states that a mold type c can be active in caster j′ in the SP, only if the same was true in the MP. MP wj′,l,c e y j′,c
The parameter “chngs” during any iteration i is obtained from the total number of changeovers in the previous feasible SP (in iteration i′). Equation 66 provides a lower bound (LB) to the objective function.
leL( j′)
The decomposition approach outlined in section 6 was applied to three problems: P1, P2, and P3. The problem data are provided in Appendix A. Problem P1 is the illustrative problem presented in section 5, whereas problems P2 and P3 may be
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7613 Table 1. Model Statistics Problem P1 binary variables continuous variables constraints
Problem P2
FS
MP
SP
FS
MP
SP
FS
MP
SP
378 2385 2029
290 556 746
158 545 1028
578 3834 2973
420 814 1043
212 742 1325
1236 10509 6153
1660 2478 2934
416 1683 2758
Table 2. Full Mixed-Integer Programming (MIP) Model Results
a
Problem P3
problem
best solution
% gap
time (s)
P1 P2 P3
20 20 NSa
0 50 NSa
10 246.79 20 000.00 20 000.00
No solution
Table 3. Decomposition Algorithm Results Computational Requirements (CPU s) problem
solution
iterations
total
MP
SP
last MP
P1 P2 P3
20 20 20
2 6 18
2.78 78.71 8160.78
2.54 29.06 6772.21
0.25 49.65 1388.57
0.91 7.11 2744.19
characterized as medium-sized problems. The iterative algorithm was implemented on the same computer and modeling software as that described in section 5 for the FS model. Integer solution feasibility emphasis and an objdif of 9.9999 were used as algorithmic options for the MP as well as the SP. In all of the examples, the weight coefficients in the objective function were
Figure 7. Gantt chart for example problem P3 using the iterative algorithm.
considered as follows: ω1 ) ω3 ) ω4 ) 10 and ω2 ) 100. The model statistics are summarized in Table 1. The model statistics for the MP and SP change with the iterations and, therefore, are shown here for the first iteration only. Table 2 shows the computational results with the FS model. Note that only problem P1 is solved to optimality within the resource limit of 20 000 s, and that no integer solution is determined for problem P3. Table 3 shows that the iterative scheme shows significant improvement in computational time over the FS model. Furthermore, all the solutions shown in Table 3 are optimal solutions. To prove optimality, it is required that all the SPs are solved to optimality and the algorithm stops when an MP is proven infeasible. Table 3 shows that a major portion of the solution time in the iterative process is spent over the MPs, which were also solved to optimality in this case. However, strict optimality is often not the requirement in practice and a good feasible solution is sufficient. In such cases, the decomposition algorithm can still be used to generate good solutions relatively quickly. For problem P3, the algorithm finds a solution with no violation of due times in 1570 s; however,
7614
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
this solution contains one changeover more than the optimal solution. The optimal solution itself is subsequently attained in 5417 s, but an additional 2744 s are required to prove optimality (by proving the final MP infeasible). The single largest time requirement among the iterations is required for proving the final MP infeasible. Therefore, the decomposition algorithm is effective in solving medium-sized problems to optimality. To solve scheduling problems for large casthouses (that often consist of ∼10 furnaces, 5 casters, and more orders), we can use the proposed algorithm as a part of a broader decomposition strategy that simplifies the large problem. The decomposition can be spatial, where a high-level planning model assigns orders to casthouses, or temporal, where orders over months are assigned to small weekly scheduling problems. Also, not all furnaces and casters are bottlenecks and often a high impact on productivity can be made just by scheduling the units that are more tightly constrained. Finally, as mentioned previously, the goal often is not to find optimal solutions but just good feasible solutions. In such cases, the algorithm can be very effective for larger problems by using stopping criteria other than optimality. For example, it was determined that, in a majority of the iterations, the MP reaches the optimal solution relatively quickly but takes a longer time to prove optimality. Hence, we could use the MP with a smaller resource limit instead of solving to optimality. Figure 7 shows the Gantt chart for an optimal solution to problem P3, which consists of 4 furnaces, 4 casters, 3 feed pots, and 16 orders for 4 alloys and 3 mold types. Each furnace and caster was allowed a maximum of 6 slots in the MP. The time window (between release and due times) within which the 3 feed pots can be used are shown at the top of the figure. The dark triangles on the feed lines are the times at which the feed is withdrawn to start a furnace operation. The numbers in parentheses are the numbers of such withdrawals at the same time. Figure 7 also indicates the release and due times for each order. The dark inverted triangles on the order lines show the times at which each order is satisfied. Several orders require multiple batches to be satisfied. Orders 14-16 are dummy orders with zero minimum demand for the washcast product. They are used to facilitate the modeling of changeovers in the furnace and, therefore, are not required to be satisfied. The changeover in Furnace 2 leads to one batch of the washcast product that goes toward order 16. 8. Conclusions In this paper, we present mixed-integer programming methods for the optimization of aluminum smelter casthouse operations. The casting process is rich with several features such as multiple stages with parallel units, variable unit capacities, variable processing times, and unique connectivity between units at any given time. The common approach in scheduling models for multistage plants is to assign and sequence batches (corresponding to orders) through the stages and units. However, the presence of nonuniform parallel units with variable capacity ranges in this case precludes a unique breakdown of the orders into a specific number of batches. To address this problem, a novel slot-based mixed-integer linear programming (MILP) model is developed that tracks material flows through the units, thus providing flexibility over the actual number of batches processed for meeting the demand. Also, to model the processing requirement in the furnaces, material balances for the furnace stage were included in the model. The connectivity requirement between furnaces and casters is modeled by introducing binary variables for the connectivity between the furnace and caster slots.
To address instances of industrial importance effectively, we developed an iterative decomposition scheme that consists of a master problem (MP) that provides an assignment of batches to orders and a solution for the furnace section of the process, followed by a sub-problem (SP) that provides a feasible solution for the entire problem. The solution of the MP provides an upper bound and is used to simplify the full model and obtain the SP, which provides a lower bound. The MP and SP are combined in an iterative algorithm, and the decomposition method is shown to solve medium- to large-sized problems to optimality in a reasonable time. Acknowledgment P.P. and C.T.M. would like to acknowledge financial support from the Graduate School and the College of Engineering of the University of Wisconsin-Madison. Appendix A. Example Data Table A1 gives process data for example problems, labeled P1, P2, and P3, whereas Table A2 gives process data that are common for all of the example problems. Table A3 describes the feed-pot availability. Table A1. Process Data for Example Problems problem
P1
number of furnaces number of casters number of feedpots number of alloys number of mold types number of ordersa horizon (h) furnace capacity (ton) Fur1 (FCMax/FTF/DTE) Fur2 (FCMax/FTF/DTE) caster capacity (ton) Cast1 C1 (CCMin/CCMax) C2 (CCMin/CCMax) C3 (CCMin/CCMax) Cast2 C1 (CCMin/CCMax) C2 (CCMin/CCMax) C3 (CCMin/CCMax) furnace slots caster slots a
P2
P3
2 2 2 3 3 10 30
3 2 2 3 3 12 30
4 4 3 4 3 16 36
55/35/10 70/40/10
70/30/5
70/30/5
35/45 35/45 60/70
25/40
25/40
5 8
6 6
35/45 20/25 45/50 6 6
Includes the orders defined for allowing washcast processing.
Table A2. Process Data Common for All Examples parameter
value
hot metal quantity per feed pot, Φ max j minimum hold time in furnace, τhold minimum wash duration in furnace, τwash cycle time for casters, τcast changeover time for casters, τch yield factors for all casters, ηj,c (C1/C2/C3)
400 tons 3h 2h 3h 6h 0.94/0.95/0.95
Table A3. Feed-Pot Availability Feed1
Feed2
Feed3
problem rFeed1 (h) dFeed1(h) rFeed2 (h) dFeed2 (h) rFeed3 (h) dFeed3 (h) P1 P2 P3
0 0 0
15 15 12
15 15 12
30 30 24
24
36
Table A4 gives order data for example problem P1, whereas Table A5 gives order data for example problems P2 and P3.
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7615 Table A4. Order Dataa for Example Problem P1 order number
alloy
mold type
Fm (h)
δm (h)
Demm (ton)
order number
alloy
mold type
Fm (h)
δm (h)
Demm (ton)
1 2 3 4 5
A1 A2 A2 A3 A1
C1 C1 C2 C2 C2
0 0 6 12 6
16 16 24 24 30
70 35 40 40 85
6 7 8 9 10
A3 A1 φ φ φ
C3 C3 C1 C2 C3
16 16 0 0 0
30 30 30 30 30
47 47 0 0 0
a
The last three orders are defined for allowing washcast processing.
Table A5. Order Dataa for Example Problems P2 and P3 Problem P2 order number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 a
Problem P3
alloy
mold type
Fm (h)
δm (h)
Demm (ton)
alloy
mold type
Fm (h)
δm (h)
Demm (ton)
A1 A2 A1 A3 A1 A3 A1 A2 A2 φ φ φ
C1 C1 C2 C2 C2 C3 C3 C2 C1 C1 C2 C3
0 0 6 12 6 16 16 6 16 0 0 0
16 16 24 24 30 30 30 16 30 30 30 30
70 35 35 35 70 35 35 35 35 0 0 0
A1 A2 A2 A3 A1 A3 A2 A4 A4 A2 A3 A1 A3 φ φ φ
C1 C1 C2 C2 C2 C3 C3 C2 C1 C2 C1 C3 C1 C1 C2 C3
0 0 0 12 6 0 16 6 16 12 0 20 12 0 0 0
16 16 16 24 30 16 30 16 30 36 30 36 36 36 36 36
70 35 35 35 70 35 35 35 35 70 105 70 105 0 0 0
The last three orders in each example are defined for allowing washcast processing.
Appendix B. Simplification of the Full (FS) Model in the Iterative Algorithm The following constraints are modified from the full mixedinteger program to form the sub-problem (SP) in the iterative algorithm. (The equation numbering scheme indicates the equation given in the main text. For example, eq B22 corresponds to eq 22 in the main text.) In eqs 21, 22, 24, 25, 34, and 35, the left-hand side summations are modified.
∑
uj,k, j′,l )
j∈JF,k∈K( j), ( j,k, j′)∈FC( j,k, j′)
∑
∀ j′ ∈ JC, l ∈ L( j′)
xj,k,s ∑ s∈S
∀ j ∈ JF, k ∈ K( j) (B22)
Ffcs, j,k,c, j′,l ) Bcj′,l,c
s∈S, j∈JF,k∈K( j) ( j,k, j′)∈FC( j,k, j′)
∀ j′ ∈ JC, l ∈ L( j′), c ∈ C (B24)
∑
Ffcs, j,k,c, j′,l ) Bfoutj,k,s
c∈C, j′∈JC,l∈L( j′) ( j,k, j′)∈FC( j,k, j′)
∀ j ∈ JF, k ∈ K( j), s ∈ S (B25)
∑
Ams,c, j′,l,m ) ηj,c
m∈MS(s)∩MC(c) ( j,m)∈CO( j,m)
∑ ∑ Ffcs, j,k,c, j′,l j∈JF keK( j)
∀ j′ ∈ JC, s ∈ S, c ∈ C, l e L( j′) (B34)
∑
min(CC min c, j )uj,k, j′,l e c∈C
(CC max ∑ ∑Ffcs, j,k,c, j′,l e max c, j )uj,k, j′,l c∈C s∈S c∈C
∀ j ∈ JF, k ∈ K(j), j′ ∈ JC, l ∈ L(j′): (j,k,j′) ∈ FC(j,k,j′) (B23, B40) Tcsj′,l - H(1 - uj,k, j′,l) e Tffj,k e Tcsj′,l + H(1 - uj,k, j′,l) ∀j ∈ JF, k ∈ K( j), j′ ∈ JC, l ∈ L( j′):( j,k, j′) ∈ FC( j,k, j′) (B26) Equations 33 and 36 are expressed over the set CO(j,m).
(B21) uj,k, j′,l )
j′∈JC,l∈L( j′) ( j,k, j′)∈FC( j,k, j′)
∑
∑wj′,l,c
c∈C
Equations 23, 26, and 40 are expressed over the set FC(j,k,j′).
∑ Ams,c, j,l,m g Demm
j∈JC leL( j) ( j,m)∈CO( j,m)
∀ m ∈ M, s ) S(m), c ) C(m) (B35)
max )zj,l,m Ams,c, j,l,m e max(ηj,c′CC j,c′ c′∈C
∀ j ∈ JC, m ∈ M, s ) S(m), c ) C(m), l e L( j): ( j,m) ∈ CO( j,m) (B33) Fmzj,l,m - s1j,l,m e Tcfj,l e δm zj,l,m + H(1 - zj,l,m) + s2j,l,m ∀ j ∈ JC, l ∈ L( j), m ∈ M: ( j,m) ∈ CO( j,m) (B36) Nomenclature Models FS ) full-space model MP ) master problem SP ) sub-problem Indices c, c′ ) mold type i, i′ ) iteration j, j′, j′′ ) unit k ) furnace slot l ) caster slot m ) order s ) alloy φ ) wash alloy
7616
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006
Sets C ) caster mold types; C ) {C1, C2, ...} C(m) ) mold type for order m; C(m) ⊆ C JC ) casters; JC ) {Cast1, Cast2, ...} JF ) furnaces; JF ) {Fur1, Fur2, ...} JP ) feed pots; JP ) {Feed1, Feed2, ...} M ) orders; M ) {1, 2, ...} MS(s) ) orders for alloy s; MS(s) ⊆ M MC(c) ) orders for mold type c; MC(c) ⊆ M S ) alloys; S ) {A1, A2, ..., φ} S* ) S\{φ} S(m) ) alloy for order m; S(m) ⊆ S Parameters Bi ) parameter used in the integer cut in the MP for iteration i CC min j,c ) minimum batch size for a mold type c in caster j CC max j,c ) maximum batch size for a mold type c in caster j chngs ) maximum number of changeovers allowed in MP and SP Demm ) minimum amount due for order m FTFj ) fill-to-full amount in furnace j DTEj ) draw-to-empty amount in furnace j H ) big-M parameter K(j) ) maximum number of slots for furnace j in FS and SP KM(j) ) maximum number of slots for furnace j in MP L(j) ) maximum number of slots for caster j in FS and SP LB ) lower bound for the MP objective LM(j) ) maximum number of slots for caster j in MP ) maximum holding capacity in furnace j FC max j rj ) release time for feed pot j dj ) due time for feed pot j Sf0j,s ) opening inventory of alloy s in the furnace j slack ) maximum allowed sum of slack variables MP x j,k,s ) value of variable xj,k,s from the MP solution MP V j,k, j′m ) value of variable Vj,k,j′,m from the MP solution MP y j′,c ) value of variable yj,c from the MP solution i R j,k, j′m ) parameter used in the integer cut in the MP during iteration i ηj,c ) yield factor in caster j and mold type c Φ min ) minimum amount available from feed pot j j
) maximum amount available amount from feed pot j Φ max j γ ) greatest common factor (GCF) of all time-related parameters in the problem data Fm ) release times for order m δm ) due times for order m τcast ) cycle time for a batch in the casters τch ) changeover (setup) time in the casters τ hold ) minimum holding time in furnace j j wash τ ) washing (changeover) time in the furnaces ω1, ω2, ω3, ω4 ) weight factors for different terms in the objective function Binary Variables qj′,j,k ) 1 when there is flow from feed source j′ to furnace slot (j,k) uj,k,j′,l ) 1 when there is flow from furnace slot (j,k) to caster slot (j′,l) Vj,k,j′,m ) 1 when furnace slot (j,k) satisfies order m via caster j′ in the MP wj,l,c ) 1 when mold type c is assigned to caster slot (j,l)
wˆ j,l ) 1 when there is a changeover in caster j immediately before slot l xj,k,s ) 1 when alloy s is assigned to furnace slot (j,k) yj,c ) 1 when atleast one batch of mold type c is processed in caster j in the MP zj,l,m ) 1 when caster slot (j,l) satisfies an order m PositiVe Variables Ams,c,j,l,m ) amount from caster slot (j,l) toward order m for product (s ) S(m), c ) C(m)) Am′j,k,j′,m ) amount from furnace slot (j,k) toward order m via caster j′ in the MP Bcj,l,c ) batch size of mold type c in caster slot (j,l) Bfj,k,s ) batch size of alloy s in furnace slot (j,k) in the MP Bfinj,k,s ) flow of alloy s into furnace slot (j,k) Bfoutj,k,s ) flow of alloy s from furnace slot (j,k) Bwashj,k,s ) amount of alloy s incorporated into the wash state in furnace slot (j,k) Ffcs,j,k,c,j′,l ) flow quantity of alloy s in furnace slot (j,k) to mold type c in caster slot (j′,l) Fpfj′,j,k ) flow quantity from feed pot j′ to furnace slot (j,k) moldsj ) number of active mold types in caster j in the MP obj′ ) objective function in the MP s1j,l,m, s2j,l,m ) slack variables in the FS and SP s1′j,k,m, s2′j,k,m ) slack variables in the MP Sfj,k,s ) inventory of alloy s at the beginning of furnace slot (j,k) SfTj,s ) inventory of alloy s at the end of the last slot in furnace j Tcsj,l ) caster slot (j,l) start time Tcfj,l ) caster slot (j,l) finish time Tfsj,k ) furnace slot (j,k) start time Tffj,k ) furnace slot (j,k) finish time Literature Cited (1) Gravel, M.; Price, W.; Gangne, C. Scheduling Jobs in an Alcan Aluminum Foundry using a Genetic Algorithm. Int. J. Prod. Res. 2000, 38 (13), 3031-3041. (2) Gravel, M.; Price, W.; Gangne, C. Scheduling Continuous Casting of Aluminum Using a Multiple Objective Ant Colony Metaheuristic. Eur. J. Oper. Res. 2002, 143, 218-229. (3) Pacciarelli, D.; Pranzo, M. Production Scheduling in a SteelmakingContinuous Casting Plant. Comput. Chem. Eng. 2004, 28, 2823-2835. (4) Nicholls, M. G. Optimizing the Operations of an Ingot Mill in an Aluminium Smelter. J. Oper. Res. Soc. 1994, 45 (9), 987-998. (5) Nicholls, M. G. Developing an Integrated Model of an Aluminium Smelter Incorporating Sub-models with Different Time Bases and Levels of Aggregation. Eur. J. Oper. Res. 1997, 99, 477-490. (6) Moon, S. D.; Hrymak, A. N. Scheduling of the Batch Annealing ProcesssDeterministic Case. Comput. Chem. Eng. 1999, 23, 1193-1208. (7) Harjunkoski, I.; Grossmann, I. E. A Decomposition Approach for the Scheduling of a Steel Plant Production. Comput. Chem. Eng. 2001, 25, 1647-1660. (8) Naphade, K. S.; Wu, S. D.; Storer, R. H.; Doshi, B. J. Melt Scheduling to Trade Off Material Waste and Shipping Performance. Oper. Res. 2001, 49 (5), 629-645. (9) Tang, L.; Luh, P.; Liu, J.; Fang, L. Steel-making Process Scheduling Using Lagrangian Relaxation. Int. J. Prod. Res. 2002, 40 (1), 55-70. (10) Schwindt, C.; Trautmann, N. Scheduling the Production of Roling Ingots: Industrial Context, Model, and Solution Method. Int. Trans., Oper. Res. 2003, 10, 547-563. (11) Reklaitis, G. V. Overview of scheduling and planning of batch process operations, NATO Advanced Study InstitutesBatch Process Systems Engineering, Antalya, Turkey, 1992. (12) Pinto, J.; Grossmann, I. E. Assignment and Sequencing Models for the Scheduling of Chemical Processes. Ann. Oper. Res. 1998, 81, 433466. (13) Shah, N. Single- and Multisite Planning and Scheduling: Current Status and Future Challenges. AIChE Symp. Ser. 1998, 94, 75.
Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7617 (14) Mendez, C. A.; Cerda, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-Art Review of Optimization Methods for Short-Term Scheduling of Batch Processes. Comput. Chem. Eng. 2006, 30, 913-946. (15) Pinto, J.; 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-3051. (16) Me´ndez, C. A.; Cerda´, J. An MILP-Based Approach to the ShortTerm Scheduling of Make-and-Pack Continuous Production Plants. OR Spectrum 2002, 24, 403-429. (17) Gupta, S.; Karimi, I. A. An Improved MILP Formulation for Scheduling Multiproduct, Multistage Batch Plants, Ind. Eng. Chem. Res. 2003, 42, 2365-2380. (18) 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-9190. (19) Kondili, E.; Pantelides, C. C.; Sargent, R. A General Algorithm for Short-Term Scheduling of Batch OperationssI. MILP Formulation. Comput. Chem. Eng. 1993, 17, 211-227. (20) Shah, N.; Pantelides, C. C.; Sargent, R. A General Algorithm for Short-Term Scheduling of Batch OperationssII. Computational Issues. Comput. Chem. Eng. 1993, 17, 229-244. (21) Pantelides, C. C. Unified Frameworks for the Optimal Process Planning and Scheduling. In Proceedings on the Second Conference on Foundations of Computer Aided Operations; 1994; pp 253-274. (22) Schilling, G.; Pantelides, C. C. A Simple Continuous-Time Process Scheduling Formulation and a Novel Solution Algorithm. Comput. Chem. Eng. 1996, 20, S1221-S1226. (23) Zhang, X.; Sargent, R. W. H. The Optimal Operation of MixedProduction FacilitiessA General Formulation and Some Approaches for the Solution. Comput. Chem. Eng. 1996, 20, S897-S904. (24) 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-4359. (25) Mockus, L.; Reklaitis, G. V. Continuous Time Representation Approach to Batch and Continuous Process Scheduling. 1. MINLP Formulation. Ind. Eng. Chem. Res. 1999, 38, 197-203. (26) Castro, P.; Barbosa-Povoa, A. P. F. D.; Matos, H. An Improved RTN Continuous-Time Formulation for the Short-Term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2001, 40, 2059-2068. (27) Giannelos, N. F.; Georgiadis, M. C. A Simple New ContinuousTime Formulation for Short-Term Scheduling of Multipurpose Batch Processes. Ind. Eng. Chem. Res. 2002, 41, 2178-2184. (28) Maravelias, C. T.; Grossmann, I. E. A New General ContinuousTime State Task Network Formulation for the Short-Term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2003, 42 (13), 30563074.
(29) Sundaramoorthy, A.; Karimi, I. A. A Simpler Better Slot-Based Continuous-Time Formulation for Short-Term Scheduling in Multipurpose Batch Plants. Chem. Eng. Sci. 2005, 60 (10), 2679-2702. (30) Maravelias, C. T. Mixed-time Representation for State-Task Network Models. Ind. Eng. Chem. Res. 2005, 44 (24), 9129-9145. (31) Floudas, C. A.; Lix, X. Continuous-Time Versus Discrete-Time Approaches for Scheduling of Chemical Processes: A Review. Comput. Chem. Eng. 2004, 28 (11), 2109-2129. (32) Burkard, R. E.; Hatzl, J. Review: Extensions and Computational Comparison of MILP Formulations for Scheduling of Batch Processes. Comput. Chem. Eng. 2005, 29, 1752-1769. (33) Harjunkoski, I.; Grossmann, I. E. Decomposition Techniques for Multistage Scheduling Problems Using Mixed-Integer and Constraint Programming Methods. Comput. Chem. Eng. 2002, 25 (11), 15331552. (34) Kelly, J. D. Chronological Decomposition Heuristic for Scheduling: Divide and Conquer Method. AIChE J. 2002, 48 (12), 29952999. (35) Kelly, J. D.; Mann, J. L. Flowsheet Decomposition Heuristic for Scheduling: A Relax-and-Fix Method. Comput. Chem. Eng. 2004, 28, 2193-2200. (36) Engell S.; Markert A.; Sand, G.; Schultz R. Aggregated Scheduling of a Multiproduct Batch Plant by Two-Stage Stochastic Integer Programming. Optim. Eng. 2004, 5 (3), 335-359. (37) Maravelias, C. T.; Grossmann, I. E. A Hybrid MILP/CP Decomposition Approach for the Continuous Time Scheduling of Multipurpose Batch Plants. Comput. Chem. Eng. 2004, 28 (10), 1921-1949. (38) Roe, B.; Papageorgiou, L. G.; Shah, N. A hybrid MILP/CLP Algorithm for Multipurpose Batch Process Scheduling. Comput. Chem. Eng. 2005, 29 (6), 1277-1291. (39) Maravelias, C. T. A Decomposition Framework for the Scheduling of Single- and Multi-Stage Processes. Comput. Chem. Eng. 2006, 30 (3), 407-420. (40) Ryan, D. M. Optimized Cell Batching for New Zealand Aluminium Smelters Ltd. Presented at the 33rd Annual Conference of the Operational Research Society of New Zealand (ORSNZ), 1998.
ReceiVed for reView May 23, 2006 ReVised manuscript receiVed August 6, 2006 Accepted August 23, 2006 IE060652H