Mixed-Time Mixed-Integer Linear Programming ... - ACS Publications

Mar 31, 2007 - the functional reference grid for the time representation, included in discrete-time ... In ref 17, a mixed-time model using a fixed-ti...
0 downloads 0 Views 680KB Size
Ind. Eng. Chem. Res. 2007, 46, 2781-2796

2781

Mixed-Time Mixed-Integer Linear Programming Scheduling Model Joakim Westerlund,*,† Mattias Ha1 stbacka,† Sebastian Forssell,‡ and Tapio Westerlund§ Centre for Industrial Engineering and Management and Process Design Laboratory, Åbo Akademi UniVersity, FIN-20500 Åbo, Finland, and Borealis Group, P.O. Box 330 FIN-06101, PorVoo, Finland

This paper presents a novel mixed-time mixed-integer linear programming (MILP) scheduling model for industrial problems where intermediate storage handling is of particular concern. The proposed model is a crossbreed between a continuous-time and a discrete-time model where a continuous-time representation is incorporated in a discrete-time grid. The mixed-time formulation is useful in the modeling of multi-stage multi-product production processes using intermediate storages with nonlinear optimal storage profiles. The model is combining useful features from continuous-time models, flexibility, and solution efficiency, with the functional reference grid for the time representation, included in discrete-time representations. The model is able to handle both short-term and periodic problems. The main purpose of the presented model is not to challenge existing models on solution efficiency or quality on theoretical problems but to present a model that is very useful in the modeling of industrial scheduling problems where intermediate storage handling is of importance. The proposed model is able to solve industrial problems which, according to our knowledge, have not previously been solvable by existing models. The novelty with the proposed model is that an entirely continuous-time representation is incorporated in a uniform time grid. In practice this means that despite the continuous-time representation, aging profiles or highly nonlinear storage inventory profiles may, for example, be addressed with the proposed model. Introduction Scheduling is concerned with the allocation of production resources to production orders in an optimized way. Scheduling is utilized to support production, purchasing, logistics, and so forth and plays a key role in process operations where it may yield great improvements of production performance.1 The research area of batch and continuous process scheduling has received great attention from both academia and industry in the past two decades.2 Consequently, a remarkable amount of mathematical models for the solution of these challenging problems have been developed recently. Some well-known models can be found in ref 3, presenting an early discrete-time mixed-integer nonlinear programming (MINLP) model for multiproduct batch plant design,4 introducing the discrete-time state-task network,5 and presenting a mixed-integer linear programming (MILP) model for the scheduling of multipurpose batch plants operated in a cyclic mode, and the resource task network model is presented in ref 6, addressing simultaneous scheduling and design. Some more recent models can be found in refs 7 and 8, introducing the event point concept or the Resource Task Network (RTN)-based continuous-time model presented in ref 9, which was further developed into the model presented in ref 10. Other interesting models are the continuoustime State Task Network (STN)-based model presented in ref 11, the immediate batch precedence model of ref 12, and the grid based common variable continuous-time model in ref 13, to mention a few. Reviews of different models and techniques can be found, for example, in refs 1, 2, 14, and 15. Some scheduling approaches are only suitable for a specific type of problem, while others are more general in the sense that they * To whom correspondence should be addressed. Tel.: +358 2 215 4705. Fax: +358 2 215 4791. E-mail address: joakim.westerlund@ abo.fi. † Centre for Industrial Engineering and Management, Åbo Akademi University. ‡ Borealis Group. § Process Design Laboratory, Åbo Akademi University.

can handle different network structures, storage policies, and even changeover requirements.16 The model presented in this paper is of a fairly general nature and is able to handle a broad number of network structures, advanced intermediate storage policies (including aging aspects), and changeovers. In ref 17, a mixed-time model using a fixed-time grid but allowing processing tasks to have variable duration was presented. Even though the model of ref 17 and the model presented in this paper are similar in the way that they allow events to take place continuously while at the same time taking advantage of the positive features with a discrete-time grid, the two models differ significantly in many aspects. The model presented in this paper differs from the model presented in ref 17, for example, in the sense that it allows tasks to start and finish anywhere within a time interval while the model of ref 17 forced tasks to start exactly at a time point but allowed tasks to finish anywhere within a time interval. It is worth noticing that by this procedure, the model presented in this paper requires intermediate storage capacity to ensure that the changeover procedure works properly in all cases. The model is, hence, not directly applicable to problems working under strict JIT discipline. The two models also differ significantly in the actual mathematical formulation. For example, the model presented in ref 17 assumes that a unit is always assigned to a task, while this is not the case with the model presented in this paper. Another important difference between the two models is in the formulation of the assignment constraints. The model presented in ref 17 uses the following variables for the assignment constraints: (i) A binary Zjn ) 1 if a task starts on unit j at time point n (ii) A binary Wsin ) 1 if task i starts at time point n (iii) A binary Wpin ) 1 if task i is being processed at time point n (iv) A binary Wfin ) 1 if task i finishes at or before time point n The model presented in this paper, on the other hand, uses the following variables for the assignment constraints:

10.1021/ie060991a CCC: $37.00 © 2007 American Chemical Society Published on Web 03/31/2007

2782

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Chart 1

(i) A binary yi,k,t ) 1 if task i starts or continues at unit k in the beginning of interval t (ii) A continuous variable wi,k,t, providing an upper limit for the momentary production as a fraction of the maximum production for a whole time grid The proposed model was initially developed to offer advanced decision-making support for production planning and intermediate storage handling for the first stage of an industrial production process including many stages. The model of ref 17 was, in contrast, developed to address scheduling as an integrated module in supply chain management problems. When formulating an industrial scheduling problem, the time representation is typically one of the first aspects needing consideration. Time representation is, in this paper, divided into two major subgroups: continuous or discrete. The proposed mixed-time model may be characterized as a crossbreed between discrete-time and continuous-time models. Continuous-time formulations are often favored as they tend to obtain better quality solutions and require significantly less computational resources.18 The academic focus has, hence, shifted more and more to continuous-time formulations, and numerous formulations have been proposed in the literature on the basis of continuous representation of time in the last 10 years.19 Continuous-time formulations are, however, not necessarily more efficient than their discrete-time counterparts, and the selection of the most efficient formulation for a particular problem is still very much an open issue.18 Discrete-time formulations, dividing the time horizon into a number of time intervals of uniform durations, offer, for example, the advantage of having a convenient reference grid for the time representation. This reference grid does not only offer an easy overview of all process equipment units but also enables the formulation of strongly nonlinear production profiles using a “piecewise linear” representation. This feature has been utilized in the proposed mixed-time model. The main disadvantage of discrete-time models is that variable processing times can be handled only as discrete approximations and that the number of intervals may be so large that the resulting model is too hard to solve.21 This paper also presents software using the proposed mixedtime model. The software is currently in continuous industrial use. 1. Time Representation Event representation over time is usually the first aspect needed be taken under consideration in the problem formulation procedure of scheduling problems. The event representation is often divided into subgroups according to their physical features. In Chart 1, a roadmap for optimization of scheduling problems is presented. The focus in network-oriented processes is, according to the classification above, on reference points to check resources, while the focus in batch-oriented processes is on the arrangement of resource utilization.

