Optimal Periodic Scheduling of Multistage ... - ACS Publications

Ind. Eng. Chem. Res. 2007, 46, 3669-3683

3669

Optimal Periodic Scheduling of Multistage Continuous Plants with Single and Multiple Time Grid Formulations Pedro M. Castro* and Augusto Q. Novais Departamento de Modelaç a˜ o e Simulaç a˜ o de Processos, Instituto Nacional de Engenharia, Tecnologia e InoVaç a˜ o, 1649-038 Lisboa, Portugal

This paper presents a new mixed integer nonlinear program (MINLP) model for the periodic scheduling of multistage, multiproduct continuous plants featuring equipment units in parallel that are subject to sequence dependent changeovers. The formulation is based on the resource task network (RTN) process representation, features combined processing and changeover tasks, and assumes that each order is executed only once on each stage. The new multiple time grid formulation is compared to an also RTN-based, single time grid formulation, through the solution of a few example problems taken from the literature. The results show that the proposed formulation is significantly more efficient computationally, essentially because smaller problems are generated and, because the number of tasks to execute can easily be predicted, no iterative procedure over the total number of event points in the time grid is required to find the optimal solution. On the other hand, the single time grid formulation is able to sometimes find better solutions to the problem because of its ability to consider, for each order, the execution of more than a single instance of the processing task per stage. 1. Introduction Extensive reviews of optimization approaches for scheduling have recently appeared in the literature.1,2 These focus mostly on batch processes, which have received considerable attention, whereas much less work has been reported in the scheduling of continuous plants despite their practical importance. These arise frequently in chemical processing industries and petroleum refineries and polymer and specialty chemical plants as well as paper machines. Furthermore, one of the trends in the chemical industry is to move toward continuous flexible multiproduct plants that can respond more quickly to demand changes and that can handle a large product portfolio. The short-term mode of operation, which is appropriate in plants where product demands are subject to rapid change, has also been preferred over a cyclic mode. Periodic scheduling is suitable wherever product demands are stable over extended periods of time. It can also be used to address noncyclic scheduling problems whenever a long time horizon is involved. In such problems, a production pattern typically exists that, when periodically repeated, achieves production very close to the optimal one. The advantage of periodic scheduling is that the size of the mathematical problem is reduced to a smaller one as a result of the consideration of a smaller time horizon and, hence, is solved much faster. The drawback is that periodic problems are generally formulated as mixed integer nonlinear programs (MINLPs) whereas short-term problems are typically of the mixed integer linear program (MILP) type which, for similar sizes, are easier to solve. Meńdez et al.1 have identified time representation as the most important issue for the classification of optimization models. While no substantial developments have occurred with discretetime formulations in recent years, continuous-time formulations have been in the spotlight in the last 10 years. A few research groups have presented their models, and the process systems engineering community is still working on finding out which * To whom correspondence should be addressed. Tel.: +351210924643. Fax: +351-217167016. E-mail: [email protected].

model is the best approach for a particular type of problem. Comparative studies have recently appeared where several of the formulations have been implemented and solved using the same software and hardware. Shaik et al.3 have considered five different continuous-time formulations for the short-term scheduling of multipurpose batch plants, while Castro and coworkers4,5 have compared six conceptually different approaches (one discrete-time and three continuous-time formulations and one constraint programming (CP) and one hybrid MILP/CP model) for the short-term scheduling of single/multistage batch plants. Time grid based continuous-time models can be divided into single and multiple time grid formulations. Shaik et al.3 have identified the formulations of Janak et al.6 as the most computationally efficient continuous-time formulation and the one of Castro et al.7 as the best single time grid formulation for the short-term scheduling of multipurpose batch plants. Both of these formulations have the drawback of requiring a few iterations over the total number of event points in the grid(s) before the optimal solution is found. The latter has the additional advantage of being easily adapted to handle periodic scheduling problems with continuous tasks.8 The periodic formulation that is closest to the former is the one by Wu and Ierapetritou9 but it is restricted to batch plants. For problems with a simpler network structure, multistage, multiproduct plants, Castro et al.4 have recently presented two very efficient, different multiple time grid formulations for problems involving sequence dependent changeovers, for which specifying the total number of event points that lead to global optimal solutions is not an issue. However, the models did not take into consideration the amount of material processed by the several tasks nor inventory profiles. The main motivation for the work included in this paper was to address these issues and to adapt the formulation to periodic problems, a process that was in the past accomplished with relative ease for the single time grid formulation.7 The definition of the periodic scheduling problem involved in multistage multiproduct plants is taken from Pinto and Grossmann10 together with most of the data for the example problems.

10.1021/ie0613570 CCC: $37.00 © 2007 American Chemical Society Published on Web 04/10/2007

3670

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

Figure 1. Multiproduct continuous plant with parallel equipment units and dedicated intermediate storage units.

The problem described by Pinto and Grossmann10 consisted of a single equipment unit per stage and was addressed by a multiple time grid continuous-time formulation. The new formulation that is proposed in this paper considers the more general case of parallel units within stages and, thus, can handle more flexible processing plants. It also has the additional advantages of not being restricted to the same product sequence on every stage and allowing continuous tasks to be performed below their maximum processing rates, whenever storage costs can be reduced. These three features are common to the single time grid formulation,8 which will be used for comparison. The new formulation considers the cycle time as a model variable and gives rise to a non-convex MINLP, like the one of Pinto and Grossmann.10 The MINLPs will be solved mostly by the local optimization solver GAMS/DICOPT, which does not guarantee global optimality of the solutions returned. Global optimization solvers are computationally more demanding, so comparatively, their practical use is restricted to mathematical problems that are lower in size. In this work, GAMS/BARON has been used as the global optimization solver and proved its usefulness in a few problem instances. Global optimization has been applied by Alle and Pinto,11 for a formulation similar to that of Pinto and Grossmann,10 and by Wu and Ierapetirou.9 2. Problem Definition In this paper, the periodic scheduling problem of multistage, multiproduct continuous plants is considered. Given are a set of orders i ∈ I that must follow a sequence of processing stages k ∈ K, to reach the condition of final products; a set of available equipment units m ∈ M, each allocated to a single stage, with set Mk including the units belonging to stage k. The plant features intermediate storage tanks for each product and for each stage of production (see Figure 1). The continuous processing tasks are characterized by maximax mum processing rates Fi,m (kg/h), which will normally be met unless operating at lower rates will reduce storage costs. Transition times arise between the processing of two successive products on the same unit and are sequence dependent, cli,i′,m (h). The price of the final products, which are subject to a minimum demand rate di (kg/h), are given by ci ($/kg). Intermediate storage costs cii,k ($/kg) are proportional to the maximum amount of material that needs to be stored on that particular stage, which is directly related to the capacity of that particular storage tank. Thus, with respect to the intermediate storage tanks, the scheduling model is also incorporating decisions related to equipment design. Final (with respect to the last stage) storage costs are assumed to be dependent both on the amount stored and the storage time ($/(kg‚h)) and are given by cfi. Finally, the cost for transitions is represented by ctk,i,i′ ($/h). The objective will be to find a cyclic schedule that maximizes the profit (product revenue minus storage and transition costs). More specifically, the problem consists in determining the (i)

cycle time; (ii) assignment of orders to equipment units; (iii) sequence of orders allocated to the same equipment unit; (iv) length of processing and cleaning tasks; (v) processing rate of continuous tasks; (vi) amounts of products to be produced; and (vii) levels of intermediate storage and final product inventories. 3. Single Time Grid Formulation (STG) The single time grid periodic formulation used in this paper is conceptually the one presented in Castro et al.8 There are, however, some differences that result from (i) the disaggregation of the different resource types from the general concept of resource given by the Resource Task Network process representation;12 (ii) the disaggregation of the general concept of task into processing and cleaning tasks; and (iii) the need to consider a continuous consumption (for accurate calculation of storage costs), at an unknown rate, of the final products. The former changes are essentially nomenclature adjustments caused by the consideration of a specific multistage network structure instead of the more general multipurpose type. As a consequence, the structural parameters are no longer required and the formulation becomes easier to comprehend. Time grid based formulations rely on one or more time grids with a pre-specified number of time points (also known as event points). This model uses a single time grid with a number of nonuniform time intervals that are specified through set T. The quality of the solution returned by the model is highly dependent on the cardinality of T, and normally, optimality can only be achieved with a sufficiently high value of |T|, which in practice involves an iterative search procedure.3,8 The maximum number of time intervals that tasks can span has also been shown13 to have a significant influence on computational effort. However, this only affects batch tasks because we can assume without loss of generality that continuous tasks only span one time interval. In this work, the continuous tasks are the processing tasks whereas the cleaning tasks are of the batch type. Because the former will last longer than the latter, it is assumed that the cleaning tasks can only span one time interval. As a consequence, all model variables characterizing the extent of tasks feature a single time index t, which give their starting event point. 3.1. Variables. The formulation features two sets of binary variables: Ni,m,t and N h i,i′,m,t. The former identify the processing of order i in unit m starting at event point t, while the latter identify the execution of the changeover task between orders i and i′. The amount of material processed is given by ξi,m,t. The time of each event point, relative to the start of the cycle, is given by Tt, while H gives us the cycle time. Other positive continuous variables include the excess resource variables. In our recent work,4 these variables have been disaggregated into three different types, although only two are strictly necessary. Variables Ci,m,t represent the excess availability of equipment unit m at a condition that enables it to process order i at event

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3671

