
Research Article


Development of Heuristics for MultiProduct MultiLevel Capacitated Lotsizing Problem with SequenceDependent Setups


M. Mohammadi,
S. M.T. Fatemi Ghomi,
B. Karimi
and
S. A.Torabi



ABSTRACT

This study considers the problem of multiproduct multilevel capacitated lotsizing and sequencing problem with sequencedependent setups. A Mixed Integer Programming (MIP) formulation of the problem is proposed which is impractical to solve in reasonable computing time for nonsmall instances. Reducing the dimensionality of the problem and allowing to solve larger instances, a modified mathematical model is developed which ignores majority of combinations. The ability to quickly find integerfeasible solutions for nonsmall instances is another aspect of this paper. Hybrid methods that mixes rollinghorizon approach and heuristic are developed. Heuristic is used to determine binary variables of current period. To test the accuracy of hybrid methods, a procedure for obtaining a lower bound on the optimal solution is developed. The tradeoffs between objective values and computing times are also provided.







INTRODUCTION
Lotsizing and scheduling problems have been an area of active research starting
from the seminal study of Wagner and Whithin (1958).
Since then there has been a considerable amount of investigation in order to
incorporate other important 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
(Fandel and StammenHegene, 2006; Merece
and Fonton, 2003; Karimi et al., 2003).
Models of lotsizing and scheduling are divided in the literature into small
bucket and big bucket problems (Eppen and Martin, 1987).
Small bucket problems break the planning horizon in small time periods such
that at most one product can be manufactured in a single period. 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.
Small bucket problems for multilevel multiproduct production are the MultiLevel
Discrete Lotsizing and Scheduling Problem (MLDLSP) and the MultiLevel Proportional
Lotsizing and Scheduling Problem (MLPLSP) (Kimms, 1996).
Both models enable simultaneous lotsizing and scheduling, but limit the number
of products to be manufactured in a period.
The MultiLevel Capacitated Lotsizing Problem (MLCLSP), a big bucket problem,
does not have this disadvantage, but it can not fix lotsizes and schedules simultaneously.
To attempt to unite the advantages of the MLPLSP and MLCLSP, Fandel
and StammenHegene (2006) have made a model formulation based on twolevel
time structure (Fleischmann and Meyr, 1997) which enables
simultaneous lotsizing and scheduling for multiproduct multilevel job shop
production.
This study deals with the deterministic dynamic models with a finite planning horizon, where the production of several different products on seriallyarranged capacitated machines is concerned. The challenging problem of efficient lotsizing on a flow shop with sequencedependent setups is modeled using a new Mixed Integer Programming (MIP) formulation. Simultaneous lotsizing and scheduling is essential if sequencedependent setup costs and setup times occur during production.
Solving the singlelevel multiperiod CLSP with sequencedependent setups is
equivalent to solving multiple dependent TSPs (Gupta and
Magnusson, 2005).
Hence, like the TSP, the CLSP with sequencedependent setups also belongs to a set of problems that are called NPhard. That means it is very difficult to optimally solve large instances of the problem. In fact, the solution time rises exponentially as either the number of variables and constraints increase. The introduction of multilevel production 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.
In this study, mathematical model and solution approaches for multiproduct multilevel capacitated lotsizing problem with sequencedependent setups have been proposed.
PROBLEM DEFINITION AND MATHEMATICAL MODEL
Assumptions: Following assumptions are made for the problem of multiproduct
multilevel capacitated lotsizing with sequencedependent setups:
• 
Several products are produced on seriallyarranged machines in flow shop
structure. There are M stages of processing which occur in a linear fashion
from 1,...,M 
• 
Each machine is constrained in capacity 
• 
When the machines are setup, sequencedependent setup costs and times
accrue 
• 
The settingup 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. 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 positivetime (and cost) setups
but rather up to N such setups. The remaining zerotime (cost) setups are
modeling phantoms and do not exist in reality (Clark
and Clark, 2000; Clark, 2003). 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 can not be produced earlier in a period than the production
of its required component is finished. In other words, production on a production
level can only start if a sufficient amount of the product from the previous
production level is available; this is called vertical interaction 
• 
To guarantee the vertical interaction, idleness before each production
is defined with the help of shadow product (Fandel and
StammenHegene, 2006) 
• 
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. The starting setup configuration on all machines are the same 
• 
The triangle inequality holds, i.e., it is never faster to change over
from one product to another by means of a third product. In other words,
a direct changeover is at least as capacity efficient as going via another
product 
• 
Setup cost has the form w_{ijm}= f_{w}.S_{ijm}
where, f_{w} is opportunity cost per unit of setup time 
• 
Infinite buffers exist between stages and before the first stage and after
the last stage 
Mathematical model
Indices
i, j, k 
: 
Product type 
n, n`, n" 
: 
Designation for a specific setup No. 
m 
: 
Level of production 
t 
: 
Period 
Parameters
T 
: 
Planning horizon 
N 
: 
No. of different products 
M 
: 
No. of production levels 
bigM 
: 
A large real No. 
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_{jm} 
: 
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_{ijm} 
: 
Sequencedependent 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_{ijm}≥0
and for i = j, S_{ijm} = 0 
W_{ijm} 
: 
Sequencedependent 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_{ijm} ≥0 and for i = j, W_{ijm} = 0 
j_{0m} 
: 
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 nth setup on machine m in
period t is from product i to product j or
not 

: 
Quantity of product j produced between nth and (n+1)th setups on machine
m in period t 

: 
Shadow product: idle time (in quantity units) before production of product
j on machine m in period t in order to ensure that direct predecessor of
this product (production of product j on machine m1 in period t) has been
completed. In other words, the gap between nth setup and its production
(in quantity units), in order to guarantee vertical interaction 
Subject to:
j = 1,...,N, t = 1,...,T 
(2) 
j = 1,...,N, m = 1,...,M1, t = 1,...,T

(3) 
j = 1,..,N, n` = 1,...,N, n`` = 1,...,N,
m = 1,...,M1, t = 1,...,T 
(4) 
n = 1,...,N, j = 1,...,N, m = 1,...,M, t = 1,...,T