As stated in the introduction section, scheduling is concerned with the allocation of jobs to tasks in time, and the time representation is frequently divided into two main groups: discrete and continuous. As a result of the massive combinatorial complexity of large scale discrete-time formulations and the fact that a discrete approximation of the time horizon leads to suboptimal solutions by definition,2 continuous-time formulations have in recent times received more attention. There are, nevertheless, also some very useful features linked to discretetime formulations. A discrete-time formulation is, for example, very practical when modeling nonlinear production profiles in multi-stage scheduling problems with intermediate storages. If aging profiles, additionally, are required for the storages, discrete-time models are preferable. A good example of a practical application where aging profiles need be taken into account in the scheduling procedure is the tissue manufacturing industry. Tissue mills using both virgin fiber and de-inked pulp in their end-products may store de-inked intermediate products in storage towers for a couple of days, using preservative chemicals. The pulp is, however, going bad after a given time in storage.20 In such cases, the optimal production plan has a highly nonlinear optimal storage profile which is essential for the production. It is, hence, a topic of great concern to develop a model that is able to take aging profiles for intermediate storages into consideration in the actual scheduling procedure. Even though continuous-time models generally outperform discrete-time models, concerning solution efficiency and quality, discrete-time models remain popular for industrial problems, keeping the debate over the effectiveness of discrete versus continuous-time representation open.21 The model presented in this paper was originally developed to handle a single very large industrial application and was gradually developed from an early purely discrete-time model to its final mixed-time form. In the following sections, a brief introduction to two basic discrete- and continuous-time models is presented. It is worth noticing that the models presented below are two specific, very uncomplicated, versions of discrete and continuous-time models. These models might not be the most well chosen ones for illustrating discrete and continuous-time models in a general way. However, the two models were selected because of the fact that they provide an easy way to illustrate how the mixedtime model utilizes features from continuous and discrete-time models. For example, the mathematical formulation of the mixed-time model is based on the discrete-time model presented below, and it is fairly easy to see how the mixed-time model is an evolution of this model. The mixed-time model is, subsequently, presented and compared to the discrete-time model, presented below, on a number of illustrative example problems. The mixed-time model is not compared to any continuous-time models as it is developed to handle a class of industrial problems which are very difficult, if not impossible, to model using a continuous-time representation. In section 1.1 the discrete-time model which was used as a starting point for the proposed mixed-time formulation is used, and in section 1.2 a simple continuous-time representation is presented. These formulations are not representative for all discrete- and continuous-time models but are presented to make it easier to discover how the mixed-time model was developed. For a more thorough review of existing scheduling models the reader is referred to ref 12 or 21. 1.1. Discrete-Time Representation. A traditional way of assigning processing units to production tasks, in discrete-time formulations, is to introduce binaries for every task starting at each unit and each discrete-time interval. The corresponding

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2783

Figure 1. Discrete-time production plan.

Figure 2. Discrete-time Gantt chart.

binary for task i starting processing at unit k at time interval t is, thus, yi,k,t. This binary is equal to 1 if task i is taking place at unit k during time interval t and 0 otherwise. If the unit to task binary yi,k,t is equal to 1, the production pi,k,t in task i at unit k within time grid t may be up to the maximum value, pi,k,max, in task i at unit k according to eq 1.1.1. A similar procedure was used in the state task network presented in ref 4.

pi,k,t e yi,k,tpi,k,max

∀i ∈ I, k ∈ K, t ∈ T

(1.1.1)

Furthermore, to manage tasks competing for shared processing equipment, eq 1.1.2 is used to ensure that not more than one task takes place at each unit at any time. I

yi,k,t e 1 ∑ i)1

∀k ∈ K, t ∈ T

(1.1.2)

In Figure 1 a discrete-time production plan for a processing unit, k, with production capacity for two tasks, A and B, is presented. The momentary production, batch size, is displayed on the y-axis and the time horizon on the x-axis. Figure 1 presents a discrete-time production plan where the fact that batches, in discrete-time models, may be of variable size, as long as the length, in the time dimension, is a multiple of the discretization grid. Figure 2 presents the same discrete-time production plan in a more familiar form as a Gantt chart. 1.2. Continuous-Time Representation. Continuous-time scheduling models have received a significant amount of attention during the last decades. In continuous-time models events or tasks are not restricted to specific grid points. In other words, events can take place at any time unlike in discretetime models. The authors of ref 19 are broadly classifying continuous-time models into three distinct categories: slot based, global-event based, and unit-specific-event based formulations. The mathematical formulation in the most basic continuoustime models may generally be classified as models using either one set of variables for assignment and sequencing, for example, yii′k )1 if task i precedes task i′ on unit k, or separate variables for assignment and sequencing, for example, yik ) 1 if task/job i is processed on unit k and yii′ ) 1 if task i precedes task i′. The constraints defining that a task may not overlap another and that due and release dates may not be violated are essential. The non-overlapping constraints may be written as in eqs 1.2.1 and 1.2.2. It is clear that if yik and yi′k are equal to 1, either yii′ or yi′i must be equal to zero, implying that either i is preceding i′ or i′ is preceding i. tsi is the starting time for task i, and pik is

the processing time for task i at unit k. Big-M formulations are, typically, used to write the disjunctive non-overlapping constraints in conjunctive form where M is an appropriate upper bound.

yii′ + yi′i g yik + yi′k - 1, tsi′ g tsi +

∀i, i′ ∈ I, i > i′, k ∈ K (1.2.1)

∑pikyik - M(1 - yii′)

k∈K

∀i, i′ ∈ I, i * i′ (1.2.2)

Equation 1.2.2 is transforming the disjunction describing that two jobs produced on the same processing unit cannot overlap in time into conjunctive form stating that task i′ should start no earlier than at the completion time of task i at unit k, if both tasks are processed at the same unit. Figure 3 is presenting a possible production plan using a continuous-time model. The batch size is displayed on the y-axis, and time is displayed on the x-axis. Figure 3 displays the momentary production in units, u, versus relative time unit, rtu. In Figure 3 the fact that a continuous-time production plan may have both variable batch size and length (in time) is displayed, u being units produced and rtu being relative time units. In Figure 4, the same production plan is presented as a Gantt chart. The main advantage with a continuous-time model is, obviously, that the time representation is by definition continuous and is not confined to discrete event points. A great disadvantage is, however, that the number of binary variables used in the non-overlapping constraints increases quadratically with the number of tasks while the binaries in a discrete-time model increase linearly with the discretization. However, the more recent continuous-time models require significantly less binary variables in the problem formulation than the illustrative one presented above. Slot-based continuous-time formulations are based on a prespecified number of time slots with unknown duration to be allocated to batches. The batches are defined a priori, and no mixing or splitting operations may take place. Batches may, further, start and finish at any time during the scheduling horizon. The main advantage with slot-based models is the significant reduction in model size when a minimum number of time slots is predefined, generally resulting in good computational performance. Slot-based models, furthermore, offer a simple model and easy representation for sequencing and allocation scheduling problems. The main disadvantage is the modeling of resource and inventory constraints and the fact that predefining the number of time slots may result in sub-optimal

2784

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 3. Continuous-time production plan.

Figure 4. Continuous-time Gantt chart. Table 1. Main Features of the Most Common Approaches to Scheduling Problems time representation: discrete event representation

time representation: continuous

global time intervals

global time points

unit-specific time events

main decisions

lot-sizing, allocation, sequencing, timing

type of process

general network

time slots

unit-specific immediate precedence

immediate precedence