Figure 2. Possible solution from single time grid (STG) formulation (|I| ) 3, |M| ) 2, |K| ) 2, and |T| ) 7).

point t. Basically, Ci,m,t is equal to 1 whenever no processing task (i,m) or cleaning task (i,i′,m) is being performed at interval t and is equal to 0 otherwise. It is important to emphasize that, although not necessary, because the model constraints ensure that Ci,m,t ) {0,1}, they could also be defined as binary variables. Also, the sum over i is linked to the availability of the equipment resource, for which variables Rm,t could be used explicitly.4 The other set of excess resource variables, Si,k,t, give the amount of material from order i at a state corresponding to the output from stage k, at event point t. As a result of the assumption of dedicated storage units, these amounts will be available in storage, so variables Si,k,t will be used to generate the storage profiles. Then, the maximum amount over time of order i in stage k will set the required capacity of the corresponding storage vessel, Vi,k. Finally, we will be using ∆i for the total production over the cycle of product i. 3.2. Conceptual Overview. The most important features of the single time grid formulation are illustrated in Figure 2, with a simple periodic schedule. The time grid consists of seven event points, where the first, due to the cyclic nature of the problem, acts both as its origin (T1 ) 0) and end (H, represented as a dashed line). The minimum number of event points that achieves feasibility is equal to the maximum number of tasks that need to be performed over all stages, which in this case is equal to 6 (3 processing + 3 cleaning tasks). This property derives from the fact that the duration of both types of tasks can be adjusted to match the duration of the several time intervals, set by other, more limiting tasks. For example, the duration of the first time interval is set by order I2 in M2, which is performed at its maximum processing rate (in the figure, tasks are performed at their maximum rates unless otherwise specified). Order I1 in M1 is executed over the same interval but is performed at a rate below its maximum value. Another example occurs in the sixth time interval. Order I3 is processed at the same rate in both stages (FI3,K1 ) FI3,K2) but while it is processed at the maximum rate in unit M1 (stage K1) the same is not true for max unit M2 (note this is the same as saying that Fmax I3,M2 > FI3,M1 because an amount of I3 is processed in the two stages). This allows for a possible reduction in storage costs while keeping the required number of event points at a minimum. The corresponding behavior for cleaning tasks can be seen in intervals 5 and 7, where the required changeover time, respectively clI1,I3,M2 and clI3,I2,M2, is actually lower than their duration, which is set by the higher clI2,I3,M1 and clI3,I1,M1 values. Any given order can be processed more than once in a particular stage (e.g., I2 in M1 and M2). No changeover is required in such cases and therefore changeover tasks featuring i ) i′ can be eliminated from the formulation. It is important to emphasize that whenever there are two instances of the same task being executed in the same unit we cannot relate their processing rates based solely on the length of the tasks. For instance, the two occurrences of I2 in M2 are executed both at

their maximum rates, but the first occurrence lasts much longer simply because a larger amount of material is being processed. Unlike short-term scheduling,4,5 where the execution of an order in stage k must be preceded by its completion in stage k - 1, in periodic scheduling the beginning and end of the time grid coincide. Thus, an order in stage k can also start at a lower event point than it does in stage k - 1 (e.g., I2 is processed twice in M2, at t ) 1 and t ) 2, while in M1 it starts at t ) 3 and t ) 4). This possibility simply means that intermediate longer storage will be required. Finally, the model allows for the use of different order sequences on different stages, as can be seen for M1 (I1-I2-I3) and M2 (I2-I1-I3), which allows us to achieve higher profits whenever the least total changeovers sequence changes from one stage to the other. 3.3. Constraints. We already mentioned that the first time point of the periodic time grid defines the boundary between consecutive cycles, which by definition have the same properties. This means that tasks starting at the last event point can be viewed simply as wrapping around to the start of the cycle when considering the mass balances over the several resources. The following wrap-around operator14 greatly facilitates the formulation of some of the model constraints.

{

t - |T|, t > |T| Ω(t) ) t, 1 e t e |T| t + |T|, t < 1

(1)

The excess resource balances given in eqs 2 and 3 are multiperiod balances where resource availability at event point t is related to that at the previous event point and increases by the occurrence of the processing/cleaning tasks starting at t 1 and decreases by those starting at t. To get an idea of the relation between the excess resource variables Ci,m,t and the binary extent variables Ni,m,t, it can be said that whenever Ni,m,t ) 1, Ci,m,t ) 0. Note in eq 2 that processing tasks consume and produce the same equipment condition while changeover tasks involve dissimilar conditions. In eq 3, the material at a state associated to stage k increases by processing tasks performed in units belonging to stage k and decreases by those performed in units belonging to stage k + 1. Note that, for the last stage, we are assuming that the final products leave the dedicated storage units at a constant rate equal to the delivery/demand rate. Thus, to determine the amount of product i that leaves the system during interval t - 1, we need to know its total production over the cycle, ∆i, the cycle time, H, and the interval length, through the times corresponding to event points t and t - 1. Beause of the third term on the righthand side (RHS), eq 3 is a nonlinear constraint.

Ci,m,t ) Ci,m,Ω(t-1) + Ni,m,Ω(t-1) - Ni,m,t +

∑ (Nh i′,i,m,Ω(t-1) -

i′ ∈ I i′*i

N h i,i′,m,t) ∀ i ∈ I, m ∈ M, t ∈ T (2)

3672


Si,k,t ) Si,k,Ω(t-1) +

∑

m∈M

( [

ξi,m,Ω(t-1)

|

m∈Mk

- ξi,m,Ω(t-1)

|

-

m∈Mk+1

]|

∆i(Tt|t*1 + H|t)1 - TΩ(t-1)) H

)

be used for other variables, have some influence in the search procedure and consequently on the quality of the solution returned, because we are dealing with a non-convex MINLP that will be solved by local optimization solvers.

k)|K|

T1 ) 0; Tt e Hmax ∀ t ∈ T; Hmin e H e Hmax

∀ i, k, t (3)

To ensure that we are not generating a schedule that requires a start-up procedure that uses more equipment units than those available at the plant, eq 4 is used. It breaks the cycle by considering the excess resource balance at the boundary with terms related to tasks starting at t ) 1 only and not to those related to tasks ending at H.

Ci,m,1 ) 1 - ∑Ni,m,1 - ∑ ∑N h i,i′,m,1 ∀ m ∈ M ∑ i∈I i∈I i∈I i′∈I

(4)

i′*i

The timing constraints given in eq 5 are written for every equipment resource m to take advantage of the fact that at most a single task can be allocated to a particular unit at any given time. It is assumed that, for any order, the equipment units can be operated below their maximum processing rates and that the changeover tasks can be made to last longer than the prespecified durations, although such durations can easily be enforced by an additional set of constraints.8

(ξi,m,t/Fmax ∑ i,m + i∈I N h i,i′,m,tcli,i′,m) ∀ m ∈ M, t ∈ T ∑ i′∈I

Tt+1|t*|T| + H|t)|T| - Tt g

The next set of constraints is not mandatory but was found to improve the computational performance of the model. Equation 10 is a classic14 anchoring of the schedule constraint aimed at reducing solution degeneracy arising from cyclic schedule permutations. We have selected to assign the execution of order one in the first stage, to the first event point. Bear in mind, however, that in some problems featuring parallel units, it is theoretically possible that the optimal solution involves the production of the same order in two or more units belonging to the first stage and starting at the first event point. In such a case, eq 10 would remove such solutions from the feasible space and an alternative task would have to be selected to anchor the schedule. Equation 11 states that the maximum number of changeovers per stage must be lower than the number of orders, which is reasonable because increasing the number of changeovers beyond the minimum value decreases the productivity of the plant. Finally, eq 12 forces all equipment units to be occupied at all times, through the excess resource variables Ci,m,t. This constraint is not limiting in any way because equipment idle time can be achieved with Ni,m,t ) 1 and ξi,m,t ) 0 (see eqs 2, 5, and 6).

∑

(5)

