- Email: [email protected]
A new heuristic method for distribution networks considering service level constraint and coverage radius
Expert Systems with Applications 36 (2009) 5620–5629
Contents lists available at ScienceDirect
Expert Systems with Applications journal homepage: www.elsevier.com/locate/eswa
A new heuristic method for distribution networks considering service level constraint and coverage radius V.R. Ghezavati *, M.S. Jabal-Ameli, A. Makui Department of Industrial Engineering, Iran University of Science and Technology, P.C. 16844, Narmak, Tehran, Iran
a r t i c l e
i n f o
Keywords: Facility location Distribution network Genetic algorithm Uncertainty modeling Inventory Heuristic .
a b s t r a c t This paper presents a new mathematical model for designing distribution networks in a supply chain system considering service level constraint optimizing strategic decisions (location), tactical decisions (inventory), and assigning decisions. In real-world cases, demand, traveling time or any parameters in classical models may change over the period of time. So, considering uncertainty yields more flexibility for the results and the proposed model. In our study, environmental uncertainty is described by discrete scenarios. In this model, we have service level constraint in order to prevent inventory lost in distribution centers (DCs). Also, we assume that customer’s demand is stochastic with Poisson distribution function and DCs have coverage radius constraints thus any DC cannot service all the customers. In this model, location of DCs is selected and optimized and the best flow of products from supplier to DCs also from DCs to customers is determined. In this way, the customers’ demand should be satisfied at least service level. To solve this nonlinear integer programming model we first present a new and robust solution based on a genetic search framework and then based on genetic algorithm results and some optimizer rules we propose a new heuristic method. Finally, some numerical examples are presented to illustrate the effectiveness of the proposed algorithms. 2008 Elsevier Ltd. All rights reserved.
1. Introduction One of the most important aspects of designing distribution networks is decision making about where to locate new facilities such as retailers, distribution centers (DCs) or plants. These strategic decisions are crucial determinant of whether materials will flow efficiently through the distribution system. Because of the fact that each parameter of this problem can change sensitively in the period of time – such as time, demands and distances – so, decisions related to the design can be made in effect. Recently, development of the distribution networks models under uncertainty has high priority for researchers (Snyder, 2006). Because inventory costs depend on the number of distribution centers, if these two decisions (strategic and tactical decisions) have made in a single model, we have cost saving significant in a distribution network. Snyder, Daskin, and Chung-Piaw (2007) developed a stochastic version of location model with risk pooling (SLMRP) that considers expected inventory costs with strategic costs into a single model. The SLMRP allows the parameters to be described by discrete scenarios.
Abbreviations: DC, distribution network; SLMRP, stochastic location model with risk pooling; GA, genetic algorithm; PDSDP, production–distribution system design problem. * Corresponding author. Tel.: +98 9126237336; fax: +98 217893177. E-mail address: [email protected] (V.R. Ghezavati). 0957-4174/$ - see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.eswa.2008.06.130
In this paper, we propose SLMRP considering service level constraint and coverage radius restrictions that handles parameters to be uncertain which described by discrete scenarios. The aim of this model is to select the best sites for distribution centers, assign customers to DCs according to coverage radius and determine inventory levels considering service level constraint. Also, inventory costs are computed applying new approach which is more realistic and practical than previous techniques. So, we introduced a nonlinear integer programming which minimizes total expected cost of the network and propose a new and robust genetic algorithm and a heuristic procedure to solve the model. 2. Literature review Many researchers have widely developed distribution network models in supply chain systems. More recently, Daskin, Coullard, and Shen (2002) studied inventory – location models that considered expected cost of inventory with cost of location and allocation concurrently. Aikens (1985); Francis, McGinnis, and White (1983) have studied previous location and demand allocation models. Some researches in uncertainty are studied in the area of multi product and multi echelon in supply chain by Tsiakis, Shah, and Pantelides (2001). In this model, objective is to find capacities, location of new facilities and the best flow of material in the network in order to minimize total expected costs.
5621
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Gabor and van Ommeren (2006) presented an approximation algorithm for a facility location problem with stochastic demand. They presented an expected value of a constraint that the probability that an arbitrary request lost was at most a. Chan and Chung (2005) applied combination of genetic algorithms with analytic hierarchy process to obtain the advantages of multi-criterion decision for solving distribution networks problem in a supply chain system. In their algorithm decision maker can consider weights for each criterion by using a pair wise comparison approach. Zhiqiang and Bostel (2007) studied a reverse supply chain system in which in their model there were three types of facilities. They considered both forward and reverse flows in the network and introduced a 0–1 mixed integer model. Amiri (2006) studied about aggregating production planning and distribution network design. He considered multi level capacity for each distribution center to achieve more utilization for capacity of DCs. He proposed a mixed integer model and developed a heuristic solution procedure. Mirchandani, Oudjit, and Wong (1985) proposed P-median model considering uncertainty which described by set of scenarios. They proposed the model as a deterministic p-median problem with jIj jSj customers instead of jIj, where I is number of customers and S is number of scenarios. MirHassani, Lucas, Mitra, Messina, and Poojari (2000) proposed a model in a distribution networks considering fixed resources. They modeled the problem as a twostages model where discrete variables are in first stage and continues are in the second stages. They used Benders decomposition algorithm to solve the model. Harpera, Shahania, Gallagherb, and Bowiec (2005) proposed location–allocation problem model for planning health services based on geographical parameters. In they model location of servers and customers who need services is determined. Lieckens and Vandaele (2007) designed the reverse logistic of a network by combining queuing theories with the classical models. They believed that by this combination they can get two advantages. The first is that the model can been proposed dynamically and the second is getting more uncertainty in reverse logistic modeling. They modeled this problem as a mixed integer nonlinear program (MINLP-model). Zhou, Min, and Gen (2002) introduced a location–allocation model based on naive balanced star spanning forest algorithm. Their purpose is to keep the best balanced of transportation costs and customer service. They developed a genetic algorithm to solve their model. Chan, Carter, and Burnes (2001) formulated location–routing problem considering multi-depot and multi vehicle. They used uncertainty of demand in the model where stochastic demands were generated by a queuing network at each service region. Dasci and Verter (2001) proposed the model in production–distribution system design problem (PDSDP) in a supply chain system in which in their new structure to model, costs and customers’ demand has continues distribution function. They show that discrete and continuous modeling approaches complement each other. Zuo-Jun, Shen, and Qi (2007) developed a model in a supply chain system where customers’ demand has uncertainty so they considered safety stock to respond fluctuation of demand. In their model decision maker should determine the number and locations of the distribution centers. The objective function of their model was minimizing total expected costs consist of fixed costs for locations, inventory costs (ordering and maintaining costs) and transportation costs. They formulated the problem as a nonlinear integer programming model and introduced Lagrangian relaxation algorithm to solve the model. Conde (2007) developed a robust location–allocation model in a supply chain. In their model the customers’ demand was uncertain and it only had an interval which was estimated. The problem was selecting locations for new facilities in order to satisfy a specific portion of demands. The objective function was minimax regret, i.e. the locations that minimize the worst-case loss in the objective function.
3. Model formulation In this section, we describe in more detail the designing distribution networks under uncertainty which we will present. There is a set of customers (N), at which demands are generated, and a set of potential sites for DCs (M), where facilities may be opened. We assume that the demands for each customer i are generated according to a Poisson process, independent of the demand at other customers in (N). At each open facility an inventory is kept such that an arriving request finds a zero inventory (and is lost), with probability at most a. We then say that (1 a) is the fill rate of the system. In this model, we have fixed location cost for opening and operating a DC, and a fixed cost to order at any DC and a holding cost for inventory, respectively. We assume that the DCs follow a (Q, R) inventory policy and total inventory cost consists of ordering cost, holding cost, inventory lost cost and not safety stock cost. 3.1. Order quantity and service level constraint In the previous location–inventory models the required order, for each DC’s, is the sum of two parts. First part is computed based on the sum of mean demand, assigned to the DC (based on EOQ model):
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2Sh ðY ij li Þ
ð1Þ
i
and the second part is the safety stock:
hZ a
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi X LT Y ij r2i
ð2Þ
i
A new approach in this article can be proposed without dividing the order to two parts with adding service level constraint. By combining the two parts, let Yj be the total amount of the order for DC j. So, if the probability of Pr(Xj 6 Yj) is at least (1 a), we can get three advantages where Xj is random variable denoting demand distribution function for DC j. (a) The required inventory can be computed automatically based on proposed constraint, (b) we can consider service level constraint in order to prevent inventory lost and (c), in the previous opinion for safety stock inventory all types of costs such as ordering cost and transportation cost did not consider in the model and just holding cost was computed, although in practice we have to pay all costs for safety stock and our idea consider all of them simultaneously in the model. For this purpose, this service level constraint must be considered. Indeed, we do not need to compute the safety stock.
0
an arbitrary request arriving
B P@ at DC j with inventory level Y j is not lost
1 C AP1a
ð3Þ
According to this assumption, order quantity is greater than customers’ mean demand and model determines optimum quantity so that customers’ demand satisfies at least (1 a)%, so there is no necessity to hold safety stock. (If we want to hold safety stock, then we should order exactly equal customers’ demand mean and by using that safety stock we can respond to demand fluctuation. but in our study we order more than customers’ demand mean.) The service level constraint considered for this is
MinðY j jPrðX j 6 Y j Þ P 1 aÞ8j 2 M
ð4Þ
Xj: random variable denoting demand distribution function for DC j. Yj: decision variable denoting total product which DC j should supply to satisfy customers’ demand with confidence level (1 a)%.
5622
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
(In this case, Xj has Poisson distribution function.) Because customers’ demands have Poisson distribution, total product that any DC needs to order has Poisson distribution. (The summation of some Poisson distribution follows Poisson distribution.) For example, if a special DC satisfies the demand of customers 1, 2 and 3 the distribution of product that this DC needs has Poisson distribution with the rate k = k1 + k2 + k3. We have another purpose from mentioned constraint and it refers to consider ‘‘Min” for the constraint. If we use the constraint as Pr(Xj 6 Yj) P 1 a, the model in the optimal solution can satisfy the constraint with any fill rate greater than (1 a)%. So, fill rate of all DCs may be different in the optimal solution, and this can lead to unsmooth for the network. But if we use the constraint as (Min(YjjPr(Xj 6 Yj) P 1 a)) the fill rate for DCs will be the same and performance of the network can be more smooth.
Decision variables 8 < 1 If customer’s demand i satisfies X ijs ¼ by DC at site j in scenario s : 0 Otherwise Yjs: variable denoting the total amount of product from supplier to DC j in scenario s probability of occurring inventory lost in DC at site j in bjs: scenario s
3.2. Coverage radius
3.4. Model formulation
In practice, based on some restrictions such as long distances, it may be occurred that each DC cannot cover all customers. In other words, all distribution centers cannot service all of the customers because each of the DCs has special critical coverage radius and if a customer does not be in this critical coverage, the DC cannot service that customer. So, we consider this parameter in the model which yields more flexibility for the proposed model. This parameter is named coverage radius and we show it by Zij. Hwang (2002) in his article used similarity from this notice. Fig. 1 illustrates this concept by deleting some possible links between DCs and customers based on coverage radius restrictions.
The cost in a specific scenario s is computed like that in Snyder et al. (2007) and is given by:
3.3. Notation We use Sets N: M: S:
the following notation index set of customer zones index set of potential distribution sites set of possible scenarios
Parameters Cij: cost of supplying one unit of demand to customer zone i from DC at site j C js : cost of supplying one unit of demand to DC at site j from supplier in scenario s fixed cost for opening and operating DC at site j Fj: mean of demand per time unit for customer i in scenario s kis: holding cost per unit in DC at site j hj: fixed cost per order placed to the supplier by DC at site j Sj: l: cost for unit lost 1 If a DC at site j can cover customer i Z ij ¼ 0 Otherwise; ps:
probability that scenarios s occurs.
Fig. 1. The figure of supply chain network with coverage radius.
We consider location variable (U) to be independent of scenario because it must be made before it is known which scenario will be occurred. Also, allocation variables should be dependent to the scenario, so we indicate the X, Y and b variables with index of scenario. Our goal is to choose locations of DCs in order to minimize the total expected cost of the network.
Min Z ¼
X
XX
fj U j þ
j2M
j2M
X
s2S
ps C js Y js þ
!
!
X
C ij kis X ijs þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2hj Sj Y js
i22N
l
ð5Þ
X ijs Z ij ¼ 1 8i 2 N; 8s 2 S
ð6Þ
þbjs
kis X ijs
i2N
Subject to:
X j
X ijs 6 Z ij U j
8i 2 N; 8j 2 M; 8s 2 S
Y js X eð Min Y js r¼0
P
k X Z i is ijs ij
Þ P k X Z r i is ijs ij r!
ð7Þ ! P1a
8j 2 M; 8s 2 S bjs þ
Y js X eð
P
k X Z i is ijs ij
r¼0
ð8Þ Þ P k X Z r i is ijs ij ¼ 1 8j 2 M; 8s 2 S r!
ð9Þ
U j 2 f0; 1g 8j 2 M
ð10Þ
X ijs 2 f0; 1g 8j 2 M; 8i 2 N; 8s 2 S
ð11Þ
Yjs P 0 8j 2 M
ð12Þ
0 6 bjs 6 1 8j 2 M; 8s 2 S
ð13Þ
The proposed model minimizes total expected cost made of: the costs associated with opening and operating the DCs, the costs of shipments from the supplier to the DCs, the costs to serve demand of customers from the DCs, total costs of ordering and holding inventory the expected costs of inventory lost in any opened DC (the expected cost = probability of inventory lost total expected demand unit inventory lost cost). Constraint set (6) guarantees that all customers allocated to exactly one DC in each scenario with considering coverage radius. Constraint set (7) says that a customer can be allocated to a DC if that customer is in the coverage radius of DC and DC is opened. Set constraint (8) forces that service level for DCs in each scenario would be at least (1 a)%. (as we said before, because customers’ demand is Poisson distribution, so total product that DC needs has Poisson distribution and sum of the demand that DC satisfies should be less than the inventory that DC receives (Yjs) and this probability should be at least (1 a)% and this constraint
5623
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
return minimum feasible value of Yjs variable). Constraint set (9) presents the probability of occurring loss of inventory in each scenario for any opened DC
ðbjs ¼ 1 PðX j 6 Y js Þ 8j 2 M; 8s 2 SÞ Constraint sets (10)–(13) determine type of variables. 4. Solution approach
No. of DC which Cover each Customer
No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 3 2
2 2 4 6
2 2 1 7
… … … …
m m 9 2
Scenario n
4
5
2
…
6
. . .
Fig. 2. Sample of chromosome structure.
It is known that uncapacitated location problems belong to the class of NP-hard problems that cannot be solve in reasonable computational time. In this paper, we propose a nonlinear and stochastic programming model for the uncapacitated network design problem which is extremely hard to optimally solve in large sizes. Thus, we use an efficient meta-heuristic method based on genetic algorithm (GA). Also, we design a heuristic algorithm based on optimizer rules which can give us sub-optimal solutions, to measure effectiveness of the algorithm. For these reasons we use this genetic algorithm similar to Zhou and Liu (2003). 4.1. Genetic algorithm GA is a stochastic search and heuristic optimization technique, which has been widely adopted by many researchers to solve various problems. This algorithm was first developed by Holland (1975). It mimics the mechanism of genetic evolution in the biological nature and consists of a population of chromosomes (strings or individuals) that are composed of genes. These genes represent a number of values, called alleles. Each chromosome (genotype) represents one potential solution (phenotype). The process of genetic operators (i.e., crossover and mutation) is carried out in the pool; after that, an evolution is completed by creating new chromosomes (offspring). This offspring is expected to be stronger than the parents, but this may not always be true. GA does not rely on analytical properties of the function to be optimized (Goldberg (1989)). In short, GA has two major processes: (1) GA is iteratively and randomly generating new solutions; and (2) these solutions are checked for the optimality according to predefined fitness functions. This becomes the most powerful principle of GA. It makes them widely suitable for finding an optimal solution in many complex problems, such as traveling salesman problem (TSP) and any forms of scheduling problems.
4.1.3. Decoding In this section, we compute the other variables (Yjs, Uj, bjs) of the model based on definition of the feasible chromosome. To compute Yjs and bjs we do the procedure described in Fig. 3: 4.1.4. Compute the fitness of all chromosomes The rank-based evaluation function is defined as the objective function for chromosome k = 1, 2, . . . , Pop_Size. 4.1.5. Selection strategy By generating a random real number, r, from the interval [0, 1], chromosome Vk is selected as a parent provided that r < P; where the parameter P is the probability of crossover or mutation operator. 4.1.6. Genetic algorithm operators As we know, one of the important parts in GA is generating new chromosomes from current chromosomes (called parents). This process can be done by genetic operators. Operators of this algorithm are crossover and mutation. 4.1.6.1. Crossover operator. We combine the chromosomes Vk, k = 1, 2, . . . , Pop_Size by crossover operation. In order to determine the parents for crossover operation, we repeat the following process from k = 1 to pop_size: generating a random real number r from the interval [0, 1], the chromosome Vk will be selected as a
4.1.1. Representation of chromosome In our GA-based approach, each chromosome or bit string (i.e., an example solution) consists of a matrix: Assign matrix. The length of this matrix is equal to: number of scenarios by the number of customers in the network. The Assign matrix represents the assignments of the customers to DCs in any scenarios. If customer i is assigned to DC k in scenario s, its entry in the Assign matrix has value k in row s and column i. Fig. 2 shows a sample of the chromosome structure. 4.1.2. Initialize chromosome population We initialize pop_size chromosomes X k ¼ ðxki;s Þ ¼ ðxk1;1 ; xk2;1 ; . . . ; k xn;1 ; . . . ; xk1;s ; xk2;s ; . . . ; xkn;s Þ k = 1, 2, . . . , Pop_Size from the feasible region {(X)jgis(xis) 6 0, xis = 1, 2, . . . , n} randomly. (xki;s Denote number of DC which must cover customer i in scenario s in chromosome k.) First for each customer, number of DCs which can cover that customer is determined. Then a random number from those DCs is generated and customer in scenario s is allocated to that DC. (In other words Xijs = 1.)
Fig. 3. Pseudo code of the decoding procedure.
5624
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
parent provided that r < Pc, where the parameter Pc is the probability of crossover. Then we group the selected parents V 01 ; V 02 ; V 03 ; . . . to the pairs ðV 01 ; V 02 Þ; ðV 03 ; V 04 Þ; . . ., without loss of generality; let us describe the crossover operator on each pair by ðV 01 ; V 02 Þ. First, for all customers in all scenarios, we generate a random number k from the open interval (0, 1), and then the crossover operator on V0 and V 02 will produce one child X as follows in Fig. 4. In fact, in our chromosome each customer in each scenario is assigned to a special DC. For this reason, in crossover process we have two parents that each customer may be assigned to different DCs in each parent. To generate offspring, we assume that each customer can be assigned to one of them randomly and this is occurs with probability 0.5. In other words, in offspring each customer is assigned to one of the DCs which covered that customer in the parents randomly. We show sample of this process in Fig. 5. 4.1.6.2. Mutation operator. We update the chromosomes Vk, k = 1, 2, . . . , Pop_Size by mutation operation. Similar to the proves of selecting parents for crossover operation, we repeat the following steps from k = 1 to pop_size: generating a random real number r from interval [0, 1], then chromosome Vk will be selected as a parent provided that r < Pm; where the parameter Pm is the probability of mutation. For each selected parent:
V k ¼ ðxk1;1 ; xk2;1 ; . . . ; xkn;1 ; . . . ; xk1;s ; xk2;s ; . . . ; xkn;s Þ we mutate it with one of the following ways. Because we have three purposes in mutation (reducing number of opened DCs, substitution an opened DC with a closed DC to check all conditions in the solution, reassigning customers’ allocation). So, we considered three types for mutation operator and any of them may be occurred with probability 1/3 for each chromosome selected for mutation. (a) Select randomly from opened DCs (Uj = 1).If selected DC can be closed, then (a DC can be closed if all customers assigned
to it, can be assigned to another opened DC):For each customer in each scenario, allocated to that old DC, select randomly from new DC which has two below conditions and assign the customer to the new DC: – New DC j is opened in this chromosome. – New DC j can cover the customer. By this way, we can reduce the number of opened DC. We show sample of this process in Fig. 6. (b) Select randomly from opened DCs and name it j, then select randomly from closed DCs and name it j0 . With below procedure we close DC j and instead of it open DC j0 .First, for customer allocated to DC j in any scenario, if the customer can be allocated to DC j0 , then allocate it in the same scenario to DC j0 .Otherwise, allocate it to another opened DC which can cover it. Sample of mutation type 3 has been illustrated in Fig. 8.By this way, we can replace an opened DC with a closed DC for more searching in feasibility space. Fig. 7 has shown an example of this process. (c) Select a customer in a chromosome randomly. Reassign its demands to another opened DC which can cover it. Sample of mutation type 3 has been illustrated in Fig. 8.
4.2. Heuristic procedure Our purpose from solving mathematical modeling is to study the behavior of solution procedure and finally use it in practice and real world. For these reasons, we present a heuristic method based on some rules obtained from optimal model and our proposed GA that lead to optimality for mentioned model in the article. Because the model of the article cannot be solved with optimizer algorithms, we use from this heuristic algorithm to get sub optimal solutions. So, by using this heuristic and genetic algorithm we can use this model in real world and practice (Miranda & Garrido, 2004). 4.2.1. The rules Note: according to coverage radius which is one of our assumptions in modeling; our problem is a set-covering problem.
∀i ∈ N , ∀s ∈ S : ( i is index of customers) If <0.5 then
(a) As there is no capacity restriction on the distribution centers, at optimality each customer will be supplied entirely from the most cost-effective located facility. This may be fortuitous because many firms require single sourcing of customers to gain efficiencies and ensure good customer service. In other words, each customer should be assigned to the nearest opened DC (Love, Morris, & Weslowsky, 1988a).
X i , s = V1′(i, s )
Else
X i , s = V2′(i, s )
End if
Fig. 4. Pseudo code of the crossover operator.
Parent 1 No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
4 4 1 2
Parent 2
5 5 1 2
No. of Gens No. of Customers Scenario1 Scenario 2
1 1 3 1
2 2 1 1
4 4 1 2
5 5 1 2
Crossover Offspring No. of Gens No. of Customers Scenario1 Scenario 2
1 1 3 1
2 2 2 1
3 3 2 3
Fig. 5. Sample of crossover operator.
3 3 3 3
4 4 1 1
5 5 3 3
5625
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Parent No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 2 2
2 2 2 3
1 1 3 2
2 2 2 2
1 1 1 2
2 2 2 1
3 3 2 3
4 4 2 2
5 5 3 2
Fig. 6. Sample of mutation type 1.
Parent No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
Mutation type 2
No. of Gens No. of Customers Scenario1 Scenario 2
3 3 2 3
4 4 3 2
5 5 3 2
Fig. 7. Sample of mutation type 2.
Parent No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
No. of Gens No. of Customers Scenario1 Scenario 2
3 3 1 1
4 4 1 2
5 5 1 2
Mutation type 3 Fig. 8. Sample of mutation type 3.
(b) Each customer must be assigned to the same DC in each scenario and not different. (According to the first rule). For example, if customer 1 is assigned to DC 4 in scenario 1, this customer must be assigned to DC 4 in any scenario too because DC 4 is the nearest opened DC to customer 1. (c) Minimum number of opened DC is obtained in the way that we can assign all customers to DCs considering coverage radius. Base on Zij matrix we open minimum number of DCs that each customer can be covered by at least one of them. This is the result of set-covering problem in Love, Morris, and Wesolowsky (1988b). (d) Because set-up costs and difficulties of installing DCs is too higher than transportation and inventory costs and also regarding to unlimited capacity for DCs, any opened DC can cover any number of customers which are in its coverage radius despite their demands so, in these problems (Set covering) it’s suitable to work with the minimum number of DCs that leads to less costs and cost saving for the network. (In uncapacitated set covering location problems number of opened DCs is less than capacitated set covering location problems. In basic model for set-covering problem the aim of modeling is to minimize number of opened DCs (Love et al., 1988b) so, we follow this objective too. This result is fit with logic of designing distribution networks. Because, if n DCs are opened, capacity of network to response customers’ demand is 1. In this situation if n + 1 DCs are opened, the capacity of network to response customers is 1 too. In the second case, we just increase strategic costs and saving other costs is little. Our GA confirms this result too. In Table 1, we solve large size problems by GA and as we can see, number of opened DCs in all of them is minimum. 4.2.2. Heuristic stages We divide this method into two stages. In the first stage, we find which DCs should be opened and in the second stage allocation
Table 1 Solving large size problem with GA # of Customers
# of DCs
# of Scenarios
Number of opened DCs
Interaction
70 70 70 75 75 75 80 80 80 85 85 85 90 90 90 95 95 95 100 100 *100
75 75 75 50 50 50 50 50 50 60 60 60 55 55 55 65 65 65 60 60 60
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
12–22–41 12–45 25–37–39 23–24–40 23–30–33 6–30–41 6–15–44 29–35–45 2–9 39–46–59 17–52–54 6–14–27 11–34–52 14–16–27 13–15–19 3–11–54 30–33–62 21–27–40 32–48–59 20–36–28 9–17–25
250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250
decisions for customers and quantity of ordering for DCs considering service level are made.
4.2.2.1. First stage. In this stage, we should find minimum number of DCs that must be opened so that all customers can be covered. For this purpose, we check opening 1 DC, then 2 DCs, 3 DCs and etc. until we find feasible solutions to assign all customers. In fact, in this section we solve set-covering problem which is mentioned in rule c. To check whether 1 DC can cover all customers, we use from below equation:
5626
Sj ¼
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
X
Z ij
j ¼ 1; 2; . . . ; n
ð14Þ
i
Above equation computes number of customers that DC j can cover. If Sk = m (m number of customers), we will finish our computations in this stage and open DC k. Sj – m for j = 1, . . . , n means that by opening just 1 DC, we cannot cover all customers, so we should find 2 DCs j and j0 that Sj;j0 ¼ m and by opening DC j and j0 we can cover all customers
0 1 XB C 0 0 Sj;j0 ¼ @Z ij þ Z ij0 Z ij Z ij0 A j; j ¼ 1; . . . ; n; j–j |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} i Ai 8 < 1 If customer i can cover 0 Ai ¼ by at least DC j or DC j : : 0 Otherwise:
ð15Þ
ð16Þ
If $j1, j2:Sj1,j2 = m, our computations in this stage are finished. If "j1, j2:Sj1,j2 – m, we should check opening more number of DCs. In this situation we present general equation as below:
1 jk P P ðZ ij1 Z ij2 Þ þ þ C B ðZ ij Þ C B j¼j1 j1
Sj1;j2;...;jk
j1
ð17Þ In general, if Sj1,j2,. . .,jk = m, we open DCs with number j1, j2, . . ., jk. For example, for k = 3:
0 Sj1;j2;j3 ¼
ðZ i;j1 þ Z i;j2 þ Z i;j3 Þ
1
C XB B ðZ i;j1 Z i;j2 þ Z i;j1 Z i;j3 þ C B C @ Z i;j2 Z i;j3 Þþ A i
ð18Þ
ðZ i;j1 Z i;j2 Z i;j3 Þ
Fig. 9. Pseudo code of computing variables.
If we want to solve the model with full enumeration for a small problem with 20 customers, 2 scenarios and 5 DCs total points that we should calculate objective function of them is We have a variable with 3 indices (Xijs) and Uj and they are (0, 1) so: Total points:
ð220 25 22 Þ 25 ¼ 220þ5þ2þ5 ¼ 232 As we can see total points are 232 that is too difficult and time consuming to enumerate all of them. These results are for a small problem. For a large size problem with 100 customers, 60 DCs and 2 scenarios total points is
When above equation for K = k is equal to m, we stop our computations in this stage. Our purpose from this stage is to find minimum number of DCs that we should open them to cover all customers and for this we first check opening 1 DC and if it is infeasible we open 2 DCs. . . At the end of this stage we may get multi solutions. (E.g., DC 1, 3 can cover all customers and DC 2, 3 can cover them too.) So we calculate next stage for all of solutions obtained from first stage.
In fact enumeration 2222 points is infeasible while our genetic algorithm get us sub optimal solution in under 3 min (problem signed with * in Table 1). So, computation time in our genetic algorithm is so little and this is one of the advantages our work that leads to validation for GA.
4.2.3. Second stage
5. Computational results
(Step-1)
We divide this section to two parts. In first part, we will present effectiveness and robustness of the genetic algorithm and in the second part we will compare genetic algorithm with heuristic method to measure optimality of genetic algorithm. In order to illustrate the effectiveness of the genetic algorithm, we give some numerical examples that are performed on a personal computer. All algorithms considered in this study were coded in Visual Basic 6 and run on a Pentium IV PC with 1.5 GHz CPU and 256MB RAM. For this purpose we solve three examples.
(Step-2)
We found the DCs which must be opened in the previous stage. In this stage, we make allocation decisions. Because DCs have unlimited capacity, so each customer must assign to the nearest opened DC in each scenario in order to have minimum transportation costs in the network. (Rules 1 and 2). For each opened DC quantity of ordering considering service level is computed as below and compute these variables:To compute Yjs we do the procedure described in Fig. 9.
ð2100 260 22 Þ 260 ¼ 2222
5.1. Measuring robustness of the GA By doing mentioned heuristic algorithm which is based on optimizer rules we can obtain solutions that are near optimal solution.
Consider a company who wants to locate new distribution centers. Assume that there are 20, 30 and 50 customers whose de-
5627
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629 Table 2 Comparison solution of example 1
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.9 0.8 0.75 0.65 1 0.95 0.8 0.75 0.75
0.33 0.30 0.27 0.23 0.18 0.33 0.28 0.23 0.18 0.22
300 300 300 300 300 300 300 300 300 300
2.3 2.3 2.5 2.3 2.3 2.3 2.3 2.3 2.3 2.3
53905.64 53851.17 53661.01 53921.09 54200.21 53789.8 53862.59 53615.76 53789.8 54029.47
851.02 817.82 817.42 745.09 806.45 860.42 785.41 740.82 860.42 806.99
0.005 0.004 0.001 0.006 0.011 0.003 0.005 0.000 0.003 0.008
Average CPU time: 92 s.
mands are Poisson distributed variable with mean kis. Also, there are 2 scenarios for kis. Suppose that a decision maker needs to select distribution centers from 8, 10 and 12 potential distribution centers to serve customers. The capacities are unlimited. In these examples service level is (1 a) = 99%. In Table 2, we compare solutions when different parameters are taken with the same generations as a stopping rule. It appears that all the minimal costs differ little from each other. In order to account for it, we present a parameter, called the percent error, i.e. (objective value – the best objective value)/the best objective value, where the best objective value is the least one of all the ten minimal costs obtained above. The last column named by ‘error’ in Table 2 is just this parameter. It follows from Table 2 that the percent error does not exceed 1.1% when different parameters are selected, which implies that the genetic algorithm is robust to the parameter settings and effective to solve model. Therefore, the variation of parameters in the algorithm has slight impact on the optimal objective, which implies that the algorithm designed in this paper is much robust. Solutions with different parameters of GA are compared on the basis of equivalent generations. In order to measure the differentia of each other, the percent error is proposed and showed in Table 3 as ‘error’. It follows from Table 3 that the percent error does not exceed 2.3% when different parameters are selected. Thus, the genetic algorithm is also robust to the parameter settings and effective to solve model. Similarly, we run the genetic algorithm for ten times with different parameters of GA on the basis of equivalent generations. In order to measure the differentia between these results, ‘error’, i.e. the percent error, is calculated and given in Table 4. From these computational results, we see that the maximal percent error does not exceed 2.3% when different parameters are chosen. Therefore, the genetic algorithm is also robust to the parameter settings and effective to solve model.
5.2. Comparing heuristic method and GA In this section, we compare solutions obtained from our GA with heuristic method and finally compute error. We solve seven examples and we describe one of them in detail and present final results for the others. Assume we have 20 customers, 6 distribution centers and 2 scenarios. We fallow the steps of the heuristic as below:
0
1 B1 B B B1 B B1 B B B1 B B1 B B1 B B B1 B B1 B B0 B Z ij ¼ B B0 B B1 B B0 B B B1 B B1 B B B1 B B1 B B1 B B @1 1
1 0 0 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1
1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 0 1
1 1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 0 1 0 1
0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
1 1 1C C C 1C C 1C C C 1C C 1C C 1C C C 1C C 1C C 1C C C 1C C 1C C 1C C C 0C C 1C C C 1C C 1C C 1C C C 0A 1
According to Zij matrix we check opening 1 DC. So for 6 DC we have S1 = 17, S2 = 15, S3 = 16, S4 = 14, S5 = 17 and S6 = 18. (m = No. of customers = 20).
Table 3 Comparison solution of example 2
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.95 0.75 0.65 0.7 0.9 0.95 0.8 0.8 0.75
0.33 0.30 0.23 0.23 0.23 0.30 0.28 0.25 0.28 0.28
300 300 300 300 300 300 300 300 300 300
3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10
60583.07 60664.29 61475.17 60974.89 61127.22 60400.95 60304.03 60073.22 60328.1 60701.14
1095 1056.3 998.8 1128.97 1154.79 1172.57 1067.19 1085.16 1093.51 1086.65
0.008 0.010 0.023 0.015 0.018 0.005 0.004 0.000 0.004 0.010
Average CPU time: 151 s.
5628
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Table 4 Comparison solution of example 3
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.95 0.8 0.7 0.95 0.9 0.8 0.85 0.95 0.85
0.33 0.30 0.27 0.25 0.28 0.30 0.30 0.30 0.30 0.28
300 300 300 300 300 300 300 300 300 300
3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10
96851.03 97252.06 98411.66 98591.62 96575.56 96558.15 97077.29 96825.47 96331.58 96754.13
1982.28 2029.1 2096.08 1898.52 1954.14 1938.79 1817.6 1952.29 1819.63 1959.11
0.005 0.010 0.022 0.023 0.003 0.002 0.008 0.005 0.000 0.004
Average CPU time: 270 s.
As we can see, none of DCs have Sj = m. So by opening 1 DC we cannot cover all customers. We must check opening 2 DCs. By using equation mentioned in the algorithm (17) for DC j and j0 we have S(1, S(2, S(3, S(4, S(5,
2) = 20 3) = 18 4) = 18 5) = 19 6) = 20
S(1, S(2, S(3, S(4,
3) = 19 S(1, 4) = 19 S(1, 5) = 19 S(1, 6) = 20 4) = 19 S(2, 5) = 20 S(2, 6) = 20 5) = 19 S(3, 6) = 19 6) = 19
As we see DCs (1, 2)–(1, 6)–(2, 5)–(2, 6) and (5, 6) can cover all customers. So we compute stage two for all of them and select the best as sub-optimal solution. 1, 1, 2, 2, 5,
2 ? total 6 ? total 5 ? total 6 ? total 6 ? total
cost = 70,576 cost = 101,558 cost = 72,556 cost = 82,976 cost = 100,633
Total cost is made of fixed cost, transportation cost, inventory cost and loss of inventory cost. As we see the best objective is related to open DCs 1 and 2. We solve this problem with our GA. Total cost of this solution is = 72,072. We compute percentage error as below:
Error is ¼
72072 70576 ¼ 0:021 ¼ 2:1% 70576
Final results for the other problems are presented in Table 5. In Table 5, we compare solutions in GA and heuristic method when different problems are taken. As we see behavior of GA is similar to heuristic method. In other words, difference between them may be appeared just in allocation while location decisions are the same. In order to count error we use from this ratio: (objective function in GA – objective function in heuristic)/objective function in heuristic. It follows from Table 5 that the percent error
does not exceed 2.3%, which implies that the genetic algorithm is effective to solve model. So, solutions obtained from GA are near sub-optimal solutions obtained from heuristic method and we can use them in practice.
6. Conclusion and future directions In this paper, we introduced notation of distribution network in a supply chain system while we have service level constraint and customer’s demand is stochastic with Poisson distribution function. In the proposed model, we assumed that DCs have coverage radius restriction. We formulated the problem as a nonlinear integer program. The advantages of the proposed model are as follows: considering coverage radius which yields more flexibility for the model, considering service level constraint that by using it the required inventory can be computed automatically, it prevents inventory lost and all types of cost for inventory are consider in the model. To solve this nonlinear integer programming model we first introduced a new and robust genetic algorithm solution and then based on genetic algorithm results and some optimizer rules we presented a new heuristic method. Our computational results showed the algorithm to be quite effective. For future research we suggest three directions: (a) We can study this model while Zij is probability and this parameter does not have deterministic value. For example, it has a distribution probability function and based on it, Zij gets value 0 or 1. (b) Setting inventory capacity constraint is another development for this model. (c) Applying the other ordering methods to calculate inventory costs in designing distribution networks. (d) Aggregating this model with other assumption such as routing and considering time window can be an interesting development of the proposed model.
Table 5 Comparing heuristic with GA solution No. of customers
No. of DCs
No. of scenarios
Opened DCs in GA
Opened DCs in heuristic
Objective function in GA
Objective function in heuristic
Error
8 10 12 15 20 25 30
5 5 6 6 6 8 12
2 2 2 2 2 2 2
1 1 4 3.6 1.2 2.7 2.9
1 1 4 3.6 1.2 2.7 2.9
26745 29105.9 28566 43210.3 72072 78511 68375.12
26745 29105.9 28566 43210.3 70576 77903 66819.5
0 0 0 0 0.021 0.008 0.023
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
References Aikens, C. H. (1985). Facility location models for distribution planning. European Journal of Operational Research, 22(3), 263–279. Amiri, Ali (2006). Designing a distribution network in a supply chain system: Formulation and efficient solution procedure. European Journal of Operational Research, 171(2), 567–576. Chan, Y., Carter, W. B., & Burnes, M. D. (2001). A multiple-depot, multiple-vehicle, location-routing problem with stochastically processed demands. Computers & Operations Research, 28(8), 803–826. Conde, E. (2007). Minmax regret location–allocation problem on a network under uncertainty. European Journal of Operational Research, 179(3), 1025– 1039. Dasci, A., & Verter, V. (2001). A continuous model for production–distribution system design. European Journal of Operational Research, 129(2), 287– 298. Daskin, M. S., Coullard, C. R., & Shen, Z.-J. M. (2002). An inventory-location model: Formulation, solution algorithm and computational results. Annals of Operations Research, 110, 83–106. Chan, Felix, T. S., & Chung, S. H. (2005). Multicriterion genetic optimization for due date assigned distribution network problems. Decision Support Systems, 39(4), 661–675. Francis, R. L., McGinnis, L. F., & White, J. A. (1983). Locational analysis. European Journal of Operations Research, 12(3), 220–252. Gabor, A. F., & van Ommeren, J. C. W. (2006). An approximation algorithm for a facility location problem with stochastic demands and inventories. Operations Research Letters, 34(3), 257–263. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning. MA: Addison-Wiley. Harpera, P. R., Shahania, A. K., Gallagherb, J. E., & Bowiec, C. (2005). Planning health services with explicit geographical considerations: A stochastic location– allocation approach. Omega, 33(2), 141–152. Holland, J. H. (1975). Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence in 1992 (2nd ed.). Cambridge: MIT Press.
5629
Hwang, H. S. (2002). Design supply-chain logistics system considering service level. Computers & Industrial Engineering, 43(1–2), 283–297. Lieckens, K., & Vandaele, N. (2007). Reverse logistics network design with stochastic lead times. Computers & Operations Research, 34(2), 395–416. Love, R. F., Morris, J. G., & Wesolowsky, G. O. (1988a). Facilities location: Models and methods (p. 186). New York: North Holland Publishing Co. Love, R. F., Morris, J. G., & Wesolowsky, G. O. (1988b). Facilities location: Models and methods (p. 174). New York: North Holland Publishing Co. Lu, Zhiqiang, & Bostel, Nathalie (2007). A facility location model for logistics systems including reverse flows: The case of remanufacturing activities. Computers & Operations Research, 34(2), 299–323. Miranda, Pablo A., & Garrido, Rodrigo A. (2004). Incorporating inventory control decisions into a strategic distribution network design model with stochastic demand. Transportation Research Part E, 40(3), 183–207. Mirchandani, P. B., Oudjit, A., & Wong, R. T. (1985). Multidimensional extensions and a nested dual approach for the m-median problem. European Journal of Operational Research, 21(1), 121–137. MirHassani, S. A., Lucas, C., Mitra, G., Messina, E., & Poojari, C. A. (2000). Computational solution of capacity planning models under uncertainty. Parallel Computing, 26(5), 511–538. Snyder, L. V. (2006). Facility location under uncertainty: A review. IIE Transactions, 38, 537–554. Snyder, L. V., Daskin, M. S., & Chung-Piaw, Teo (2007). The stochastic location model with risk pooling. European Journal of Operational Research, 179(3), 1221–1238. Tsiakis, P., Shah, N., & Pantelides, C. C. (2001). Design of multi-echelon supply chain networks under demand uncertainty. Industrial and Engineering Chemistry Research, 40, 3585–3604. Zhou, J., & Liu, B. (2003). New stochastic models for capacitated location–allocation problem. Computers & Industrial Engineering, 45(1), 111–125. Zhou, G., Min, H., & Gen, M. (2002). The balanced allocation of customers to multiple distribution centers in a supply chain network: A genetic algorithm approach. Computers and Industrial Engineering, 43(1–2), 251–261. Zuo-Jun, M., Shen & Qi, L. (2007). Incorporating inventory and routing costs in strategic location models. European Journal of Operational Research, 179(2), 372–389.
Contents lists available at ScienceDirect
Expert Systems with Applications journal homepage: www.elsevier.com/locate/eswa
A new heuristic method for distribution networks considering service level constraint and coverage radius V.R. Ghezavati *, M.S. Jabal-Ameli, A. Makui Department of Industrial Engineering, Iran University of Science and Technology, P.C. 16844, Narmak, Tehran, Iran
a r t i c l e
i n f o
Keywords: Facility location Distribution network Genetic algorithm Uncertainty modeling Inventory Heuristic .
a b s t r a c t This paper presents a new mathematical model for designing distribution networks in a supply chain system considering service level constraint optimizing strategic decisions (location), tactical decisions (inventory), and assigning decisions. In real-world cases, demand, traveling time or any parameters in classical models may change over the period of time. So, considering uncertainty yields more flexibility for the results and the proposed model. In our study, environmental uncertainty is described by discrete scenarios. In this model, we have service level constraint in order to prevent inventory lost in distribution centers (DCs). Also, we assume that customer’s demand is stochastic with Poisson distribution function and DCs have coverage radius constraints thus any DC cannot service all the customers. In this model, location of DCs is selected and optimized and the best flow of products from supplier to DCs also from DCs to customers is determined. In this way, the customers’ demand should be satisfied at least service level. To solve this nonlinear integer programming model we first present a new and robust solution based on a genetic search framework and then based on genetic algorithm results and some optimizer rules we propose a new heuristic method. Finally, some numerical examples are presented to illustrate the effectiveness of the proposed algorithms. 2008 Elsevier Ltd. All rights reserved.
1. Introduction One of the most important aspects of designing distribution networks is decision making about where to locate new facilities such as retailers, distribution centers (DCs) or plants. These strategic decisions are crucial determinant of whether materials will flow efficiently through the distribution system. Because of the fact that each parameter of this problem can change sensitively in the period of time – such as time, demands and distances – so, decisions related to the design can be made in effect. Recently, development of the distribution networks models under uncertainty has high priority for researchers (Snyder, 2006). Because inventory costs depend on the number of distribution centers, if these two decisions (strategic and tactical decisions) have made in a single model, we have cost saving significant in a distribution network. Snyder, Daskin, and Chung-Piaw (2007) developed a stochastic version of location model with risk pooling (SLMRP) that considers expected inventory costs with strategic costs into a single model. The SLMRP allows the parameters to be described by discrete scenarios.
Abbreviations: DC, distribution network; SLMRP, stochastic location model with risk pooling; GA, genetic algorithm; PDSDP, production–distribution system design problem. * Corresponding author. Tel.: +98 9126237336; fax: +98 217893177. E-mail address: [email protected] (V.R. Ghezavati). 0957-4174/$ - see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.eswa.2008.06.130
In this paper, we propose SLMRP considering service level constraint and coverage radius restrictions that handles parameters to be uncertain which described by discrete scenarios. The aim of this model is to select the best sites for distribution centers, assign customers to DCs according to coverage radius and determine inventory levels considering service level constraint. Also, inventory costs are computed applying new approach which is more realistic and practical than previous techniques. So, we introduced a nonlinear integer programming which minimizes total expected cost of the network and propose a new and robust genetic algorithm and a heuristic procedure to solve the model. 2. Literature review Many researchers have widely developed distribution network models in supply chain systems. More recently, Daskin, Coullard, and Shen (2002) studied inventory – location models that considered expected cost of inventory with cost of location and allocation concurrently. Aikens (1985); Francis, McGinnis, and White (1983) have studied previous location and demand allocation models. Some researches in uncertainty are studied in the area of multi product and multi echelon in supply chain by Tsiakis, Shah, and Pantelides (2001). In this model, objective is to find capacities, location of new facilities and the best flow of material in the network in order to minimize total expected costs.
5621
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Gabor and van Ommeren (2006) presented an approximation algorithm for a facility location problem with stochastic demand. They presented an expected value of a constraint that the probability that an arbitrary request lost was at most a. Chan and Chung (2005) applied combination of genetic algorithms with analytic hierarchy process to obtain the advantages of multi-criterion decision for solving distribution networks problem in a supply chain system. In their algorithm decision maker can consider weights for each criterion by using a pair wise comparison approach. Zhiqiang and Bostel (2007) studied a reverse supply chain system in which in their model there were three types of facilities. They considered both forward and reverse flows in the network and introduced a 0–1 mixed integer model. Amiri (2006) studied about aggregating production planning and distribution network design. He considered multi level capacity for each distribution center to achieve more utilization for capacity of DCs. He proposed a mixed integer model and developed a heuristic solution procedure. Mirchandani, Oudjit, and Wong (1985) proposed P-median model considering uncertainty which described by set of scenarios. They proposed the model as a deterministic p-median problem with jIj jSj customers instead of jIj, where I is number of customers and S is number of scenarios. MirHassani, Lucas, Mitra, Messina, and Poojari (2000) proposed a model in a distribution networks considering fixed resources. They modeled the problem as a twostages model where discrete variables are in first stage and continues are in the second stages. They used Benders decomposition algorithm to solve the model. Harpera, Shahania, Gallagherb, and Bowiec (2005) proposed location–allocation problem model for planning health services based on geographical parameters. In they model location of servers and customers who need services is determined. Lieckens and Vandaele (2007) designed the reverse logistic of a network by combining queuing theories with the classical models. They believed that by this combination they can get two advantages. The first is that the model can been proposed dynamically and the second is getting more uncertainty in reverse logistic modeling. They modeled this problem as a mixed integer nonlinear program (MINLP-model). Zhou, Min, and Gen (2002) introduced a location–allocation model based on naive balanced star spanning forest algorithm. Their purpose is to keep the best balanced of transportation costs and customer service. They developed a genetic algorithm to solve their model. Chan, Carter, and Burnes (2001) formulated location–routing problem considering multi-depot and multi vehicle. They used uncertainty of demand in the model where stochastic demands were generated by a queuing network at each service region. Dasci and Verter (2001) proposed the model in production–distribution system design problem (PDSDP) in a supply chain system in which in their new structure to model, costs and customers’ demand has continues distribution function. They show that discrete and continuous modeling approaches complement each other. Zuo-Jun, Shen, and Qi (2007) developed a model in a supply chain system where customers’ demand has uncertainty so they considered safety stock to respond fluctuation of demand. In their model decision maker should determine the number and locations of the distribution centers. The objective function of their model was minimizing total expected costs consist of fixed costs for locations, inventory costs (ordering and maintaining costs) and transportation costs. They formulated the problem as a nonlinear integer programming model and introduced Lagrangian relaxation algorithm to solve the model. Conde (2007) developed a robust location–allocation model in a supply chain. In their model the customers’ demand was uncertain and it only had an interval which was estimated. The problem was selecting locations for new facilities in order to satisfy a specific portion of demands. The objective function was minimax regret, i.e. the locations that minimize the worst-case loss in the objective function.
3. Model formulation In this section, we describe in more detail the designing distribution networks under uncertainty which we will present. There is a set of customers (N), at which demands are generated, and a set of potential sites for DCs (M), where facilities may be opened. We assume that the demands for each customer i are generated according to a Poisson process, independent of the demand at other customers in (N). At each open facility an inventory is kept such that an arriving request finds a zero inventory (and is lost), with probability at most a. We then say that (1 a) is the fill rate of the system. In this model, we have fixed location cost for opening and operating a DC, and a fixed cost to order at any DC and a holding cost for inventory, respectively. We assume that the DCs follow a (Q, R) inventory policy and total inventory cost consists of ordering cost, holding cost, inventory lost cost and not safety stock cost. 3.1. Order quantity and service level constraint In the previous location–inventory models the required order, for each DC’s, is the sum of two parts. First part is computed based on the sum of mean demand, assigned to the DC (based on EOQ model):
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2Sh ðY ij li Þ
ð1Þ
i
and the second part is the safety stock:
hZ a
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi X LT Y ij r2i
ð2Þ
i
A new approach in this article can be proposed without dividing the order to two parts with adding service level constraint. By combining the two parts, let Yj be the total amount of the order for DC j. So, if the probability of Pr(Xj 6 Yj) is at least (1 a), we can get three advantages where Xj is random variable denoting demand distribution function for DC j. (a) The required inventory can be computed automatically based on proposed constraint, (b) we can consider service level constraint in order to prevent inventory lost and (c), in the previous opinion for safety stock inventory all types of costs such as ordering cost and transportation cost did not consider in the model and just holding cost was computed, although in practice we have to pay all costs for safety stock and our idea consider all of them simultaneously in the model. For this purpose, this service level constraint must be considered. Indeed, we do not need to compute the safety stock.
0
an arbitrary request arriving
B P@ at DC j with inventory level Y j is not lost
1 C AP1a
ð3Þ
According to this assumption, order quantity is greater than customers’ mean demand and model determines optimum quantity so that customers’ demand satisfies at least (1 a)%, so there is no necessity to hold safety stock. (If we want to hold safety stock, then we should order exactly equal customers’ demand mean and by using that safety stock we can respond to demand fluctuation. but in our study we order more than customers’ demand mean.) The service level constraint considered for this is
MinðY j jPrðX j 6 Y j Þ P 1 aÞ8j 2 M
ð4Þ
Xj: random variable denoting demand distribution function for DC j. Yj: decision variable denoting total product which DC j should supply to satisfy customers’ demand with confidence level (1 a)%.
5622
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
(In this case, Xj has Poisson distribution function.) Because customers’ demands have Poisson distribution, total product that any DC needs to order has Poisson distribution. (The summation of some Poisson distribution follows Poisson distribution.) For example, if a special DC satisfies the demand of customers 1, 2 and 3 the distribution of product that this DC needs has Poisson distribution with the rate k = k1 + k2 + k3. We have another purpose from mentioned constraint and it refers to consider ‘‘Min” for the constraint. If we use the constraint as Pr(Xj 6 Yj) P 1 a, the model in the optimal solution can satisfy the constraint with any fill rate greater than (1 a)%. So, fill rate of all DCs may be different in the optimal solution, and this can lead to unsmooth for the network. But if we use the constraint as (Min(YjjPr(Xj 6 Yj) P 1 a)) the fill rate for DCs will be the same and performance of the network can be more smooth.
Decision variables 8 < 1 If customer’s demand i satisfies X ijs ¼ by DC at site j in scenario s : 0 Otherwise Yjs: variable denoting the total amount of product from supplier to DC j in scenario s probability of occurring inventory lost in DC at site j in bjs: scenario s
3.2. Coverage radius
3.4. Model formulation
In practice, based on some restrictions such as long distances, it may be occurred that each DC cannot cover all customers. In other words, all distribution centers cannot service all of the customers because each of the DCs has special critical coverage radius and if a customer does not be in this critical coverage, the DC cannot service that customer. So, we consider this parameter in the model which yields more flexibility for the proposed model. This parameter is named coverage radius and we show it by Zij. Hwang (2002) in his article used similarity from this notice. Fig. 1 illustrates this concept by deleting some possible links between DCs and customers based on coverage radius restrictions.
The cost in a specific scenario s is computed like that in Snyder et al. (2007) and is given by:
3.3. Notation We use Sets N: M: S:
the following notation index set of customer zones index set of potential distribution sites set of possible scenarios
Parameters Cij: cost of supplying one unit of demand to customer zone i from DC at site j C js : cost of supplying one unit of demand to DC at site j from supplier in scenario s fixed cost for opening and operating DC at site j Fj: mean of demand per time unit for customer i in scenario s kis: holding cost per unit in DC at site j hj: fixed cost per order placed to the supplier by DC at site j Sj: l: cost for unit lost 1 If a DC at site j can cover customer i Z ij ¼ 0 Otherwise; ps:
probability that scenarios s occurs.
Fig. 1. The figure of supply chain network with coverage radius.
We consider location variable (U) to be independent of scenario because it must be made before it is known which scenario will be occurred. Also, allocation variables should be dependent to the scenario, so we indicate the X, Y and b variables with index of scenario. Our goal is to choose locations of DCs in order to minimize the total expected cost of the network.
Min Z ¼
X
XX
fj U j þ
j2M
j2M
X
s2S
ps C js Y js þ
!
!
X
C ij kis X ijs þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2hj Sj Y js
i22N
l
ð5Þ
X ijs Z ij ¼ 1 8i 2 N; 8s 2 S
ð6Þ
þbjs
kis X ijs
i2N
Subject to:
X j
X ijs 6 Z ij U j
8i 2 N; 8j 2 M; 8s 2 S
Y js X eð Min Y js r¼0
P
k X Z i is ijs ij
Þ P k X Z r i is ijs ij r!
ð7Þ ! P1a
8j 2 M; 8s 2 S bjs þ
Y js X eð
P
k X Z i is ijs ij
r¼0
ð8Þ Þ P k X Z r i is ijs ij ¼ 1 8j 2 M; 8s 2 S r!
ð9Þ
U j 2 f0; 1g 8j 2 M
ð10Þ
X ijs 2 f0; 1g 8j 2 M; 8i 2 N; 8s 2 S
ð11Þ
Yjs P 0 8j 2 M
ð12Þ
0 6 bjs 6 1 8j 2 M; 8s 2 S
ð13Þ
The proposed model minimizes total expected cost made of: the costs associated with opening and operating the DCs, the costs of shipments from the supplier to the DCs, the costs to serve demand of customers from the DCs, total costs of ordering and holding inventory the expected costs of inventory lost in any opened DC (the expected cost = probability of inventory lost total expected demand unit inventory lost cost). Constraint set (6) guarantees that all customers allocated to exactly one DC in each scenario with considering coverage radius. Constraint set (7) says that a customer can be allocated to a DC if that customer is in the coverage radius of DC and DC is opened. Set constraint (8) forces that service level for DCs in each scenario would be at least (1 a)%. (as we said before, because customers’ demand is Poisson distribution, so total product that DC needs has Poisson distribution and sum of the demand that DC satisfies should be less than the inventory that DC receives (Yjs) and this probability should be at least (1 a)% and this constraint
5623
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
return minimum feasible value of Yjs variable). Constraint set (9) presents the probability of occurring loss of inventory in each scenario for any opened DC
ðbjs ¼ 1 PðX j 6 Y js Þ 8j 2 M; 8s 2 SÞ Constraint sets (10)–(13) determine type of variables. 4. Solution approach
No. of DC which Cover each Customer
No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 3 2
2 2 4 6
2 2 1 7
… … … …
m m 9 2
Scenario n
4
5
2
…
6
. . .
Fig. 2. Sample of chromosome structure.
It is known that uncapacitated location problems belong to the class of NP-hard problems that cannot be solve in reasonable computational time. In this paper, we propose a nonlinear and stochastic programming model for the uncapacitated network design problem which is extremely hard to optimally solve in large sizes. Thus, we use an efficient meta-heuristic method based on genetic algorithm (GA). Also, we design a heuristic algorithm based on optimizer rules which can give us sub-optimal solutions, to measure effectiveness of the algorithm. For these reasons we use this genetic algorithm similar to Zhou and Liu (2003). 4.1. Genetic algorithm GA is a stochastic search and heuristic optimization technique, which has been widely adopted by many researchers to solve various problems. This algorithm was first developed by Holland (1975). It mimics the mechanism of genetic evolution in the biological nature and consists of a population of chromosomes (strings or individuals) that are composed of genes. These genes represent a number of values, called alleles. Each chromosome (genotype) represents one potential solution (phenotype). The process of genetic operators (i.e., crossover and mutation) is carried out in the pool; after that, an evolution is completed by creating new chromosomes (offspring). This offspring is expected to be stronger than the parents, but this may not always be true. GA does not rely on analytical properties of the function to be optimized (Goldberg (1989)). In short, GA has two major processes: (1) GA is iteratively and randomly generating new solutions; and (2) these solutions are checked for the optimality according to predefined fitness functions. This becomes the most powerful principle of GA. It makes them widely suitable for finding an optimal solution in many complex problems, such as traveling salesman problem (TSP) and any forms of scheduling problems.
4.1.3. Decoding In this section, we compute the other variables (Yjs, Uj, bjs) of the model based on definition of the feasible chromosome. To compute Yjs and bjs we do the procedure described in Fig. 3: 4.1.4. Compute the fitness of all chromosomes The rank-based evaluation function is defined as the objective function for chromosome k = 1, 2, . . . , Pop_Size. 4.1.5. Selection strategy By generating a random real number, r, from the interval [0, 1], chromosome Vk is selected as a parent provided that r < P; where the parameter P is the probability of crossover or mutation operator. 4.1.6. Genetic algorithm operators As we know, one of the important parts in GA is generating new chromosomes from current chromosomes (called parents). This process can be done by genetic operators. Operators of this algorithm are crossover and mutation. 4.1.6.1. Crossover operator. We combine the chromosomes Vk, k = 1, 2, . . . , Pop_Size by crossover operation. In order to determine the parents for crossover operation, we repeat the following process from k = 1 to pop_size: generating a random real number r from the interval [0, 1], the chromosome Vk will be selected as a
4.1.1. Representation of chromosome In our GA-based approach, each chromosome or bit string (i.e., an example solution) consists of a matrix: Assign matrix. The length of this matrix is equal to: number of scenarios by the number of customers in the network. The Assign matrix represents the assignments of the customers to DCs in any scenarios. If customer i is assigned to DC k in scenario s, its entry in the Assign matrix has value k in row s and column i. Fig. 2 shows a sample of the chromosome structure. 4.1.2. Initialize chromosome population We initialize pop_size chromosomes X k ¼ ðxki;s Þ ¼ ðxk1;1 ; xk2;1 ; . . . ; k xn;1 ; . . . ; xk1;s ; xk2;s ; . . . ; xkn;s Þ k = 1, 2, . . . , Pop_Size from the feasible region {(X)jgis(xis) 6 0, xis = 1, 2, . . . , n} randomly. (xki;s Denote number of DC which must cover customer i in scenario s in chromosome k.) First for each customer, number of DCs which can cover that customer is determined. Then a random number from those DCs is generated and customer in scenario s is allocated to that DC. (In other words Xijs = 1.)
Fig. 3. Pseudo code of the decoding procedure.
5624
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
parent provided that r < Pc, where the parameter Pc is the probability of crossover. Then we group the selected parents V 01 ; V 02 ; V 03 ; . . . to the pairs ðV 01 ; V 02 Þ; ðV 03 ; V 04 Þ; . . ., without loss of generality; let us describe the crossover operator on each pair by ðV 01 ; V 02 Þ. First, for all customers in all scenarios, we generate a random number k from the open interval (0, 1), and then the crossover operator on V0 and V 02 will produce one child X as follows in Fig. 4. In fact, in our chromosome each customer in each scenario is assigned to a special DC. For this reason, in crossover process we have two parents that each customer may be assigned to different DCs in each parent. To generate offspring, we assume that each customer can be assigned to one of them randomly and this is occurs with probability 0.5. In other words, in offspring each customer is assigned to one of the DCs which covered that customer in the parents randomly. We show sample of this process in Fig. 5. 4.1.6.2. Mutation operator. We update the chromosomes Vk, k = 1, 2, . . . , Pop_Size by mutation operation. Similar to the proves of selecting parents for crossover operation, we repeat the following steps from k = 1 to pop_size: generating a random real number r from interval [0, 1], then chromosome Vk will be selected as a parent provided that r < Pm; where the parameter Pm is the probability of mutation. For each selected parent:
V k ¼ ðxk1;1 ; xk2;1 ; . . . ; xkn;1 ; . . . ; xk1;s ; xk2;s ; . . . ; xkn;s Þ we mutate it with one of the following ways. Because we have three purposes in mutation (reducing number of opened DCs, substitution an opened DC with a closed DC to check all conditions in the solution, reassigning customers’ allocation). So, we considered three types for mutation operator and any of them may be occurred with probability 1/3 for each chromosome selected for mutation. (a) Select randomly from opened DCs (Uj = 1).If selected DC can be closed, then (a DC can be closed if all customers assigned
to it, can be assigned to another opened DC):For each customer in each scenario, allocated to that old DC, select randomly from new DC which has two below conditions and assign the customer to the new DC: – New DC j is opened in this chromosome. – New DC j can cover the customer. By this way, we can reduce the number of opened DC. We show sample of this process in Fig. 6. (b) Select randomly from opened DCs and name it j, then select randomly from closed DCs and name it j0 . With below procedure we close DC j and instead of it open DC j0 .First, for customer allocated to DC j in any scenario, if the customer can be allocated to DC j0 , then allocate it in the same scenario to DC j0 .Otherwise, allocate it to another opened DC which can cover it. Sample of mutation type 3 has been illustrated in Fig. 8.By this way, we can replace an opened DC with a closed DC for more searching in feasibility space. Fig. 7 has shown an example of this process. (c) Select a customer in a chromosome randomly. Reassign its demands to another opened DC which can cover it. Sample of mutation type 3 has been illustrated in Fig. 8.
4.2. Heuristic procedure Our purpose from solving mathematical modeling is to study the behavior of solution procedure and finally use it in practice and real world. For these reasons, we present a heuristic method based on some rules obtained from optimal model and our proposed GA that lead to optimality for mentioned model in the article. Because the model of the article cannot be solved with optimizer algorithms, we use from this heuristic algorithm to get sub optimal solutions. So, by using this heuristic and genetic algorithm we can use this model in real world and practice (Miranda & Garrido, 2004). 4.2.1. The rules Note: according to coverage radius which is one of our assumptions in modeling; our problem is a set-covering problem.
∀i ∈ N , ∀s ∈ S : ( i is index of customers) If <0.5 then
(a) As there is no capacity restriction on the distribution centers, at optimality each customer will be supplied entirely from the most cost-effective located facility. This may be fortuitous because many firms require single sourcing of customers to gain efficiencies and ensure good customer service. In other words, each customer should be assigned to the nearest opened DC (Love, Morris, & Weslowsky, 1988a).
X i , s = V1′(i, s )
Else
X i , s = V2′(i, s )
End if
Fig. 4. Pseudo code of the crossover operator.
Parent 1 No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
4 4 1 2
Parent 2
5 5 1 2
No. of Gens No. of Customers Scenario1 Scenario 2
1 1 3 1
2 2 1 1
4 4 1 2
5 5 1 2
Crossover Offspring No. of Gens No. of Customers Scenario1 Scenario 2
1 1 3 1
2 2 2 1
3 3 2 3
Fig. 5. Sample of crossover operator.
3 3 3 3
4 4 1 1
5 5 3 3
5625
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Parent No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
No. of Gens No. of Customers Scenario 1 Scenario 2
1 1 2 2
2 2 2 3
1 1 3 2
2 2 2 2
1 1 1 2
2 2 2 1
3 3 2 3
4 4 2 2
5 5 3 2
Fig. 6. Sample of mutation type 1.
Parent No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
Mutation type 2
No. of Gens No. of Customers Scenario1 Scenario 2
3 3 2 3
4 4 3 2
5 5 3 2
Fig. 7. Sample of mutation type 2.
Parent No. of Gens No. of Customers Scenario1 Scenario 2
1 1 1 2
2 2 2 1
3 3 2 1
Offspring
4 4 1 2
5 5 1 2
No. of Gens No. of Customers Scenario1 Scenario 2
3 3 1 1
4 4 1 2
5 5 1 2
Mutation type 3 Fig. 8. Sample of mutation type 3.
(b) Each customer must be assigned to the same DC in each scenario and not different. (According to the first rule). For example, if customer 1 is assigned to DC 4 in scenario 1, this customer must be assigned to DC 4 in any scenario too because DC 4 is the nearest opened DC to customer 1. (c) Minimum number of opened DC is obtained in the way that we can assign all customers to DCs considering coverage radius. Base on Zij matrix we open minimum number of DCs that each customer can be covered by at least one of them. This is the result of set-covering problem in Love, Morris, and Wesolowsky (1988b). (d) Because set-up costs and difficulties of installing DCs is too higher than transportation and inventory costs and also regarding to unlimited capacity for DCs, any opened DC can cover any number of customers which are in its coverage radius despite their demands so, in these problems (Set covering) it’s suitable to work with the minimum number of DCs that leads to less costs and cost saving for the network. (In uncapacitated set covering location problems number of opened DCs is less than capacitated set covering location problems. In basic model for set-covering problem the aim of modeling is to minimize number of opened DCs (Love et al., 1988b) so, we follow this objective too. This result is fit with logic of designing distribution networks. Because, if n DCs are opened, capacity of network to response customers’ demand is 1. In this situation if n + 1 DCs are opened, the capacity of network to response customers is 1 too. In the second case, we just increase strategic costs and saving other costs is little. Our GA confirms this result too. In Table 1, we solve large size problems by GA and as we can see, number of opened DCs in all of them is minimum. 4.2.2. Heuristic stages We divide this method into two stages. In the first stage, we find which DCs should be opened and in the second stage allocation
Table 1 Solving large size problem with GA # of Customers
# of DCs
# of Scenarios
Number of opened DCs
Interaction
70 70 70 75 75 75 80 80 80 85 85 85 90 90 90 95 95 95 100 100 *100
75 75 75 50 50 50 50 50 50 60 60 60 55 55 55 65 65 65 60 60 60
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
12–22–41 12–45 25–37–39 23–24–40 23–30–33 6–30–41 6–15–44 29–35–45 2–9 39–46–59 17–52–54 6–14–27 11–34–52 14–16–27 13–15–19 3–11–54 30–33–62 21–27–40 32–48–59 20–36–28 9–17–25
250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250 250
decisions for customers and quantity of ordering for DCs considering service level are made.
4.2.2.1. First stage. In this stage, we should find minimum number of DCs that must be opened so that all customers can be covered. For this purpose, we check opening 1 DC, then 2 DCs, 3 DCs and etc. until we find feasible solutions to assign all customers. In fact, in this section we solve set-covering problem which is mentioned in rule c. To check whether 1 DC can cover all customers, we use from below equation:
5626
Sj ¼
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
X
Z ij
j ¼ 1; 2; . . . ; n
ð14Þ
i
Above equation computes number of customers that DC j can cover. If Sk = m (m number of customers), we will finish our computations in this stage and open DC k. Sj – m for j = 1, . . . , n means that by opening just 1 DC, we cannot cover all customers, so we should find 2 DCs j and j0 that Sj;j0 ¼ m and by opening DC j and j0 we can cover all customers
0 1 XB C 0 0 Sj;j0 ¼ @Z ij þ Z ij0 Z ij Z ij0 A j; j ¼ 1; . . . ; n; j–j |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} i Ai 8 < 1 If customer i can cover 0 Ai ¼ by at least DC j or DC j : : 0 Otherwise:
ð15Þ
ð16Þ
If $j1, j2:Sj1,j2 = m, our computations in this stage are finished. If "j1, j2:Sj1,j2 – m, we should check opening more number of DCs. In this situation we present general equation as below:
1 jk P P ðZ ij1 Z ij2 Þ þ þ C B ðZ ij Þ C B j¼j1 j1
Sj1;j2;...;jk
j1
ð17Þ In general, if Sj1,j2,. . .,jk = m, we open DCs with number j1, j2, . . ., jk. For example, for k = 3:
0 Sj1;j2;j3 ¼
ðZ i;j1 þ Z i;j2 þ Z i;j3 Þ
1
C XB B ðZ i;j1 Z i;j2 þ Z i;j1 Z i;j3 þ C B C @ Z i;j2 Z i;j3 Þþ A i
ð18Þ
ðZ i;j1 Z i;j2 Z i;j3 Þ
Fig. 9. Pseudo code of computing variables.
If we want to solve the model with full enumeration for a small problem with 20 customers, 2 scenarios and 5 DCs total points that we should calculate objective function of them is We have a variable with 3 indices (Xijs) and Uj and they are (0, 1) so: Total points:
ð220 25 22 Þ 25 ¼ 220þ5þ2þ5 ¼ 232 As we can see total points are 232 that is too difficult and time consuming to enumerate all of them. These results are for a small problem. For a large size problem with 100 customers, 60 DCs and 2 scenarios total points is
When above equation for K = k is equal to m, we stop our computations in this stage. Our purpose from this stage is to find minimum number of DCs that we should open them to cover all customers and for this we first check opening 1 DC and if it is infeasible we open 2 DCs. . . At the end of this stage we may get multi solutions. (E.g., DC 1, 3 can cover all customers and DC 2, 3 can cover them too.) So we calculate next stage for all of solutions obtained from first stage.
In fact enumeration 2222 points is infeasible while our genetic algorithm get us sub optimal solution in under 3 min (problem signed with * in Table 1). So, computation time in our genetic algorithm is so little and this is one of the advantages our work that leads to validation for GA.
4.2.3. Second stage
5. Computational results
(Step-1)
We divide this section to two parts. In first part, we will present effectiveness and robustness of the genetic algorithm and in the second part we will compare genetic algorithm with heuristic method to measure optimality of genetic algorithm. In order to illustrate the effectiveness of the genetic algorithm, we give some numerical examples that are performed on a personal computer. All algorithms considered in this study were coded in Visual Basic 6 and run on a Pentium IV PC with 1.5 GHz CPU and 256MB RAM. For this purpose we solve three examples.
(Step-2)
We found the DCs which must be opened in the previous stage. In this stage, we make allocation decisions. Because DCs have unlimited capacity, so each customer must assign to the nearest opened DC in each scenario in order to have minimum transportation costs in the network. (Rules 1 and 2). For each opened DC quantity of ordering considering service level is computed as below and compute these variables:To compute Yjs we do the procedure described in Fig. 9.
ð2100 260 22 Þ 260 ¼ 2222
5.1. Measuring robustness of the GA By doing mentioned heuristic algorithm which is based on optimizer rules we can obtain solutions that are near optimal solution.
Consider a company who wants to locate new distribution centers. Assume that there are 20, 30 and 50 customers whose de-
5627
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629 Table 2 Comparison solution of example 1
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.9 0.8 0.75 0.65 1 0.95 0.8 0.75 0.75
0.33 0.30 0.27 0.23 0.18 0.33 0.28 0.23 0.18 0.22
300 300 300 300 300 300 300 300 300 300
2.3 2.3 2.5 2.3 2.3 2.3 2.3 2.3 2.3 2.3
53905.64 53851.17 53661.01 53921.09 54200.21 53789.8 53862.59 53615.76 53789.8 54029.47
851.02 817.82 817.42 745.09 806.45 860.42 785.41 740.82 860.42 806.99
0.005 0.004 0.001 0.006 0.011 0.003 0.005 0.000 0.003 0.008
Average CPU time: 92 s.
mands are Poisson distributed variable with mean kis. Also, there are 2 scenarios for kis. Suppose that a decision maker needs to select distribution centers from 8, 10 and 12 potential distribution centers to serve customers. The capacities are unlimited. In these examples service level is (1 a) = 99%. In Table 2, we compare solutions when different parameters are taken with the same generations as a stopping rule. It appears that all the minimal costs differ little from each other. In order to account for it, we present a parameter, called the percent error, i.e. (objective value – the best objective value)/the best objective value, where the best objective value is the least one of all the ten minimal costs obtained above. The last column named by ‘error’ in Table 2 is just this parameter. It follows from Table 2 that the percent error does not exceed 1.1% when different parameters are selected, which implies that the genetic algorithm is robust to the parameter settings and effective to solve model. Therefore, the variation of parameters in the algorithm has slight impact on the optimal objective, which implies that the algorithm designed in this paper is much robust. Solutions with different parameters of GA are compared on the basis of equivalent generations. In order to measure the differentia of each other, the percent error is proposed and showed in Table 3 as ‘error’. It follows from Table 3 that the percent error does not exceed 2.3% when different parameters are selected. Thus, the genetic algorithm is also robust to the parameter settings and effective to solve model. Similarly, we run the genetic algorithm for ten times with different parameters of GA on the basis of equivalent generations. In order to measure the differentia between these results, ‘error’, i.e. the percent error, is calculated and given in Table 4. From these computational results, we see that the maximal percent error does not exceed 2.3% when different parameters are chosen. Therefore, the genetic algorithm is also robust to the parameter settings and effective to solve model.
5.2. Comparing heuristic method and GA In this section, we compare solutions obtained from our GA with heuristic method and finally compute error. We solve seven examples and we describe one of them in detail and present final results for the others. Assume we have 20 customers, 6 distribution centers and 2 scenarios. We fallow the steps of the heuristic as below:
0
1 B1 B B B1 B B1 B B B1 B B1 B B1 B B B1 B B1 B B0 B Z ij ¼ B B0 B B1 B B0 B B B1 B B1 B B B1 B B1 B B1 B B @1 1
1 0 0 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1
1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 0 1
1 1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 0 1 0 1
0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
1 1 1C C C 1C C 1C C C 1C C 1C C 1C C C 1C C 1C C 1C C C 1C C 1C C 1C C C 0C C 1C C C 1C C 1C C 1C C C 0A 1
According to Zij matrix we check opening 1 DC. So for 6 DC we have S1 = 17, S2 = 15, S3 = 16, S4 = 14, S5 = 17 and S6 = 18. (m = No. of customers = 20).
Table 3 Comparison solution of example 2
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.95 0.75 0.65 0.7 0.9 0.95 0.8 0.8 0.75
0.33 0.30 0.23 0.23 0.23 0.30 0.28 0.25 0.28 0.28
300 300 300 300 300 300 300 300 300 300
3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10
60583.07 60664.29 61475.17 60974.89 61127.22 60400.95 60304.03 60073.22 60328.1 60701.14
1095 1056.3 998.8 1128.97 1154.79 1172.57 1067.19 1085.16 1093.51 1086.65
0.008 0.010 0.023 0.015 0.018 0.005 0.004 0.000 0.004 0.010
Average CPU time: 151 s.
5628
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
Table 4 Comparison solution of example 3
1 2 3 4 5 6 7 8 9 10
Pop_size
Pc
Pm
Gen
No. of opened DC
Total cost
Cost of inventory lost
Error
20 20 20 20 20 30 30 30 30 30
1 0.95 0.8 0.7 0.95 0.9 0.8 0.85 0.95 0.85
0.33 0.30 0.27 0.25 0.28 0.30 0.30 0.30 0.30 0.28
300 300 300 300 300 300 300 300 300 300
3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10 3.10
96851.03 97252.06 98411.66 98591.62 96575.56 96558.15 97077.29 96825.47 96331.58 96754.13
1982.28 2029.1 2096.08 1898.52 1954.14 1938.79 1817.6 1952.29 1819.63 1959.11
0.005 0.010 0.022 0.023 0.003 0.002 0.008 0.005 0.000 0.004
Average CPU time: 270 s.
As we can see, none of DCs have Sj = m. So by opening 1 DC we cannot cover all customers. We must check opening 2 DCs. By using equation mentioned in the algorithm (17) for DC j and j0 we have S(1, S(2, S(3, S(4, S(5,
2) = 20 3) = 18 4) = 18 5) = 19 6) = 20
S(1, S(2, S(3, S(4,
3) = 19 S(1, 4) = 19 S(1, 5) = 19 S(1, 6) = 20 4) = 19 S(2, 5) = 20 S(2, 6) = 20 5) = 19 S(3, 6) = 19 6) = 19
As we see DCs (1, 2)–(1, 6)–(2, 5)–(2, 6) and (5, 6) can cover all customers. So we compute stage two for all of them and select the best as sub-optimal solution. 1, 1, 2, 2, 5,
2 ? total 6 ? total 5 ? total 6 ? total 6 ? total
cost = 70,576 cost = 101,558 cost = 72,556 cost = 82,976 cost = 100,633
Total cost is made of fixed cost, transportation cost, inventory cost and loss of inventory cost. As we see the best objective is related to open DCs 1 and 2. We solve this problem with our GA. Total cost of this solution is = 72,072. We compute percentage error as below:
Error is ¼
72072 70576 ¼ 0:021 ¼ 2:1% 70576
Final results for the other problems are presented in Table 5. In Table 5, we compare solutions in GA and heuristic method when different problems are taken. As we see behavior of GA is similar to heuristic method. In other words, difference between them may be appeared just in allocation while location decisions are the same. In order to count error we use from this ratio: (objective function in GA – objective function in heuristic)/objective function in heuristic. It follows from Table 5 that the percent error
does not exceed 2.3%, which implies that the genetic algorithm is effective to solve model. So, solutions obtained from GA are near sub-optimal solutions obtained from heuristic method and we can use them in practice.
6. Conclusion and future directions In this paper, we introduced notation of distribution network in a supply chain system while we have service level constraint and customer’s demand is stochastic with Poisson distribution function. In the proposed model, we assumed that DCs have coverage radius restriction. We formulated the problem as a nonlinear integer program. The advantages of the proposed model are as follows: considering coverage radius which yields more flexibility for the model, considering service level constraint that by using it the required inventory can be computed automatically, it prevents inventory lost and all types of cost for inventory are consider in the model. To solve this nonlinear integer programming model we first introduced a new and robust genetic algorithm solution and then based on genetic algorithm results and some optimizer rules we presented a new heuristic method. Our computational results showed the algorithm to be quite effective. For future research we suggest three directions: (a) We can study this model while Zij is probability and this parameter does not have deterministic value. For example, it has a distribution probability function and based on it, Zij gets value 0 or 1. (b) Setting inventory capacity constraint is another development for this model. (c) Applying the other ordering methods to calculate inventory costs in designing distribution networks. (d) Aggregating this model with other assumption such as routing and considering time window can be an interesting development of the proposed model.
Table 5 Comparing heuristic with GA solution No. of customers
No. of DCs
No. of scenarios
Opened DCs in GA
Opened DCs in heuristic
Objective function in GA
Objective function in heuristic
Error
8 10 12 15 20 25 30
5 5 6 6 6 8 12
2 2 2 2 2 2 2
1 1 4 3.6 1.2 2.7 2.9
1 1 4 3.6 1.2 2.7 2.9
26745 29105.9 28566 43210.3 72072 78511 68375.12
26745 29105.9 28566 43210.3 70576 77903 66819.5
0 0 0 0 0.021 0.008 0.023
V.R. Ghezavati et al. / Expert Systems with Applications 36 (2009) 5620–5629
References Aikens, C. H. (1985). Facility location models for distribution planning. European Journal of Operational Research, 22(3), 263–279. Amiri, Ali (2006). Designing a distribution network in a supply chain system: Formulation and efficient solution procedure. European Journal of Operational Research, 171(2), 567–576. Chan, Y., Carter, W. B., & Burnes, M. D. (2001). A multiple-depot, multiple-vehicle, location-routing problem with stochastically processed demands. Computers & Operations Research, 28(8), 803–826. Conde, E. (2007). Minmax regret location–allocation problem on a network under uncertainty. European Journal of Operational Research, 179(3), 1025– 1039. Dasci, A., & Verter, V. (2001). A continuous model for production–distribution system design. European Journal of Operational Research, 129(2), 287– 298. Daskin, M. S., Coullard, C. R., & Shen, Z.-J. M. (2002). An inventory-location model: Formulation, solution algorithm and computational results. Annals of Operations Research, 110, 83–106. Chan, Felix, T. S., & Chung, S. H. (2005). Multicriterion genetic optimization for due date assigned distribution network problems. Decision Support Systems, 39(4), 661–675. Francis, R. L., McGinnis, L. F., & White, J. A. (1983). Locational analysis. European Journal of Operations Research, 12(3), 220–252. Gabor, A. F., & van Ommeren, J. C. W. (2006). An approximation algorithm for a facility location problem with stochastic demands and inventories. Operations Research Letters, 34(3), 257–263. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning. MA: Addison-Wiley. Harpera, P. R., Shahania, A. K., Gallagherb, J. E., & Bowiec, C. (2005). Planning health services with explicit geographical considerations: A stochastic location– allocation approach. Omega, 33(2), 141–152. Holland, J. H. (1975). Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence in 1992 (2nd ed.). Cambridge: MIT Press.
5629
Hwang, H. S. (2002). Design supply-chain logistics system considering service level. Computers & Industrial Engineering, 43(1–2), 283–297. Lieckens, K., & Vandaele, N. (2007). Reverse logistics network design with stochastic lead times. Computers & Operations Research, 34(2), 395–416. Love, R. F., Morris, J. G., & Wesolowsky, G. O. (1988a). Facilities location: Models and methods (p. 186). New York: North Holland Publishing Co. Love, R. F., Morris, J. G., & Wesolowsky, G. O. (1988b). Facilities location: Models and methods (p. 174). New York: North Holland Publishing Co. Lu, Zhiqiang, & Bostel, Nathalie (2007). A facility location model for logistics systems including reverse flows: The case of remanufacturing activities. Computers & Operations Research, 34(2), 299–323. Miranda, Pablo A., & Garrido, Rodrigo A. (2004). Incorporating inventory control decisions into a strategic distribution network design model with stochastic demand. Transportation Research Part E, 40(3), 183–207. Mirchandani, P. B., Oudjit, A., & Wong, R. T. (1985). Multidimensional extensions and a nested dual approach for the m-median problem. European Journal of Operational Research, 21(1), 121–137. MirHassani, S. A., Lucas, C., Mitra, G., Messina, E., & Poojari, C. A. (2000). Computational solution of capacity planning models under uncertainty. Parallel Computing, 26(5), 511–538. Snyder, L. V. (2006). Facility location under uncertainty: A review. IIE Transactions, 38, 537–554. Snyder, L. V., Daskin, M. S., & Chung-Piaw, Teo (2007). The stochastic location model with risk pooling. European Journal of Operational Research, 179(3), 1221–1238. Tsiakis, P., Shah, N., & Pantelides, C. C. (2001). Design of multi-echelon supply chain networks under demand uncertainty. Industrial and Engineering Chemistry Research, 40, 3585–3604. Zhou, J., & Liu, B. (2003). New stochastic models for capacitated location–allocation problem. Computers & Industrial Engineering, 45(1), 111–125. Zhou, G., Min, H., & Gen, M. (2002). The balanced allocation of customers to multiple distribution centers in a supply chain network: A genetic algorithm approach. Computers and Industrial Engineering, 43(1–2), 251–261. Zuo-Jun, M., Shen & Qi, L. (2007). Incorporating inventory and routing costs in strategic location models. European Journal of Operational Research, 179(2), 372–389.