INTRODUCTION
Since its formulation in 1915, the squarerootformula for the Economic Order
Quantity (EOQ) has been used in the inventory literature for a long time. This
formula is based on the assumption of a constant demand. The discrete case of
the dynamic version of EOQ was first discussed by Wagner
and Whitin (1958). Regarding the continuoustime dynamic EOQ models, Silver
and Meal (1969) were the first to suggest a simple modification of the classical
squarerootformula in the case of timevarying demand. However, Donaldson
(1977) discussed, for the first time, the classical noshortage inventory
policy for the case of a linear, timedependent demand. His treatment was fully
analytical and needed extensive computational effort to obtain the optimal solution.
The question of inventory shortages and backlogging were not considered at
all by the aforementioned researchers. Deb and Chaudhuri
(1987) were the first to incorporate shortages into the inventory lotsizing
problem with a linearly increasing timevarying demand. EOQ models for deteriorating
items with a trended demand were considered by Goswami and
Chaudhuri (1991, 1992), Xu and
Wang (1990), Chung and Ting (1993, 1994),
Kim (1995), Hariga (1995, 1996),
Benkherouf (1995), Jalan et al.
(1996), Jalan and Chaudhuri (1999), Giri
and Chaudhuri (1997) and Lin et al. (2000),
etc. In the model of Deb and Chaudhuri (1987), shortages are allowed in all
cycles except the last one. Each of the cycles in which shortages are permitted
starts with replenishment and ends with a shortage which is backlogged in the
next cycle. Numerous studies have been carried out to address the problems of
imperfect quality EMQ model with rework (Hayek and Salameh,
2001; Chiu, 2003; Chiu et al.,
2004; Jamal et al., 2004). Chiu
and Chiu (2006) studied optimal replenishment policy for an imperfect quality
EMQ model with backlogging and failure in repair using conventional approach.
That is to derive the optimal lot size by utilizing differential calculus on
the expected cost function with the need to prove optimality first. Grubbstrom
and Erdem (1999) first introduced an algebraic method to solve the classic
EOQ and Economic Production Quantity (EPQ) model without the use of derivatives.
A few researches have then been carried out using the same method (CardenasBarron,
2001; Wee and Chung, 2007). In these researches that
extend the algebraic approach to an imperfect quality EMQ model examined by
Chiu and Chiu (2006), the use of differential calculus
is replaced with an algebraic method and the optimal lot size solution is obtained
under the expected cost minimization.
Taleizadeh et al. (2008c) introduced an Economic
Order Quantity (EOQ) model with payment in advance to purchase highprice raw
materials. In this study we relax and change some of their assumptions. At first
we considered transportation cost as a linear function. Total discount policy
is considered instead of incremental discount one. Also we developed their model
based on two warehouses strategy to store raw material in which holding cost
is different for each of warehouses.
PROBLEM DEFINITION
Inventory holding as a tacticallevel decision against nonsecure situations
to increase confidence level of responding to customer demands and also
to resist against all kind of nondeterministic situations in receiving
raw materials made inventory control and management an important concept
in supply chain management. The model that is proposed in this research
is an applied model developed based on the real constraints and environments
of manufacturing companies. In these companies, the annual demands of
different products are first estimated at the start of the year. Then,
based on these estimates, the inventory and planning departments proceed
to material planning. In some cases, improper material planning and control
policies and loss of sufficient raw materials in proper times and quantities,
result to customer complaints and loss of customer and market share. In
some instances gaining the lostmarketshares needs more promotion and
advertising costs and causes increased production cost. To confront with
these instances and to minimize rawmaterial shortages, in this research
a model for material planning and control is developed.
In a manufacturing company, the steps involved in the ordering process
of the materials are:
• 
Forecast the number of finished products 
• 
Forecast the required rawmaterials 
• 
Order the required materials that consist of: 
• 
Releasing orders of the materials to a supplier 
• 
Paying a percentage (α%) of the purchasing cost at t_{0
} 
• 
Paying the remaining payment (1α)% at t_{Ic } 
• 
Transportating materials and receiving them at T 
Figure 1 shows the inventory control cycle of a material.
There are two payment methods in a material ordering process; (1) the
credit transaction with a specific lead time and (2) cash method that
has smaller lead time. In the cash method the purchasing cost of the total
ordered material is paid to the supplier at the ordering time. However,
in the credit transaction method α% of the material purchasing cost
is paid at the order release time (t_{0}) and the remaining (1α)%
is paid at the starting time of the material transportation (t_{lc}).
The ordered materials are received at T.

