Two-echelon, multi-commodity supply chain network design with mode selection, lead-times and inventory costs

Computers & Operations Research 39 (2012) 1345–1354

Contents lists available at SciVerse ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

Two-echelon, multi-commodity supply chain network design with mode selection, lead-times and inventory costs Hannan Sadjady, Hamid Davoudpour n Department of Industrial Engineering, Amirkabir University of Technology, Tehran, Iran

a r t i c l e i n f o

abstract

Available online 16 August 2011

Designing distribution networks – as one of the most important strategic issues in supply chain management – has become the focus of research attention in recent years. This paper deals with a twoechelon supply chain network design problem in deterministic, single-period, multi-commodity contexts. The problem involves both strategic and tactical levels of supply chain planning including locating and sizing manufacturing plants and distribution warehouses, assigning the retailers’ demands to the warehouses, and the warehouses to the plants, as well as selecting transportation modes. We have formulated the problem as a mixed integer programming model, which integrates the above mentioned decisions and intends to minimize total costs of the network including transportation, lead-times, and inventory holding costs for products, as well as opening and operating costs for facilities. Moreover, we have developed an efficient Lagrangian based heuristic solution algorithm for solving the real-sized problems in reasonable computational time. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Distribution Supply chain network design Facility location problems Mixed-integer programming Lagrangian heuristic

1. Introduction Facility location problems have been studied for almost a century. One could consider Weber’s studies in 1909, which led to his industrial location theory, to be the first studies in this research area. Weber’s studies were later extended by Hakimi’s in 1964. Various forms of location models have been developed since then, which can be categorized into four groups, based on ReVelle et al. [15], namely: analytic models, continuous models, network models, and discrete location models. Despite the differences between these models, they all have some main features in common. They all involve a set of customers whose locations are spatially distributed but known, and a set of facilities to serve their demands, whose locations should be determined. After the quality revolution in the late 1970s, the global competition between companies increased [3]. Context of supply chain management (SCM) was introduced in 1982 by Webber and Oliver, as a response to these competitive pressures on enterprises. In SCM, integration and coordination between a series of suppliers, manufacturers, and distributors are made to increase the competitive advantages of the entire supply chain. During recent decades, the SCM concept has been developed and attracted the attention of many researchers from different research areas. As a result, facility location models have been gradually applied in designing distribution networks—as one of the most important strategic issues in SCM [11].

n

Corresponding author. Tel.: þ98 2166466497. E-mail address: [email protected] (H. Davoudpour).

0305-0548/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2011.08.003

Among the four introduced categories of facility location models, discrete models are recognized as being more suitable for designing physical distribution networks. Discrete facility location problems could also be divided into six different groups: median problems, covering problems, center problems, uncapacitated facility location problems (UFLP), capacitated facility location problems (CFLP), and a group of problems which could be called supply chain network design (SCND) problems based on Melo et al. [11]. The first five groups are well studied in the literature (see, e.g. [9,13,15]), and have some similar characteristics, such as single facility type, deterministic and static parameters, single commodity, and location-allocation decisions. As Melo et al. mentioned in their review article, these characteristics are so far from practical settings of network design. Consequently, many extensions to models of these five groups have been developed and formed the SCND group. These extensions include considering different types of facilities and multi-echelon networks, stochastic and dynamic parameters of demands and costs, along with inclusion of multiple types of products. Since considering all of the above extensions and/or other realistic issues together leads to practically insolvable models, researchers include selective features in their studies. Pirkul and Jayaraman [14] considered two-echelon, multi-commodity, capacitated facility location problem, in which the objective was locating a given number of plants and warehouses to satisfy the customer demands with the minimum total network costs. Later, Jayaraman and Pirkul [8] extended their work by considering the raw material vendors for supplying the manufacturing plants, although no location decision was made in this level. They presented a mixed integer programming (MIP) formulation of

1346

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

the problem, and developed a heuristic solution procedure based on Lagrangian relaxation of the model. Eskigun et al. [4] studied design of an outbound supply chain network for an automotive company. The location of assembly plants and their demand areas were known, and manufactured vehicles could either be delivered by trucks or through distribution centers. Decisions to be made included determining the locations of distribution centers and choice of transportation mode, considering total network costs and lead-times. They proposed a Lagrangian heuristic to solve the presented MIP model. Wu et al. [17] considered a single-echelon network with general facility setup cost function, in which several facilities could be located in each site. The goal was to determine the location of sites and the size (number) of facilities in each open site. They formulated the problem as a MIP model and proposed a Lagrangian heuristic algorithm to solve it. Amiri [2] studied a two-echelon supply chain distribution network design problem. He proposed an integrated model to determine the locations, numbers, and also capacities of the facilities in both layers. The presented MIP model was solved by a heuristic procedure based on Lagrangian relaxation of the problem. Amiri resolved a common drawback in most of the previous studies by considering multiple capacity levels for each facility. In this paper, designing a multi-commodity distribution network including multiple manufacturing plants, warehouses, and retailers is investigated. The location of retailers and their demand for each kind of products are known in advance. These demands should be satisfied through the distribution warehouses, which receive their specific product requirements from manufacturing plants. The locations and capacities of plants and warehouses are not predetermined and should be chosen from candidate options. Also, as is the case in the real world, we have considered different transportation modes between facilities, with different characteristics including their costs, traveling times, and service frequencies. The parameters showing these specifications for each transportation mode could be varied for different combinations of origin-destination. Moreover, they could be varied for different products as we have considered multiple products with different volumes/weights. Hence, the problem is locating and sizing the plants and the warehouses, along with determining the best distribution strategy throughout the network, considering the choice of transportation mode. We have formulated the problem as a MIP model, in which the objective function is to minimize total cost of the distribution network, including the fixed costs of opening and operating the facilities, transportation costs, lead-time costs and inventory holding costs. An efficient Lagrangian based heuristic algorithm is presented for solving real-sized problems in reasonable computational time. The proposed model in this article is an extension of the model presented by Amiri [2], which resolves the major drawback in most of previous research studies by considering multiple levels of capacities available for each facility. The distinction between our research and the work of Amiri [2] is basically in the transportation mode selection and multi-commodity formulation, as well as considering lead-time and inventory holding costs in the objective function. Moreover, besides the inherent and general differences between our solution procedure and the one presented by Amiri, in particular, the heuristic algorithm and the method of updating Lagrangian multipliers in our methodology differs from Amiri’s. Considering lead-time and inventory holding costs independently in the objective function enables managers and decision makers to examine different scenarios. Managers can obtain better insight into the quantitative aspects of different scenarios, through validating and analyzing the solution results by performing sensitivity analysis or employing other analytical methods. The relatively small running time of the proposed solution procedure enables managers to rule out non-desirable scenarios and identify good solutions, in a short time.