general precedence

allocation, sequencing, timing sequential

material balances

network flow equations

network flow equations

critical modeling issues

time interval duration, scheduling period (data dependent)

number of time points (iteratively estimated)

number of time points (iteratively estimated)

number of time slots (estimated)

number of batch tasks sharing units (lot sizing)

number of batch tasks sharing units (lot sizing)

number of batch tasks sharing resources (lot sizing)

critical problem features

variable processing time, sequencedependent changeovers

intermediate due dates and raw material supplies

intermediate due dates and raw material supplies

resource limitations

inventory, resource limitations

inventory, resource limitations

inventory

solutions if the number of slots is not well chosen. Furthermore, the model size and complexity is greatly dependent on the number of predefined time slots. For more details on slot-based models, readers are referred to ref 22 or 23. Unit-specified direct precedence based models, on the other hand, are assuming a priori defined batches, no mixing and splitting operations, and that batches may start and finish at any time during the scheduling horizon. The advantages being that sequencing is explicitly considered in model variables and changeover times and costs are easy to implement. Disadvantages with unitspecific direct precedence models are the large number of sequencing variables and the fact that resource and material balances are hard to model. For more details on unit-specific direct precedence based models, readers are referred to ref 24. The unit-specific direct precedence based models were followed by the global direct precedence; readers are referred to refs 25 and 12 for more details on global direct precedence based models. The global general precedence based model was the next step in the development of precedence based models and offers the additional advantage of cutting down the number of sequencing variables, resulting in faster solution. The sequencing decisions may in global general precedence models be extrapolated to other resources. For more details on global general precedence models, readers are referred to refs 26-28. Table 1 presents the most common optimization approaches, as well as their main features. The table is a simplified version of the one presented in ref 1. In this paper we present a mixed-time formulation, utilizing helpful features from discrete and continuous-time models alike.

batch oriented

The result is a model where changeovers may occur continuously at any point in time at all process units. A discrete reference grid, using global time intervals to keep track of the time dimension, is used to approximate possible nonlinear production profiles and aging in storage. By using a discretetime grid nonlinear production profiles are modeled by piecewise linear representations. 1.3. Mixed-Time Representation. The formulation presented in this paper is a mixed-time formulation using a global uniform time grid but at the same time enabling events, such as changeovers, to take place at any given time, within any of the time grids, during the scheduling horizon. This procedure makes the mixed-time model significantly more flexible than a traditional discrete-time model where changeovers may only take place at the beginning/end of a time grid. The mixed-time formulation, as presented in this paper, may easily be incorporated into other discrete-time formulations, allowing for events to take place continuously within the discrete-time intervals. In this paper the mathematical formulation is based on the discretetime model presented in section 1.1, making it more easy to illustrate how the model is formulated. The mixed-time formulation is basically an improvement of a traditional discrete-time formulation, such as the one mentioned in section 1.1. The main novelty in the considered formulation is the ability to make task/ product changeovers at any time within the scheduling horizon, not only at discrete time points. Tasks may not only stop but also start at any time within the time intervals. The unit to task constraints used in the presented formulation are presented in eq 1.3.1.

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2785

pi,k,t e wi,k,tpi,k,max

∀i ∈ I, k ∈ K, t ∈ T (1.3.1)

The presented formulation uses two continuous variables, wi,k,t and pi,k,t, to determine the momentary production of each task, i, at each unit, k, during the time sequence, t, where a traditional discrete-time model uses one discrete variable and one continuous. pi,k,t represents the momentary production during time grid t and wi,k,t is providing a fractional upper limit for the momentary production and is thereby bounded to take values between 0 and 1, while pi,k,max is the global maximum production rate for task i at unit k during one grid. The wi,k,t variable is determined by the binaries yi,k,t and yi,k,t+1 as shown in eq 1.3.2. The wi,k,t variable is constrained to assume continuous values between 0 and 1. This procedure implies that instead of, as in a strict discrete-time formulation, using the binary yi,k,t to determine whether or not a given task is performed during a specific time grid, the binary yi,k,t is now used to highlight if a task, i, starts or continues at unit k in the beginning of the time grid t. yi,k,t may, hence, be equal to zero although the task i at unit k is performed during the time grid t in case it will not start in the beginning of the grid. A changeover may, thus, take place anywhere (any time) within a grid and not only at the beginning or end of a grid. The conditions for wi,k,t are given below. The yi,k,t variable is defined in the same way as in the purely discrete formulation presented in section 1.1. If required, a “smoothing procedure” may be applied where fluctuation of the production between adjacent time intervals is penalized, resulting in a more smooth production, for example, in multi-stage problems where the production capacity of the subsequent production stages is less than the capacity of the earlier stages.

wi,k,t e yi,k,t + yi,k,t+1

∀i ∈ I, k ∈ K, t ∈ {1, 2, 3, ..., N - 1} (1.3.2)

At the end of the time horizon, during time grid N, the wi,k,t variable is, obviously, determined solely by the final yi,k,t variable because no yi,k,t+1 exists.

wi,k,N e yi,k,N

∀i ∈ I, k ∈ K

(1.3.3)

To make sure that not more than one task is performed at any unit at any time, except for the case when a changeover is taking place, eq 1.3.4 is used. Furthermore, to make sure that no unit is producing anything at a level exceeding the maximum throughput, eq 1.3.5 is used. I

yi,k,t e 1 ∑ i)1

∀k ∈ K, t ∈ T

(1.3.4)

∀k ∈ K, t ∈ T

(1.3.5)

I

wi,k,t e 1 ∑ i)1

To make sure that task i may be performed at unit k during time-sequence t, if yi,k,t is equal to 1, eq 1.3.6 is used.

yi,k,t e wi,k,t

∀k ∈ K, t ∈ T

(1.3.6)

In Figure 5, a short hypothetical production plan using a mixed-time formulation is presented. The underlying table presents potential values on the yi,k,t and wi,k,t variables. Figure 5 presents the fact that both batch size and length are variable using the mixed-time model. Figure 6 is presenting the same production plan as a Gantt chart. As can be observed from Figure 5, the wi,k,t variable gives the fractional production capacity limit for each time grid. For

example, if wA,k,t is equal to 0.8 and the maximum production capacity, pA,k,max for A at unit k for the time grid t is equal to 10 tonnes, the momentary production, pA,k,t, of task A at unit k during time t should be no more than 8 tonnes. It is worth noticing that like in the discrete-time case, the wi,k,t variable (yi,k,t in the discrete-time case) may be 1 although the actual momentary production is not at 100%. wi,k,t is only giving a limit to the maximum allowed momentary production as seen in eq 1.3.1. Assuming constant production rates within the grids, the exact time of a changeover may, retrospectively to the actual optimization, be calculated using the solution values of the wi,k,t variables according to eq 1.3.7.

∆ti,k,t )

wi,k,t I

∆t

(1.3.7)

