Planning and Scheduling with Multiple Intermediate Due DatesAn

8314

Ind. Eng. Chem. Res. 2005, 44, 8314-8322

Planning and Scheduling with Multiple Intermediate Due DatessAn Effective Discrete Time Model Formulation Christopher Suerie* FG Produktion & Supply Chain Management, Institut fu¨ r Betriebswirtschaftslehre, FB Rechts- und Wirtschaftswissenschaften, Technische Universita¨ t Darmstadt, Hochschulstrasse 1, 64289 Darmstadt, Germany

The planning and scheduling of chemical plants still leaves several opportunities for optimization. Therefore, many authors have made contributions to this field, resulting in generally two modeling paradigms: continuous time and discrete time model formulations. As both modeling approaches have their advantages, the model formulation presented here is based on a discrete time scale while, at the same time, taking aspects of time continuity into account. It is shown that the proposed modeling approach compares favorably with three model formulations from the literature. Furthermore, a closer evaluation shows that all these benchmarks suffer from minor flaws. 1. Introduction Although much research has been done in the area of planning and scheduling of chemical plants, this topic remains important, as the best possible utilization helps to recoup the large sums invested there. Consequently, a large stream of research both in economics literature (e.g., Blo¨mer and Gu¨nther1 or Kallrath2) and in chemical engineering literature (e.g., Shah3) has evolved. Both are unified by the aim to either schedule as much product as possible in a given time (time-oriented goals) or to maximize or minimize some profit or cost function, respectively (cost-oriented goals). Even though decomposition approaches and heuristics have been proposed (e.g., Lamba and Karimi4) to solve larger problem instances from industry, mathematical programming model formulations form the basis of almost all solution approaches. Mathematical programming-based approaches are attractive, because they allow for the addition of additional restrictions into the model easily as they come up in the underlying production system. Furthermore, models can be set up in a modular fashion such that several modules (sets of constraints) model different aspects of the production system and may be combined freely. Further advantages of mathematical programming-based approaches are that an optimal solution is guaranteed to be found (if enough time can be spent), a feasible solution will be found (if one exists), and bounds on the objective are available throughout the solution process, which altogether give a kind of performance guarantee. Moreover, standard software (CPLEX, LINDO, and Xpress-MP) can be used to solve mathematical programming models, and these software packages have seen remarkable improvements concerning their capabilities in recent years.5,6 However, it is still important to develop and use “good” model formulations, because bad formulations may require much bigger computational efforts to solve a specific problem instance. This paper is organized as follows. First, a description of the underlying planning problem is given. Then, a * Author to whom correspondence should be addressed. Tel.: +49 6151 163563. Fax: +49 6151 165162. E-mail: [email protected], [email protected].

new model formulation based on a discrete time scale is proposed in Section 3. In the fourth section, a test instance frequently used in the literature is introduced. The optimal solution of the proposed model formulation for this test instance is compared to solutions of the same test instance taken from literature. Discrepancies are analyzed, and computational performance is assessed. Finally, some additional computational results are presented, and Section 5 provides a summary and some conclusions. 2. Problem Description The problem tackled here consists of the planning and scheduling of a chemical production site. On this site, several products can be produced on different resources. Not all products may be produced on each resource, but there does not exist a unique product-resource assignment either. Production capacity is scarce, and in addition to production, maintenance activities have to be scheduled in certain predefined time intervals. The production rate is different for each product on each resource and is allowed to vary between certain (resourcedependent) bounds. Changeover activities are sequencedependent. This means the associated costs depend on the order of the products. Changeovers will only be considered in the (cost) objective function, because resources require a minimum run length for production, which is much longer than the changeover time. Demand occurs at various points in time, not only at the end of the planning interval. However, at these intermediate due dates, demand may be either fulfilled or backlogged. Thus, orders that cannot be met immediately need to be fulfilled later (or have to be backlogged until the end of the planning interval). Backlogging is not desirable. Therefore, a shortage penalty is associated with backlogged demand. Furthermore, as the planning data usually contains uncertainty with respect to demands, a certain safety stock should be held for each product. Consequently, there is also a shortage penalty for missing this safety stock target. The aim is to derive a production schedule for the near future (weeks/months) minimizing a cost function consisting of shortage penalties (backlog and safety stock

10.1021/ie048771p CCC: $30.25 © 2005 American Chemical Society Published on Web 09/30/2005

Ind. Eng. Chem. Res., Vol. 44, No. 22, 2005 8315

Figure 1. Different time scales.

violations), changeover costs, and holding costs while, at the same time, respecting the restrictions imposed by the production system (production rates, minimum run-lengths, and capacities). 3. Modeling Approach The modeling of time, either continuous or discrete (see Figure 1), constitutes an important cornerstone of the modeling approach. An overview of production planning models based on continuous and discrete time scales is given, for example, by Suerie.7 Moreover, Sivanandam et al.8 have conducted a comparison of various continuous time approaches. Continuous time scales incorporate the advantage that events (e.g., changeovers) may be scheduled at any point in time. Continuous time scales can be used if time-oriented goals prevail, for example, production is to be maximized in the planning interval. However, if customer demand at certain due dates needs to be met, some sort of time discretization is needed to align the production system’s output with the customer’s requests. The new model formulation presented here rests on a discrete time scale. The model’s periods are defined by the due dates unless a finer resolution is deemed necessary. A finer resolution needs to be applied if the periods defined by due dates are “too big” with respect to the problem-defining data (see Figure 1). As “too big” is not an exact measure, this term will subsequently be defined more precisely. The model formulation that will be proposed is based on a lot-sizing model formulation called PLSP (PLSP ) proportional lot-sizing and scheduling problem; Drexl and Haase9). PLSP is a planning model which decides upon lot sizes and the production sequence simultaneously. However, as it is a period-oriented model, it does not allow for the limitation of production quantities which overlap several periods but belong to one production run. It is based on the assumption that setup states can be carried over period boundaries and that in each period, at most, one changeover is allowed. Thus, if two products are to be produced in one period, the sequence is given, because the product that has been set up in the period is the second product. If at most one changeover is allowed in each period, at most two different products can be produced in each period. However, the difference between two due dates may be long enough to schedule three or more products, and this may also be optimal. Thus, a finer resolution is needed. On the other hand, defining the resolution too fine will result in a model formulation with too many variables and constraints to be solvable efficiently. Consequently, a time resolution needs to be looked for that allows all feasible solutions (including the optimal solution) to be found while, at the same time, using as few periods as possible. Some data defining the problem (e.g., minimal run lengths) can be used to derive a valid time discretization. For example, if the difference between two due dates is 1 month (30 days) and the