The remainder of this paper is organized as follows. Section 2 is dedicated to mathematical formulation of the problem. In Section 3, a Lagrangian relaxation of the presented model is introduced. The heuristic solution procedure is explained in Section 4. The computational experiments are reported in Section 5, followed by the conclusions in Section 6.

2. Mathematical formulation We use the following notations for mathematical formulation of our model. Note that parameters are considered on annual basis in this article; however, any other units of time can be applied. Index sets: I J K L T E G

set of retailer locations, indexed by iAI set of potential warehouse sites, indexed by jAJ set of potential plant sites, indexed by kAK set of different product types, indexed by lAL set of available transportation modes, indexed by tAT set of possible capacity levels for warehouses, indexed by eAE set of possible capacity levels for plants, indexed by gAG

Parameters fwej fpgk dli wcje pckg vl cwrijlt

cpwltjk

twrijlt tpwltjk mlt whlj phlk wsfijlt psfjklt

fixed annual cost of opening and operating a warehouse at site j, with capacity level of e fixed annual cost of opening and operating a plant at site k, with capacity level of g annual demand of retailer i for product l capacity level e for warehouse j capacity level g for plant k unit volume (weight) of product l unit cost of transporting product l by mode t from warehouse j to retailer i (including the warehouse operational costs for a unit of product l) unit cost of transporting product l by mode t from plant k to warehouse j (including the unit manufacturing cost of product l at plant k) delivery lead-time of product l from warehouse j to retailer i via mode t delivery lead-time of product l from plant k to warehouse j via mode t monetary value per unit of lead-time for product l in mode t annual inventory holding cost for one unit of product l at warehouse j annual inventory holding cost for one unit of product l at plant k service frequency of mode t for product l from warehouse j to retailer i service frequency of mode t for product l from plant k to warehouse j

Decision variables Xijlt

lte Yjk

fraction (with respect to dli ) of retailer i’s demand for product l delivered from warehouse j via transportation mode t fraction (with respect to wcje ) of product l delivered from plant k to warehouse j with capacity level of e via transportation mode t

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

Wje

binary variable equals to 1 if a warehouse with capacity level e is opened at site j, otherwise to 0 binary variable equals to 1 if a plant with capacity level g is opened at site k, otherwise to 0

Pkg

Considering the above notation, the model can be stated as follows: Problem P. XX X X g g XXXX Min ZP ¼ fwej Wje þ fpk Pk þ cwrijlt dli Xijlt jAJ eAE

þ

X X XXX j A J k A K l A L t A Te A E

þ

X X XXX

iAI jAJ lALt AT

k A Kg A G lte cpwltjk wcje Yjk þ

XXXX mlt twrijlt dli Xijlt iAI jAJ lALt AT

lt

m

lte tpwltjk wcje Yjk þ

j A J k A K l A L t A Te A E

þ

X X XXX 1 j A J k A K l A L t A Te A E

S.t. XX Xijlt ¼ 1

2

XXXX 1 2 iA I jA J l A Lt A T

whlj wsfijlt dli Xijlt

lte phlk psfjklt wcje Yjk

8 iA I, l A L

ð1Þ

jAJ tAT

X XXX dli vl Xijlt r wcje Wje iAI lALt AT

X

Wje r1

8 jAJ

ð2Þ

eAE

8j A J

ð3Þ

eAE

XX X XX dli vl Xijlt r wcje Yjklte iAI tAT

XXXX

lte wcje Yjk r

j A J l A L t A Te A E

X

8j A J, l A L

ð4Þ

8k A K

ð5Þ

k A Kt A Te A E

Pkg r 1

X

pckg Pkg

gAG

8kA K

ð6Þ

gAG

Xijlt Z 0

8i A I, j A J, l A L, t A T

lte Z0 Yjk

8j A J, k A K, l A L, t A T, eA E

ð7Þ ð8Þ

Wje A f0,1g

8j A J, eA E

ð9Þ

Pkg A f0,1g

8k A K, g A G

ð10Þ

1347

