MIP-based heuristics for lotsizing in capacitated pure flow shop with sequence-dependent setups

Lotsizing in capacitated pure flow shop with sequence-dependent setups has been considered in this paper. An exact formulation of the problem is provided as a mixed-integer program. It is well known that the capacitated lotsizing and scheduling problem (CLSP) is NP-hard. The introduction of serially arranged machines and sequence-dependent setups makes the problem even more complicated. Five MIP-based heuristics based on iterative procedures are provided. The first three heuristics are based on the original model but to solve non-small instances of problem, the last two heuristics are based on permutation flow shop problem which ignores the majority of combinations. To test the accuracy of heuristics, two lower bounds are developed and compared against the optimal solution. The trade-offs between solution quality and computational times of heuristics are also provided.

Keywords: lotsizing and scheduling; pure flow shop; sequence-dependent; MIP-based algorithm; relax-and-fix

Smooth and cost-efficient running of a factory often depends on its manager's ability to select appropriate lot sizes and production schedules (Gupta and Magnusson [14]). For this reason, such problems have been the focus of a large number of articles in the operations research literature.

Simultaneous lotsizing and scheduling is essential if sequence-dependent setup costs and setup times occur during production (Fandel and Stammen-Hegene [12], Almada-Lobo et al.[1]). Various models can be used to solve the problem of simultaneous lotsizing and scheduling. These differ with regard to their specific features. Among the characteristic features of the models for lotsizing and scheduling are the segmentation of the planning horizon, the time dependence of the model parameters, the information degree of the model parameters, the number of products and production stages, the production structure and the capacity restrictions (Karimi et al.[16], Merece and Fonton [20], Fandel and Stammen-Hegene [12]).

Models of lotsizing and scheduling are divided in the literature into small bucket and big bucket problems (Eppen and Martin [11]). Small bucket problems break the planning horizon in small time periods which limits the number of products manufactured in a single period. Small bucket problems are divided in the literature into discrete lotsizing and scheduling problem (DLSP), proportional lotsizing and scheduling problem (PLSP) and continuous setup lotsizing problem (CSLP) (Karimi et al.[16]). In DLSP at most one product can be manufactured in a single period, but PLSP allows up to two products (Drexl and Kimms [10]). Consequently, if a setup is performed, the entire time interval must be devoted to the setup. That is, setups and production runs comprise of an integer number of time intervals (see, for example, Tempelmeier and Buschkuhl [21]). Small bucket problems for multi-level multi-product production are the multi-level discrete lotsizing and scheduling problem (MLDLSP) (Kimms [18]) and the multi-level proportional lotsizing and scheduling problem (MLPLSP) (Kimms [17], Kimms and Drexl [19]). Both models enable simultaneous lotsizing and scheduling, but limit the number of products to be manufactured in a period.

The multi-level capacitated lotsizing problem (MLCLSP), a big bucket problem, does not have this disadvantage, but it cannot fix lotsizes and schedules simultaneously. To attempt to unite the advantages of the MLPLSP and MLCLSP, Fandel and Stammen-Hegene ([12]) have made a model formulation based on two-level time structure (Fleischmann and Meyr [13]), which enables simultaneous lotsizing and scheduling for multi-product multi-level job shop production. Because of the model's complexity caused by the great number of variables, only problems with a low number of products and periods can be solved optimally. This great number of variables is caused by the differentiation between three types of binary variables.

The challenging problem of efficient lotsizing on the pure flow shop with sequence-dependent setups is modelled using a new mixed integer programming (MIP) formulation in this paper.

Solving the single-level multi-period CLSP with sequence-dependent setups is equivalent to solving multiple dependent TSPs. Hence, like the TSP, the CLSP also belongs to a set of problems that are called NP-hard. That means it is very difficult to optimally solve large instances of the problem. The introduction of serially arranged machines makes the problem even more complicated. Therefore, it is necessary to find reasonable heuristic solutions for medium and large instances. Also it is important to develop a computable lower bound in order to test the accuracy of the heuristics.

The paper has the following structure. Section 2 introduces a detailed description of the problem and its underlying assumptions, Section 3 deals with the development and comparison of lower bounds in detail and Section 4 provides the heuristics. Section 5 reports the numerical experiments and finally Section 6 discusses the concluding remarks and recommendations for future studies.