(6) 
n = 1,...,N, j = 1,...,N, m = 1,...,M,
t = 1,...,T 
(7) 
j ≠ j_{0m}, i = 1,...,N, m
= 1,...,M 
(8) 
n =1,...,N1, i = 1,...,N, m = 1,...,M, t = 1,...,T

(10) 
i = 1,...,N, m = 1,...,M, t = 2,...,T 
(11) 
In this model, Eq. 1 represents the objective function which
minimizes the sum of the sequencedependent 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 inflows to each node (j, m, t) is equal to outflows from that node. In
other words, the sum of stock on hand of product j on machine m at the beginning
of period t and production volume of product j on machine m during period t
is equal to the sum of the stock on hand of product j on machine m at the end
of period t and production volume of product j on machine m+1 during
period t.
Equation 4 guarantees within one period that product j on
machine m is produced before its direct successor (product j on machine m+1).
Left side of Eq. 4 is equal to the time between the beginning
of period t and the end of production of product j on machine m if production
of product j on machine m can take place between n`th and (n`+1)th setups in
period t, else it is a negative number. Right side of Eq. 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 production of product j on machine
m+1 can take place between n”th and (n”+1)th setups in period t, else it is
a big number.
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.
Equation 10 and 11 represent the relationship
between successive setups.
Equation 12 and 13 represent the type
of variables.
Equation 14 indicates that at the beginning of planning
horizon there is no onhand inventory.
DEVELOPMENT OF LOWER BOUNDS
So far, we have successfully formulated the problem. However, the formulation presented in the earlier section is not a practical approach to solve large instances of the problem. Present experiments show that computation time grows exponentially with the number of products, the number of machines and the number of periods in the planning horizon. Therefore, it is necessary to develop a computable lower bound in order to test the accuracy of the heuristics. Note that the heuristic solution is by definition an upper bound on the optimal solution.
Here, we obtain two lower bounds on the problem. First lower bound is achieved
by relaxing all binary variables and relaxing Eq. 4. The
latter is made by adding the following equation to the first lower bound.
Table 1: 
Comparison of lower bounds and exact optimal solutions.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 and the
original model 