results in the presented model is more reasonable and better reflects the trade-off between the speed and cost of transportation modes. Finally, the last two summations represent the average annual inventory holding costs of warehouses and plants. The service frequency of modes is related to the number of times that products are transported between facilities in one year. For example, a product may be shipped once, twice, or more frequently in a year. If mode t transports product l from plant k to warehouse j, 4 times in a year, then its service frequency for that product would be psfjklt ¼ 1=4. Therefore, the average inventory holding cost for such product, would be one forth of a product being shipped only once in a year. The service frequency of modes may be determined from exogenous factors, such as restricted time tables or schedules for road, rail, air and water modes, or dictated by 3PL companies in the contract. They may also be determined by endogenous factors, including management decisions to restrict the number of shipments due to economies of scale or other considerations. The reason of multiplying these terms by ½ is to estimate the inventory holding cots for the average amount of inventory. The presented model can be used for evaluating different scenarios (in this case different service frequencies), and determining the best distribution strategy. Constraint set (1) guarantees that the demand for each commodity of each retailer is totally satisfied through warehouses. Constraint set (2) ensures that the distributed products from each of the open warehouses do not exceed its capacity limit. Constraint set (5) ensures capacity limits for plants, similar to constraint set (2) for warehouses. Constraint sets (3) and (6) guarantee that in each potential site, at most one facility (warehouse and plant, respectively) with only one capacity level could be located. Constraint set (4) implies that for each commodity, total retailers’ demand satisfied by a warehouse must be in balance with the total available units of that product at the warehouse, which have been transported from the open plants. The non-negativity restrictions on the corresponding decision variables are enforced by constraint sets (7) and (8). Finally, constraint sets (9) and (10) are due to the binary nature of decision variables. We can extend the introduced model in several ways to cope with different situations in reality. For example, in some studies (e.g. [8,14]) there is a limitation on the total number of facilities. This could be for a result of budget limitations or because of the pre-located facilities. In general, if the total number of facilities to be located is supposed to be between fmin and fmax, additional constraint sets could be introduced. Constraint sets (11) and (12) impose these limits on warehouses and plants, respectively, XX wmin r Wje rwmax ð11Þ jAJ eAE

The objective function minimizes total variable and fixed costs of the network. The first two summations are annual fixed costs of opening and operating the warehouses and plants, respectively. The third term represents the variable costs of storing and transferring the retailers’ demands from warehouses. The forth term is associated with the variable production costs of products and their transportation costs from plants to warehouses. The subsequent two terms are delivery lead-time costs from warehouses to retailers, and from plants to warehouses, respectively. Since our model allows for different transportation modes, the delivery lead-time of products may vary from one mode to another. The lead-time cost estimates the monitory value of products from the decision maker’s point of view. The costs imposed by possible tardiness penalties, back orders, and lost sales, as well as the deterioration and obsolescence of products can be taken into account in these terms. Taking these issues into consideration may lead the decision maker to choose a transportation mode, which is not necessarily the one with the smallest transportation cost. It is possible to consider these terms as parts of transportation costs; however, producing different scenarios and analyzing the

pmin r

XX

Pkg r pmax

ð12Þ

k A Kg A G

Another important extension of the basic model could be imposing capacity constraints on transportation modes. In the basic model, we assumed infinite capacity for transportation modes. This assumption could be justified when we use the public transportation systems or third party logistics (3PL) companies for distributing our goods. This is not far from the real world scenarios, especially in using transportation modes with high investment and maintenance costs (like air, rail, and water modes). However, whenever it is applicable, we can consider the capacity limits on individual transportation modes, by introducing the following constraint sets: X dli vl Xijlt r captij Zijt 8i A I, j A J, t A T ð13Þ lAL

XX

lte t wcje Yjk rcaptjk Zjk

t A Te A E

8kA K, j A J, t A T

ð14Þ

1348

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

In above equations, captab is the capacity of mode t between t location sites a and b. The term Zab is an integer decision variable which determines the required number of mode t for delivery of goods between a and b. Considering constraint sets (13) and (14), the fixed costs of using transportation modes should also be included in the objective function XXX XXX t t Z 0 ¼ ZP þ vcijt Zijt þ vcjk Zjk ð15Þ iAI jAJ tAT

Problem LR. MinZLR ¼ ZP þ

þ

XX

XX X XX lte dli vl Xijlt  wcje Yjk

jAJ lAL

iAI tAT

blj

XX

0

gli @

iAI lAL

k A Kt A Te A E

XX Xijlt 1A

ð18Þ

jAJ tAT

Subject to (2), (3), (5)–(10), (16) and (17).

j A J k A Kt A T

t in the new summations represents the unit fixed The term vcab cost of using transportation mode t between a and b. Hence, the model would determine the required number of transportation modes between each pair of facilities by replacing the objective function of the basic model with Eq. (15), and including the constraint sets (13) and (14) in the new model. The mixed integer programming formulation of Problem P presented in this section is an extension of the capacitated facility location problem (CFLP), which includes the classical uncapacitated facility location problem (UFLP), as a special case. The UFLP is proved to be NP-hard (see [7,10]); therefore, the proposed model in this article also belongs to NP-hard class. The computational times required to solve reasonably sized instances of the proposed model can grow prohibitively long, as there is no polynomial time algorithm for solving NP-hard problems. Therefore, we will present a heuristic solution procedure based on Lagrangian relaxation of the model in the following sections.

Problem LR could be further decomposed into two distinct subproblems LR1 and LR2: Problem LR1. XX XXXX l MinZLR1 ¼ fwej Wje þ ðcwrijlt dli þ gli þ bj dli vl ÞXijlt jAJ eAE

iAI jAJ lALt AT

XXXX XXXX 1 whlj wsfijlt dli Xijlt þ mlt twrijlt dli Xijlt þ 2 iAI jAJ lALt AT iA I jA J l A Lt A T ð19Þ Subject to (2), (3), (7), (9) and (16). Problem LR2. X X g g X X XXX l lte fpk Pk þ ðcpwltjk bj Þwcje Yjk MinZLR2 ¼ j A J k A K l A L t A Te A E

k A Kg A G

þ

X X XXX

lte mlt tpwltjk wcje Yjk

j A J k A K l A L t A Te A E

þ

X X XXX 1 j A J k A K l A L t A Te A E