wi,k,t ∑ i)1 ∆ti,k,t represents the changeover time calculated from the beginning of the time grid, and ∆t is the length of the specified discretization grid. If the discretization grid, in the example above, is specified to 4 h in Figure 5, the first changeover from task A to task B takes place exactly at 7.2 h from the beginning of the time horizon. Using the ∆ti,k,t variable, one is, for example, able to produce continuous-time Gantt charts to represent production plans. The drawback connected to the utilization of the ∆ti,k,t variable is that no idling will be possible in the illustration of the production plan. In practice, idling is, however, possible (as can be seen in eq 1.3.3) even though the ∆ti,k,t variable is indicating that a certain task is taking place. The ∆ti,k,t variable is, consequently, only used to enable a more easy way of illustrating the exact changeover time for each changeover. One critical requirement for the model to work properly is the existence of intermediate storages. If no intermediate storage is available, the mixed-time model may yield solutions where the changeover procedure is not working properly. 1.3.1. Changeovers. If additional sequence dependent changeover variables (defined in eq 1.3.8) or sequence independent variables (defined in eq 1.3.9) are introduced, the number of changeovers at each unit, k, can be minimized. xCO,ii′,k,t is a continuous variable equal to 1 if a changeover from task i to task i′ takes place in the beginning of time grid t, and it is 0 otherwise. The changeover variable does not need be specified as a binary variable because its integer relaxation will always be 0 or 1 by using eq 1.3.8 if the changeovers are penalized and, consequently, minimized, in the objective function. Because the summarized changeovers are penalized in the objective function, the xCO,ii′,k,t variable will always be equal to 0 if a changeover is not taking place, thereby preventing unnecessary penalizing in cases where all changeovers are penalized equally. xCO,k,t is a sequence independent variable determining whether or not a changeover from or to task i is taking place in the beginning of time-grid t. This variable is also a continuous variable assuming a value of either 0 or 1. It is also possible to incorporate sequence-dependent or sequence-independent changeover times in the formulation. This is done by adding an additional constraining term to the momentary production constraint in eq 1.3.1. The momentary production is then limited to be a certain amount less than the maximum production if a changeover takes place, thereby introducing a changeover time. If changeover times are used, eq 1.3.1 is replaced by eq 1.3.10 where pii′,k,CO is the sequence-dependent changeover time related decrease in production during time interval t at machine k. Using

2786

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 5. Mixed-time production plan.

Figure 6. Mixed-time Gantt chart.

Figure 7. Actual continuous nonlinear storage profile.

this procedure, the changeover time from task i to task i′ is incorporated in the processing time of task i.

xCO,ii′,k,t g yi,k,t-1 + yi′,k,t - 1 xCO,k,t g yi,k,t-1 - yi,k,t xCO,k,t g yi,k,t - yi,k,t-1

∀i, i′ ∈ I, i * i′ (1.3.8)

} }

pi,k,t e wi,k,tpi,k,max - xCO,ii′,k,tpii′,k,CO i * i′

∀i ∈ I

(1.3.9)

∀i, i′ ∈ I, k ∈ K, t ∈ T (1.3.10)

In the proposed model, changeover times or penalties are used only for events happening at consecutive time intervals. The objective function is, consequently, not penalized in cases where “production stops” take place for a whole time grid. If the objective function would be penalized for such events, manually added service or maintenance stops would penalize the objective function with costs or changeover times which, in practice, do not appear. 1.3.2. Aging in Intermediate Storage. Using the mixed-time formulation one is allows easy incorporation of aging profiles for intermediates in storages. This is an important aspect in many industrial applications where the quality of intermediates may be degraded with the storage time. In the worst case, intermediates with limited storage properties may totally be destroyed. The general formulation for the intermediate storages is, in the presented model, formulated as in most discrete-time models, and the storage balance for each storage s may be formulated according to eq 1.3.11.

ms,t+1 ) ms,t + ms,in,t+1 - ms,out,t+1

the time interval. The total flow in or out for a given time period ∆t is given in eq 1.3.12.

ms,in,t+1 )

(1.3.12)

In Figure 8, the piecewise linear approximation of the highly nonlinear storage profile in Figure 7 is presented using either a bar chart or a piecewise linear function graph. A piecewise approximation of the storage profiles is easily achieved from the intermediate storage material balances (in eq 1.3.11) in any discrete-time model. This simple, yet useful, representation may be of great help in industrial problems where an approximation of the storage profile is required to ensure an optimal production scheme. In many cases it is, however, not enough to know how the storage profiles develop over time, and it may also be necessary to know how old the substance in storage is. In the mixed-time formulation, aging of specific intermediates may be penalized in the objective function using the constraints in either eqs 1.3.13 and 1.3.15 or 1.3.14 and 1.3.15, where the amount of “old intermediates” is presented. tkeep represents the time, in grids, after which the intermediate in storage is getting old. tkeep is a function of the grid length ∆t and the actual intermediate specific aging time limit told. t

mold,s,t g ms,t -



z)t-told+1

ms,in,z

told ∈ Z

(1.3.13)

t

mold,s,t g ms,t-told -

(1.3.11)

ms,t+1 is the amount of intermediate product in storage s at time t+1. ms,in,t+1 and ms,out,t+1 are the total batch sizes in and out of storage under a given time interval. According to eq 1.3.11, the storage content is always determined at the end of

∫t t+∆tm˘ s,in,t dt

told )



z)t-told+1

tkeep ∆t

ms,out,z

told ∈ Z (1.3.14)

told, ∆t ∈ Z

(1.3.15)

The aging procedure is demonstrated in Figure 9 and Table 2 below. In Figure 9 and Table 2 the discretization grid is 12 h,

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2787

Figure 8. Piecewise linear representation of nonlinear storage profile.

Figure 9. Illustration of aging profiles in intermediate storage. Table 2. Storage Inventory and Flows into and out of Storage

storage inventory old substance flux in flux out

1

2

3

4

5

6

7

8

9

10

20 0 20 0

45 0 25 0

35 0 10 20

45 15 20 10

50 5 25 20

45 5 15 20

45 30 0 0

25 25 0 20

15 5 10 20

25 0 20 10

and storage of intermediates for more than 24 h is penalized in the objective function. The summarized amount of “old substance”, s, during each time grid with the length ∆t is penalized with a substance-specific aging penalty, pold,s, in the objective function, where Cold,tot is the total aging penalty/cost. The aging penalty is presented in eq 1.3.16. T

Cold,tot ) pold,s

∑ t)1

mold,s,t

(1.3.16)

The amount of “old substance”, substance having been in storage for more than 24 h, at time 5 is consequently, according to eq 1.3.13, ms,5 - (ms,in,5 + ms,in,4) ) 50 - (20 + 25) ) 5. Using eq 1.3.14, the corresponding amount is ms,3 - (ms,out,5 + ms,out,4) ) 35 - (10 + 20) ) 5. In Figure 9 the flows into and out of storage are demonstrated using bars, and the storage profile (inventory) over time as well as the amount of “old substance” are represented as areas. It is important to notice that the discrete-time grid enables a piecewise linearization of highly nonlinear storage profiles: the finer the grid, the better the linearization is describing the actual profile. There is, however, a tradeoff between more accurate linearization (finer grid) and solution time. 2. Illustrative Example Problems The mixed-time model has been tested and proven efficient on a number of test problems. The illustrative example problems are concerned with the profit maximization of a multi-product facility of the flow-shop type. Three different problems are presented to show how the mixed-time formulation works in

Table 3. Production Capacities and Raw Material Price for Illustrative Examples 1 and 2 production capacity (raw material) a

b

c

production capacity (end product) d

A

machine A 5 t/h 6 t/h 4 t/h machine B 6 t/h 3 t/h 4 t/h machine C 6 t/h machine D raw material 90 110 100 120 price euro/t euro/t euro/t euro/t