The following assumptions are made for the problem of multi-product multi-level capacitated lotsizing with sequence-dependent setups:

• • Several products are produced on serially arranged machines in a pure flow shop structure.
• • Each machine can only produce one product type. Machine m can only process the mth level intermediary products that are required for production of the product type of level m + 1.
• • Each machine is constrained in capacity.
• • When the machines are setup, sequence-dependent setup costs and times accrue.
• • Setup costs have the form Wijm = fw · Sijm where fw is opportunity cost per unit of setup time.
• • The setting-up of a machine must be completed in a period.
• • There must be precisely N (number of products) setups in each period on each machine, even if a setup is just from a product to itself. Setups follow the so-called triangular inequalities, i.e. it is never faster to change over from one product to another by means of a third product: otherwise over N setups in a time period could be obtained. Since a setup time (and cost) from a product to itself is zero, note that the model does not force a machine to have exactly N positive-time (and cost) setups but rather up to N such setups. The remaining zero-time (cost) setups are modelling phantoms and do not exist in reality (Clark and Clark [6], Clark [5]). This feature makes possible for a lotsize, or production run, to continue over consecutive time periods without incurring real setup for later period (setup carry over).
• • The required resources and parts must be ready for production.
• • External demand exists for final products and is satisfied at the end of each period.
• • There are no lead times between the different production levels for transportation or cooling the products.
• • Shortages are not permitted.
• • A component cannot be produced earlier in a period than the production of its required component is finished, this assumption remains true even if there is inventory of the required component. In other words, production on a production level can only start if all products from the previous production level are available; this is called vertical interaction.
• • To guarantee the vertical interaction, idleness between each setup and its production is defined with the help of shadow product (Fandel and Stammen-Hegene [12]).
• • There are no demand and no storage costs for shadow products.
• • At the beginning of the planning horizon each machine is setup for a defined product.

In order to formulate the model, we rely on the following data:

Indices

• i , j , k Product type.
• Designation for a specific setup number.
• m Level of production.
• t Period.

Parameters

• T Planning horizon.
• N Number of different products.
• M Number of production levels/number of machines.
• BigM A large real number.
• C m , t Available capacity of machine m in period t (in time units).
• d j , t External demand for product j at the end of period t (in units of quantity).
• h j , m Storage costs unit rate for product j in level m.
• b j , m Capacity of machine m required to produce a unit of product (or shadow product) j (in time units per quantity units).
• p j , m , t Production costs to produce one unit of product j on machine m in period t (in money unit per quantity unit).
• S i , j , m Sequence-dependent setup time for the setup of the machine m from production of product i to production of product j (in time units); for i ≠ j , S i , j , m ≥ 0 and for i = j , S i , j , m = 0.
• W i , j , m Sequence-dependent setup cost for the setup of the machine m from production of product i to production of product j (in money units); for i ≠ j , W i , j , m ≥ 0 and for i = j , W i , j , m = 0.
• j 0 m The starting setup configuration on machine m.

Decision variables

• I j , m , t Stock of product j at level m at the end of period t.
• Binary variable, which indicates whether the n th setup on machine m in period t is from product i to product j (= 1) or not (= 0).
• Quantity of product j produced after n th setup on machine m in period t.
• Shadow product: the gap (in quantity units) between n th setup (to product j) on machine m in period t and its related production in order to ensure that direct predecessor of this product (production of product j on machine m − 1 in period t) has been completed. In other words, idle time (in quantity units) before production of product j on machine m in period t in order to guarantee vertical interaction.

Graph

Subject to

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

In this model, Equation (1) represents the objective function which minimises the sum of the sequence-dependent setup costs, the storage costs and the production costs. Equation (2) ensures the demand supply in each period. Equation (3) shows that in a network, total of in-flows to each node is equal to out-flows from that node.

Equation (4) guarantees within one period each typical product j on machine m is produced before its direct successor (product j on machine m + 1).