3. Lagrangian relaxation and decomposition of Problem P Many problems could be considered as easy problems, which a set of side constraints has made them complicated. Based on this observation, Lagrangian relaxation scheme simplifies the original problem by dualizing the complicated constraints and adding their penalty functions into the objective function. Solving the relaxed problem in minimization models provides a lower bound for the main problem. Different methods, such as subgradient optimization algorithm, may then be used to find the best Lagrangian multipliers in order to achieve the highest possible lower bound. Also, feasible solutions to the original problem can be constructed from the solutions of Lagrangian relaxation, using heuristic algorithms. Any feasible solution to the main problem is an upper bound for the optimal objective function value. The availability of both lower and upper bounds is an advantage of Lagrangian heuristic procedure over other heuristics and metaheuristics approaches, as the relative gap between these two bounds can be used as a measure of the maximum error of the solution. The Lagrangian relaxation scheme has been used successfully for solving various distribution network design models, as well as other optimization problems. Detailed discussion on Lagrangian relaxation methodology is beyond the scope of this paper and the interested reader is referred to Ref. [1,5,6]. We obtain the Lagrangian relaxation of the Problem P by relaxing l the constraint sets (1) and (4), using multipliers gli A R and bj A R þ associated to these constraint sets, respectively. Also, the following redundant constraint sets are imposed into the model, in order to tighten the lower bounds gained from the relaxed model: X Xijlt r 1 8i A I, j A J, l A L ð16Þ tAT

XXX

1

!

2

lte phlk psfjklt wcje Yjk

ð20Þ

Subject to (5), (6), (8), (10) and (17). Problem LR1 and LR2 can be solved independently and their solution could be used to provide a solution for the Lagrangian Problem LR. The solution procedure for LR1 and LR2 are described in the following subsections. Knowing the optimal value of these problems, we can calculate the optimal value of Problem LR based on the following equation, the proof of which is obvious: XX MinZLR ¼ MinZLR1 þ Min ZLR2  gli ð21Þ iAI lAL

3.1. Solution procedure of LR1 Problem LR1 can be separated into 9J9 sub-problems (one for each potential warehouse site jAJ), which is denoted here by LR1j. Leaving constraint set (3) aside, each of these sub-problems could be further decomposed into 9E9 independent sub-problems, denoted by LR1ej . Problem LR1ej . XXX XXX l Min ZLR1ej ¼ ðcwrijlt dli þ gli þ bj dli vl ÞXijlt þ mlt twrijlt dli Xijlt iA I l ALt A T

þ

XXX 1 2 iAI lALt AT

iAI lALt AT

whlj wsfijlt dli Xijlt þ fwej Wje

Subject to XXX dli vl Xijlt r wcje Wje iAI l A Lt A T

X Xijlt r 1

8iA I, l A L

tAT

Xijlt Z0

8i A I, l A L, t A T

Wje A f0,1g Yjklte r 1

8j A J, k A K

ð17Þ

l A L t A Te A E

The solution to LR1j is given by the following equation. If we denote the optimal value of problem (  ) by O(  ), then for all jAJ:

Thus, we get a Lagrangian relaxation of Problem P, denoted by LR, given below.

OðLR1j Þ ¼ MinfMine A E ðOðLR1ej ÞÞ,0g

ð22Þ

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

Proof. For all jAJ, the binary decision variable Wje is either equal to 0 or 1. If Wje is equal to 0, it is obvious that O(LR1j) is equal to 0, as the decision variable Xijlt would be equal to zero for all iAI, lAL, and tAT because of the constraint set (2). Also, constraint (3) in problem LR1j ensures that at most one capacity level will be selected at site j, so that Wje could be equal to 1 only for one of the eAE, and that would be the e with the smallest objective value ðMinðOðLR1ej ÞÞÞ. & eAE

Considering the above observation, for solving LR1j, we only need to solve LR1ej for all eAE, with preset Wje ¼ 1. The solution corresponding to capacity level e with the most negative objective function is retained as the solution for LR1j, and Wje ¼ 0 for the other capacity levels. If there is no such capacity level, decision variables Xijlt and Wje , and the optimal objective value of LR1j would be zero vectors. It should be mentioned that when we set the value of Wje to 1, Problem LR1ej becomes a linear programming problem which could be solved by any LP solver; however, as it has a similar structure to the continuous knapsack problem, it is more efficient to use the enumerative approach presented in the appendix. The optimal solution to Problem LR1ej can be computed extremely rapidly using the described procedure, which is based on the greedy algorithm of continuous knapsack problems. The greedy algorithm for the continuous knapsack problem and the proof of its optimality is described in Ref. [12, p. 450].

1349

lower bound on the optimal value of Problem P. The quality of this bound is depended on the value of Lagrangian multipliers ((g,b) vector). The Lagrangian multiplier vector which leads to the best lower bound is derived by solving the below problem n

ZLR ðgn , b Þ ¼

Max

g A R, b A R þ

fðZLR ðg, bÞg

As computation of the optimal set of multipliers is difficult, here we use an iterative solution procedure based on the subgradient optimization method (see, e.g. [1,6]). Each iteration of this solution procedure consists of three major steps: 1. Solving the LR problem using current Lagrangian multipliers (lower bound estimation). 2. Finding a feasible solution based on LR results (upper bound calculation). 3. Updating the Lagrangian multipliers. The overall solution procedure is described in Fig. 1. In the previous section, we explained the first step of this iterative procedure. In this section, the Lagrangian based heuristic algorithm for finding feasible solutions, and the multipliers updating technique are described.

3.2. Solution procedure of LR2

4.1. Heuristic feasibility algorithm

Similar to the solution approach described for Problem LR1, Problem LR2 can be decomposed into 9K9 sub-problems (one for each potential plant site kAK), denoted here by LR2k. Leaving constraint set (6) aside, each of these sub-problems could be further decomposed into 9G9 independent sub-problems, denoted by LR2gk .