minimal run length on a resource is 10 days, then, at most, three changeovers will take place within this interval. Thus, a time discretization fulfilling the requirements will be to split the month into three periods, each with a length of 10 days. In general, the following time discretization is proposed: Each period s (with s denoting the time structure given by due dates) is divided into r periods of equal length, such that r is the smallest integer value equal to or greater than the capacity (duration) of s divided by the minimum of the minimal run lengths of all products that can be produced in period s. This time scale will be denoted by the index t. Of course, this time discretization may lead to a very fine resolution of time if the minimum of the minimal run lengths is short compared to the length of each period. Such a fine resolution of time resulting in models with many binary variables has been known to be a major shortcoming of all models based on a discrete time scale (e.g., Zentner et al.10). However, it is obvious to use these kinds of models only in those circumstances in which the nature of the planning problem (in Section 4, mid-term planning) allows for such a discretization of time. The notation given in the Nomenclature Section is introduced to describe the mathematical model. The mathematical model consists of an objective function (eq 1) and constraints (eqs 2-16).

Min

∑ ∑blpjtIBjt + ∑ ∑sspjtISSVjt + j∈J t∈u j∈J t∈u ∑ ∑hjsIj,[last(s-1)] + ∑ ∑ ∑ ∑ (hjs/2)Xmjt + j∈J s∈S j∈J s∈S t∈u m∈M s

∧j∈Jm

∑ ∑ ∑ ∑scmijYmijt i∈J j∈J m∈M t∈u

(1)

∧j∈Jm ∧i∈Jm

IBjt + Ijt-1 +

∑ Xmjt ) Ijt + IBjt-1 + djt

m∈M ∧j∈Jm

∀ j∈J, t∈u (2)

ISSVjt g sstj - Ijt XTmjt e cmt(Zmjt + Zmjt-1)

∀ j∈J, t∈u

(3)

∀ m∈M, j∈Jm, t∈u\{0} (4)

Ymijt g Zmjt + Zmit-1 - 1 ∀ m∈M; i, j∈Jm; i*j; t∈u\{0} (5)

∑ Zmjt ) 1

∀ m∈M, t∈u

(6)

j∈Jm

∑ Ymijt e 1 - Zmjt-1

i∈Jm i*j

∀ m∈M, j∈Jm, t∈u (7)

8316

Ind. Eng. Chem. Res., Vol. 44, No. 22, 2005

∑ XTmjt ) cmt

∀ m∈M , t∈u

(8)

j∈Jm

∀ m∈M, j∈Jm, t∈u (9)

minratemjXTmjt e Xmjt

∀ m∈M, j∈Jm, t∈u

Xmjt e maxratemjXTmjt

∀ m∈M, j∈Jm\{0}, t∈u (11)

Kmjt e Kmjt-1 + XTmjt Kmjt g Kmjt-1 + XTmjt -

(10)

∑ maxrlmjYmijt+1

i∈Jm i*j

∀ m∈M, j∈Jm\{0}, t∈u\{T} (12) Kmjt e maxrlmj(1 -

∑ Ymijt+1)

i∈Jm i*j

∀ m∈M, j∈Jm\{0}, t∈u\{T} (13) Kmjt-1 + XTmjt g minrlmj

∑ ∑ Ymikt

i∈Jmk∈Jm i*j k*i

∀ m∈M, j∈Jm\{0}, t∈u (14) Ijt g 0, IBjt g 0, IBj0 ) 0, ISSVjt g 0, Kmjt g 0, Xmjt g 0, XTmjt g 0, Ymijt g 0 ∀ j∈J, t∈u (15) Zmjt ∈ {0;1}