B

C

4 t/h 5 t/h

4 t/h

Table 4. Product Recipe, Market Demand, and Price for Illustrative Examples 1 and 2 raw material a b c d

product recipe A

B

100%

30% 70%

C 30% 60% 10%

A

B

C

demand (tonnes/week)

150

200

550

product price (RMU/ton)

175

200

225

comparison to a traditional discrete-time formulation using one binary variable for each task starting processing at each machine at each time grid. In practice this means that the same model as in the mixed-time version is used but the production equation (eq 1.3.1) using the continuous unit-to-task variable wi,k,t is replaced by eq 1.1.1 using the discrete variable yi,k,t. A similar discrete-time formulation is used in ref 4. The illustrative examples are purely fictional and consist of a chemical process in two steps. In the first step the raw material enters machine A or B. The raw material is processed at machines A or B and thereafter at machine C or D in the second production step. In illustrative example 1, intermediate storage capacity is available for one intermediate product, while in examples 2 and 3, intermediate storage capacity is available for all intermediate products. In both examples the objective function presented in eq 2.1 is used, and the cost, capacity, and demand data are displayed in Tables 3 and 4. In section 3, the mixed-time model, applied to an actual industrial application, is demonstrated. The

2788

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 10. Flow sheet for illustrative production process 1.

Figure 11. Schematic flow sheet for production process in examples 2-4.

industrial application in section 3 demonstrates, for example, how the aging penalties work in practice. The objective function, presented in eq 2.1, is concerned with the total profit maximization, taking into account raw material price, end product market price, and changeover costs. Pip is the market price for product i, Pir is the raw material price for raw material r, pip,kn,t is the production of product p at machine kn during time grid t, and pir,kn,t is the production of raw material ir at machine kn during time grid t and is further presented in section 1.3. PCO is the changeover cost/penalty and xCO,kn,t is a variable indicating whether or not a changeover is taking place at machine kn during time grid t and is further presented in section 1.3.1. Kr represents the set of machines A and B, and Kp is the set of machines C and D while K represents the whole set of machines A, B, C, and D. Ir represents the set of raw materials a, b, c, and d, and Ip is a set of end products; end product A, B, and C. If additional penalizing of intermediates getting old in storage is needed, the aging penalty of eq 1.3.16 should further be included in the objective function. T

T

T

pi ,k,t - ∑ Pi ∑ ∑ pi ,k,t - PCO∑ ∑ xCO,k,t ∑ Pi k∈K ∑∑ i ∈I t)1 i ∈I k∈K t)1 t)1 k∈K p

p

p

p

p

r

r

r

r

p

(2.1)

2.1. Illustrative Example 1. In this example, the process is assumed to operate under just-in-time (JIT) discipline when it comes to intermediate storages for all but one intermediate. All process equipment is, however, equipped with small feeding containers to prevent raw material shortages due to shorter interruptions in the process at earlier stages and to sustain an even quality of the raw material feed. The feeding container may be utilized as a kind of intermediate storage but not for time periods over one shift. Similar arrangements may be found in many industrial applications and is a requirement for the

mixed-time model to work properly if no intermediate storage is provided for the intermediate products. A schematic flow sheet of the process is given in Figure 10, and the production capacities of the machines are in Table 2. The product recipes and market demand (minimum production requirements) are, additionally, given in Table 2. The process is assumed to operate 24 h a day, 7 days a week. Furthermore, the process is working periodically using a 7 days horizon. The objective function is taking into consideration the maximization of profit subject to raw material costs and end product price. The objective function is, further, penalized by a sequence-independent changeover cost of 1000 RMU (relative money unit) for all changeovers at all machines. It may be worth mentioning that example 1 is a hard problem to solve even using a fairly coarse discretization grid because the product recipes for products B and C require mixing of intermediate products, and intermediate storage capacity is provided only for one intermediate product. 2.2. Illustrative Examples 2 and 3. These illustrative examples consist of the same chemical process as in section 2.1, but intermediate storage capacity is included for all intermediates. The flow sheet of the considered process is displayed in Figure 11. The maximum capacity of the intermediate storages in illustrative example 2 are 100 tonnes and 50 tonnes in example 3. 2.2.1. Computational Results. Illustrative example 1 is used to demonstrate the huge influence the mixed-time formulation has on the objective function value for cases with limited intermediate storage capacity. Illustrative examples 2 and 3 are used to demonstrate the difference in performance using a traditional discrete-time representation and a mixed-time representation in cases where a more flexible intermediate storage capacity is available. Illustrative examples 2 and 3 may be seen as more fair comparisons because the results using a discretetime model are highly dependent on intermediate storage

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2789 Table 5. Computational Results for Illustrative Example 1a time model [time discretization horizon representation] grid [h] [days]

CPU time [s]

objective function value [RMU/week]

MIP gap [%]

discrete discrete discrete

24 12 12

7 7 7

5575.9 1000a

infeasible 71322.9 71284.3

opt 42.13

mixed mixed mixed

24 12 12

7 7 7

46.52 21108.9b 1000a

137891.4 126093.4 125846.2

opt 13.1 16.96

a

CPU time limit reached. b Memory limit reached.

Table 6. Computational Results for Illustrative Example 2 time model [time discretization horizon representation] grid [h] [days]

a

CPU time [s]

objective function MIP value gap [RMU/week] [%]

discrete discrete discrete discrete

12 6 6 6

7 7 7 7

530.7 120a 1000a 10000a

149835.4 148360.4 150414.1 151160.0

opt 4.45 2.35 0.82

mixed mixed mixed mixed

12 6 6 6

7 7 7 7

487.8 120a 1000a 10000a

152235.1 152142.9 152142.9 152142.9

opt 3.10 2.22 1.50

CPU time limit reached.

Table 7. Computational Results for Illustrative Example 3 time model [time discretization horizon representation] grid [h] [days]

CPU time [s]

objective function value [RMU/week]

MIP gap [%]

discrete discrete discrete discrete discrete

24 12 6 6 6

5 7 7 7 7

44.38 10000a 120a 1000a 10000a

138052.6 143972.0 140964.0 143418.0 144858.0

opt 2.79 10.59 7.81 5.70

mixed mixed mixed mixed mixed

24 12 6 6 6

7 7 7 7 7

46.47 10000a 120a 1000a 10000a

149596.8 148491.5 144348.9 146010.3 147980.3

opt 2.19 9.02 7.09 4.88

a

CPU time limit reached.

capacity and time discretization especially in processes where mixing is included. It is worth noticing that no comparison to any continuous-time model is done because the proposed mixedtime model is developed to handle problems where continuoustime models are not applicable. It is, however, fairly clear that most continuous-time models would outperform the mixed-time model, regarding computational performance and solution quality, in cases where, for example, aging profiles are not of interest and continuous-time models are applicable.