The solution obtained from Problem LR may violate some of the constraints which have been relaxed. In this subsection, we introduce a heuristic procedure, in which these infeasible results are converted into feasible solutions. This heuristic algorithm consists of two steps. In the first step, it determines the open warehouses and their capacities based on the LR1 results. Also, the distribution strategy for satisfying the retailers’ demands is decided in this step. The second step decides about the open plants and their capacities, using the results of the first step and the solution to Problem LR2. Then, it assigns the open warehouses to open plants. Three forms of violation could occur in the first step, regarding the solution to LR1. First of all, the total capacity of open warehouses may be insufficient for satisfying the demands of all retailers. Second, there may be some retailers whose demands are allocated to more than one warehouse. Also, some retailers may not be assigned to any of the open warehouses. The first step of heuristic algorithm for feasible assignment of retailer demands to open warehouses is as follows:

Problem LR2gk . XXXX l lte ðcpwltjk bj Þwcje Yjk Min ZLR2g ¼ k

j A J l A L t A Te A E

þ

XXXX

lte mlt tpwltjk wcje Yjk

j A J l A L t A Te A E

þ

XXXX 1 j A J l A L t A Te A E

Subject to XXXX

2

lte phlk psfjklt wcje Yjk þ fpgk Pkg

lte wcje Yjk r pckg Pkg

j A J l A L t A Te A E

XXX

lte Yjk r1

8j A J

l A L t A Te A E lte Yjk Z0

8j A J, l A L, t A T, eA E

Pkg A f0,1g The solution to LR2k is given by the following equation, whose proof is similar to Eq. (22); therefore, we omit it. OðLR2k Þ ¼ MinfMing A G ðOðLR2gk ÞÞ,0g

ð23Þ

Following the same strategy, we need to solve 9G9 linear problems to find the optimal solution for LR2k. Again, any LP solver could be used for solving these problems, but we include an efficient solution procedure based on the special structure of these problems in the appendix. 4. Heuristic solution procedure The solution to Lagrangian relaxation sub-problems LR1 and LR2, which has been described in the previous section, provides a

1. If the total capacity of open warehouses is insufficient for satisfying the retailer demands, do the following sub-steps, else go to step 2. a. For each of the open warehouses, calculate cost of increasing the current capacity to the next higher level, and find the cheapest choice. b. Find the cheapest potential warehouse and its capacity to be opened. c. Choose one of the above strategies for increasing total capacity, which leads to the smallest increase in LR1 objective value; then go back to step 1. 2. For each of the retailers and commodities, calculate the unit l assignment cost (cij ) to each of the open warehouses, as follows:   1 clij ¼ Min cwrijlt þmlt twrijlt þ whlj wsfijlt tAT 2 3. For each of the retailers and commodities that are not assigned yet, determine two open warehouses with available capacities (jli1 and jli2 ), which are associated to the two least unit assignl l l ment costs (cij Z cij2 Z cij1 ).

1350

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

Fig. 1. Overall solution procedure.

4. Calculate unit penalty rate (UPRli ) for each of the retailers and commodities that are not assigned yet, as follows: l

UPRli ¼

6. Repeat steps 2 through 5 until all the retailers’ demands for all products are satisfied.

l

ðcij2 cij1 Þ vl

This term shows the penalty of assigning demand of retailer i for product l to the second cheapest warehouse, per unit of warehouse capacity usage. 5. Assign the retailer demand for associated product with the highest unit penalty rate to the cheapest warehouse: if it is possible, satisfy the total retailer demand from the warehouse and update the available capacity of the warehouse; otherwise, assign the remaining available capacity of warehouse for satisfying the retailer demand; remove the warehouse from the list of available warehouses and update the retailer demand to its new value.

Here, the first step of the heuristic algorithm is completed. Since the open warehouses and the assignment of retailers to them are determined in the first step, the second step of the algorithm would have similar structure to its first step. In the second step, the open plants and their capacities, as well as assignment of open warehouses to them are determined based on the solution to LR2 and the results of the first step. At this stage, the heuristic procedure provides us a feasible solution, which is an upper bound for the optimal solution of Problem P. 4.2. Updating the Lagrangian multipliers Adjusting the Lagrangian multipliers is the final step in each iteration of overall solution procedure. In this step, the multipliers

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

are updated based on the results of two previous steps, using the following equations. We tried some other methods for updating the multipliers (e.g. [2,16]); however, the gap percentage between the bounds of those methods was much bigger in comparison with the method presented in this article 1

ðgli Þq þ 1 ¼ ðgli Þq þ yq ðdev1li Þ l

l

8i A I, l A L

2

ðbj Þq þ 1 ¼ Maxð0,ðbj Þq þ yq ðdev2lj ÞÞ

ð24Þ

8j A J, p A P

ð25Þ

where in above equations: XX dli ðXijlt Þq dli 8i A I, l A L dev1li ¼

ð26Þ

jAJ tAT

dev2lj ¼

XX X XX lte dli vl ðXijlt Þq  wcje ðYjk Þq iAI tAT

8j A J, l A L

ð27Þ

k A Kt A Te A E

y1q ¼ lq P

BUBðZLR Þq P l 2 iAI l A L ðdev1i Þ

ð28Þ

y2q ¼ lq P

BUBðZLR Þq P l 2 jAJ l A L ðdev2i Þ

ð29Þ

lte The terms ðXijlt Þq and ðYjk Þq are decision variables of Problem LR at iteration q. The best upper bound value among the first q feasible solutions is denoted by BUB and (ZLR)q is the lower bound value at iteration q. lq is a scalar with the initial value of 2. Whenever the best lower bound found so far has failed to increase in a specified number of successive iterations (here denoted by Nmax), lq would be reduced by a factor of 2, and the Lagrangian multipliers would be reset to their values which had previously led to finding the best lower bound. The iterative solution procedure described in this paper terminates after a predefined number of iterations (here denoted by Max_Iteration), or earlier if the gap percentage between the bounds shrinks to a value less than %GAPmin (Gap ¼ ðBUBBLBÞ=BLB  100).