Fig. 1: 
Inventory control picture 
Let us assume that the credit transaction method is used for payments
in which the leadtime is deterministic and constant. Materials ordering
and their transportation are done in batch form and the amount of raw
materials in each batch and the number of batches in each vehicle is deterministic.
Based on the product groups’ consistency, the ordered materials can be
carried together by finitecapacity vehicles with maximum capacity of
F (in volume).
In this study, we plan to determine the optimal order quantity of each
material in a joint replenishment policy such that the total cost is minimized.
The costs are: (1) purchasing cost under total discount for each order,
(2) holding cost for on hand inventory (including capital, warehouse and
insurance costs), (3) capital cost of the next order, (4) transportation
cost, (5) clearance cost and (6) fixedorder costs. Furthermore, we assume
that the lead times are less than the cycle times of the materials.
All parameters of the problem are crisp and the quantity of the orders
must be integer multiples of packets, each containing more than one product.
Also we assume that the lead time is less than a cycle length. The constraints
of the problem are space, budget and upper limit for the number of orders
per year. Space constraint is different by common one incase if the total
required space in each cycle exceeds than capacity of warehouse, the company
should rent another warehouse.
PROBLEM MODELING
The parameters and the variables: For i = 1,2,…,n and j = 1,2,…,P
the parameters and the variables of the model are:
P 
= 
No. of the products 
D_{j} 
= 
Annual constant demand of the jth product 
n_{j} 
= 
No. of items in the packets of the jth product 
h_{1} 
= 
Holding cost per unit of on hand inventory during a period in first
warehouse (is same for all products) 

= 
Holding cost per unit of on hand inventory during a period in rented
warehouse (is same for all products) 
h2_{j} 
= 
Capital cost per unit of the jth ordered product during the period
before the payment of the remaining purchasing cost (10% of the total
purchasing cost) 
h3_{j} 
= 
Capital cost per unit of the jth ordered product during the period
after the payment of the remaining purchasing cost 

= 
Transportation cost for each unit of the jth product 

= 
Clearance cost for each unit of the jth product 

= 
Purchasing cost of the jth product in the ith discount break point 
q_{ij} 
= 
ith discount break point of the jth product 
ROP_{j} 
= 
Reorder point of the jth product 
f_{j} 
= 
Space required for each packet of the jth product 
L 
= 
Constant joint lead time for each order 
A 
= 
Fixed order cost per each order 
Q_{j} 
= 
Decision variable representing the order quantity of the jth product 
m_{j} 
= 
Decision variable representing the number of packets that have
been ordered for the jth product 
TB 
= 
Total available budget 
T 
= 
Decision variable representing the joint cycle length 
N 
= 
No. of orders in each year (N = 1/T) 
N_{T} 
= 
Upper limit for number of orders 
C_{H} 
= 
Annual total holding cost of the products 
C_{c} 
= 
Annual total clearance cost 
C_{P} 
= 
Annual total purchasing cost of the products 
C_{T} 
= 
Annual total transportation cost of the products 
Z 
= 
Annual total costs 
Here, the inventory model of the problem is developed.
The objective function: The objective function of the model is
to minimize the total cost of the joint replenishment problem given in
Eq. 1.
The terms in the right hand side of Eq. 1 are derived
as follows.
Transportation cost (C_{T}): The total transportation
cost is calculated based on Eq. 2, in which
is the transportation cost for each unit of the jth product and Q_{j}
is the number of order of jth product.
Purchasing cost under total discount (C_{P}): The purchasing
cost of the company for the jth product at the beginning of a period can
be calculated using the total discount policy. Let the total discount
policy be as below in which q_{ij} and C_{ij}; i = 1,2,…,n
are the discount points and the purchasing costs for each unit of the
jth product that corresponds to its lth discount break point, respectively.
In order to include the discount policy in the inventory model, we use
Eq. 4 to model the total discount policy. Figure
2 shows this policy.
Clearance cost (C_{C}): The clearance cost of the company
for the jth product at each period is jth, the order quantity is and
the number of order per year is Q_{j}. Hence, the total annual
clearance cost of the company will be:

Fig. 2: 
Total discount policy to purchase the jth product 
Fixed order cost (C_{A}): The fixed order cost of each
order per period is A and the number of orders per year is N. Hence and
the total annual fixed order cost in a disjoint ordering policy will be
NA. However, since we are using the joint replenishment policy, the fixed
order cost will change to:
Holding cost (C_{H}): The total available first warehouse
space is F and the space required for each unit of the jth product is
f_{j}. If the total required space for all products be more than
F, the company should rent another warehouse. So the holding cost during
a cycle will be:
By introducing the binary variables Y_{1}, Y_{2} Eq.
8 will change to Eq. 9 as below:
The first part of the capital cost occurs between the α% and (1α)%
payment times and is derived as:
The other part of the capital cost occurs between [t_{lc}, T]
and is calculated as:
So, the total holding cost during a year will be:
Also, according to (Fig. 1) we have:
Knowing that
Equation 12 can be written as Eq. 14
for a joint order case:
The constraints: Furthermore, since the total available budget
is TB and the purchasing cost that is calculated is:
the budget constraint will be:
For the sake of convenience in ordering, clearance, transportation and
some other activities, we assume an upper limit for the number of orders
per year. In other words:
Moreover, the quantities of the orders must be integer multiples of packets,
each containing more than one product. That is:
Knowing that:
We have:
Finally the model of the problem becomes:
A SOLUTION ALGORITHM
Since the model in Eq. 21 is mixedintegernonlinear in
nature, reaching an analytical solution (if any) to the problem is difficult
(Gen and Cheng, 1997). Hence, we need to employ a metaheuristic
search algorithm to solve it. Furthermore, as efficient treatment of mixed integer
nonlinear optimization is one of the most difficult problems in practical optimization
(ElSharkawi, 2008), we need to employ a metaheuristic
search algorithm to solve it.
New ways have been found to optimize problems for less than a century, but
nature has used various ways of optimization for millions of million years.
Recently scientists mimicked nature to solve different kinds of complex optimization
problems. Most of these problems are so complicated and time consuming that
we cannot use an exact algorithm to solve them. Thus, typically some nonprecise
algorithms are used to find a near optimum solution in a shorter period. We
call these algorithms metaheuristic (Dorigo and Stutzle, 2004).
Many researchers have successfully used metaheuristic methods to solve complicated
optimization problems in different fields of scientific and engineering disciplines.
Some of these metaheuristic algorithms are simulating annealing (Aarts
and Korst, 1989; Taleizadeh et al., 2008a),
threshold accepting (Dueck and Scheuer, 1990), Tabu search
(Joo and Bong, 1996), genetic algorithms (AlTabtabai
and Alex, 1999; Taleizadeh et al., 2008b),
neural networks (Gaiduk et al., 2002), ant colony
optimization (Dorigo and Stutzle, 2004), fuzzy simulation
(Taleizadeh et al., 2009), evolutionary algorithm
(Laumanns et al., 2002) and harmony search (Lee
and Geem, 2004; Geem et al., 2001; Taleizadeh
et al., 2008c).
Simulated Annealing (SA) is a local search method that finds its inspiration
in the physical annealing process studied in statistical mechanics (Aarts
and Korst, 1989). It is an effective and efficient method that produces
good suboptimal solutions and has diffused rapidly into many combinatorial optimization
areas (Kirkpatrick et al., 1994). An SA algorithm
repeats an iterative neighbor generation procedure and follows search directions
that improve the objective function value. While exploring solution space, the
SA method offers the possibility of accepting worse neighbor solutions in a
controlled manner in order to escape from local minima. More precisely, in each
iteration, for a current solution x characterized by an objective function value
f(x’), a neighbor x’ is selected from the neighborhood of x denoted N(x) and
defined as the set of all its immediate neighbors. For each move, the objective
difference Δ = f(x’)f(x) is evaluated. For minimization problems, x’ replaces
x whenever Δ≤0. Otherwise, x’ could also be accepted with a probability
The acceptance probability is compared to a number RN ε Uni(0,1)
generated randomly and x’ is accepted whenever P>RN.
The factors that influence the acceptance probability are the degree
of objective function value degradation Δ (smaller degradations induce
greater acceptance probabilities) and the parameter T called temperature
(higher values of T give higher acceptance probability). The temperature
can be controlled by a cooling scheme specifying how it should be progressively
reduced to make the procedure more selective as the search progresses
to neighborhoods of good solutions. There exist theoretical schedules,
which require infinite computing time guaranteeing asymptotic convergence
toward the optimal solution. In practice, much simpler and finite computing
time schedules are preferred even if they do not guarantee an optimal
solution.
A typical finite time implementation of SA consists of decreasing the
temperature T in S steps, starting from an initial value T_{0}
and using an attenuation factor α, (0<α<1). The initial
temperature T_{0} is supposed to be high enough to allow acceptance
of any new neighbor proposed in the first step. In each step s, the procedure
generates a fixed number of neighbor solutions N_{sol} and evaluates
them using the current temperature value Ts = α^{s}T_{0}.
The whole process is commonly called cooling chain or also markov chain.
Adaptation of SA to an optimization problem consists in defining its specific
components: a neighbor generation procedure, objective function evaluation,
a method for assigning the initial temperature, a procedure to change
the temperature and a cooling scheme including stopping criteria. These
adaptation steps for our new adaptations of SA for the inventory model
of Eq. 21 are described here.
A neighbor generation procedure: The neighbor generation is an
important component of SA. It has to be designed to allow easy neighbor
generation and fast calculation of Δ and guarantee accessibility
for the entire solution space.
In this study, the initial solutions are generated in two groups. In
the first group, the procedure randomly chooses the initial solution from
the possible solution space. For the other groups of initial solutions
the procedure employs some solutions of a Genetic algorithm described
here.
Evaluating the objective function: After creating each solution,
the objective function should be evaluated to compare the jth and the
ith solutions. In a maximization problem, if the objective function of
the new solution (j) is bigger than the previous solution (i), then (j)
will be accepted. Otherwise, by generating a random number the better
solution will be selected.
Assigning the initial temperature: Temperature is one of the important
parameters of the SA algorithm that affects the acceptance or nonacceptance
of the objective function changes. The initial temperature should be chosen
such that at the first level a big number of inappropriate solutions are
accepted. In this way, more changes are possible and hence more solutions
are explored. Furthermore, the time required for the process repetition
along with the annealing process depends on the initial temperature. In
this research, the large values of 1000, 2000 and 5000 are chosen for
initial temperatures.
Changing the temperature: One of the primary aspects of the annealing
process in a SA algorithm is the range of temperature change. The temperature
plays an important role in the possibility of selecting a bad solution.
On one hand, when the temperature assumes a high value, a big number of
bad solutions are accepted, leading to selecting a local optimum point.
On the other hand, when the temperature is low, the probability of the
solution to be a local optimum is high. In this study, we change the temperature
of the SA algorithm based on a geometric function given in Eq.
22 with α = 0.9, 0.95 and 0.99.
Ts = αT_{s1}; s =
1,2,… and 0<α<1 
(22) 
Cooling scheme and stopping criterion: In a specific temperature
of a SA algorithm, it is necessary to analyze the equilibrium state after
a couple of renitence to see if the annealing process will be continued
in that temperature or it will be stopped and moved to the next temperature.
An epoch is the number changes that should occur in a specific temperature.
Different types of stopping criteria have been presented for different
SA algorithms in the literature. Some of them are:
• 
A lower bound on the final temperature 
• 
An upper bound on the total number of states 
• 
Reaching a solidification point based on the value of a function
evaluated on the value of the objective function and the temperature 
• 
An upper bound on the total number of exchanges accepted along the
annealing process 
• 
An upper bound on the total number of refused changes 
In this research, we select the first criterion and stop the algorithm
after it reaches the final temperature T_{F}. Furthermore, we
use 100, 200 and 500 for different values of N(t).
In short, the steps involved in the SA algorithm used in this research
are:
1 
Choosing an initial solution i from the group of the
feasible solutions S 
2 
Choosing the initial temperature T_{0}>0 
3 
Selecting the number of iterations N(t) at each temperature 
4 
Selecting the final temperature T_{F} 
5 
Determining the process of the temperature reduction until it reaches
T_{F} 
6 
Setting the temperature exchange counter n to zero for each temperature.
(Balancing process) 
7 
Creating the j solution at the neighborhood of the i solution 
8 
Evaluating the objective function f = Z at any temperature. 
9 
Calculating Δ = f(j)f(i) 
10 
Accepting the solution j, if Δ<0. Otherwise,
generating a random number RN~Uniform (0,1). If
then selecting the j solution