∀ m∈M, j∈Jm, t∈u

(16)

The objective is to minimize a cost function consisting of penalty costs for backlog (term 1) and safety stock violations (term 2), inventory holding costs (term 3 and 4), and changeover costs (term 5). As the penalty costs should not be charged more often if a finer time resolution is chosen, the associated cost coefficients (blpjt and sspjt) may be set to zero if the end of period t does not coincide with a due date. Inventory costs are assessed such that a full holding cost coefficient (hjt) is charged if products reside in the inventory for the complete period, but half the cost is charged if the product has only entered the inventory in that period. This scenario assumes that products leave the inventory only at the due dates (end of a period) and the inflow in the inventory is spread evenly over the period. This calculation takes into account only the time discretization imposed by due dates (periods s). Although a more precise calculation of holding costs is possible taking into account the time discretization (periods t) used in this model, the former calculation is employed to allow for a fair comparison to other model formulations from the literature, because the evaluation of holding costs based on a time discretization based on due dates has been used there. Finally, changeover costs are defined as sequence-dependent. Constraint 2 defines inventory flow constraints. Inflows are production in period t, inventories at the beginning of t, and the backlog at the end of period t, while the outflows comprise demand, inventory at the end of t, and the backlog at the beginning of t. Constraint 3 calculates the safety stock violations on the basis of the current inventory position Ijt and the safety stock target sstj. Constraints 4-7 are typical PLSP constraints. At the end of each period on each resource, a distinct setup state is defined by constraint 6. If the resource has been set up for product i at the end of period t-1 (Zmit-1 ) 1)

and it is set up for product j at the end of period t (Zmjt ) 1), then a changeover from product i to j must have taken place in period t on resource m (Ymijt ) 1) (constraint 5). Constraint 7 forbids changeovers to products that are already set up, while constraint 4 allows only the production of products that are either set up at the beginning of a period or at the end of a period. Constraint 8 allocates the available capacity (measured in units of time) of each resource in each period to the products that are allowed to be produced in the period on that resource. However, at most, two of the production variables XTmjt on the left-hand side can differ from zero because of setup constraint 4, in accordance with the basic assumption of the underlying PLSP model. As idling of resources is possible, a dummy product (j ) 0) is included in the set of products to allow for idle time. Constraints 9 and 10 transform the production time variables XTmjt into mass units Xmjt respecting minimal and maximal production rates. Finally, constraints 11-14 are used to observe additional restrictions of the production system. Here, each resource requires a minimal run length for each production run. Therefore, the production time of product j of subsequent periods (XTmjt) is accumulated into variables Kmjt unless a changeover to product j takes place in period t (constraints 11 and 12). If any changeover takes place in t, the minimal run length (minrlmj) has to be met (constraint 14). Constraint 13 ensures that Kmjt is reinitialized to zero before a new production run starts. Constraints 15 and 16 state the domain of the variables. In this model formulation, only setup state variables Zmjt need to be defined as binary variables. Variables Ymijt take only binary values in feasible solutions due to the objective function (eq 1) and constraints 5 and 7. But note that because of constraint 5, variables Ymijt are only bounded from below. This means they must take the value 1 whenever a changeover occurs in a period. However, they may take a value greater than zero also on other occasions, for example, if there is no cost associated with the corresponding variable (e.g., for a changeover to idle in our test set). In fact, these changeovers do not happen, because they can be filtered out by applying constraint 17 to all changeovers. However, for the construction of a feasible (optimal) solution, variables Ymijt are not necessary, because variables Zmjt suffice to define the setup pattern.

Ymijt e Zmjt

∀ m∈M; i, j∈Jm; i*j; t∈u

(17)

4. Computational Studies 4.1. Description of the Test Instance and Customization of the Modeling Approach. A test instance introduced into the literature by Karimi and McDonald11 will be used to assess the performance of the model formulation presented. This test instance has also been solved by the models of Ierapetritou et al.12 and Lee et al.13 in their papers, which will serve as benchmarks here. The test instance describes a multiproduct plant with 7 resources and 14 products working 7d/24h. Starting from an initial plant state, the aim is to derive a detailed schedule for the next 3 months, which needs to be more detailed in the first month. Customer demand (due dates) awaits fulfillment at the end of the first four weeks and at the end of months two and three. A