5. Computational results In this section we will describe the computational experiments which have been designed for evaluating the performance of our proposed solution approach. The algorithm was coded in Matlab 7.0.0 and all computational studies have been implemented on a Pentium 4 PC with 2.4 GHz CPU and 1 GB of RAM. Eight problem sets with different sizes and structures are generated as it has summarized in Table 1. For small-sized problems, the Lagrangian heuristic (LH) solution is compared with the optimal solution resulting from LINGO 8.00 optimization software. Since LINGO or other commercial software cannot solve the real-world largesized problems, test problems of this kind are only solved by the

1351

proposed LH. Near-optimal solutions for small-sized problems, and small gaps between the solution bounds in large-sized problems, show that we can trust the presented algorithm for both small and large-sized problems. For each problem size, we have considered 4 different problem types according to their capacity and cost structures. The variable costs (transportation, lead-time and inventory holding costs) represent 10–20% of total network cost in the first two problem types, and 40–50% in the 3rd and 4th types. Also, total available capacities for facilities in the 2nd and 4th problem types are 1.5 to 2 times greater than the 1st and 3rd types. For each of these problem types, a series of problems with different parameter values are solved in order to achieve a reasonable level of confidence in the performance of the solution procedure on that problem structure. In each test problem, the values of parameters are generated based on the following orders. The distances among the locations of retailers and warehouses, and warehouses and plants are generated from a uniform distribution between 150 and 1500. Retailers’ demands for each product type are generated by choosing a random integer from the set {10, y, 100}. The volume (weight) of each product type is then selected randomly from the set {1, y, 10}. In practice, products are categorized into different groups based on their characteristics, e.g. volume (weight), value, hazard, fragility, security, etc. Here, the products are divided into 4 different groups based on their volumes (weights), as described in Table 2. For generating cost coefficient cwrijlt , we first compute an average unit transportation cost cwrijt for every transportation mode t, every retailer i, and every warehouse j, by multiplying the distance between i and j by a random integer chosen from the set {10, y, 20}. Then, for every product l, the cost coefficient cwrijlt is drawn from a uniform distribution in the interval ½y1 cwrijt , y2 cwrijt , where y1 and y2 are the coefficients determined in Table 2. The cost coefficient cpwltjk is generated by a similar method. Note that the unit variable cost of manufacturing and warehousing could be included in cost coefficients cwrijlt and cpwltjk , respectively. For generating twrijlt and tpwltjk , an average delivery lead-time Llt for each product l and transportation mode t is first selected randomly from the set {1, y, 15}. Then, for every retailer– warehouse and warehouse–plant pair, twrijlt and tpwltjk are computed, Table 2 Product categories and their related coefficients. Coefficients

1 2 3 4

Volume range

[1,2] [3,5] [6,8] [9,10]

Coefficients

h1

h2

0.5 0.5 0.75 1

0.75 1 1.25 1.25

Table 1 Size of test problems. Problem code

Small sizes 5-3-2-3-2-5 10-5-3-5-3-5 15-8-5-3-5-5 20-10-5-3-3-5 Large sizes 20-15-10-3-3-5 40-20-10-3-3-5 50-30-15-3-3-5 100-30-20-3-3-5

No. of customers

No. of warehouses

No. of plants

No. of products

No. of transportation modes

No. of capacities (9E9¼9G9)

5 10 15 20

3 5 8 10

2 3 5 5

3 5 3 3

2 3 5 3

5 5 5 5

5 20 40 50 100

3 15 20 30 30

2 10 10 15 20

2 3 3 3 3

2 3 3 3 3

5 5 5 5 5

1352

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

respectively, from the uniform distribution in the interval ½y3 Llt , y4 Llt , where y3 and y4 are the coefficients based on the origin-destination distances, as described in Table 3. For each product l and transportation mode t, mlt is chosen randomly in the interval ½y5 m , y6 m , in which m  is an average monetary value generated randomly between 500 and 2000, y5 ¼0.75 and y6 ¼ 1.25. For each warehouse j, an average annual unit inventory holding cost wh j is first chosen randomly from the interval [18,200, 25,500], then for each product l, whlj is drawn from a  uniform distribution in the interval ½y1 wh j , y2 whj . Similarly, for every plant k and product l, ph is chosen from the interval k [27,300, 38,300], and phlk is generated from the interval  ½y1 ph k , y2 phk . To generate the service frequencies, average values  of wsfij and psfjk are first selected randomly from the interval [0.04, 1], for each retailer–warehouse and warehouse–plant pairs, respectively. Then, for every mode t and product l, wsfijlt and psfjklt values are computed uniformly from the ½y5 wsfij , y6 wsfij  and ½y5 psfjk , y6 psfjk  intervals, respectively. For low capacity instances (1st and 3rd problem types), the capacity levels of facilities are determined as follows. First an average capacity level for each facility is computed based on the following equation, where 9F9 is the number of potential facilities (9J9 for warehouses and 9K9 for plants): P P dl vl cap ¼ i A I l A L i ð30Þ 9F9