11 
Setting n = n+1. If n is equal to N(t) then go to 12. Otherwise,
go to 7 
12 
Reducing the temperature. If it reaches T_{F} then stop.
Otherwise, go to 6 
NUMERICAL EXAMPLE
Consider a multiproduct inventory control problem with twenty products
and general data given in Table 1. The total available
warehouse space for first store, the total budget and the maximum number
of order are F = 500, TB = 170,000,000 and N_{T} = 6, respectively.
Also, L = 2.5 months, A = 2000, (t_{0}) = 0, h_{l}= 40,
=
50 and (t_{lc}) = 1.5 months. Table 2 shows
different values of the SA parameters used to obtain the solution. In
this research, all of the possible combinations of the SA parameters are
employed and using the min criterion the best combination of the parameters
has been selected. Table 3 shows the best combination
of the SA algorithm and the best result is shown in Table
4. Furthermore, the convergence path of the best result of the objective
function is shown in Fig. 3.
Table 1: 
General data of numerical example 

Table 2: 
The parameters of the SA algorithm 

Table 3: 
The best combination of the SA parameters 


Fig. 3: 
Convergence graph of SA algorithm 
CONCLUSION AND RECOMMENDATIONS FOR FUTURE RESEARCH
In this study, a joint replenishment multiproduct multi constraint inventory
model with two warehouse space with different holding cost based on EOQ
model to purchase high price raw materials was investigated. A mathematical
modeling by total discount, clearance, fixed order, holding, shortage
and transportation costs was introduced and shown to be integernonlinear
programming problems. Then, a metaheuristic algorithm (Simulated Annealing)
has been proposed to solve the integer nonlinear problem.
• 
Some recommendations for future works follow: 
• 
Shortage can be included in the model 
• 
Stochastic demands or lead time can be taken into account 
• 
It can be considered that lead time may exceed than cycle length 
• 
Deterioration rate for inventory in the warehouse can be considered 