minimal run length of 10 days is required for each product on each resource. Furthermore, test and maintenance operations (called “planned outages”) have to be scheduled. These have a defined duration and must take place in certain time windows. Other problem data include initial inventory levels, safety stock targets, production rates, production suitability, and cost coefficients. All data are taken directly from the paper of Karimi and McDonald.11 However, their data include changeover costs for changeovers from product 7 (8) to product 8 (14), but product 8 cannot be produced on any resource that is able to produce products 7 or 14. Thus, these changeovers are not possible. Consequently, these cost coefficients will be ignored here. As Karimi and McDonald11 and the other authors12,13 who used this test instance pointed out, the scheduling of this plant can be done in five subproblems, because the assignment of products to resources allows this kind of decomposition. These subproblems will be called examples 1-5 in accordance with refs 11-13. With respect to the time structure, the time discretization as proposed above is applied. For most subproblems, this yields 10 periods with the first month (28 days) consisting of four periods (7 days each) and the second and third months (30 days) consisting of three periods (10 days each). The planned outages can be dealt with easily. In the new model, planned outages are simply products with some special characteristics, that is, products, for which backlog is not allowed (IBjt ) 0), minimal and maximal run length (minrlmj and maxrlmj) equal their duration, and production outside the assigned time windows is not allowed (corresponding variables Zmjt ) 0). Furthermore, as their time windows may only partially lie in a period, their associated run-length variables (XTmjt) may be (upper-) bounded. Furthermore, the time structure chosen may be affected by the duration of planned outages. For example, test operation 1 has a duration of 4.166 days and needs to be scheduled between days 15 and 30. Thus, two changeovers may be possible in weeks 3 (days 1521) or 4 (days 22-28). In line with the proposed time discretization algorithm, periods are split into r periods of equal length with r being the smallest integer equal to or greater than the period length (7 days) divided by the smallest minimal run length (4.166). Consequently, these two weeks are split into r ) 7/4.166 ) 1.68 ) 2 periods, each with a length of 3.5 days. Finally, inventory variables Ij0, setup state variables Zmj0, and accumulation variables Kmj0 need to be set according to the initial plant state. 4.2. Analysis of Solutions. In this section, the optimal solution of the test instance of Karimi and McDonald11 is analyzed in detail and compared to solutions presented by other authors,11-13 while in the next section, computational performance is assessed and compared to the same benchmarks. This is a necessary exercise, because different authors presented different optimal solutions (objective function values, tables, and Gantt charts), sometimes without even commenting on this issue, and of course all of them used different platforms (hardware and software) to obtain their results. Thus, in our view, a fair comparison of these different approaches (models) has not been done yet. In Karimi et al.14 and Sundaramoorthy and Karimi,15 the prerequisites of a fair comparison have been extensively discussed. It is shown that different hardware and

Table 1. Optimal Objective Function Values and Their Components (Suerie7) example 1 2 3 4 5

objective holding setup backlog safety stock ($) costs ($) costs ($) penalty ($) penalty ($) 128003 16704 350816 797129 43272

4889 2731 4713 1495 2301

8400 8600 16600 4800 12200

114713 5373 319303 760348 19052

0 0 10200 30486 9720

software strongly influence computational results even if varied independently. Here, the term “fair” is based on the computational study and reflects that all model formulations have been tested on the same data set with the same solver on the same hardware. Consequently, we coded all model formulations to run them on the same hardware (Pentium 4, 1.7 GHz, 256 MB RAM) and with the same commercial solver (Xpress-MP, release 2003G) to obtain our results. However, we detected that all published model formulations11-13 needed minor modifications to deliver the optimal solution to the test instance. These modifications are explained with the discussion of the optimal solutions as well as in the Supporting Information (A.III, B.III, and C.III), which presents all benchmark models in their entirety. However, we could not reproduce all solutions provided in the literature, but this may also be due to the fact that the different authors might have used slightly different data sets. The optimal solutions (objective function values and Gantt charts) are presented in Table 1 and Figure 2. But also note that schedules other than that depicted in Figure 2 might be optimal, if they have the same cost. This is especially true with respect to different resource utilization combinations. For example, on resource 1 in month 3, any other combination of production of products 3, 11, and idle is optimal as long as the minimal run length for product 11 is fulfilled. Thus, the production run of 11 might be shortened or extended within a certain range. Moreover, if the production rate of 3 was reduced (halved), the resource would not be idle at all. From this point on, we regard solutions as being the same if they have the same objective function value (cost) and produce the same amount of material (mass) in the same periods (defined by due dates). Going into more detail, regarding example 1 (resource 1), all benchmarks report slightly different optimal objective function values (125 602,11 125 603,12 and 128 00213). However, all report a Gantt chart which is compatible with our solution as well as our implementations of the benchmarks given the objective function values presented in Table 1. Thus, these differences might be due to rounding differences or a slightly different test data set being used. With respect to example 2 (resource 2), again, all authors provide different objective function values (16 201,11 16 138,12 and 16 73713). However, this time, the Gantt charts also differ. The schedule of Karimi and McDonald11 starts the second production run of 10 in week 4 instead of month 2, as proposed by the other benchmarks12,13 and our solution. However, this has only a small cost effect and might be due to the fact that they stopped the optimization run with a tolerance of $100 of the optimal solution. More striking is the fact that none of the authors found out that also postponing the first production run of product 10 can save holding costs. Postponing this production run is possible, be-

8318


Figure 2. Optimal Solution to Test Instance (Suerie7).