In Table 5, the computational results from a couple of test runs of illustrative example 1 are presented. One may observe that the mixed-time formulation is, in some cases, able to come up with feasible solutions to problems which are infeasible using a traditional discrete-time model. This is, obviously, the case when a very coarse discretization grid is used and the production constraints are fairly hard to meet. An interesting fact, shown in Tables 5-7, is that it in many cases appears to be unprofitable to use a fine discretization for the mixed-time formulation. In these cases the computational complexity is increased, resulting in poor computational performance, without resulting in any significant increase of solution quality. Note that a number of computational runs using the same discetization grid and time horizon but different CPU time limits are displayed in Tables 5-7. In cases where the CPU time is highlighted with a superscript a, the CPU time limit was reached. In cases where the CPU time is highlighted with a superscript b, the memory limit was reached and the optimization terminated. Production plans for example 1, using both time representations, are presented in Figures 12 and 13. The discretization (12 h, 7 days) is presented using dotted lines. The Gantt charts are presenting production plans for 14 days of operation for the four different machines included in the production process. In all Gantt charts, a, b, c, and d correspond to processing of raw materials a, b, c, and d while A, B, and C correspond to processing of end products A, B, and C. In Figure 13 one may observe how the mixed-time model works in practice. Product changeovers are taking place continuously at any time, regardless of the time discretization. It is clear from the Gantt charts that the comparison in example 1 is not very fair to the discrete-time model as it does not take advantage of the feeding containers and is, hence, producing highly suboptimal solutions. Example 1, however, illustrates the huge influence on solution quality that a continuous-time formulation may have in many cases. In example 2, three different CPU time limits (120 s, 1000 s, and 10 000 s) have been used to illustrate the convergence speed. It is worth mentioning that the proposed model is intended for industrial applications, where it is of utmost importance to come up with good quality feasible solutions within a reasonable CPU time; global optimality is of less importance. The mixed-time model is clearly very efficient in finding good quality solutions compared to the discrete-time model. The mixed-time model is already in 487.8 CPUs, using a discretization grid of 12 h, reaching a profit that is higher than what is reached in 10 000 CPUs using a discrete-time model and a discretization grid of 6 h. A discrete-time Gantt chart for example 2, using a discretization grid of 24 h, is presented in Figure 14. In the Gantt charts in Figures 13, 16, and 20, one may observe how the mixed-

Figure 12. Gantt chart for illustrative example 1 using a discrete-time model with a discretization grid of 12 h and a time horizon of 7 days.

Figure 13. Gantt chart for illustrative example 1 using the mixed-time model with a discretization grid of 12 h and a time horizon of 7 days.

2790

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 14. Gantt chart for illustrative example 2 using a discrete-time model with a discretization grid of 24 h and a time horizon of 7 days.

Figure 15. Intermediate storage profiles for illustrative example 2 using a discrete-time model with a discretization grid of 24 h and a time horizon of 7 days.

Figure 16. Gantt chart for illustrative example 2 using a mixed-time model with a discretization grid of 24 h and a time horizon of 7 days.

Figure 17. Intermediate storage profiles for illustrative example 2 using a mixed-time model with a discretization grid of 24 h and a time horizon of 7 days.

Figure 18. Gantt chart for illustrative example 3 using a discrete-time model with a discretization grid of 6 h and a time horizon of 7 days.

time model is allowing changeovers to happen continuously within all time intervals. It is worth noticing, though, that the mixed-time model does not allow more than one changeover within each grid. The discretization grid may, hence, be of some importance in cases where many shorter batches are needed. In many cases a finer discretization grid will, however, only deteriorate the computational performance without affecting the solution quality positively. This may be seen, for example, in Table 5. It is, thus, recommended to choose a discretization grid with appropriate length, suiting the specific process under consideration, preferably starting with a coarse grid. In cases where the operation is performed in shifts, setting the discretization grid equal to the shift length is usually convenient if more frequent changeovers are not required. Storage profiles from illustrative example 2 are displayed in Figures 15 and 17. The storage inventory, in tonnes, is displayed on the y-axis, and the time discretization is displayed on the x-axis. Because the problem is assumed to consider a periodic type of operation, the intermediate storage inventories for time

sequence 0 and 7 are the same in Figures 15 and 17. By decreasing the intermediate storage size from 100 tonnes (example 2) to 50 tonnes (example 3), the problem is made more complex, as seen in Table 7. One may observe the difference in intermediate storage utilization between solutions gained with a discrete-time model and the solutions of the mixed-time model. Generally, the pressure on intermediate storage utilization is lower using the mixed-time model because of the fact that the mixed-time model offers a more flexible form of operation. The production schemes (Gantt charts) and intermediate storage profiles for illustrative example 3 are presented in Figures 18-21. The figures are presenting the results using a discretization grid of 6 h and a production cycle of 7 days. The runs, using both time representations, were terminated at a CPU time limit of 10 000 s. In Table 7, one may observe that the mixed-time model already in 46.47 CPUs reached a clearly more profitable solution than what the discrete-time model reached in 10 000 CPUs. It is clear from all runs that

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2791

Figure 19. Intermediate storage profiles for illustrative example 2 using a discrete-time model with a discretization grid of 6 h and a time horizon of 7 days.

Figure 20. Gantt chart for illustrative example 3 using a mixed-time model with a discretization grid of 6 h and a time horizon of 7 days.

Figure 21. Intermediate storage profiles for illustrative example 2 using a mixed-time model with a discretization grid of 6 h and a time horizon of 7 days.

the mixed-time model completely outperforms a traditional discrete-time model concerning computational efficiency. This is due to the fact that the mixed-time model is able to reach good quality solutions using a very coarse discretization grid whereas a traditional discrete-time model is forced to use a much finer grid to even come close to similar profits as the mixedtime model. The illustrative example problems were solved using CPLEX 10 (default settings) on a 3.02 GHz Pentium 4 Windows XP computer. In the following section industrial prototype software, based on the mixed-time formulation, is briefly presented. A brief background of the industrial application in question is presented in a concealed way to avoid leakage of sensitive data. 3. Industrial Software Prototype The mixed-time model was applied to a large industrial application, and a prototype software, mixed-integer strategic planning tool, MISPT, was built based on the proposed model. The MISPT tool was applied to an industrial production plant at the raw material preparation end of the production. The process in question consists of three main production steps 1-3. In the first production step the raw material is processed and fed into intermediate storages. In the second production step the actual product is produced, and in the third production step the product is refined to form different end products. A schematic flow-diagram of the process is presented in Figure 22. All production steps consist of a number of non-identical parallel machines. The MISPT tool was incorporated in the production planning procedure of the factory to handle the production planning of the first production step, and the intermediate storage handling between production steps 1 and 2. Production steps 2 and 3 were already incorporated in the corporate ERP system and

Figure 22. Schematic flow diagram of the industrial production process.

were, hence, not incorporated in the MISPT tool. The production plans of production step 2 were, consequently, used as input data to the MISPT tool and the production scheduling of production step 1 and the intermediate storage between step 1 and step 2. The actual production planning optimization was, hence, concerned with cost minimization related to changeovers at production step 1 and minimization of preservative chemical costs related to the aging in intermediate storage. The input data from production step 2 contains batch size, duration, and product recipe. The production at production step 2 contains around 100 different products with different product recipes, containing one to five raw materials each. Certain raw materials, processes in production step 1, may be substituted with pre-processed “market raw material” which require significantly less processing. Market raw material is processed exclusively at machine 1d. Machines 1a and 1b are able to handle the same raw materials, RM1, RM2, and RM3, while machine 1c handles RM4 and RM5. A more detailed flow diagram of production steps 1 and 2 is presented in Figure 23. Intermediate storage IS1 is exclusively used for RM1, IS2 is exclusively used for RM3 while IS3 is a shared tank for RM1 or RM3. IS4 and IS5 are exclusively used for RM2, IS6 is used for RM4, and IS7 is used for RM5. All raw materials have limited preservation properties and will get downgraded, or even completely spoiled, after a certain time in storage. Because of this fact, preservatives are added to the intermediate storages, adding costs to the production process. Aging penalties/costs