max max H Ni,m,t ∀ i ∈ I, m ∈ M, t ∈ T ξi,m,t e Fi,m

(6)

The amount of material that is available at a particular event point must be lower than the capacity of the corresponding storage vessel; see eq 7. This set of variables will be present in the objective function but are not required for the last processing stage.

Si,k,t e Vi,k ∀ i ∈ I, k ∈ K, k * |K|, t ∈ T

(7)

An important problem constraint is that the minimum demand rates must be met for all product orders. Equation 8 is a linear constraint that ensures that the total production of product i over the cycle, ∆i (a model variable), is greater than the minimum amount that needs to be produced, which in turn is given by the product of the minimum demand rate di (a parameter of the model) and the cycle time H. Recall that ∆i is determined through eq 3.

∆i g diH ∀ i ∈ I

(8)

The following set of constraints (see eq 9) define the time associated to the first event point (see also Figure 2) and place appropriate bounds on the timing variables and on the cycle time. Note that the upper and lower bounds on the cycle time are defined a priori and, like all other bounds that may eventually

N1,m,1 ) 1

(10)

m∈Mk)1

i′*i

Equation 6 states that order i can only be processed in unit m at event point t if the corresponding task is active. For the upper bound we have used the maximum amount of material that can be processed over the cycle, which is given by the product of the maximum processing rate by the maximum cycle time.

(9)

∑ ∑ ∑ ∑ Nh i,i′,m,t e |I| ∀ k ∈ K

(11)

Ci,m,t ) 0 ∀ i ∈ I, m ∈ M, t ∈ T

(12)

m∈Mk i∈I i′∈I t∈T i′*i

The MINLP model is complete with the objective function (eq 13), which maximizes the hourly profit of the plant. The first term on the numerator reflects the revenue due to product sales, the second gives the total changeover cost, and the last two give the total storage costs for the intermediates and for the final products, respectively.

[

∑c ∆ - ∑ ∑ ∑ ∑ ∑ct

max

i

i

i∈I

k∈K m∈Mk i∈I

h i,i′,m,t k,i,i′cli,i′,mN

i′∈I t∈T i′*i

-

∑ ∑ ci

i,kVi,k

k∈K i∈I k*|K|

H

∑ ∑ cf /2(S i

i∈I

i,|K|,t

]

-

+ Si,|K|,Ω(t+1))(Tt+1|t*|T| + H|t)|T| - Tt)

t∈T

H

(13)

4. New Multiple Time Grid (MTG) Formulation The new continuous-time formulation is conceptually similar to formulation CT4I developed by Castro et al.4 It uses multiple time grids, one for each equipment resource, four-index binary variables to identify the execution of combined processing and changeover tasks, and is suitable for multistage plants. However, CT4I is a short-term scheduling formulation that does not handle the material resources involved in order processing over the stages. Furthermore, it is suitable for batch plants while the new formulation deals with continuous plants. The current problem features force us to rigorously model inventory profiles for the intermediate and final product resources. Intermediate material


Figure 3. Possible multiple time grid for MTG (|T| ) 3, |M| ) 3).

profiles result from inputs from equipment units belonging to the stage where the material is produced and outputs from the subsequent stage. Thus, they are resources that are shared by a few equipment units. Handling of shared resources by multiple time grid formulations is particularly tricky,1,6,10,15 and a large set of constraints is usually involved. The proposed method has little resemblance to the one embedded in the other multiple time grid formulation10 that has tackled problems of the type considered here. 4.1. Handling of Time. The new formulation uses |M| time grids with a number of pre-specified time points, |T|. Contrary to the single time grid formulation, where an iterative procedure is required to find the optimal value, the cardinality of T can easily be set. Similarly to CT4I,4 all tasks only need to span a single time interval, which means that the value that ensures optimality is |T| ) |I|. Again we will be using the same number of time points in all grids, although it is straightforward to do otherwise. Figure 3 shows a simple example of how the multiple time grids may look like for a problem involving three orders and three equipment units. Its purpose serves to illustrate the features related to the relation of time points of one grid to those of different grids. The most important characteristic is that each time grid is totally independent from the next. This means that it is irrelevant if two grids belong to the same stage or not. This concept, when applied to multistage problems, is new to the best of our knowledge, because there are always constraints relating the timing variables of time grids belonging to consecutive stages of production.6,10 As a direct consequence of the novel modeling technique, all time grids can be placed inside the same boundaries, [0, H], whereas in the model by Pinto and Grossmann,10 both upper and lower bounds normally increase from one stage to the next. The other important aspect is that no task is forced to end exactly at the upper bound, H. Tasks starting at the last time point may continue past the end of the cycle, which is the same as saying that they are wrapped around to the beginning of the same cycle and end at the first time point, which, unlike in the single time grid (see Figure 2), does not need to be located exactly at the boundary (see grid linked to M2 in Figure 3). However, for one time grid, it is convenient to set the time value of the first time point to zero to reduce the number of degenerate solutions without reducing the feasible space. Any time grid can be chosen, and we have selected the one associated to the first equipment unit, M1. In the case that there are fewer than the maximum number of tasks to be performed on a given machine, one or more time points will not be required, meaning that in the optimal solution there may be two points from that grid with the same time value, as illustrated for M3 in Figure 3, time points 2 and 3 (3 is superimposed on 2), and in Figure 4, time points 1 and 2. 4.2. Variables. While STG uses one set of binary variables for processing tasks and another for cleaning tasks, the new MTG formulation requires a single set because it considers

combined processing and changeover tasks. Variables N h i,i′,m,t now represent the execution of order i at time point t, together with the required changeover task for order i′ to immediately follow in unit m. In cases where a single order is allocated to a particular unit, no changeover will be required and cases with i ) i′ must be part of the domain, not like in the STG. Four other sets of binary variables, Y1i,k-Y4i,k, are required to identify the location of task i in stage k (k * |K|) with respect to the processing of the same order in stage k + 1. These will be required to accurately determine the capacity of the intermediate storage vessel, Vi,k. Other sets of variables common to both formulations are ξi,m,t and Ci,m,t. Note that the excess variables linked to material states have been dropped from the MTG because these are only appropriate when using a single time grid.4 The timing variables Tt,m now feature an additional index to account for the fact that |M| time grids are used (see Figure 3) and H remains the cycle time. As a result of the use of combined processing and changeover tasks, the time points will be associated to the start of the processing tasks but normally not to their end. For inventory purposes, the exact ending time of the processing part of the task must be determined, and four other sets of positive continuous timing variables are needed. Combined tasks starting at time point t unit m end their processing at Tpt,m. The start and finishing time of order i in stage k is given by Tsi,k and Tfi,k, respectively. To avoid numerical problems with the nonlinear problems, we will be using ∆Ti,k for the processing time (k ) |K|). Finally, the processing rate of order i in stage k is given by Fi,k, while Di will give the delivery/demand rate of product i. 4.3. Conceptual Overview. Figure 4 illustrates how MTG works with a simple example comprising three orders, three units, and two stages, the last one featuring two units in parallel. Each order is processed just once on every stage resulting from the execution of a particular combined task. In equipment units involving two or more orders (e.g., M1 and M3), the combined tasks will always feature processing and cleaning tasks and will thus have different order indexes (e.g., N h I2,I1,M1,1 ) 1). In contrast, unit M2 handles only I3, meaning that the combined task consists of just the processing part and the two order indexes are the same (N h I3,I3,M2,1 ) 1). In this case, the task length is equal to the time horizon, H, which means that the order processing rate is equal to its delivery rate, Di. Note that, like in STG, the processing rate of a task (now stage rather than unit related), Fi,k, can be lower than its maximum processing rate, to reduce inventory costs. On the other hand, MTG does not need to be concerned with relaxing the changeover times when necessary because, due to the multiple time grids, the number of time points that is required to solve the problem is not affected by them. Thus, the cleaning part of the combined task lasts exactly the pre-specified value, cli,i′,m. As mentioned in section 4.1, all time grids are placed inside the interval [0,H] that comprises the cycle under consideration, meaning that the timing variables Tt,m are subject to those bounds. These time points identify the starting points of tasks so it is natural that variables Tsi,k are subject to the same limits. However, because of the possibility of task wrap-around, it is likely that a few tasks end at t ) 1, with T1,m > 0 (e.g., task (I2,I1) in M3). In such a case, if one makes the corresponding Tfi,k ) T1,m, then we would have Tsi,k ) T|T|,m > T1,m, which is not possible because by definition ∆Ti,k ) Tfi,k - Tsi,k g 0. This can be easily overcome by adding H to such Tsi,k and corresponding Tpt,m variables and by setting their upper bounds to 2H. Note that the round dotted line located outside the cycle

3674