cause the resource is idle at the beginning of the planning interval anyway. The original model formulation by Karimi and McDonald11 is not able to find the “true” optimal solution, because one of their constraints forbids idle time in a period, if in a later period backlogging occurs. From first sight, this seems to make sense economically: Why should a machine run idle if a customer does not receive his order? On the other hand, if this order is too small to establish an additional production run, which possibly blocks the resource for a minimal run length, it might be cheaper to backlog this order for a short amount of time, as is the case here. Also, for example 3 (resources 3 and 4), all authors report different objective function values (350 257,11 350 216,12 and 345 94413) and Gantt charts. Again, Karimi and McDonald11 suffer from their additional constraints, which do not allow idle time on a resource if any other product is backlogged. Here, they schedule product 14 at the end of the third month to satisfy these constraints. Ierapetritou et al.12 present the same Gantt chart as shown in Figure 2. Thus, they might have used a slightly different data set. The solution presented by Lee et al.13 obviously suffers from two failures. First, scheduling test 2 directly after test 1 in contrast to the schedule presented in Figure 2 saves $5000 in transition costs but incurs $5360 in backlog penalties. Second, they report an ending inventory for product 7 well above the safety stock target (7348 vs 1080). Producing less, which is allowed as the minimal run length is already fulfilled, would reduce holding costs by $130. For example 4 (resources 5 and 6), again, three different objective function values are reported in the literature (794 385,11 794 386,12 and 800 27813). The original model formulation by Karimi and McDonald11 switches the first lot of product 12 from resource 5 to resource 6. As both resources are identical (with respect to production rates), this solution is an alternative optimal solution and, therefore, is considered to match ours, as depicted in Figure 2. The solution of Ierapet-

ritou et al.12 shows production of product 12 on both resources in the first week. This is not optimal, because holding cost coefficients are lower for product 5 than for product 12. However, the savings associated with this difference is only $1.03 and might explain the different optimal objective value compared to that of Karimi and McDonald.11 Finally, the Gantt chart by Lee et al.,13 again, does not show the optimal solution and can be improved by inspection. First, they also schedule product 12 on both resources in the first week, and second, they propose to schedule the maintenance operation at the beginning of the third month, which incurs additional transition costs of $5000 but does not yield any benefit. Presumably, the reason is initialization constraint 23 of their paper, which forbids outages at the end of the planning horizon. On the other hand, they should have obtained the correct optimal solution if they had modeled a product for idle time in their tests, because, then, a transition from outage to idle would have occurred at the end of the planning horizon. From the solution statistics (number of binary variables) they report, we suspect that they only include an idle product in their test data if this is necessary in the optimal solution. Thus, they omitted the idle product, if possible, to reduce the count of binary variables. This is most obvious if examples 1 and 2 are compared but also holds true for the other examples. There, they report the same number of binary variables for both examples, although example 1 involves the scheduling of three products (and no idle time necessary), while example 2 involves only the scheduling of two products (plus idle time). However, as one usually does not know the optimal solution beforehand, the dummy product for idle time should be modeled in each example. For the last example, two different objective function values have been reported in the literature (42 07211,12 and 43 27213). However, the optimal solutions (Gantt charts and tables) are all compatible with the one presented in Figure 2.

Ind. Eng. Chem. Res., Vol. 44, No. 22, 2005 8319 Table 2. Solution Statistics (Example 1) Karimi and McDonald constraints continuous variables binary variables nonzero elements nodes LP after cut generation time (s)

Ierapetritou et al.

Lee et al.

new

initial

presolve

initial

presolve

initial

presolve

initial

presolve

556 432 64 1985

468 404 64 1870 147 20487 3

1171 455 48 4624

813 225 144 2729 128 15318 2

489 177 120 1533

430 144 116 1360 85 24381 1

460 323 40 1695

406 267 40 1542 250 20263 1

Table 3. Solution Statistics (Example 2) Karimi and McDonald constraints continuous variables binary variables nonzero elements nodes LP after cut generation time (s)

Ierapetritou et al.

Lee et al.

new

initial

presolve

initial

presolve

initial

presolve

initial

presolve

423 273 48 1353

342 238 48 1263 33 13452 0

657 284 32 2719

419 150 71 1487 43 7555 0

325 125 72 957

277 99 69 819 7 16000 0

297 199 30 939

264 164 30 849 31 14938 0


Ierapetritou et al.

Lee et al.

new

initial

presolve

initial

presolve

initial

presolve

initial

presolve

1473 1459 158 6038

1023 1001 152 4458 193 327987 17

3225 1007 126 11992

2221 432 374 7124 7353 264369 102

1327 443 362 4519

1097 273 349 3682 155 333854 7

1918 1488 125 9689

1436 970 125 7360 1692 327016 23