2792

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 23. Flow diagram of the industrial production process steps 1 and 2.

Figure 24. Screen shot of the main window of MISPT.

are added to the objective function using the chemical costs for preservatives as input. The aging penalties are discussed in section 1.3.2. In the following section, an illustrative example using the MISPT tool is presented. 3.1. Illustrative Example. The prototype software, MISPT, was customized to the needs of the actual end-users, which in practice means that the focus was more on usability than on the theoretical aspects of the mathematical model. The previous chapters have shown that the mixed-time formulation utilizes some features from a traditional discrete-time model, such as ability to model storage profiles, making it applicable to problems were continuous models are not applicable. These features were of great relevance when formulating the models

used in the MISPT tool. A screenshot of the MISPT graphical user interface is displayed in Figure 24. The demand in production step 2 is loaded from the local ERP system, and MISPT converts the production program/ demand automatically into discrete-time grids. The discrete demand will always be an approximation of the real world, because of delays in the industry reporting systems. The length of the discrete-time grids can be set to 2, 4, or 8 h. It is worth noticing that the problems may differ with different length of the grids. Problem formulations with more dense discretization are more sensitive to rapid changes in demands which may affect the intermediate storage profiles but also the production plan in the production step 1. That is why the problems with different

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2793 Table 8. Solution for Machines 1a and 1b and Storages IS1-IS5 hours

days

time (s)

MIP gap

value (RMU)

8 4 2

7 7 7

7 28 635

opt opt opt

221689.36 229585.16 232659.46

Table 9. Solution for Machine 1c, IS6 and IS7 hours

days

time (s)

MIP gap

value (RMU)

8 4 2

7 7 7

5 122 900a

opt opt 0.72

411160.97 423584.34 430204.06

a

CPU time limit.

grid length often have to be considered as separate problems and not better approximations of the same production program. Further discussion about the optimal length of the grids in this application is beyond the scope of this paper. The real-world problem in this example is divided into two separate problems. The first contains machines 1a and 1b and IS1 to IS5 (Table 8). The other one contains machine 1c and IS6 and IS7 (Table 9). This is done automatically in the MISPT tool. The two first machines are independent from the third, and both groups of storages are separated from each other; see Figure 23. Each problem is solved separately, which gives shorter optimization times, because of smaller problems. Fast solution is more desirable than global optimal solutions in this particular industrial application. The time limit was set to 15 min in this example. The problems were solved using a AMD Athlon 64 3200+, 1024 MB, and ILOG CPLEX 10.0. Input data, such as storage levels at the start and the discrete approximation of the demand, have a tendency to have major influence on solution times. In the situations where the demand is more evenly spread over the whole optimization horizon, the optimization is significantly more difficult. It is, however, possible to modify the model parameters in the application to emphasize some special desired features. For example, the properties for machines 1a and 1b may differ from those of machine 1c, and it may be desirable

Figure 25. Production schedule using an 8 h grid.

Figure 26. Production schedule using a 4 h grid.

to treat them differently. The production at machine 1a and 1b includes lower production capacities, and changeovers are also more expensive in real life compared to in machine 1c where production capacity is not a bottleneck. Normally it is not necessary to change the parameters, because the default settings are good enough for normal daily situations. The objective function value includes, for example, raw material costs, changeover penalties, and preservation chemical costs caused by “too old” of a substance in storage. The user can change the length of the grids, as mentioned before, to 2, 4, or 8 h. The schedules at production step 1 often appear to be the same, even though the strategies for the storages may differ slightly, especially when several storages can be used for the same raw material and the aging horizon is very tight. In the rare case of problems with a very dense discretization grid and where changeovers are acceptable in intervals similar to the length of the discretization grid, the mixed-time formulation does not offer any greater benefits than a normal discretetime formulation. In cases where a very fine discretization grid is used, the optimization may leave unused grids to circumvent the changeover penalty due to the fact that the objective function is penalized by changeover penalties only for adjacent time intervals. The changeover penalty procedure could easily be changed to also penalize cases where one grid is left “empty” to prevent changeover penalty. In this application, however, it was found more practical to use the changeover procedure described in section 1.3.1. In the MISPT tool, a grid length of 8 h is used as the default as it is very handy in most processes using continuous work shifts of 8 h length. The first grid will in most cases be shorter than the following ones, because the grids are automatically synchronized to the start time of the next predefined work shift. Figures 25-27 illustrate production schedules for machine 1b, using the same approximated demand at production step 2 (Figure 22) with different grid lengths. Another useful feature in the MISPT tool is the piecewise linear representation of highly nonlinear storage profiles and especially the aging profile of the storages. The user can directly

2794

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

Figure 27. Production schedule using a 2 h grid.

Figure 28. Storage profile and aging profile (darker bars), 4 h grid, aging horizon 1 day. Table 10. Storage Inventory and Flow into and out of Storage. grid

flow in

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

4.9 19.4 19.4 19.4 19.4

flow out

4.8 11.0 11.0 10.0 1.4 1.4 1.4 16.1 16.1 14.9 11.1 11.1

storage inventory

old substance

32.7 52.1 47.3 36.3 44.8 54.2 72.2 70.8 69.3 69.3 69.3 53.2 37.1 22.2 22.2 11.1 0.0 0.0

7.0 7.0 2.2

11.0 11.0 11.0 14.3 17.7 22.2 22.2 11.1

see where and when there is “too old” of a substance in the intermediate storages. The user can define, in reference to, for example, the time of year, the time before the old substance will be penalized in the objective function discussed in chapter 1.3.2. The penalty for each product can be set separately. It is worth mentioning that the time before substances are defined as “too old” is defined in days, and the number of grids taken into account is translated according to the grid length. More dense grids have a tendency to be more sensitive to old substances, because the storage levels are always updated at the end of each grid. An example of a storage profile and the aging profile are shown in Figure 28. The underlying data of Figure 28 is represented in Table 10. In the actual industrial software, the user is also provided with

report tables where all flows of materials and amounts of old substances are given in numbers for each grid. 4. Conclusions In this paper a novel mixed-time formulation for large scale industrial scheduling problems is presented. The novelty with the model is that it is introducing a continuous-time representation into a discrete-time grid, keeping track of the time dimension. In practice this means that changeovers, at all machines, do not have to happen at discrete time points and that the model is using a discrete reference grid for the time dimension, practicable for storage profiles over time and so forth. The model is proven efficient in multi-stage applications using intermediate storages but is, however, not applicable to cases where the process is assumed to operate under a strict JIT discipline. Another important feature with the proposed mixed-time formulation is the ability to handle strongly nonlinear storage profiles and aging aspects of intermediates in storage. In addition to the ability to handle the time dimension in various practicable ways, the mixed-time model requires a significantly less dense discretization grid to come up with good quality solutions compared to a traditional discrete-time model, in addition to the fact that changeovers may happen continuously. This is a major advantage with the model as it enables solution of large scale problems, to similar objective function values, significantly faster than a traditional discrete-time model. The model is, however, limited in the same way as a traditional discrete-time model by the computational complexity combined with a rising number of binaries. The model should, consequently, not be applied to problems where a continuous-time model is applicable because the solution efficiency of the mixed-

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007 2795