Figure 4. Possible solution from multiple time grid formulation MTG (|I| ) 3, |M| ) 3, |K| ) 2, and |T| ) 3).

Figure 5. Four different possibilities of task interaction between consecutive stages.

boundaries in Figure 4 has exactly the same length as the continuation of task (I2,I1) in M3, in the early periods of the cycle (Tp3,M3 ) TfI2,K2 ) T1,M3 - clI2,I1,M3 + H). The model assumes that there is no net accumulation of materials over the cycle, meaning that (i) the amount of material processed in stage k must equal the amount processed in stage k + 1 for a particular order and (ii) the total amount of product produced over the cycle exits the plant at a constant rate Di. If processing tasks of the same order, in consecutive stages, are not processed simultaneously, it is clear that the capacity of the storage vessel must equal the total amount produced. On the other hand, if they overlap at least partially, then the worst case scenario for intermediate storage will not occur. The exact formula to use to determine Vi,k depends on the tasks’ relative position, and four different possibilities exist (one for each binary variable), as can be seen in Figure 5. For Y3i,k ) 1 or Y4i,k ) 1, one task completely overlaps the other, meaning that the one that lasts longer has a lower processing rate (because they both process the same amount of material). For Y3i,k ) 1, the period of highest inventory occurs at Tsi,k+1 because, after that point, stage k + 1 consumes the intermediate at a higher rate than its rate of production in stage k, and the material in storage decreases steadily until zero at Tfi,k+1. Then, the inventory goes up again until Tfi,k and remains constant until the end of the cycle and, because we are dealing with a periodic problem, up to Tsi,k where it starts to increase again. For Y4i,k ) 1, the period of maximum inventory occurs for Tfi,k and that of minimum inventory for Tsi,k. In both cases, Vi,k is equal to the duration of the highest rate task times the difference between the stage processing rates. When the overlap between tasks of consecutive stages is just partial, for example, Y1i,k ) 1 or Y2i,k ) 1, the formula to use to determine Vi,k depends on the relation between the stage processing rates, which we do not know a priori. Considering

case one as an example, if Fi,k g Fi,k+1, the generation of the intermediate material occurs at a higher rate than its depletion, so the point of maximum inventory occurs at Tfi,k. After that, the inventory decreases until Tfi,k+1, which must necessarily be the point of minimum inventory ()0). Consequently, we would have Vi,k ) (Tfi,k+1 - Tfi,k)Fi,k+1. If, on the other hand, Fi,k < Fi,k+1, the maximum inventory point occurs at Tsi,k+1 and then goes down to zero at Tfi,k+1. It remains zero up to end of the cycle, back to the beginning of the same cycle, and until order i starts to be processed in stage k, at Tsi,k. Thus, Vi,k ) (Tsi,k+1 - Tsi,k)Fi,k. 4.4. Constraints. A few model constraints of the new formulation are like those presented for the STG as a result of the fact that they both rely on the RTN representation and use similar sets of variables. For this reason, we will be focusing on the main adjustments resulting from the alternative continuous-time treatment. The MTG formulation considers a single type of tasks, so the excess resource balances corresponding to equipment states feature a single set of binary variables. The other difference in eqs 14 and 15 with respect to eqs 2 and 4 is that the summation is done for all i′ ∈ I, because as discussed in section 4.3, it is possible that no changeovers are required (in cases where a particular equipment unit is processing just one order). Nevertheless, it is important to emphasize that, whenever we are dealing with problems featuring |M| ) |K|, we can remove variables N h i,i,m,t from consideration.

Ci,m,t ) Ci,m,Ω(t-1) +

(N h i′,i,m,Ω(t-1) ∑ i′∈I

N h i,i′,m,t) ∀ i ∈ I, m ∈ M, t ∈ T (14)

Ci,m,1 ) 1 - ∑∑N h i,i′,m,1 ∀ m ∈ M ∑ i∈I i∈I i′∈I

(15)

With multiple time grids we cannot keep track of intermediates like in single time grid formulations. The best that can be done is to use the overall mass balance, which states that the amount of order i handled in stage k is equal to that processed in stage k + 1; see eq 16. Although we are assuming that no mass is lost from one stage to the next, other mass proportions are straightforward to implement.

∑ ∑(ξi,m,t|m∈M - ξi,m,t|m∈M

m∈M t∈T

k

) ) 0 ∀ i ∈ I, k ∈ K, k * |K|

k+1

(16)

Equation 17 is equivalent to eq 6, while eq 18a is related to eq 8. Now, however, we are using variables Di for the delivery rate of product i so the exact definition is required. The left-


hand side (LHS) consists of the product of two continuous variables, and as a result eq 18a is a nonlinear constraint. Naturally, the delivery rate must be greater than the minimum demand rate (eq 18b). max max ξi,m,t e Fi,m H

N h i,i′,m,t ∀ i ∈ I, m ∈ M, t ∈ T ∑ i′∈I

(17)

∑ ∑ ξi,m,t ∀ i ∈ I m∈M t∈T

(18a)

Di g d i ∀ i ∈ I

(18b)

DiH )

Every order is processed only once on each stage. MTG has as an important model assumption: the execution of a single combined task per stage for each order. Thus, eq 19 is in order.

∑ ∑ ∑ Nh i,i′,m,t ) 1 ∀ i ∈ I, k ∈ K

(19)

m∈Mk t∈T i′∈I

The duration of a particular time interval must be greater than the duration of the combined tasks taking place, which in turn is greater than the minimum processing time for a given amount of material handled (processed at maximum rate), plus the appropriate changeover time. Equation 20 corresponds to eq 5. In STG and as a result of the fact that the time of the first event point is equal to zero (eq 9), the time of the last event point cannot exceed H. In MTG, the time of the first time point can be set to zero just for one equipment unit (see eq 49), so eq 21 is needed.

∑ i′∈I

∑ i∈I

Hmax(1 -

∑ ∑Nh i,i′,m,t′) ∀ i ∈ I, k ∈ K, m ∈ Mk, t ∈ T i′∈I t′∈T

(25)

t′et

Tfi,k g Tpt,m Hmax(1 -

∑ ∑Nh i,i′,m,t′) ∀ i ∈ I, k ∈ K, m ∈ Mk, t ∈ T i′∈I t′∈T

(26)

t′gt

|K|

TΩ(t+1),m + H|t)|T| - Tt,m g

Tfi,k e Tpt,m +

max (ξi,m,t/Fi,m +

N h i,i′,m,tcli,i′,m) ∀ m ∈ M, t ∈ T (20)

T|T|,m e H ∀ m ∈ M

(21)

Equation 27 defines the processing rate of order i in stage k and is a nonlinear constraint. The RHS gives the amount of i produced in all units m belonging to stage k. Note that eqs 1719 ensure that there will be exactly one unit m and time interval t for which ξi,m,t * 0. The total amount produced of order i in stage k is also given by the product of the processing rate by the difference between the ending and the starting times of order i in stage k (LHS). The processing rate on a particular stage must be lower than the order maximum processing rate, eq 28. The difference between the ending and starting event points of a task is equal to its processing time, which needs to be defined explicitly only for the last processing stage; see eq 29. To avoid numerical problems with the solvers, those variables need to be properly bounded. Equation 30 can be replaced by a tighter constraint, eq 31, whenever |M| ) |K|.

Fi,k(Tfi,k - Tsi,k) )

∑ ∑ ξi,m,t ∀ i ∈ I, k ∈ K

(27)

m∈Mk t∈T

max N h i,i′,m,tFi,m ∀ i ∈ I, k ∈ K ∑ ∑ ∑ m∈M i′∈I t∈T

(28)

∆Ti,|K| ) Tfi,|K| - Tsi,|K| ∀ i ∈ I

(29)

∆Ti,|K| e H ∀ i ∈ I

(30)

∆Ti,|K| + ∑ ∑ ∑ ∑ N h i,i′,m,tcli,i′,m ) H ∑ i∈I m∈ M i∈I i′∈I t∈T

(31)

Fi,k e

k

|K|

i′*i

Besides Tt,m, MTG features four other sets of timing variables that are required to calculate the inventory costs. Combined tasks starting at Tt,m will end their processing part at Tpt,m and their changeover part at TΩ(t+1),m; see eq 22 and Figure 4. Those two sets of variables will correspond to starting, Tsi,k, or ending, Tfi,k, times of processing tasks, respectively, whenever ∑i′∈IN h i,i′,m,t ) 1 (m ∈ Mk). Equations 23-26 are the required big-M constraints, which are conceptually similar to those used in Castro et al.4 to relate timing variables Tt,m with transfer time variables. As explained in that paper, for performance reasons we include inside the big-M term other binary variables than those expected (i.e., with t′ ) t).