Finally, we would like to point out that the model formulation by Lee et al.13 contains another flaw, although this one does not affect any of the solutions to examples 1-5. In their model, constraints for subtour elimination are missing (see the Supporting Information). The impact of this flaw can be explained by means of an example: Assume the changeover from product 13 to 9 was easy and had cost only $600 instead of $3000, then the optimal solution of the model by Lee et al.13 would propose changeovers from 13 to 9, 8 to test 3, test 3 to idle, and idle to 8 in month 3 on resource 7. These changeovers cost only $5600, but as they contain a cycle (8 f test 3 f idle f 8), this schedule is not feasible and cannot be implemented. The optimal sequence costs $8000 in changeover costs and involves transitions from 13 to 8, 8 to test 3, and test 3 to 9. The model will yield the correct solution, if standard subtour elimination constraints known, for example, from the traveling salesman problem (e.g., Miller et al.16) are added. We would like to point out that our implementation of all four model formulations [the new one and the three (slightly modified) benchmarks, see the Supporting Information] finds the same optimal solution to the data set used in our tests. 4.3. Computational Performance. In this section, the computational performance of the new model formulation and the three benchmarks will be assessed. Note, that here our implementations of the four models are to be compared. Thus, the model sizes (number of constraints and the variables) differ from what the other authors reported in their papers. However, the benchmark model formulations have only been altered if the published original version has been shown to be incorrect or incomplete (see Section 4.2 and the Supporting Information). The most remarkable difference (from the

literature) is the huge number of binary variables in the model by Lee et al.13 However, these are necessary if subtour elimination constraints are added. Otherwise, the model by Lee et al.13 would come with the second lowest number of binary variables of the models tested but would not guarantee the finding of an implementable solution in every case (see counterexample in Section 4.2). Furthermore, the solution times of the model by Ierapetritou et al.12 are much longer here than reported in their paper. They have made use of priorities to guide the branch-and-bound process. Nothing similar has been done here, because we found the standard settings of the MIP solver to be rather robust. Tables 2-6 provide solution statistics to examples 1-5. Each table provides data for one example for all four models tested. Given are the size of each individual model (number of constraints, binary and continuous variables, and nonzero coefficients in the matrix) as well as the number of nodes and the time (in seconds) to obtain the optimal solution and prove optimality. As standard solvers are able to improve the initial (given) model formulation considerably, reduced model sizes (after presolve) and the linear programming (LP) relaxation (after cut generation) are also given. For the models by Karimi and McDonald11 and Ierapetritou et al.,12 the number of time slots (event points) needs to be determined before the model is solved. If the number is chosen to be too big, the problem size is inflated unnecessarily, while the optimal solution might not be found if this number is chosen to be too low. One can choose an iterative procedure to determine the optimal number of time slots per period, but this would require the solving of each model several times. Different procedures on how to determine the “best” number of slots have been discussed extensively in Sundaramoorthy and Karimi.15 Here, we chose one

8320



Ierapetritou et al.

Lee et al.

new

initial

presolve

initial

presolve

initial

presolve

initial

presolve

946 641 107 3237

684 508 107 2719 21 795040 2

1333 480 75 5459

826 235 143 2885 3432552 542364 22424

694 226 176 2131

581 164 169 1785 11 796700 0

718 461 67 2515

598 356 67 2155 330 791618 2


Ierapetritou et al.

Lee et al.

new

initial

presolve

initial

presolve

initial

presolve

initial

presolve

703 608 67 2689

487 414 67 1957 453 18262 4

1323 500 52 4960

865 235 133 2825 2602 2434 12

567 236 130 1850

442 140 124 1409 125 24300 2

655 475 42 2743

465 321 42 2024 374 15740 1

setting that provided the optimal solution to all examples, but generally, this may understate the true computational effort associated with these two models. In the model formulation by Karimi and McDonald,11 two time slots have been assigned to each week and four time slots to each of the monthly periods. In the model formulation by Ierapetritou et al.,12 two event points have been assigned to each week and three (four) event points to the second (third) month. As these examples are rather small problems, solution times are (with the model by Ierapetritou et al.12 being an exception) rather short. The model by Ierapetritou et al.12 exhibits the smallest lower bounds after cut generation, which leads to the largest search trees (number of nodes) and longest solution times, and leaves some space for improvement by applying a scheme to prioritize variables within the branch-and-bound search. The other model formulations show an approximately equal performance regarding the time to prove optimality of the solutions, with the new model being the quickest if only one resource has to be scheduled (examples 1, 2, and 5), while the model of Lee et al.13 solves the quickest if there are two parallel resources available (examples 1-4). With respect to the model size, the model by Lee et al.13 exhibits the smallest number of variables (binary plus continuous) as well as the smallest number of constraints if there are parallel resources, while the new model formulation exhibits the smallest number of constraints if only one resource needs to be planned for. Presolve is most effective for the model by Ierapetritou et al.,12 with model sizes reduced by roughly one-third. Interestingly, the number of binary variables increases in the presolving stage for this model. According to Dash Optimization, the reason for this is that the presolve procedure finds variables in the model which are implied integers (binary). Promoting these variables to integers is then beneficial for branching and cutting. However, for the other models, remarkable reductions are achieved also. Most of the binary variables of the model by Lee et al.13 are induced by subtour elimination constraints. These do not seem to be used for branching. This explains why their model solves very fast despite the huge amount of binary variables present in their model formulation. The number of binary variables, which is an indicator for the scalability (worst case

performance for branching) of the model formulation, is lowest for the new model formulation. 4.4. Example Modifications. To judge the robustness of the proposed model formulation, some further tests have been carried out. The first test was aimed at a more precise production plan. In this test, the second month is planned for in weekly buckets. Consequently, the orders of the second month have been split into four weekly periods while retaining the proportional pattern of orders of the first month. To retain the initial time structure, each week lasts 7.5 days in month 2 (30 days divided by 4). In a second (third) test, the flexibility of the production equipment has been reduced additionally; in the second test, the minimal production rate has been set to be 80% of the maximal production rate, while in the third test, a constant production rate matching the maximal production rate of each product has been assumed. The resulting optimal solution to the third test is depicted in Figure 3 and Table 7. We would like to point out that, again, all models (the new model and our implementation of the three benchmarks) found the same optimal solutions to all modified examples. As demands occur earlier in the modified test instance compared to those of the original data set (some demand of month 2 must be met already at the end of week 5), production is shifted toward the beginning of the planning interval. Furthermore, production runs tend to be shorter, because this allows more different products to be produced early in the schedule. With respect to the cost components in the optimal solutions, penalties for safety stock violation and backlog tend to be equal or higher than before, because backlogging is penalized at the end of each period. While the same is true for setup costs, holding costs are sometimes reduced (examples 4 and 5). This is the result of more frequent changeovers caused by the new demand pattern. Computational results of these tests are provided in Table 8. The different models are ranked according to the time required to prove the optimality for each (modified) example (1 ) fastest; 4 ) slowest). The mean rank is chosen as a measurement of computational speed to account for the difference in the structure of the test problems. If, for example, solution times were averaged directly, this would have heavily biased the