time model is significantly worse compared to any recent continuous-time model, applied to problems solveable using continuous-time models. It is, additionally, important to notice that the mixed-time model does not offer a purely continuoustime representation. Critical production events, such as product changeovers, may happen only once within every discrete-time sequence. Overall the proposed model offers a simple, yet effective, way of modeling scheduling problems where aging aspects in intermediate storages is of concern. The model, additionally, challenges traditional discrete-time models on computational efficiency also in problems where no aging profiles in intermediate storage are required. The benefits of the proposed model are demonstrated through a brief introduction to actual industrially used software based on the proposed model. Work still remains to make the proposed model more general; the model as it is presented in this paper offers, however, a practicable platform for industrial production planning and scheduling problems with advanced intermediate storage requirements. Nomenclature Indices i ) process task k ) unit t ) time interval Parameters pi,k,max ) maximum momentary production rate i at unit k for any time interval pii′,k,CO ) changeover time from i to i′ at unit k pold,s ) substance-specific aging penalty told ) age limit before substance getsis getting old in intermediate storage ∆t ) length of time interval t Variables Cold,tot ) total aging penalty mold,s,t ) storage inventory of “old substance” in storage s in the beginning of interval t ms,in,t ) sum of substances into storage s during time interval t ms,out,t ) sum of substances out of storage s during time interval t ms,t ) storage inventory for storage s in the beginning of time interval t pi,k,t ) production of i at unit k during interval t wi,k,t ) fractional upper limit of maximum production for momentary production i at unit k at interval t xCO,ii′,k,t ) continuous sequence dependent changeover variable equal to 1 if a change over from i to i′ takes place at k during interval t xCO,k,t ) continuous sequence independent changeover variable equal to 1 if a changeover takes place at k during interval t yi,k,t ) binary variable equal to 1 if task i starts/continues at unit k in the beginning of interval t ∆ti,k,t ) changeover time for I at unit k during time interval t Literature Cited (1) Me´ndez, 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. (2) Floudas, C. A.; Lin, X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Comput. Chem. Eng. 2004, 28, 2109-2129.

(3) Grossmann, I. E.; Sargent, R. Optimal Design of Multipurpose Chemical Plants. Ind. Eng. Chem. Process Des. DeV. 1979, 18, 343348. (4) Kondili, E.; Pantelides, C.; Sargent, R. W. H. A general algorithm for short-term scheduling of batch operations - I. MILP formulation. Comput. Chem. Eng. 1993, 17, 211-227. (5) Shah, N.; Pantelides, C. C.; Sargent, R. Optimal Periodic Scheduling of Multipurpose Batch Plants. Ann. Oper. Res. 1993, 42, 193-228. (6) Pantelides, C. C. Unified Frameworks for the optimal process planning and scheduling. Proceedings of the second conference on foundations of computer aided process operations; Cache Publications: New York, 1994; p 253. (7) 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. (8) Ierapetritou, M. G.; Floudas, C. A. Effective continuous time formulation for short-term scheduling. 3. Multiple intermediate due dates. Ind. Eng. Chem. Res. 1999, 38, 3446-3461. (9) Castro, P. M.; Barbosa-Po´voa, A. P.; Matos, H. A. An Improved RTN Continuous-Time Formulation for the Short-term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2001, 40, 20592068. (10) Castro, P. M.; Barbosa-Po´voa, A. P.; Matos, H. A.; Novais, A. Q. Simple Continuous-time Formulations for Short-term Scheduling of Batch and Continuous Processes. Ind. Eng. Chem. Res. 2004, 43, 105-118. (11) Giannelos, N. F.; Georgiadis, M. C. A Novel Event Driven Formulation for Short-term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2002, 41, 2431-2439. (12) Gupta, S.; Karimi, I. A. An improved MILP formulation for scheduling multiproduct, multistage batch plants. Ind. Eng. Chem. Res. 2003, 42, 2365-2380. (13) Maravelias, C. T.; Grossmann, I. E. New General Continuous-time State Task Network Formulation for Short-term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2003, 42, 3056-3074. (14) Pinedo, M. Scheduling: Theory, Applications and Systems; Prentice Hall: Englewood Cliffs, NJ, 1995. (15) Shah, N. Single- and multisite planning and scheduling: Current status and future challenges. Presented at the second conference on foundations of computer aided process operations, 1998, Snowbird, UT, U.S.A. (16) Castro, P. M.; Grossmann, I. E. Multiple Time Grid ContinuousTime Formulation for the Short Term Scheduling of Multiproduct Batch Plants. Proceedings of the 16th European Symposium on Computer Aided Process Engineering and the 9th International Symposium on Process Systems Engineering; Elsevier: Amsterdam, 2006; pp 20932098. (17) Maravelias, C. T. Mixed-time Representation for State-Task Network Models. Ind. Eng. Chem. Res. 2005, 44, 9129-9145. (18) Castro, P. M.; Westerlund, J.; Forssell, S. Optimal Periodic Scheduling of a Paper Mill with Recycling of Byproducts. Comput. Chem. Eng. 2006, submitted. (19) Shaik, M. A.; Janak, S. L.; Floudas, C. A. Slot-Based vs. Global Event-Based vs. Unit-Specific Event-Based Models in Scheduling of Batch Plants. Proceedings of the 16th European Symposium on Computer Aided Process Engineering and the 9th International Symposium on Process Systems Engineering; Elsevier: Amsterdam, 2006; pp 1923-1928. (20) Westerlund, J.; Castro, P.; Forssell, S. Strategic planning and design using MILP: an industrial application from the tissue manufacturing industry. Proceedings of the 16th European Symposium on Computer Aided Process Engineering and the 9th International Symposium on Process Systems Engineering; Elsevier: Amsterdam, 2006, pp 2087-2092. (21) Maravelias, C. T.; Grossmann, I. E. On the Relation of Continuousand Discrete-Time State-Task Network Formulations. AIChE J. 2006, 52, 843-849. (22) 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-3051. (23) 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, 2679-2702. (24) Cerda´, J.; Henning, G. P.; Grossmann, I. E. A mixed-integer linear programming model for short-term scheduling of single-stage multiproduct batch plants with parallel lines. Ind. Eng. Chem. Res. 1997, 36, 16951707.

2796

Ind. Eng. Chem. Res., Vol. 46, No. 9, 2007

(25) Me´ndez, C. A.; Henning, G. P.; Cerda´, J. Optimal scheduling of batch plants satisfying multiple product orders with different due-dates. Comput. Chem. Eng. 2000, 24, 2223-2245. (26) Me´ndez, C. A.; Henning, G. P.; Cerda´, J. An MILP continuoustime approach to short-term scheduling of resource-constrained multi-stage flowshop batch facilities. Comput. Chem. Eng. 2001, 25, 701-711. (27) Me´ndez, C. A.; Cerda´, J. An MILP continuous-time framework for short-term scheduling of multipurpose batch processes under different operation strategies. Optimization and Engineering 2003, 4, 7-22.

(28) Me´ndez, C. A.; Cerda´, J. An MILP framework for batch reactive scheduling with limited discrete resources. Comput. Chem. Eng. 2004, 28, 1059-1068.

ReceiVed for reView July 28, 2006 ReVised manuscript receiVed January 8, 2007 Accepted February 16, 2007 IE060991A