The next set of constraints activates the proper binary variable (see Figure 5) according to the existing interaction between processing tasks of the same order in consecutive stages. Because a given relation between any two start/ending time variables is valid for at least two alternative interaction cases, the big-M term features two or four binary variables. By doing this, we ensure that at most a single binary variable becomes active for order i, stage k. In eqs 32, 33, 36, and 37, the big-M value is equal to the upper bound on the cycle time, because the starting times (first term on the RHS) are bounded by Hmax; see eq 46. In contrast, the ending times of orders are bounded by twice the maximum cycle time, so we have used this value in eqs 34 and 35.

TΩ(t+1),m + H|t)|T| - Tpt,m )

Tsi,k e Tsi,k+1 +

∑ ∑Nh i,i′,m,tcli,i′,m ∀ m ∈ M, t ∈ T i∈I i′∈I

(22)

Tsi,k g Tsi,k+1 -

Tsi,k e Tt,m + Hmax(1 -

∑ ∑ i′∈I t′∈T

N h i,i′,m,t′) ∀ i ∈ I, k ∈ K, m ∈ Mk, t ∈ T (23)

t′et

∑ ∑Nh i,i′,m,t′) ∀ i ∈ I, k ∈ K, m ∈ Mk, t ∈ T i′∈I t′∈T t′et

Hmax(1 - Y2i,k - Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (33) Tfi,k e Tfi,k+1 + 2Hmax(1 - Y1i,k - Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (34)

Tsi,k g Tt,m Hmax(1 -

Hmax(1 - Y1i,k - Y2i,k) ∀ i ∈ I, k ∈ K, k * |K| (32)

(24)

Tfi,k g Tfi,k+1 - 2Hmax(1 - Y2i,k Y3i,k) ∀ i ∈ I, k ∈ K, k * |K| (35)

3676

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 max Fi,k.l ) max Fi,m ∀ i ∈ I, k ∈ K

Tfi,k g Tsi,k+1 - Hmax(1 - Y1i,k - Y2i,k - Y3i,k Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (36) Tfi,k+1 g Tsi,k - Hmax(1 - Y1i,k - Y2i,k - Y3i,k Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (37) As can also be seen in Figure 5, the formula to calculate the capacities of the intermediate storage vessels depends on the type of task interaction and, for some cases, even on the relation between the processing rates of consecutive stages. Seven big-M constraints are needed, eqs 38-44, where eq 44 gives the worst case scenario corresponding to no interaction, where the storage capacity is equal to the total amount of material processed over the cycle. The big-M value was made order dependent, equal to the maximum amount of material that can be processed over the cycle given by the product of maximum processing rate of the limiting stage by the maximum cycle time; see eq 45.

Vi,k g (Tsi,k+1 - Tsi,k)Fi,k Mi(1 - Y1i,k) ∀ i ∈ I, k ∈ K, k * |K| (38) Vi,k g (Tfi,k+1 - Tfi,k)Fi,k+1 Mi(1 - Y1i,k) ∀ i ∈ I, k ∈ K, k * |K| (39) Vi,k g (Tsi,k - Tsi,k+1)Fi,k+1 Mi(1 - Y2i,k) ∀ i ∈ I, k ∈ K, k * |K| (40) Vi,k g (Tfi,k - Tfi,k+1)Fi,k Mi(1 - Y2i,k) ∀ i ∈ I, k ∈ K, k * |K| (41) Vi,k g (Tfi,k+1 - Tsi,k+1)(Fi,k+1 - Fi,k) Mi(1 - Y3i,k) ∀ i ∈ I, k ∈ K, k * |K| (42) Vi,k g (Tfi,k - Tsi,k)(Fi,k - Fi,k+1) Mi(1 - Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (43) Vi,k g

∑ ∑ξi,m,t - Mi(Y1i,k + Y2i,k + Y3i,k +

m∈Mk t∈T

Y4i,k) ∀ i ∈ I, k ∈ K, k * |K| (44) max ∀i∈I Mi ) (min max Fmax i,m )H k∈K m∈Mk

(45)

Equation 46 places upper and lower bounds on the timing variables. The explanation for some of them has already been given in section 4.3. The only exception is the lower bound on ∆Ti,k, which is determined as the minimum amount that can be processed over the cycle, divided by the maximum processing rate in stage k. Equation 47 initializes (attribute .l) the processing rates to their maximum values to help with convergence.

Hmin e H e Hmax Tt,m e Hmax; Tpt,m e 2Hmax ∀ t ∈ T, m ∈ M Tsi,k e Hmax; Tfi,k e 2Hmax; max e ∆Ti,|K| e Hmax ∀ i ∈ I (46) diHmin/ max Fi,m m∈M|K|

(47)

m∈Mk

For the anchoring of the schedule constraint, we assign the processing of order one in the first stage to the first time point (eq 48). Additionally, it is desirable to set the time value of the first time point to zero, but this can only be done for a single equipment unit to avoid reducing the feasible space. Any unit can be selected, so we have selected the first one, as can be seen in eq 49. The last set of optional but computational efficient constraints is eq 12, presented for STG, and that can be used by MTG wherever there is a single equipment unit per stage.

∑ ∑ N1,i′,m,1 ) 1

(48)

T1,1 ) 0

(49)

m∈Mk)1 i′∈I

Finally, the new MINLP, multiple time grid continuous-time formulation is complete with the profit maximization objective function, eq 50. Notice that when compared to eq 13, both product revenue and final inventory cost terms have changed as a result of the use of different sets of variables. The latter is calculated by the triangle area as can be seen in Figure 6 (where the start of the task coincides with time zero for simplicity).

[

∑∑∑

max

i∈ I

m∈M|K|

ciξi,m,t -

t∈T

k∈K m∈Mk

∑∑ k∈K k*|K|

∑∑∑∑∑

i∈I

i∈I

i′∈I

t∈T

H

cii,kVi,k

∑

ctk,i,i′cli,i′,mN h i,i′,m,t -

]

cfi/2[(Fi,|K| - Di)(∆Ti,|K|)2 + Di(H - ∆Ti,|K|)2]

i∈I

H

(50)

5. Computational Results The performance of the continuous-time formulations is going to be illustrated through the solution of a few example problems, where most of the data was taken from Pinto and Grossmann.10 According to the authors, such examples arise, for instance, in continuous multiproduct plants for manufacturing lubricants. Some changes in the data were required because a more general problem is being handled. More specifically, we consider a stage specific transition cost that is proportional to the total changeover time in that stage. Additional equipment units were added to a two stage, three products example to illustrate that the formulations are perfectly capable of handling cases with more than a single unit per stage. Three examples featuring parallel units together with the base case will be thoroughly analyzed in section 5.1, while larger examples with a single unit per stage are left for section 5.2, where the focus will be more in computational performance and less on practical insights. Data for the former set of problems are given in the text, whereas for the latter set it will be provided as Supporting Information instead. Their main characteristics are given in Table 5. The mathematical formulations were implemented and solved in GAMS build 22.2. DICOPT was used as the MINLP solver, which in turn relied on CONOPT as the NLP solver and CPLEX as the MILP solver. BARON was used as an alternative MINLP solver in some problem instances. The computer consisted of a Pentium-4 3.4 GHz processor with 2 GB of RAM, running Windows XP Professional, and the computational effort required to solve the 10 example problems by the two formulations can be found in Tables 6 and 7. It is important to note at this point that both formulations rely on a single starting point for the solution of the MINLP. For STG, the starting point features all variables at zero except


Figure 6. Final inventory profiles. Table 1. Product Related Data for EX1 product i A B C

price, ci ($/kg)

minimum demand, di (kg/h)

final inventory cost, cfi ($/(kg‚h))

0.15 0.4 0.65

50 100 250

0.00406 0.00406 0.00406

Table 2. Product and Unit/Stage Related Data for EX1 maximum processing max (kg/h) rate, Fi,m product i

intermediate inventory cost, cii,k ($/kg)

unit m ) 1

unit m ) 2

stage k ) 1

800 1200 1000

900 600 1100

0.1406 0.1406 0.1406

A B C

Table 3. Sequence Dependent Changeover Times, cli,i’,m (h), and Costs ctk,i,i′ ($/h) for EX1 unit m ) 1 product i A B C a

Aa 10 4

Product

Ba

Ca

5

8 3

7

unit m ) 2 Aa 2 6

stage k ) 1

Ca

Aa

11

1 5

760 760 760 760 750 750 750 750 770 770 770 770

1

Ba

Ca

stage k ) 2

Ba

Aa

Ba

Ca

i′.