Figure 3. Optimal Solution to Modified Test Instance. Table 7. Optimal Objective Function Values and Their Components of Modified Example example

objective holding setup backlog safety stock ($) costs ($) costs ($) penalty ($) penalty ($)

1 (mod) 2 (mod) 3 (mod) 4 (mod) 5 (mod)

138754 17661 418042 1236208 60320

5597 3688 5455 1132 1883

10800 8600 19600 8400 12800

120412 5373 381091 1196190 35918

1944 11896 30486 9720

Table 8. Mean Ranks (Time to Prove Optimality) of Different Models (1 ) Fastest) McDonald11

Karimi and Ierapetritou et al.12 Lee et al.13 new

examples 1-5

test 1

test 2

test 3

2.4 3.2 1.2 1.6

2.2 3.2 2.4 1.6

2.4 3.4 2.0 1.6

3.2 3.0 2.0 1.2

outcome based on the results of the most difficult test instance (example 3). While the model by Lee et al.13 solved fastest for the original data set, it is only placed third if a more detailed plan needs to be derived (test 1). The model by Karimi and McDonald11 solves relatively fast if the second month is divided into four weekly periods and if the production rate is restricted to be at least 80% of the maximal rate. However, if the production rate is fixed (test 3), this model performs the worst on the basis of the mean rank. Nonetheless, absolute running times of this model do not exceed 55 s, whereas the results of the model by Ierapetritou et al. are always distorted by an outlier for example 4 (>1800 s). The new model formulation solves fastest for all three tests on the basis of the mean rank.

of chemical plants. It has been shown that this model formulation compares favorably with model formulations presented in the literature on the basis of a test instance broadly used in the literature.11-13 However, the selection of a model formulation for a specific problem in practice has to be done with great care. Thus, the practical problem must meet the assumptions of the model and vice versa. Otherwise, the resulting “optimal” production plan might not be the optimal solution to the problem at hand. With respect to the new model formulation, basically two assumptions have to be met in order to achieve a reasonable result: First, time discretization must be chosen such that the minimal run length exceeds a period’s length. In this case, the fundamental assumption of the PLSP that there is, at most, one changeover in each period is always fulfilled. Second, as idle time is modeled as a product, this dummy product is also scheduled according to the fundamental PLSP assumption. This means that it is not possible to have a changeover from product A to idle and afterwards to product B in one period. However, as long as production rates are allowed to vary in rather wide ranges, such as in the original test instance, this will not affect solutions too much. Finally, to give some directions for further research, all four examined models still suffer from time discretization. First, all models allow for changes in production rates at due dates (e.g., see the end of week 1 on resource 1 in Figure 2), although for a majority of production environments, assuming that the parameters of running processes are changed to “save some holding costs” is not realistic. Moreover, one could think of a more continuous representation of holding costs.

5. Summary and Conclusions

Acknowledgment

A new model formulation has been proposed to deal with the short-term production planning and scheduling

The author thanks Bjo¨rn-Ragnar Weber for helping to implement the benchmark model formulations.

8322