levels pckg for each potential plant k is computed as pck1 ¼ 0:5cap, pck2 ¼ 0:75cap, pck3 ¼ cap, pck4 ¼ 1:25cap and pck5 ¼ 1:5cap. For high capacity instances (2nd and 4th problem types), these values are simply multiplied by a random number between 1.5 and 2. Finally, for generating the fixed setup cost of opening and operating a warehouse at site j, first an average fixed cost fw j for unit of capacity at site j is selected randomly between 15,000 and 20,000. Then, for each capacity level e, the annual fixed setup cost fwej is computed by multiplying fw j by the available warehouse e capacity (fwej ¼ fw j wcj ). Similarly, the average fixed setup cost fp k for unit capacity of a plant at site k, is randomly selected from g the interval [75,000, 100,000], and fpgk is computed as fpgk ¼ fp k pck for gAG. For high variable cost instances (3rd and 4th problem types), the facility fixed costs are multiplied by 0.2. Several instances of all problem sizes and types have been created using the above definitions for parameter generation. The computational results for small-sized problems are reported in Table 4. We have also summarized the results for large-sized problems in Table 5. In the both tables, the results are described by providing the problem code and type, as well as the number of variables and constraints of the model. Also, the average and maximum (worst) CPU times for both methods are reported in seconds. For small-sized problems, %Error is introduced as a quality criterion of the LH solution in comparison with the optimal solution resulting from LINGO:

Five different capacity levels wcje for each potential warehouse j are then computed aswcj1 ¼ 0:4cap,wcj2 ¼ 0:6cap, wcj3 ¼ 0:8cap, wcj4 ¼ cap and wcj5 ¼ 1:2cap. Similarly, the available capacity

%Error ¼

Table 3 Origin-destination distance ranges and their related coefficients. Category no.

1 2 3

Distance range

[150,500] [501,1000] [1001,1500]

Coefficients

h4

h3

1.25 1.75 2.25

1 1.25 1.75

BUBLINGO optimal value  100 LINGO optimal value

ð31Þ

The maximum (worst) and average values of this criterion are reported for each type and size of the small-sized problems in Table 4. It could be observed from Table 4 that the percentage error of proposed solution approach ranges from 0.00% to 2.62% with an average value of 0.46%, which shows that the solution procedure does well in finding near optimal solutions. Also, we could observe that unlike LH computational time, the running time of LINGO optimization software is completely related to the parameter configurations and problem types. As an example, Table 4 shows that for each problem size, the LINGO computational time of the fourth problem type, is obviously less than the other types. Also, despite the Lagrangian heuristic method, the average CPU time for LINGO is

Table 4 Computational results for small-sized problems. Problem code

Type

%Error

LH Time (s)

LINGO Time (s)

Variables

Avg.

Total

Max

Avg.

Max

Avg.

Max

Constraints Integer

5-3-2-3-2-5

1 2 3 4 Avg.

0.00 0.00 0.08 0.00 –

0.00 0.00 0.02 0.00 0.00

6 5 5 4 –

5 5 4 3 4

5 19 8 3 –

4 13 6 2 6

295

25

305

10-5-3-5-3-5

1 2 3 4 Avg.

0.31 0.86 1.37 0.14 –

0.11 0.30 0.45 0.05 0.22

27 28 27 25 –

25 27 26 23 25

18 24 21 10 –

12 20 15 8 13

1915

40

1967

15-8-5-3-5-5

1 2 3 4 Avg.

0.09 0.90 0.76 2.62 –

0.06 0.53 0.41 2.16 0.79

30 32 29 62 –

25 27 28 35 28

229 888 171 22 –

158 342 42 11 138

4865

65

4896

20-10-5-3-3-5

1 2 3 4 Avg.

1.04 1.3 0.94 1.88 –

0.91 0.69 0.79 0.98 0.84

41 43 43 40 –

40 42 41 38 40

2227 1053 54 25 –

1285 736 34 15 517

4125

75

4171

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

1353

Table 5 Computational results for large-sized problems. Problem code

Type

%GAP

LH time (s)

Max

Avg.

Max

20-15-10-3-3-5

1 2 3 4 Avg.

3.15 2.1 1.84 6.03 –

2.75 1.94 1.63 2.56 2.22

58 59 59 57 –

40-20-10-3-3-5

1 2 3 4 Avg.

4.1 1.88 2.49 3.83 –

2.86 1.45 1.77 1.62 1.92

50-30-15-3-3-5

1 2 3 4 Avg.

4.82 2.57 2.69 8.4 –

100-30-20-3-3-5

1 2 3 4 Avg.

3.85 2.48 4.23 4.09 –

Variables Avg.

Total

Constraints Integer

57 57 56 53 55

9575

125

9606

149 148 148 143 –

138 146 147 141 143

16,350

150

16,441

2.73 2.53 1.98 4.52 2.94

276 274 270 258 –

272 271 266 252 265

33,975

225

34,081

2.95 1.86 3.24 2.61 2.66

734 728 730 778 –

729 724 723 751 731

54,250

250

54,491

increasing exponentially with rise in the number of binary variables. This makes it impossible to find the optimal solution to real-sized problems in a reasonable time (as an example LINGO failed to find the optimal solution to the first (smallest) instance of the large-sized problems, after 24 h of running). To evaluate the performance of Lagrangian heuristic approach for large-sized problems, we have included the maximum and average percentage gaps, as well as the LH times, in Table 5. It is evident from Table 5 that the percentage gap is correlated neither to the size, nor to the type of problems. As Table 5 indicates, the average percentage gap for large-sized problems is 2.43%, with the maximum of 8.4%. The maximum running time of Lagrangian heuristic for the largest test problem is 778 s, which shows the time efficiency of the proposed algorithm in comparison with LINGO. It is evident from the computational results of Tables 4 and 5 that the proposed heuristic solution is an effective method for both small and large-sized problems with near-optimal solutions in reasonable computing times. Therefore, the presented model and solution approach can be used for obtaining a better insight into supply chain network design.

cope with different situations in reality. In fact, the presented solution procedure can be used for solving the model in presence of constraint sets (11) and (12), with small modifications (by setting the fixed cost of opening the pre-located facilities to zero, and choosing the fmax ‘‘best’’ (minimum cost) solutions obtained after solving the LR1ej and LR1gk sub-problems). However, the presented solution procedure needs major modifications for solving the model in presence of constraint sets (13) and (14), which impose capacity limitation on transportation modes. Presenting proper solution procedure for solving such problem will be the subject of subsequent research. We also note that the proposed model and solution procedure can be easily extended to cope with multi-echelon network design problems. In this case, the decomposition of the main problem into sub-problems will remain the same, after relaxing the flow constraints between each consecutive level of facilities. The feasible solutions can then be obtained, using the same heuristic approach. Application of our model in practical situations and problems, as well as investigation of some further extensions is currently under research.