the cycle time H, which is at its lower bound. For MTG, the initialization consists of having the delivery rate Di, the processing times ∆Ti,k, and H at their lower bounds, the processing rates Fi,k at their upper bounds (see eqs 18b, 46, and 47), and all other variables at zero. From these points, DICOPT could always find very good solutions to the problems so the whole solution process is robust. Another relevant issue is that STG is a more general formulation than MTG and so can theoretically find better solutions to the problem, as will be discussed later on. Thus, different solutions to the same problem should be expected. More importantly, global optimal solutions from MTG can be worse than local optimal solutions from STG. 5.1. Motivating Examples with Parallel Units within Stages. The first example, EX1, is essentially example 2 of Pinto and Grossmann10 and has been recently solved in Castro et al.16 (example 3). It consists of a two stage, three products problem, and the data are given in Tables 1-3. Solving the problem with the single time grid formulation (STG) requires, as discussed in section 3 and in the literature1,3,5,8,13,15,16 (in Castro et al.16 for this particular example), an iterative search procedure for the optimal number of event points (|T|), the one that leads to the best solution. Although this procedure has been performed for all examples solved, only the best value of |T| will be reported to prevent losing focus on other more important issues. Recall that the proper selection of |T| is not an issue for the multiple time grid formulation (MTG), as discussed in section 4.1.

The best solution for EX1 was found by STG for 11 event points. It features an optimal profit of $355.09/h for a cycle time of 184.646 h, which is better than $352.93/h, the best value reported in Castro et al.16 Interestingly, in that work, the problem was solved by the exact same formulation but with different versions of the solvers, clearly highlighting the non-convex property of the feasible region associated to the NLP and the existence of several local optima. The solution shown in Figure 7, if not the global optimum, should be close to it because a number of sensible heuristic rules can be observed: (i) the product sequences are A-B-C and B-A-C for units M1 and M2, corresponding to a minimum total changeover of 16 h; (ii) the delivery rates of the least valuable products (A and B) are equal to the minimum demand rates; (iii) the processing rates are equal to their maximum values for all but three tasks (the first instance of the production of B in unit M2 is performed at 98.55% of the maximum rate, and the first two instances of C in unit M2 are at 71.74 and 90.91%, respectively); (iv) all but one cleaning task, which has its duration extended by 1 h, last exactly what they should; and (v) the production of intermediate C in stage 1 is accompanied most of the time by its consumption in stage 2, to get a low-capacity storage vessel, a critical issue for the product with highest demand. Concerning point four it can be argued that the cost of the extra cleaning hour ($750) should be considered (see eq 13), which would cause the profit to be reduced by $4.09/h to $351/h, making it a worse solution than the best reported. This can be avoided by turning changeover tasks into ZW tasks,8,16 but this would generally make good solutions harder to obtain because of the higher computational effort caused by the requirement of more event points. The best solution found by MTG has a slightly lower profit, $352.45/h (0.7% less), and a lower cycle time, 173.519 h. It is solved for |T| ) 3 and, in each stage, features the execution of a single processing task per product (see Figure 7). It is clear that the Gantt charts and profiles are very similar to those obtained from STG, with the particularity that maximum inventory points are usually slightly lower for MTG due to the lower cycle time, translating into lower intermediate and final inventory costs. The tradeoff, as first noted by Pinto and Grossmann,10 is that a lower cycle time also means higher changeover costs because more cycles and hence more changehfovers are required on a yearly basis. Overall, MTG still has the lowest cost solution, but STG has better product revenue as a result of a higher delivery rate of C, 789.18 vs 785.01 kg/h. Suppose now that we have an additional equipment unit available and need to decide where to allocate it to increase the profit as much as possible. In Figure 7, unit M1 is working at its maximum processing rate all the time, contrary to M2, suggesting that stage one is the limiting one. Thus, the first variant of the problem, EX1a, allocates the extra unit to stage one. Its maximum processing rate is about half that of M1 (its competitor; see Table 4) and is named M2. Unit 2 of EX1 is renamed M3. The best solution for EX1a is obtained by the multiple time grid formulation (MTG) with BARON, corresponding to a profit of $391.61/h and a cycle time of 146.18 h. Hence, the decision to allocate the extra unit to stage one increased the profit of the plant by 10.3%. Figure 8 shows the optimal schedule obtained, featuring three event points, where it can be seen that the results are not as good as expected because unit M2 is idle most of the cycle. It is interesting to see that the least valuable product, which is also the one with lowest maximum processing rate in stage one (A), is processed in the least powerful equipment, M2. As a consequence of unit M1 not handling product A, more

3678


Figure 7. Best solutions from STG (left) and MTG (right) for EX1. Table 4. Problem Data for Additional Unit(s), Examples EX1a,b,c cli,i′m (h) product i A B C a

Product

max Fi,m (kg/h)

Aa

400 500 450

5 6

Ba

Ca

9

11 7

4

i′.

C can be produced in stage one, the limiting stage for C (see Table 2), and the delivery rate can be increased to 825.46 kg/h. The optimal product sequence in M1 becomes B-C, and because of the production of a single material in M2 no changeovers are required in this unit. The product sequence in M3 remains unchanged leading to a total changeover of 14 h, 2 h less than for EX1. The contribution of the total changeover cost for the objective function is thus lower for EX1a so it is no surprise that the optimal cycle time has decreased to also decrease inventory levels (see Figure 8) and the total inventory cost. The results for STG are very similar, and the values of the cycle time are even closer than for EX1. The profit difference is almost negligible (0.16%) with STG featuring a slightly higher revenue due to the production of 825.85 kg/h of C, lower changeover and higher inventory costs caused essentially by the larger value of the cycle time because the storage profiles are almost identical (the differences are mostly caused by the different relative position of the tasks with respect to the cycle origin). The second variant of EX1, EX1b, allocates the extra unit, now named M3, to stage two. STG leads to the best solution, and the optimal profit for EX1b is equal to $468.16/h, 31.8% higher than for EX1. It is thus a better option to allocate the extra equipment unit to stage two instead of stage one. From Figure 9 one can see that the cycle time has increased substantially for STG and even more for MTG, when compared to previous cases. We have already discussed that the cycle time is influenced by the tradeoff between the total changeover and the inventory cost, the latter being dominated by the final inventory cost. With respect to changeovers, the solution from

STG features an A-B-C sequence in stage one and an A-C sequence in unit M2, corresponding to a total changeover of 19 h. Interestingly, in the optimal solution, the changeover task from C to A has been replaced by two consecutive changeovers C-B and B-A with a total changeover of just 3 h instead of 6. This behavior highlights one of the limitations of the single time grid formulation where, because changeover tasks are not combined with processing tasks (like in the MTG), it is possible that a particular changeover is replaced by a sequence of tasks (with the same initial and final states) with lower total duration. To overcome this, one could add a constraint that would limit the number of cleaning tasks that could be executed, but because we are dealing with problems with units in parallel we do not know a priori how many changeovers will be required. Contrary to STG, MTG has opted for producing a single product in M2 and for an A-B sequence in M3. The accounted total changeovers in stage two are equal to 14 h instead of 3 h, meaning that the changeover costs are higher for MTG, and thus explaining the need to use a higher cycle time. It is interesting to note that adding the equipment unit to stage one or two has more or less the same effect on total production because the delivery rate for C is for EX1b equal to 824.44 kg/h and 830.23 kg/h, for STG and MTG, respectively. The advantage of adding the unit to stage two is because a reduction in final inventory cost is more important than a reduction in intermediate inventory cost. With the additional unit in stage two, there is one equipment handling a single product for which the processing rate can be made equal to the delivery/demand rate, thus leading to a zero final inventory cost for that product (see the optimal profiles in Figure 9, for product C in MTG, and, although only partly, for product B in STG). Finally, the best schedule features the allocation of product B, the one with the lowest processing rate in unit M2 (the least powerful unit), to M3, its competitor in stage two. The third and final variant of EX1, EX1c, adds one unit to each stage, and these are named M2 and M4. When compared to the base case, the profit for EX1c is almost doubled to


Figure 8. Best solutions from STG (left) and MTG (right) for EX1a.

Figure 9. Best solutions from STG (left) and MTG (right) for EX1b.

$708.44/h for STG. MTG can only reach a profit of $668.25/h (5.7% less), when using DICOPT (BARON failed to find a solution up to 6000 s). This is the most significant difference so far between the two continuous-time formulations and serves to illustrate one of the limitations of the multiple time grid approach. As can be seen in Figure 10, the best solution features C in both units of stage one. With this, it is possible to achieve a delivery rate of 1100 kg/h, the maximum processing rate of C in M3 (see Table 2). M3 is dedicated exclusively to C, and

because this product is not produced in M4 no final storage is required, as can be seen in the profile chart. On the other hand, MTG has as a fundamental model assumption that each order can only be processed once on each stage, which is enforced by eq 19. Hence, the delivery rate of C is limited by the maximum processing rate of stage one (unit M1), 1000 kg/h. This lower revenue from the most valuable product is the main reason for the profit disparity between STG and MTG. As a result of the implicit upper bound on the delivery rate of C, the delivery rate of B exceeds, for the first time, the minimum