Supporting Information Available: Benchmark model formulations of Karimi and McDonald,11 Ierapetritou et al.,12 and Lee et al.13 used in the computational tests and comments on their modifications compared to the original sources. This material is available free of charge via the Internet at http://pubs.acs.org. Nomenclature Indices and Index Sets i, j, and k ) products: i, j, and k ∈ J J ) set of products j Jm ) set of products j producible on resource m m ) resources: m ∈ M M ) set of resources m s ) periods (time discretization based on due dates): s ∈ S S ) set of periods s based on due dates t ) periods (time discretization based on algorithm): t ∈ u u ) set of periods t us ) set of periods t belonging to period s Data blpjt ) backlog penalty cost of one unit of backlog of product j at the end of period t cmt ) capacity of resource m in period t djt ) demand for product j due at the end of period t hjt ) holding cost for one unit of product j in period t last(s) ) last period t of the set of periods us belonging to period s maxrlmj ) maximal production quantity (run length) for product j on resource m maxratemj ) maximal production rate of product j on resource m minrlmj ) minimal production quantity (run length) for product j on resource m minratemj ) minimal production rate of product j on resource m scmij ) setup cost for a changeover from product i to product j on resource m sspjt ) safety stock penalty cost associated with the violation of the safety stock target for one unit of product j at the end of period t sstj ) safety stock target for product j Variables Ijt ) inventory of product j at the end of period t IBjt ) backlog of product j at the end of period t ISSVjt ) safety stock violation for product j at the end of period t Kmjt ) campaign variable for product j on resource m in period t (current campaign quantity up to period t) Xmjt ) production quantity of product j on resource m in period t XTmjt ) production time of product j on resource m in period t Ymijt ) changeover variable ()1 if a setup operation from product i to product j is performed on resource m in period t, )0 otherwise)

Zmjt ) setup state variable ()1 if product j is set up on resource m at the end of period t, )0 otherwise)

Literature Cited (1) Blo¨mer, F.; Gu¨nther, H.-O. LP-based Heuristics for Scheduling Chemical Batch Processes. Int. J. Prod. Res. 2000, 38, 10291051. (2) Kallrath, J. Planning and Scheduling in the Process Industry. OR Spectrum 2002, 24, 219-250. (3) Shah, N. Single- and Multisite Planning and Scheduling: Current Status and Future Challenges. Foundations of ComputerAided Process Operations. AIChE Symp. Ser. 1998, 94, 75-90. (4) Lamba, N.; Karimi, I. A. Scheduling Parallel Production Lines with Resource Constraints. 2. Decomposition Algorithm. Ind. Eng. Chem. Res. 2002, 41, 790-800. (5) Atamtu¨rk, A.; Savelsbergh, M. W. P. Commercial Integer Programming Software Systems. Research Report BCOL.03.01. University of California at Berkeley: Berkeley, CA, 2003. (6) Bixby, R. E.; Fenelon, M.; Gu, Z.; Rothberg, E.; Wunderling, R. MIP: Theory and PracticesClosing the Gap. System Modelling and Optimization: Methods, Theory and Applications; Powell, M. J. D., Scholtes, S., Eds.; Kluwer: Boston, MA, 2000, 19-49. (7) Suerie, C. Time Continuity in Time-Indexed Model Formulations. New Modeling Approaches for Production Planning in the Process Industries. Ph.D. Thesis. Technische Universita¨t Darmstadt, Germany, 2004. (8) Sivanandam, S. P.; Balla, G.; Sundaramoorthy, A.; Karimi, I. A. A Comparison of Continuous-time Models for Scheduling Noncontinuous Plants. AIChE Annual Meeting, November 7-12, 2004, Austin Convention Center, Austin, TX; Paper 425m. (9) Drexl, A.; Haase, K. Proportional Lotsizing and Scheduling. Int. J. Prod. Econ. 1995, 40, 73-87. (10) Zentner, M. G.; Pekny, J. F.; Reklaitis, G. V.; Gupta, J. N. D. Practical Considerations in Using Model-based Optimization for the Scheduling and Planning of Batch/Semicontinuous Processes. J. Process Control 1994, 4, 259-280. (11) Karimi, I. A.; McDonald, C. M. Planning and Scheduling of Parallel Semicontinuous Processes. 2. Short-Term Scheduling. Ind. Eng. Chem. Res. 1997, 36, 2701-2714. (12) Ierapetritou, M. G.; Hene´, T. S.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 3. Multiple Intermediate Due Dates. Ind. Eng. Chem. Res. 1999, 38, 3446-3461. (13) Lee, K.-H.; Heo, S.-K.; Lee, H.-K.; Lee, I.-B. Scheduling of Single-Stage and Continuous Processes on Parallel Lines with Intermediate Due Dates. Ind. Eng. Chem. Res. 2002, 41, 58-66. (14) Karimi, I. A.; Tan, Y. L.; Bhushan, S. An Improved Formulation for Scheduling an Automated Wet-Etch Station. Comput. Chem. Eng. 2004, 29, 217-224. (15) Sundaramoorthy, A.; Karimi, I. A. A simpler better slotbased continuous-time formulation for short-term scheduling in multipurpose batch plants. Chem. Eng. Sci. 2005, 60, 2679-2702. (16) Miller, C. E.; Tucker, A. W.; Zemlin, R. A. Integer Programming Formulation of Traveling Salesman Problems. J. ACM 1960, 7, 326-329.

Received for review December 21, 2004 Revised manuscript received June 25, 2005 Accepted August 16, 2005 IE048771P

Planning and Scheduling with Multiple Intermediate Due DatesAn

Recommend Documents