Appendix 6. Conclusion In this paper, we have considered the problem of designing a two-echelon supply chain network, which allows multiple levels of capacities for the facilities of both stages. The problem involves locating and sizing the plants and warehouses and determining the best strategy for distributing the products, as well as the choice of transportation modes. We have formulated the problem as a mixed integer programming model which intends to minimize total costs of the network, including transportation, leadtimes, and inventory holding costs for products, as well as opening and operating costs for facilities. A heuristic solution procedure based on Lagrangian relaxation of the problem has been developed due to computational complexity of the model. The computational experiments indicate that the proposed heuristic algorithm is an effective and efficient solution approach for both small- and large-sized problems. As it is stated in the article, the constraint sets (11)–(14) are included as examples of possible extensions of the basic model to

In the following, the solution algorithm for solving LR1ej problem with a preset value of 1 for Wje is presented. In this case, the problem becomes a continuous knapsack problem with the constant term of fwej in the objective function. The greedy solution algorithm for this very problem, which offers the optimal solution, is as follows. The interested reader can find the greedy algorithm of the general form of the continuous knapsack problem, and the proof of its optimality in Ref. [12, p. 450]. 1- Calculate the term Dlti for all iAI, lAL, and tAT as follows: l

Dlti ¼

cwrijlt dli þ gli þ bj dli vl þ mlt twrijlt dli þð1=2Þwhlj wsfijlt dli dli vl

2- Let L be a regular set of non-positive values of Dlti , arranged in ascending order. 3- Let e^ be the unoccupied capacity of the warehouse, with initial value ofwcje .

1354

H. Sadjady, H. Davoudpour / Computers & Operations Research 39 (2012) 1345–1354

4- If La |, select the retailer i, product l, and the mode t associated with the first Dlti in the regular set L, whose demand has not been satisfied yet (dli a 0), otherwise go to step 7. 5- If e^ Z dli vl , do the following sub-steps, otherwise go to step 6. a. Satisfy total demand of the selected retailer i for product l, from warehouse j with capacity level e, via mode t (Xijlt ¼ 1). b. Update the available capacity of the warehouse: e^ ¼ e^ dli vl c. Update demand of retailer i for the product l (dli ¼ 0), and back to step 4. 6- Assign the remaining available capacity of warehouse for satisfying the demand of retailer i for product l, via mode t (Xijlt ¼ ðe^ =dli vl Þ). 7- Calculate the objective function, using obtained values of Xijlt and stop. A similar solution procedure is applied for solving LR2gk subproblems. The only difference is in regarding the maximum capacity level of warehouses as their demand for each product. References [1] Ahuja RK, Magnanti TL, Orlin JB. Network flows: theory, algorithms, and applications. Englewood Cliffs, NJ: Prentice Hall; 1993 (pp. 508–637). [2] Amiri A. Designing a distribution network in a supply chain system: formulation and efficient solution procedure. European Journal of Operational Research 2006;171:567–76. ¨ - SS, Simpson NC, Vakharia AJ. Integrated production/distribution [3] Erenguc planning in supply chains: an invited review. European Journal of Operational Research 1999;115:219–36.

[4] Eskigun E, Uzsoy R, Preckel PV, Beaujon G, Krishnan S, Tew JD. Outbound supply chain network design with mode selection, lead times and capacitated vehicle distribution centers. European Journal of Operational Research 2005;165:182–206. [5] Fisher M. The Lagrangian relaxation method for solving integer programming problems. Management Science 1981;27:1–18. [6] Gavish B. On obtaining the ‘best’ multipliers for a Lagrangean relaxation for integer programming. Computers & Operations Research 1979;5:55–71. [7] Gourdin E, Labbe M, Laporte G. The uncapacitated facility location problem with client matching. Operations Research 2000;48:671–85. [8] Jayaraman V, Pirkul H. Planning and coordination of production and distribution facilities for multiple commodities. European Journal of Operational Research 2001;133:394–408. [9] Klose A, Drexl A. Facility location models for distribution system design. European Journal of Operational Research 2005;162:4–29. [10] Krarup J, Pruzan PM. The simple plant location problem: survey and synthesis. European Journal of Operational Research 1983;12:36–81. [11] Melo T, Nickel S, Saldanha da Gama F. Facility location and supply chain management—a review. European Journal of Operational Research 2009;196:401–12. [12] Nemhauser GL, Wolsey LA. Integer and combinatorial optimization. New York, NY: Wiley; 1988. [13] Owen SH, Daskin MS. Strategic facility location: a review. European Journal of Operational Research 1998;111:423–47. [14] Pirkul H, Jayaraman V. A multi-commodity, multi-plant, capacitated facility location problem: formulation and efficient heuristic solution. Computers & Operations Research 1998;25:869–78. [15] ReVelle CS, Eiselt HA, Daskin MS. A bibliography for some fundamental problem categories in discrete location science. European Journal of Operational Research 2008;184:817–48. [16] Siamitros C, Mitra G, Poojari C. Revisiting Lagrange relaxation for processing large-scale mixed integer programming problems. The Centre for the Analysis of Risk and Optimisation Modelling Applications (CARISMA), Brunel University, 2004. Available online at; /http://carisma.brunel.ac.uk/papers/ CTR-04-Christos.pdfS. [17] Wu LY, Zhang XS, Zhang JL. Capacitated facility location problem with general setup cost. Computers & Operations Research 2006;33:1226–41.