3680


Figure 10. Best solutions from STG (left) and MTG (right) for EX1c. Table 5. Results Overview (DICOPT Solver) STG

MTG

problem

products

stages

units

Hmin (h)

Hmax (h)

|T|

profit ($/h)

H (h)

|T|

profit ($/h)

H (h)

EX1 EX1a EX1b EX1c EX2 EX3 EX4 EX5 EX6 EX7

3 3 3 3 4 5 7 3 4 5

2 2 2 2 2 2 2 3 3 3

2 3 3 4 2 2 2 3 3 3

100 100 200 200 400 400 400 200 200 400

250 250 600 600 800 800 800 600 600 800

11 8 9 7 9 10

355.087 390.966 468.156 708.44 7074.61 6397.78

184.646 148.117 403.692 499.103 682.283 561.828

8 8

5789.04 5166.12

407.15 503.451

3 3 3 3 4 5 7 3 4 6

352.446 382.647 450.548 668.25 7099.19 6514.49 5572.57 5923.90 5312.76 4986.59

173.519 145.674 501.315 312.725 674.304 569.69 560.763 419.186 477.037 499.383

demand, and is naturally greater for MTG, which is more constrained on the production of C, than for STG. The values are 415.12 kg/h and 306.35 kg/h, respectively. The total duration of changeover tasks in the optimal schedules is equal to 33 h for STG and 21 h for MTG. Because the cycle times are respectively 499.103 and 312.725 h, the trend that lengthier changeovers lead to higher cycle times because of the need to the reduce idle times over a fixed time span (e.g., 1 year) is observed once more. To conclude, let us indicate the values of the binary variables associated to the interaction of processing tasks of the same order in consecutive stages: Y1C,1 ) 1, Y2B,1 ) 1, and Y3A,1 ) 1. Notice, however, that the exact same solution could be obtained by MTG with Y3C,1 ) 1, Y4B,1 ) 1, and Y2A,1 ) 1 (see Figures 5 and 10). 5.2. Larger Example Problems. A few more problems were solved to provide additional results for the comparison between the two continuous-time formulations. Examples EX2-7 feature a single unit per stage and more products than the previous ones. They gradually increase in size, as a result of an increase either in the number of products or in the number of stages. Additional increments resulted in problems from MTG that were intractable by DICOPT, while EX4 and EX7 from STG could also not be solved. The limits for the BARON solver are even lower. It is useless on problems resulting from STG because the smallest one (EX1) cannot even find a feasible solution in reasonable

Table 6. Results with Solver BARON

problem

model

|T|

profit ($/h)

best possible solution ($/h)

EX1 EX1 EX1a EX1b EX1c EX2 EX5

STG MTG MTG MTG MTG MTG MTG

6 3 3 3 3 4 3

no solution 352.446 391.613 325.346 no solution 7095.12 5924.20

472.85 352.446 411.24 782.85 1271.54 7102.17 5949.39

a

H (h) 173.520 146.18 326.214 673.302 418.80

CPU time (s) 4000a 732 5000a 6000a 6000a 6000a 53700a

Maximum resource limit.

time for the minimum number of event points (6) that ensures feasibility; see Table 6 (the MINLP uses 108 binary variables, 230 continuous variables, and 147 constraints). The MINLP from EX1 and MTG features 46 binary variables, 125 continuous variables, and 196 constraints (see Table 7) and can be solved to global optimality in 732 s. Nevertheless, the solution is the same as that obtained by DICOPT in just 1.5 s. The second smallest problem is generated from EX5 with 78 binaries, 191 continuous variables, and 309 constraints, which cannot be solved to global optimality in almost 15 h. The best solution is worth $5924.2/h for a possible maximum of $5949.4/h (relative optimality gap ) 0.4%), while DICOPT can reach almost the same solution ($5923.9/h) in only 4.9 s. For EX2, BARON

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3681 Table 7. Overview of Computational Statistics (DICOPT Solver) STG problem EX1 EX1a EX1b EX1c EX2 EX3 EX4 EX5 EX6 EX7

MTG

binary variablesa

continuous variablesa

constraintsa

CPU time (s)a

total CPU (s)b

CPU time (s)

binary variables

continuous variables

constraints

198 216 243 252 288 500

415 424 476 477 523 822

262 250 280 270 280 381

33.6 53.1 32.5 74.8 239 7899

125 (12) 158 (9) 186 (10) 964 (8) 2151 (10) 7899 (10)

216 384

451 694

299 388

111 1863

1.5 1.3 1.6 2.9 3.2 287 71598 4.9 58.3 33284

48 93 93 120 112 220 616 78 176 340

125 194 194 245 230 387 905 191 350 587

196 257 257 319 307 442 784 309 481 689

430 (9) 32532 (9)

a Corresponding to the value of |T| listed in Table 5. b For all iterations starting with the minimum number of event points that achieves feasibility up to the value indicated between brackets.

found an optimal solution of $7095.12/h for a possible maximum of $7102.17/h (in 6000 s), which means that the optimal solution from DICOPT ($7099.19/h) is within 0.04% of the global optimal solution. Therefore, the practical use of BARON lies somewhere in between the size of problem EX1 (three products, two stages, and two units) and the size of problem EX2/EX5 (four products, two stages, and two units/three products, three stages, and three units). Table 5 gives an overview of the results obtained in terms of the values of the two most important model variables, the profit and the cycle time, H. Besides the main characteristics of the problem, we show the upper and lower bounds that were used for the cycle time and the number of event points leading to the best solution. It is important to emphasize that the specification of the cycle time range does not play a part as important as that reported in the work of Wu and Ierapetritou,9 where the optimal cycle time showed a significant dependence on the chosen range. In this work, we observed that given inadequate upper or lower bounds, the cycle time converged to one of the limits so we simply decreased Hmin or increased Hmax and resolved the problem. This distinct behavior is probably due to the fact that continuous processing tasks are being considered instead of batch tasks. In addition to the above procedure, a sufficiently wide time span was considered where [Hmin, Hmax] ⊃ [Hmin, 2Hmin] to avoid cutting solutions with a cycle time lower than Hmin from the feasible space. Note that such schedules can be generated by repeating the production pattern as many times as those required (e.g., two times) to get a cycle time within the given bounds (STG). For the same problem, a longer cycle time is associated to the typical requirement of a larger number of event points (schedule of higher complexity), and because of that, the formulation will be more demanding computationally.9 The most important result in Table 5 is that the single time grid formulation leads to the best solutions (in bold) only for EX1, EX1b, and EX1c, three of those discussed in the previous section. This behavior is entirely due to two performance factors because in the limit of an infinite number of event points, the feasible region from the STG includes that of the MTG. Evidence of this can be found in EX1 for which the global optimal solution from MTG (see Table 6) is worse than the local optimal solution from STG. The first performance factor is that the STG generates harder MINLPs (see Table 7), due to the use of more event points. Finding the optimal solution involves iterating on |T|, and because a single increment may involve an increase in computational effort up to 1 order of magnitude, the mathematical problems rapidly become intractable, probably before reaching the value of |T| that would allow generating the same solution as that from the MTG. The second

is that the problem is non-convex, and local optimization solvers may fail to find the global optimal solution for a given |T| value. This is particularly important because it also affects the search procedure over |T|. More specifically, we have found that a single increase in |T| sometimes led to a worse solution despite the fact that the feasible region had increased. EX1 and EX3 match examples 2 and 3 of Pinto and Grossmann,10 two of the three examples for which the optimal solution was shown. A larger example featuring eight products in three stages was also discussed,10 but this was found to be intractable by both our continuous-time formulations. EX1 is the best example to highlight the higher generality of STG and MTG over the multiple time grid formulation of Pinto and Grossmann.10 Besides their ability to handle parallel units within stages, STG and MTG have two additional advantages. The first is that they are not limited to the same product sequence in all stages. The second is that the processing rates may be lower than the maximum values. Because of this, the solution reported by Pinto and Grossmann10 corresponds to a profit of just $297/ h, well below the values of $355/h and $352.4/h found respectively by STG and MTG. The differences for EX3 are not as significant, mainly because the lowest possible duration of changeover tasks is achieved with the same sequence (AB-C-E-D) in both stages. Nevertheless, the use of lower processing rates for some tasks (production of E and D in stage 1) enabled MTG to find a slightly better solution ($6514.49/h, H ) 569.69 h) than the $6513/h (H ) 569.8 h) reported by Pinto and Grossmann.10 This is shown in Figure 11, where in contrast to the previous charts, the inventory profiles are now given for the intermediate materials only. Figure 12 shows the optimal schedule obtained by MTG for EX7, the hardest three-stage problem. Again the product sequence A-D-E-B-C is the same in all stages. The total duration of changeover tasks is equal to 90 h. Stage 2 (unit M2) is the limiting stage because all tasks are performed at the maximum rate. Notice that like in EX3, all orders are more or less being processed in the three stages simultaneously, to reduce inventory profiles. The delivery for product C is equal to 718 kg/h while all other products are at their minimum demand rates. 6. STG vs MTG The benefits from the use of a more general formulation often come at a cost and this work is no exception. The results in Table 7 clearly show that an important decrease in generality (recall that the single time grid formulation can also handle problems with other features, like more complex network structures), is accompanied by a substantial decrease in computational effort, when going from STG to MTG. Both the