The left side of Equation (4) is equal to the time between the beginning of period t and the end of production of product j on machine m if n′th setup in machine m and period t is from every product i to product j (for n′ > 1, i ≠ j), else it is a negative number. In other words, if cannot get a positive value the left side of Equation (4) would get a negative value. The right side of Equation (4) is equal to the time between the beginning of period t and the beginning of production of product j on machine m + 1 if n″th setup in machine m + 1 and period t is from every product i to product j (for n″ > 1, i ≠ j), else it is a big number. In other words, if cannot get a positive value, the right side of Equation (4) would get a big value.

Equation (5) represents the capacity constraints of machines during periods. Equation (6) indicates that setup is considered in production process. Equation (7) indicates the relationship between shadow products and setups. Equations (8) and (9) guarantee that for each machine, the first setup at the beginning of the planning horizon is from a defined product. Equations (10) and (11) represent the relationship between successive setups. Equations (8) to (11) ensure that for each triple (n, m, t) there is exactly one pair (i, j) which . Equations (12) and (13) represent the type of variables. Equation (14) indicates that at the beginning of the planning horizon there is no on-hand inventory.

To illustrate the model, a small problem has been presented. Consider a production line with two serially-arranged machines that produces three products in two periods (N = 3, M = 2, T = 2). At the beginning of the planning horizon, machines are setup to produce product 1. Product demands are shown in Table 1. Unit production times and costs are shown in Table 2. Holding cost unit rates are shown in Table 3. Setup times and costs required by products in machines are shown in Table 4. Table 5 shows the capacity of machines.

Table 1. Demand of products.

Table 2. Unit production times and costs.

Table 3. Holding cost unit rate.

Table 4. Sequence-dependent setup times and costs.

Table 5. Capacity of machines.

Based on the data shown in Tables 1 to 5, the generalised algebraic modelling system (GAMS) is utilised to find the optimal solution. The optimal sequence in both machines in period 1 is 1-2-3 and lotsizes in both machines are 100, 50 and 200. It means that demand of product 3 in period 2 is produced in period 1. The sequence of both machines in period 2 is 1-2 and lotsizes in both machines are 100 and 150. The optimal objective function is 2510. Table 6 shows all non-zero variables in the optimal solution.

Table 6. Non-zero decision variables.

So far, we have successfully formulated the problem. However, the formulation presented in the previous section is not a practical approach to solve large instances of the problem.

In this section we obtain two lower bounds on the problem. The first lower bound is achieved by solving model M1 that is obtained from the initial model by relaxing all binary variables. After relaxing the binary variables, Equation (4) does not have considerable impact on the problem because for non-integer values of , the left side of Equation (4) would be negative and the right side of Equation (4) would be a big number. It means by relaxing binary variables, Equation (4) does not guarantee vertical interaction, therefore, Equation (4) would also be relaxed. The second lower bound is obtained by solving a new model M2 that is derived from M1, adding the following equation:

Graph

a j,m,t is a binary variable.

Equation (15) is similar to the right hand side of Equation (6),

Graph

In Equation (15), we aggregate

Graph

by summing over all n.

Equation (15) is valid to M1.

Suppose that Equation (15) is not valid to M1, it means that there is at least one triple (j, m, t), where

Graph

Suppose that

Graph

and = 1 (k ≠ j). For n = n′ to N, in all setups to product j, j would be replaced by k. By this modification, in triple (j, m, t),

Graph

According to the fact that setup costs from a product to itself is zero and considering the triangle inequality:

Graph

Therefore, by assuming that Equation (15) is not valid to M1, there would be a solution better than or equal to the optimum of M1 and it is impossible.

Rolling-horizon heuristics are usually used in dynamic lotsizing and scheduling problems, where demands are gradually revealed during the planning horizon and part types have to be allocated to machines in an on-going fashion as new orders arrive. On the other hand, a rolling-horizon approach is still suitable when all parameters are perfectly known (Clark and Clark [6], Clark [5], Merce and Fontan [20], Araujo et al.[2], [3], Beraldi et al.[4]). In this case, rolling-horizon heuristics have been used to overcome computational infeasibility for large MIP problems by substituting most of the binary variables and constraints with continuous variables and constraints. The approach initially adopted decomposes the model into a succession of smaller MIPs, each with a more tractable number of binary variables.

Each iterative procedure decomposes the planning horizon into three sections (Merce and Fontan [20]). For a given iteration k:

• • The first section (beginning section) is composed of the (k − 1) first periods. Within this section, the decisions have been partially or completely frozen by previous iterations according to a selected freezing strategy.
• • The second section (central section) includes the kth period. In this section the whole problem is considered.
• • The third section (ending section) includes the last periods (from period k + 1 to period T). There, the model is simplified according to a selected simplification strategy.

At the end of iteration k, all sections roll one period and a new iteration can then be performed. The procedure stops when there is no longer any ending section. The last iteration defines decision variables over the entire horizon. Figure 1 demonstrates the iterative procedure.

Graph: Figure 1. Demonstration of the iterative procedure.

Beginning section: All decisions related to the beginning section are completely frozen.

Central section: Consists of one period, the whole problem is considered.

Ending section: Binary variables and Equation (4) are relaxed for the ending section. Equation (15) is also considered for the ending section.

Using a simplified representation for the ending section in the rolling horizon is less difficult to solve and hence permits the solution of larger problems.

The solution approach initially adopted decomposes the model into a succession of T smaller MIPs, each with a more tractable number of binary variables.

Beginning section: Only binary variables related to the beginning section are frozen.

The central section and the ending section of this heuristic are similar to the former heuristic (H1).

The only difference between H3 and H2 is in their ending section. Even more computational time is economised by eliminating the majority of variables from the ending section. (n > 1), (n > 1) and are eliminated from the ending section. Except Equations (2), (3) and (5), other constraints are ignored for the ending section.

For the ending section, bj,m and pj,m,t should be modified to estimate the capacity consumption of future setups. We assume A1 is the objective value of the second lower bound and A2 is the sum of variable production costs of the second lower bound. For the ending section we would replace bj,m and pj,m,t with and to estimate consumption of setups:

Graph

Suppose the optimal policy were to meet the demand of every period on a lot-for-lot basis with no between-periods stocks. By ignoring setup costs and using this modified unit production costs and times, the objective function would be equal to the second lower bound.

Each iteration of former heuristics (H1, H2 and H3) has N3. M binary variables in its central section. Because the time taken to solve a MIP tends to explode exponentially as the number of binary variables increases, a faster solution approach is needed to make decisions by solving a succession of smaller MIPs.

In this heuristic (H4) the search is limited to the permutation flow shop problem, a restricted problem where products have similar sequences in all machines.

Assuming that in each period products have similar sequences in all machines, the number of binary variables in each iteration reduces from N3. M to N3.

Permutation flow shop problem

is used instead of because sequencing of products in all machines are the same. Other variables and parameters are similar to those of the model explained in Section 2.2. We also assume that the starting setup configuration of all machines are the same.

Graph

Subject to

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

Graph

The iterative procedure used in this heuristic is similar to H3.

Heuristic 4 (related to the modified model) could solve some instances that former heuristics (H1, H2 and H3) are not able to solve in limited time, but by increasing the size of the problem, the time taken to solve a MIP in H4 tends to explode exponentially.

An approach to find feasible solutions for larger instances is using fix and relax method (Dillenberger et al.[8]). Dillenberger et al. ([9]), formulated a lot sequencing and sizing model with representation of sequence-independent setup times on multiple machines. The resulting mixed integer programming (MIP) model is difficult to solve optimally for industrial problems, and so the authors resorted to the fix-and-relax method (Dillenberger et al.[8]), more widely known as relax-and-fix (Wolsey [22]). This involves the solution of a series of partially relaxed MIPs, each with a number of binary variables that is small enough to be quickly and optimally solved by conventional branch-and-bound methods. As the series progresses, each set of binary variables is permanently fixed at their solution values, and the relaxed variables are reduced in number, eventually disappearing. The procedure is broadly similar to a depth-first identification of an initial integer solution for a MIP model in a branch-and-bound search. Speed is its major advantage. Araujo et al. ([2]) and Beraldi et al. ([4]) are two recent applications of relax-and-fix method to solve MIPs. The formal procedure contains the following steps.

Solve the partial linear programming relaxation, where the values of {yj0,i,1¹∀i} are constrained to be 0 or 1, while the other may vary continuously between 0 and 1.

For n = 2, ..., N, solve the partial linear programming relaxation, with the {yj0,i,1¹∀i} fixed at their 0 or 1 solution values from step 1, with the values {yj,i,1²∀i} to {} fixed at their 0 or 1 solution values from previous applications of step 2, and with the values of {} are constrained to be 0 or 1, while the remaining binary variables may vary continuously between 0 and 1.