Equation 15 is similar to the right hand side of Eq.
6, . In
Eq. 15, we aggregate by
summing over all n. Therefore second relaxation is a lower bound on the original
problem.
In order to ascertain the accuracy of these lower bounds, we performed many numerical tests. Table 1 shows the results of such tests in some instances of the problem with N = 3, M = 2 and T = 3.
To apply the Optimal Enumeration Method (OEM) on original problem and lower bounds, 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 these problems are extracted from the following uniform distributions:
b_{j,m}≈U(1.5,2), d_{j,t}≈U(0,180), h_{j,m}≈U(0.2,0.4), p_{j,m,t}≈U(1.5,2),
W_{ijm=} = S_{ijm}≈U(35,70), C_{m,t}≈U(200,300).
Table 1 confirms the advantages of the second lower bound, it is therefore used.
SOLUTION METHOD
Solving the original model
Solution procedure: Rollinghorizon 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 (Beraldi
et al., 2008; Araujo et al., 2007;
Araujo et al., 2008; Merce and Fontan, 2003; Clark,
2003; Clark and Clark, 2000). The approach initially
adopted decomposes the model into a succession of smaller MIPs, each with a
more tractable number of binary variables.
Relaxing all binary variables and determining the setup pattern of current period by a heuristic instead of solving a MIP is our approach to solve this hardtosolve MIP problem.
Present rollinghorizon approach decomposes the planning horizon into three
sections. For a given iteration k:
• 
The first section (beginning section) is composed of the (k1) first periods.
Within this section, setup pattern have been frozen by the previous iterations 
• 
The second section (central section) includes the kth period. In this
section the whole relaxed problem is considered and setup pattern of this
period is determined by heuristic 
• 
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. 
Each iteration consists of solving a linear programming problem. 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.
Setup state selection rule for current period: According to this method,
all binary variables of current period would be determined. Note that according
to Eq. 811 for each triple (n, m, t)
there is exactly one pair (i, j) for which
and for (i`, j`)≠ (i, j), .
In other words, for each triple (n, m, t) of current period, this heuristic
specifies a pair (i, j) which .
According to our method, if there is sufficient amount of inventory to satisfy demand of current period (I_{j,m,t–1}>d_{j,t}), product j in stage m would not be produced in current period.
The part of heuristic that is used to determine ordering of products in current
period is similar to (g/2,g/2) Johnson`s rullbased heuristic has been used
by Kurz and Askin (2003, 2004)
to schedule flexible flow lines with sequencedependent setups.
To order products in stages of current period, it is necessary to define job duration. Duration of a job in each stage is, by definition, the minimum time to produce demand of current period. These job durations only are used to determine sequences of products in current period and would not be used as lotsizes.
Let [i] indicates the ith job in an ordered sequence in the following. In the
following heuristic, a job duration is used. For product j in stage m and period
t it is specified by D (j, m, t) and is defined as:
D (j, m, t) = D(j, m1, t)+d (j, t).b(j, m)+min_{i}{S_{ijm}} 
Simplification strategy: More computatinal time is economized by eliminating
the majority of variables. ,
and
are
eliminated. Except Eq. 2, 3 and 5,
other constraints are ignored for ending section. All setup costs (and times)
for ending section are assumed to be 0.
b_{j,m} and p_{j,m,t} should be modified to estimate the capacity
consumption of future setups. We assume A_{1} is the objective value
of the second lower bound and A_{2} is the objective value of the original
problem by relaxing all binary variables, ignoring Eq. 4 and
all setup costs and times equal to 0. we would replace b_{j,m} and p_{j,m,t}
with b^{`}_{j,m} and p^{`}_{j,m,t} as follows:
A simplified representation for latter periods in the rolling horizon is less
difficult to solve and hence permits the solution of larger problems.
The whole algorithm (H1) All binary variables are relaxed
Begin:
t←1
while t≤T loop
for tt>t
• 
relax Eq. 4, 10 and
11 
• 
and are
equal to 0 
• 
and are
used instead of p_{j,m,tt} and b_{j,m} (for tt>t) 
In current period (t)
Create job durations D(j, m, t) as follows:
If I_{j,m,t–1}>d_{j,t}; D(j, m, t) = 0
else if ,
D(j, m, t) = D(j, m1, t)+d_{j,t}.b_{j,m} +min_{i(i≠ j)}
{S_{ijm}};
else D (j, m, t) = D (j, m1, t)+d_{j,t}.b_{j,m};
Create and
as follows:
,
Let U = and
V = .
The set U is the set of jobs that moves through stages 1 to [M/2] faster than
they move through stages [M/2]+1 to M. The set V is the set of jobs that moves
through stages [M/2]+1 to M faster than they move through stages 1 to [M/2].
Arrange jobs in U in nondecreasing order of and
arrange jobs in V in nonincreasing order of .
Append the ordered list V to the end of U, creating an ordered list to use in
the next step. Delete jobs which I_{j,l,t–1}>d_{j,t}.
Repeat the last member of the list unless the list contains N members.
If i = 1, ,
else,
m←2
while m≤M loop
Order the jobs in nondecreasing order (SPT) of D(j, m, t). Delete jobs which I_{j,l,t–1}≥d_{j,t}. Instead of deleted jobs, we replace the last true job.
if i = 1, ,
else, =1
m ←m+1
end loop
Solve LP and fix
for current period (t).
For t >1, i s
the last sequenced job in stage m and period t.
t← t+1
end loop.
Solving the modified model
The modified model: Reducing the dimensionality of the problem and
allowing to solve larger instances, a modified mathematical model is developed
which ignores the majority of combinations. In this model, in each period, products
have the same sequence and size in all stages.
Following proposes the modified model:
Subject to:
j = 1,...,N, t = 1,...,T 
(17) 
n` = 1,...,N, m = 1,...,M1, t = 1,...,T 
(18) 
m = 1,...,M, t = 1,...,T 
(19) 
n = 1,...,N, j = 1,...,N, t = 1,...,T 
(20) 
n = 1,...,N, j = 1,...,N, m = 1,...,M,
t = 1,...,T 
(21) 
j ≠j_{0}, i = 1,...,N 
(22) 
n = 1,...,N1, i = 1,...,N, t = 1,...,T 
(24) 
The whole algorithm (H2)
All binary variables are relaxed
Begin:
t←1
while t≤T loop
for tt >t
• 
relax Eq. 18, 24 and 25 
• 
for
(n>1), for
(n>1) and are
equal to 0 
• 
and are
used instead of P_{j,m,tt} and b_{j,m} (for tt>t) 
In current period (t)
Create job durations D (j, m, t) as follows:
If I_{j,t–1}>d_{j,t;} D (j, m, t) = 0
else if =
D(j, m1, t)+d_{j,t}.+b_{j,m}+min_{i(i≠ j)} {S_{ijm}}
Else D(j, m, t) = D(j, m1, t)+d_{j,t}.b_{j,m} ;
Create and
as follows:
Let U = and
V = . The set
U is the set of jobs that moves through stages 1 to [M/2] faster than they move
through stages [M/2]+1 to M. The set V is the set of jobs that moves through
stages [M/2]+1 to M faster than they move through stages 1 to [M/2].
Arrange jobs in U in nondecreasing order of and
arrange jobs in V in nonincreasing order of .
Append the ordered list V to the end of U, creating an ordered list to use in
the next step. Delete jobs which I_{j,t–1}>d_{j,t}.
Repeat the last member of the list unless the list contains N members.
If i = 1, ,
else, =1
Solve LP and fix for
current period (t).
For t>1, end_{t–1} is the last sequenced job in period t.
t←?t+1
end loop.
NUMERICAL EXPERIMENTS
To evaluate and compare the performance of developed heuristics, 20 problems
with different sizes are selected to test. For each problem size, 5 problem
instances are randomly generated and the required parameters for these problems
are extracted from the following Uniform distributions:
b_{j,m}≈U(1.5,2), d_{j,t}≈U(0,180),
h_{j,m}≈U(0.2,0.4), p_{j,m,t}≈U(1.5,2), W_{i,j,m}=
S_{i,j,m}≈U(35,70) 
C_{m,t} is calculated in accordance to satisfy demand of each period
on a justintime basis with average setups.
The heuristics are coded in Matlab programming language and are run on a personal computer with a Pentium IV, with a 3.4 GHz processor and a 4 GB of RAM.
Table 2: 
Comparison between objective functions of the second lower bound and
heuristics 

: Indicates that there is not enough memory to solve this
instance 
Table 3: 
Comparison between CPU times of the second lower bound and heuristics.
The values inside the brackets are the computational times in seconds 

*: Indicates that finding the optimum value for the second
lower bound requires more than 7200 seconds and the objective function at
this time has been considered, : Indicates that there is not enough memory
to solve this instance 
Table 2 and 3 compare the objective functions
and cpu times of heuristics and the second lower bound.
In accordance to the experiments there is not enough memory to solve instances with more than 5 products through the first algorithm (H1). The second algorithm (H2) is also able to solve instances with 10 products.
CONCLUSIONS AND RECOMMENDATION FOR FUTURE STUDIES
The contribution of the study has been to derive and test two models and two heuristics for multiproduct multilevel capacitated lotsizing problem with sequencedependent setups.
Heuristic H1 is based on the original model and is able to solve only small size problems. Heuristic H2 is based on the modified model that assumes similar sequences and sizes of products in all machines. H2 is able to solve medium size problems.
According to the experiments H1 is superior for small size problems. For larger problems that H1 is not able to solve, H2 is used.
Because of the expanding role of metaheuristic approaches to solve complicated
lotsizing problem (Jans and Degraves, 2007; Defersha
and Chen, 2008), the application of metaheuristic approaches to face this
hard to solve problem is recommended as an area for future research.

REFERENCES 
Araujo, S.A., M.N. Arenales and A.R. Clark, 2007. Joint rollinghorizon scheduling of materials processing and lotsizing with sequencedependent setups. J. Heuristics, 13: 337358. CrossRef  Direct Link 
Araujo, S.A., M.N. Arenales and A.R. Clark, 2008. Lot sizing and furnace scheduling in small foundries. Comput. Operat. Res., 35: 916932. CrossRef  Direct Link 
Beraldi, P., G. Ghiani, A. Grieco and E. Guerriero, 2008. Rollinghorizon and fixandrelax heuristics for the parallel machine lot sizing and scheduling problem with sequencedependent setup costs. Comput. Operat. Res., 35: 36443656. CrossRef  Direct Link 
Clark, A.R. and S.J. Clark, 2000. Rollinghorizon lot sizing when setup times are sequencedependent. Int. J. Prod. Res., 38: 22872307. Direct Link 
Clark, A.R., 2003. Optimization approximations for capacity constrained material requirements planning. Int. J. Prod. Econ., 84: 115131. CrossRef  Direct Link 
Defersha, F.M. and M. Chen, 2008. A linear programming embedded genetic algorithm for an integrated cell formation and lot sizing considering product quality. Eur. J. Operat. Res., 187: 4669. CrossRef  Direct Link 
Eppen, G.D. and R.K. Martin, 1987. Solving multiitem capacitated lotsizing problems using variable redefinition. Operat. Res., 35: 832848. Direct Link 
Fandel, G. and C. StammenHegene, 2006. Simultaneous lot sizing and scheduling for multiproduct multilevel production. Int. J. Prod. Econ., 104: 308316. CrossRef  Direct Link 
Fleischmann, B. and H. Meyr, 1997. The general lotsizing and scheduling problem. OR Spekt., 19: 1121. CrossRef  Direct Link 
Gupta, D. and T. Magnusson, 2005. The capacitated lot sizing and scheduling problem with sequencedependent setup costs and setup times. Comput. Operat. Res., 32: 727747. CrossRef  Direct Link 
Jans, R. and Z. Degraeve, 2007. Metaheuristics for dynamic lot sizing: A review and comparison of solution approaches. Eur. J. Operat. Res., 177: 18511875. CrossRef 
Karimi, B., S.M.T. Fatemi Ghomi and J.M. Wilson, 2003. The capacitated lot sizing problem: A review of models and algorithms. Omega, 31: 365378. CrossRef  Direct Link 
Kimms, A., 1996. Multilevel,singlemachine lotsizing and scheduling (with initial inventory). Eur. J. Operat. Res., 89: 8699. CrossRef  Direct Link 
Kurz, M.E. and R.G. Askin, 2003. Comparing scheduling rules for flexible flow lines. Int. J. Prod. Econ., 85: 371388. CrossRef  Direct Link 
Kurz, M.E. and R.G. Askin, 2004. Scheduling flexible flow lines with sequencedependent setup times. Eur. J. Operat. Res., 159: 6682. CrossRef  Direct Link 
Merece, C. and G. Fonton, 2003. MIPbased heuristics for capacitated lotsizing problems. Int. J. Prod. Econ., 85: 97111. CrossRef  Direct Link 
Wagner, H.M. and T.M. Whitin, 1958. Dynamic version of the economic lot size model. Manage. Sci., 5: 8996. CrossRef  Direct Link 