3682

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 Table 8. Strengths and Weaknesses of STG and MTG Formulations formulation single time grid (STG) strengths

weaknesses

Figure 11. Best solution from MTG for EX3.

Figure 12. Best solution from MTG for EX7.

partial (for the value of |T| that leads to the best solution) and cumulative computational efforts are given for DICOPT, following the suggestion of Shaik et al.,3 who recently noted that the latter metric is a fairer performance indicator when comparing single to multiple time grid formulations. It can be seen that the computational time for the best iteration of STG is at least 1 order of magnitude larger than for MTG, while the cumulative effort is about 2 orders of magnitude greater. STG requires a larger number of binary and continuous variables, mainly because set T has more elements, while the number of constraints is more or less the same because MTG requires a large number of constraints just to determine the capacity of the intermediate storage vessels. The number of binary variables in MTG increases rapidly with the number of products, which is not surprising because double order index variables are used. When compared to the formulation of Pinto and Grossmann,10 which requires 36 binary variables, 191 continuous variables, and 293 constraints for EX1 and 100, 527, and 727 for EX3, MTG uses more binaries but fewer continuous variables and constraints. On the basis of the computational performances listed in Pinto and Grossmann,10 their formulation is significantly faster, particularly if one takes into account that the current hardware is significantly more powerful than the computer used then.

suitable for more complex network structures, i.e., multipurpose plants can lead in principle to better solutions than MTG for |T| f ∞ any order can be processed in two or more units that belong to the same stage gives rise to mathematical problems that quickly become intractable as the problem size increases search for the optimal solution requires iterative procedure over |T| correct changeover task may be replaced by a sequence of tasks with the same initial and final equipment state and with lower total changeover

new multiple time grid (MTG) very good computational performance in multistage problems straightforward specification of |T| correct execution of changeover tasks due to the use of combined processing and changeover tasks limited to multistage problems

a particular order can be processed only once on each stage, i.e., parallel units cannot be used to their full extent smaller feasible space may lead to suboptimal solutions when compared to STG

We find it convenient to summarize the strengths and weaknesses of the single and multiple time grid formulations that were mentioned throughout the paper. These are listed in Table 8 for easier consultation. Keep in mind that STG can handle more flexible units where a particular order is processed in two or more units belonging to the same stage. As a consequence, the number of tasks to execute is not easy to predict and an iterative procedure on the number of time points making up the time grid is necessary. MTG, on the other hand, cannot handle production of orders in parallel units because it assumes the execution of just one order per stage. The advantage is that the required number of tasks can be easily predicted, and no iterative procedure is needed for finding the total number of event points to use. 7. Conclusions This paper has presented a new continuous-time formulation for the optimal periodic scheduling of multistage continuous plants with sequence dependent changeovers. It uses multiple time grids, one per equipment resource, where the number of event points that achieves optimality can be easily determined; that is, it has no effect in the quality of the solution. Although based on a conceptually similar short-term scheduling formulation for batch plants,4 it now determines not only the allocation of orders to equipment units but also the amount that needs to be processed. Furthermore, it is able to calculate the capacity that is required for the intermediate storage vessels. This has been the most substantial development due to the fact that the inputs and outputs of such storage vessels result from equipment units belonging to different time grids that, because of that, are difficult to relate. The formulation gives rise to MINLPs, where


the cycle time is an important model variable and can handle problems involving units in parallel. The performance of the new formulation has been illustrated through the solution of 10 example problems taken from the literature.10 These were also solved by a single time grid formulation8 for comparative purposes. The results have shown that the new formulation is computationally much more efficient, although a larger feasible space is associated with the single time grid formulation in the limit of an infinite number of event points. Because of this, the latter can and did find better solutions for the simpler example problems. Practical application limits for the continuous-time formulations for problems involving a single unit per stage were found to be three products in two stages and four products in two stages for the single time grid formulation (STG) and seven products in two stages and five products in three stages for the multiple time grid (MTG). Both formulations were also compared to the one of Pinto and Grossmann.10 They are significantly more general, more so the single time grid formulation, because they are not limited to a single unit per stage, to the same product sequence on every stage, or even to tasks always executed at their maximum processing rates. However, they also give rise to mathematical problems that are considerably more difficult to solve, so it is fair to conclude that the new multiple time grid formulation should be preferred over that of Pinto and Grossmann10 only in cases where such features are relevant. Supporting Information Available: Tables of data for examples EX2-EX7. This material is available free of charge via the Internet at http://pubs.acs.org. Literature Cited (1) Meńdez, C. A.; Cerda´, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art Review of Optimization Methods for Short-Term Scheduling of Batch Processes. Comput. Chem. Eng. 2006, 30, 913. (2) Floudas, C. A.; Lin, X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: A review. Comput. Chem. Eng. 2004, 28, 2109. (3) Shaik, M.; Janak, S.; Floudas, C. Continuous-Time Models for ShortTerm Scheduling of Multipurpose Batch Plants: A Comparative Study. Ind. Eng. Chem. Res. 2006, 45, 6190.

(4) Castro, P. M.; Grossmann, I. E.; Novais, A. Q. Two New ContinuousTime Models for the Scheduling of Multistage Batch Plants with Sequence Dependent Changeovers. Ind. Eng. Chem. Res. 2006, 45, 6210. (5) 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. (6) Janak, S. L.; Lin, X.; Floudas, C. A. Enhanced Continuous-Time Unit-Specific Event-Based Formulation for Short-Term Scheduling of Multipurpose Batch Processes: Resource Constraints and Mixed Storage Policies. Ind. Eng. Chem. Res. 2004, 43, 2516. (7) Castro, P. M.; Barbosa-Po´voa, A. P.; Matos, H. A.; Novais, A. Q. Simple Continuous-time Formulation for Short-Term Scheduling of Batch and Continuous Processes. Ind. Eng. Chem. Res. 2004, 43, 105. (8) Castro, P. M.; Barbosa-Po´voa, A. P.; Novais, A. Q. Simultaneous Design and Scheduling of Multipurpose Plants Using Resource Task Network Based Continuous-Time Formulations. Ind. Eng. Chem. Res. 2005, 44, 343. (9) Wu, D.; Ierapetritou, M. Cyclic Short-term Scheduling of Multiproduct Batch Plants using Continuous-time Representation. Comput. Chem. Eng. 2004, 28, 2271. (10) Pinto, J. M.; Grossmann, I. E. Optimal Cyclic Scheduling of Multistage Continuous Multiproduct Plants. Comput. Chem. Eng. 1994, 18, 797. (11) Alle, A.; Pinto, J. M. A General Framework for Simultaneous Cyclic Scheduling and Operation Optimization of Multiproduct Continuous Plants. Braz. J. Chem. Eng. 2002, 19, 457. (12) Pantelides, C. C. Unified Frameworks for the Optimal Process Planning and Scheduling. In Proceedings of the Second Conference on Foundations of Computer Aided Operations; Cache Publications: New York, 1994; p 253. (13) Castro, P. M.; Barbosa-Po´voa, A. P.; Matos, H. A. Optimal Periodic Scheduling of Batch Plants Using RTN-Based Discrete and ContinuousTime Formulations: A Case Study Approach. Ind. Eng. Chem. Res. 2003, 42, 3346. (14) Shah, N.; Pantelides, C. C.; Sargent, R. Optimal Periodic Scheduling of Multipurpose Batch Plants. Annals Oper. Res. 1993, 42, 193. (15) 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. (16) Castro, P. M.; Barbosa-Po´voa, A. P.; Novais, A. Q. A Divide and Conquer Strategy for the Scheduling of Process Plants Subject to Changeovers Using Continuous-Time Formulations. Ind. Eng. Chem. Res. 2004, 43, 7939.

ReceiVed for reView October 23, 2006 ReVised manuscript receiVed March 13, 2007 Accepted March 16, 2007 IE0613570

Optimal Periodic Scheduling of Multistage ... - ACS Publications

Recommend Documents