For t = 2, ..., T, repeat steps 1 and 2 for {yj,i,t¹∀i} to {}, fixing their values at those of the solutions in steps 1 and 2.

Note that each cycle of steps 1 and 2 involves solving N problems with just N2 binary variables each. N(N − 1) of N2 binary variables in each problem will newly have 0 due to Equation (25). Thus, the application of the cyclic fix-and-relax approach involves the solution of N.T MIPs with N2 binary variables each. As mentioned earlier, in each MIP of heuristics H1 to H3, there exists N3 · M binary variables and in each MIP of heuristic H4, there exists N3 binary variables. Therefore heuristic 5 (H5) is able to solve instances for quite sizeable values of N and T.

In order to ascertain the accuracy of the mentioned lower bounds, we performed some numerical tests. Tables 7 and 8 respectively show the results of such tests in some instances of the problem with (N = 3, M = 2, T = 3) and (N = 3, M = 3, T = 3).

Table 7. Comparison of lower bounds and exact optimal solutions in problem size N = 3, M = 2 and T = 3. The values inside the brackets are the computational time in seconds and the percentage values are the difference between the objective values of the lower bound against the original model.

Table 8. Comparison of lower bounds and exact optimal solutions in problem size N = 3, M = 3 and T = 3. The values inside the brackets are the computational time in seconds and the percentage values are the difference between the objective values of the lower bound against the original model.

Comparing the results of the second columns of Tables 7 and 8 shows that computation time grows exponentially by increasing the dimension of the problem. According to Table 7, average computational time for problems with (N = 3, M = 2, T = 3) is 168.87s. According to Table 8, average computational time for problems with (N = 3, M = 3, T = 3) is 5567.08s. It means that by increasing one level to production levels, average computational time increases more than 32 times. Tables 7 and 8 confirm the advantages of the second lower bound, therefore it has been used to compare heuristics.

To apply the exact model, lower bounds and heuristics, GAMS models are provided using GAMS IDE (ver 2.0.19.0) and solved using OSL 2. GAMS models are run on a personal computer with a Pentium 4 processor running at 3.4 GHZ. The required parameters for all numerical experiments are extracted from the following uniform distributions:

Graph

C m,t is calculated in accordance to satisfy demands of each period on a just-in-time basis with average setups.

To evaluate and compare the performance of developed heuristics, 20 problems with different sizes are selected to test. For each problem size, five problem instances are randomly generated and the required parameters for these problems are extracted from the mentioned uniform distributions.

In order to evaluate the quality of heuristics in an equal amount of time, the computational time of each MIP is limited to (7200/T) seconds for H1 to H4 and to (7200/T · N) seconds for H5.

Table 9 compares the solution times and objective values of heuristics and the second lower bound.

Table 9. Comparison of the second lower bound and heuristics. The values inside the brackets are the computational time in seconds.

In accordance to experiments, it is obvious that problems with dimensions greater than 5 * 5 * 5 can not be solved with H1 and H2. H3 is able to solve problems with dimensions of 5 * 5 * 7. H4 can solve problems with 7 * 7 * 10 and greater problems can be solved using H5.

The contribution of the paper has been to derive and test one exact formulation and five MIP-based heuristics all relied on successive resolutions of reduced size MIPs for lotsizing in capacitated pure flow shop with sequence-dependent setups.

Heuristics H1, H2 and H3 are based on the original model and are able to solve only small-size problems. Heuristics H4 and H5 are based on the restricted model that assumes similar sequences of products in all machines. H4 is able to solve medium-size problems, and H5 is recommended to solve larger problems.

In equal limited computational time, for small-size problems, the tests showed that H2 is superior. For medium-size problems, H4 is good and large problems can only be solved throuht H5.

One straightforward area for future research to extend the capabilities of H1, H2, H3 and H4 is to develop setup state selection rules to find production sequences in each period instead of solving an MIP.

Because of the expanding role of meta-heuristic approaches to solve complicated lotsizing problems (Jans and Degraves [15], Defersha and Chen [7]), the application of meta-heuristic approaches to face this hard to solve problem can be another area for future research.