# Designing a Supply Chain Network under  the Risk of Disruptions

Authors: Armin Jabbarzadeh; Seyed Gholamreza Jalali Naini; Hamid Davoudpour; Nader Azad
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper studies a supply chain design problem with the risk of disruptions at facilities. At any point of time, the facilities are subject to various types of disruptions caused by natural disasters, man-made defections, and equipment breakdowns. We formulate the problem as a mixed-integer nonlinear program which maximizes the total profit for the whole system. The model simultaneously determines the number and location of facilities, the subset of customers to serve, the assignment of customers to facilities, and the cycle-order quantities at facilities. In order to obtain near-optimal solutions with reasonable computational requirements for large problem instances, two solution methods based on Lagrangian relaxation and genetic algorithm are developed. The effectiveness of the proposed solution approaches is shown using numerical experiments. The computational results, in addition, demonstrate that the benefits of considering disruptions in the supply chain design model can be significant.


---

## Body

## 1. Introduction

Under today’s highly competitive business environment, supply chain network design is a critical and difficult decision. Many of supply chain design decisions such as facility location are strategic in nature and very expensive to change. In particular, supply chain design involves both strategic decisions of facility location and tactical decisions of inventory. Traditional supply chain design models in the literature treat location and inventory decisions separately. However, ignoring interaction between long-term decisions of location and short-term decisions of inventory can lead to suboptimality [1–3]. Thus, integrated supply chain design models incorporating location and inventory decisions have emerged in recent years [4, 5].Daskin et al. [1] and Shen et al. [6] develop a basic integrated location inventory model that explicitly considers expected inventory costs when making facility location decisions. The model simultaneously determines facility location and demand assignment decisions in order to minimize the total cost including location, inventory, and shipment costs. This basic integrated location-inventory model is extended by researchers [3, 7–12]. For instance, Snyder et al. [7] propose a stochastic version of integrated location-inventory that handles uncertainty by describing discrete scenarios. Max Shen and Qi [3], also, add routing decisions to joint location-inventory model. That is, they examine an integrated model which determines location, inventory, and routing decisions. Ozsen et al. [8] consider an integrated location inventory when capacities of facilities are limited. Yao et al. [11] consider a joint location-inventory problem in which a company allows its retailers to be sourced by more than one warehouse. Shen [4] and Melo et al. [5] present comprehensive literature reviews on integrated supply chain design models.A limitation of the most existing studies on the integrated supply chain design is that they implicitly assume that facilities are always available and perfectly reliable. However, facilities in real world are always vulnerable to partial or complete disruptions. In fact, different factors, such as natural disasters, labor strikes, power outages, parts shortages, quality rejections, poor communications of customer requirements, transportation damages, and machine breakdowns, can lead to unreliable performance of facilities and huge losses. For instance, many companies like Intel, Wal-Mart, Ford, Isuzu Motors, and Suzuki had to stop production because of cut-off of electricity and water supply at their suppliers in 2008 [13]. Disruption at a Philips Semiconductor plant in 2001 caused shortages of cell phone components for Ericsson which resulted in losing a substantial portion of Ericsson’s market to its rival, Nokia [14]. Furthermore, smaller-scale disruptions at facilities take place much more frequently. For example, Wal-Mart’s Emergency Operations Center receives a call almost every day from a store with some sort of crisis [15]. These examples show that disruptions can significantly impact firm’s operations and highlight the need to account for disruptions during design of supply chain network.Facility location models with disruptions consideration have gained much attention recently [16–21]. Snyder and Daskin [16] present two facility location models based on P-median problem and incapacitated fixed location problem, in which facilities are disrupted with the same probabilities. In both models, the objective functions include terms of expected operational cost and expected failure cost to model a trade-off between these two costs. Their models rely on the assumption that the disruption probabilities of the facilities are equal. This assumption is relaxed by [17–19].Berman et al. [17] study a P-median problem where each facility fails independently with a certain probability. The objective minimizes the expected P-median cost plus penalty cost for not being able to serve a customer. Due to intractability of the model, the authors suggest heuristic to solve the problem. Lim et al. [18] also focus on a facility location problem in presence of disruptions with the option of hardening selected facilities. Their model is formulated based on the assumption that disruption probabilities are independent. Cui et al. [19] propose a continuum approximation model to study incapacitated fixed charge location problem in which facilities are disrupted with site-dependent probabilities. Their model seeks to minimize initial setup costs and expected transportation costs in normal and failure scenarios. Snyder et al. [22] and Snyder and Daskin [23] review the broad range of facility location models under presence of disruptions. Among the above works in the area of facility location with disruptions, no model considers inventory costs.Aryanezhad et al. [24] and Chen et al. [25] propose location-inventory models when facilities are subject to disruptions. They propose integer programming models that minimize the sum of facility construction costs, expected inventory costs, and expected customer costs under normal and failure cases. Their models are based on the restrictive assumption that facilities fail independently with an equal probability. Also, they do not consider partial disruptions in their models. Qi and Shen [26] concentrate on integrated supply chain design under yield uncertainty; however, they do not consider disruptions at facilities. Qi et al. [27] formulate a joint location-inventory model in which facilities can be disrupted. They assume that all inventory at a facility is destroyed if the facility is disrupted. In fact, none of the above models consider both partial and complete disruptions at facilities.This paper investigates an integrated supply chain design problem with multiple distribution centers subject to various types of disruptions. The problem is formulated as a nonlinear integer programming model which determines the location of distribution centers and the assignments of customers to distribution centers. The proposed model in this study builds upon recent developments of integrated supply chain design models that simultaneously consider location, inventory, and shipment decisions in the same model. In order to obtain near-optimal solutions with reasonable computational requirements, two solution methods based on Lagrangian relaxation and genetic algorithm are developed. Numerical experiments are conducted to test the performance of solution approaches and draw managerial insights on the benefits of considering disruptions during design of supply chain networks.The proposed model in this paper differs from the earlier works in literature of integrated supply chain design. Unlike most of joint location-inventory models in the literature, the model takes into account the different disruptions scenarios at facilities during decision making. Also, the profit-maximizing model in this study does not restrictively assume that every potential demand has to be satisfied. Typically, cost-minimizing models in the literature require every potential demand must be met. However, for a profit-maximizing business, it may not always be optimal to serve all potential customers, especially if the additional cost is higher than the additional revenue associated with satisfying some demands [26, 28]. This article is also different from the literature on facility location problems with disruptions. First, the model does not ignore nonlinear inventory costs. In addition, it takes into account the possibility of both partial and complete facility disruptions in the same model.The remainder of the paper is organized as follows. In Section2, the problem is described and the formulation model for the problem is presented. Section 3 proposes two solution approaches based on Lagrangian relaxation and genetic algorithm. In Section 4, numerical experiments are conducted to evaluate the performance of solution methods and to study the benefits of considering disruptions in the supply chain design model. Finally, Section 5 concludes the paper along with directions for future research.


## 2. Model Formulation

### 2.1. Problem Description and Assumptions

This paper addresses a supply chain design network problem under random facility disruptions. As typically considered in the supply chain design literature, a three-tiered supply chain consisting of a single supplier, distribution centers (DCs), and customers is studied. The supplier ships one type of product to a set of customers in order to satisfy their demands. The supply chain is flexible in determining which customers to serve. In other words, if the cost of serving some customers is prohibitive, they are not served at all. DCs function as the direct intermediary between the plant and customers for shipment of the product. That is, DCs combine the orders from different customers and then order to the supplier. Similar to [17–19, 21, 25, 27], we assume that the lead time for order delivery is negligible and the demand rate is fixed. Also, the capacities of DCs are assumed to be infinite as assumed in [1–3, 6, 7, 17–19, 26–28].The key problem is that DCs may face different amounts of disruptions from time to time. If a customer is assigned to a DC but the DC is disrupted and cannot meet all of the customer’s demands, the unmet demands are lost. In this case, the system incurs a lost-sales cost for each unit of lost demand. The percentage of disruptions at each DC is uncertain. In other words, each DC may experience different amounts of disruptions with different probabilities. DCs are not required to be the same in the problem; as a result, the probabilities of disruptions for each DC can be different from those of others. To formulate the disruptions at the DCs, a scenario-based modeling approach is used, in which each scenario specifies the percentage of disruptions for each DC. For instance, if a DC experiences complete disruptions in a scenario, the percentage of the disruptions is considered hundred for that DC. Note that the scenario-based modeling framework is flexible enough to consider both complete and partial disruptions. Also, it allows us to model the complex situation in which the probability failures of DCs are dependent.The problem lies in simultaneously determining (1) how many DCs are opened, and where to locate them; (2) which subset of customers is served; (3) which DCs are assigned to which customers; (4) how much and how often to order at each DC. Figure1 demonstrates the structure of the paper for solving this problem. The problem is formulated as a mixed-integer nonlinear program which maximizes the expected total profit. That is, the objective is to maximize the expected difference between total revenue and total cost. The total cost includes three main components: (1) the fixed cost to locate DCs, (2) the working inventory cost (including order costs, shipment costs from supplier to DCs, holding costs, and lost-sales costs) at the located DCs, and (3) shipment cost from located DCs to customers.Figure 1
Strategy of the present article for determining the supply chain design decisions.


### 2.2. Notation

To develop the integrated model, the following notations are used throughout the paper. Additional notations will be given out when required.Parameters
(i)
I: set of customers indexed by i;(ii)
J: set of candidate DC locations indexed by j;(iii)
S: set of disruption scenarios indexed by s;(iv)
λi: demand rate at customer i, for each i∈I;(v)
Ri: selling price at customer i, per unit of demand for each i∈I;(vi)
fj: fixed cost of locating a DC at j, for each j∈J;(vii)
aj: fixed cost of placing an order at j, for each j∈J;(viii)
bj: fixed cost per shipment from the supplier to DC at j, for each j∈J;(ix)
cj: per-unit shipment cost from the supplier to DC at j, for each j∈J;(x)
h: inventory holding cost per unit of product per year;(xi)
dij: per-unit cost to ship from distribution center j to customer i, for each i∈I and for each j∈J;(xii)
β: weight factor associated with the shipment cost;(xiii)
θ: weight factor associated with the inventory cost;(xiv)
rsj: percentage of supply at distribution center j which is disrupted in scenario s, for each s∈S and for each j∈J;(xv)
πj: penalty cost for not being able to satisfy customers’ demands due to disruptions at j, per unit of demand, for each j∈J (it may be a lost-sales cost during disruptions, or the cost of filling the supply by purchasing product from a competitor on an emergency basis to compensate the disrupted supply at j);(xvi)
qs: probability that scenario s occurs, for each s∈S.Decision Variables
(i)
Xj=1, if j is selected as a DC location, and 0, otherwise, for each j∈J;(ii)
Yij=1, if customer i is assigned to a DC based at j, and 0 otherwise, for each i∈I and j∈J;(iii)
Zi=1, if customer i is not selected to be served, and 0 otherwise, for each i∈I.The total demand which is assigned to a distribution center is unknown in advance. However, the assigned demand to each distribution centerj can be obtained by the assignment decisions (Y) as follows:
(2.1)Dj=∑i∈IλiYij,
where Dj denotes the total demand that is assigned to the DC at j.


### 2.3. Working Inventory Cost

This subsection formulates the working inventory cost including costs of ordering, shipment from supplier to DCs, holding, and lost-sales. For the moment, letDj denote the unknown total demand that is assigned to the DC at j (it is obvious that Dj=∑i∈IλiYij). Also, let n be the unknown number of orders per year.Then, the expected shipment size per shipment from the supplier to DC atj is equal to Dj/n and the working inventory cost at distribution center j can be obtained by
(2.2)ajn+β(bj+cjDjn)n+∑sqs((1-rsj)Dj2nθh)+∑sqs(rsjDj2nπj).The first term of (2.2) is the fixed cost of placing n orders. The second term indicates the cost of shipping n orders of size Dj/n, assuming the shipment cost from the supplier to the distribution center j has a fixed cost bj and volume-dependent cost cj. The third term represents the cost of holding average of ∑sqs((1-rsj)Dj/2n) units. To explain this amount of average inventory, we note that if the DC at j was not subject to any disruptions, the average units to hold would be Dj/2n [29]. Also, recall that if scenario s was occurred with certainty, rsj% of supply at distribution center j would be disrupted and the average inventory would be (1-rsj)Dj/2n. Thus, knowing that each scenario s is occurred with probability of qs, the average units to hold will be ∑sqs((1-rsj)Dj/2n). Also, the average disrupted supply at distribution center j will be ∑sqs(rsjDj/2n). Therefore, the last term in (2.2) indicates the lost-sales cost due to disruptions at distribution center j, where the penalty cost for losing a unit of supply due to disruptions is denoted by πj. In order to determine the optimal number of orders, we take derivative of (2.2) with respect to n and set the derivative to zero:
(2.3)aj+βbj-∑sqs((1-rsj)Dj(2n)2θh)-∑sqs(rsjDj(2n)2πj)=0.Solving (2.3) for n, we obtain n=∑sqs(rsjπj+(1-rsj)θh)Dj/2(aj+βbj). Plugging this into (2.2), working inventory cost at distribution center j can be calculated as follows:
(2.4)2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)Dj+βcjDj.SinceDj=∑i∈IλiYij, (2.4) can be rewritten as follows:
(2.5)2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij.


### 2.4. Integrated Model

The problem is formulated as follows:(2.6)Max⁡∑i∈IRiλi(1-Zi)-(∑j∈JfjXj)-β∑j∈J∑i∈IdijλiYij-∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)
subject to
(2.7)∑j∈JYij+Zi=1∀i∈I,(2.8)Yij≤Xj∀i∈I,∀j∈J,(2.9)Xj∈{0,1}∀j∈J,(2.10)Yij∈{0,1}∀i∈I,∀j∈J.The objective function (2.6) is composed of four components. The first component indicates the sales revenue which is gained by serving the customers. The second component represents the fixed cost of locating DCs. The third component indicates the expected shipment cost from the DCs to customers. Finally, the fourth component represents the working inventory cost. Constraints (2.7) require each customer to be assigned to exactly one DC, or not to be served at all. Constraints (2.8) state that customers can only be assigned to candidate sites that are selected as DCs. Constraints (2.9) and (2.10) are binary constraints.


## 2.1. Problem Description and Assumptions

This paper addresses a supply chain design network problem under random facility disruptions. As typically considered in the supply chain design literature, a three-tiered supply chain consisting of a single supplier, distribution centers (DCs), and customers is studied. The supplier ships one type of product to a set of customers in order to satisfy their demands. The supply chain is flexible in determining which customers to serve. In other words, if the cost of serving some customers is prohibitive, they are not served at all. DCs function as the direct intermediary between the plant and customers for shipment of the product. That is, DCs combine the orders from different customers and then order to the supplier. Similar to [17–19, 21, 25, 27], we assume that the lead time for order delivery is negligible and the demand rate is fixed. Also, the capacities of DCs are assumed to be infinite as assumed in [1–3, 6, 7, 17–19, 26–28].The key problem is that DCs may face different amounts of disruptions from time to time. If a customer is assigned to a DC but the DC is disrupted and cannot meet all of the customer’s demands, the unmet demands are lost. In this case, the system incurs a lost-sales cost for each unit of lost demand. The percentage of disruptions at each DC is uncertain. In other words, each DC may experience different amounts of disruptions with different probabilities. DCs are not required to be the same in the problem; as a result, the probabilities of disruptions for each DC can be different from those of others. To formulate the disruptions at the DCs, a scenario-based modeling approach is used, in which each scenario specifies the percentage of disruptions for each DC. For instance, if a DC experiences complete disruptions in a scenario, the percentage of the disruptions is considered hundred for that DC. Note that the scenario-based modeling framework is flexible enough to consider both complete and partial disruptions. Also, it allows us to model the complex situation in which the probability failures of DCs are dependent.The problem lies in simultaneously determining (1) how many DCs are opened, and where to locate them; (2) which subset of customers is served; (3) which DCs are assigned to which customers; (4) how much and how often to order at each DC. Figure1 demonstrates the structure of the paper for solving this problem. The problem is formulated as a mixed-integer nonlinear program which maximizes the expected total profit. That is, the objective is to maximize the expected difference between total revenue and total cost. The total cost includes three main components: (1) the fixed cost to locate DCs, (2) the working inventory cost (including order costs, shipment costs from supplier to DCs, holding costs, and lost-sales costs) at the located DCs, and (3) shipment cost from located DCs to customers.Figure 1
Strategy of the present article for determining the supply chain design decisions.


## 2.2. Notation

To develop the integrated model, the following notations are used throughout the paper. Additional notations will be given out when required.Parameters
(i)
I: set of customers indexed by i;(ii)
J: set of candidate DC locations indexed by j;(iii)
S: set of disruption scenarios indexed by s;(iv)
λi: demand rate at customer i, for each i∈I;(v)
Ri: selling price at customer i, per unit of demand for each i∈I;(vi)
fj: fixed cost of locating a DC at j, for each j∈J;(vii)
aj: fixed cost of placing an order at j, for each j∈J;(viii)
bj: fixed cost per shipment from the supplier to DC at j, for each j∈J;(ix)
cj: per-unit shipment cost from the supplier to DC at j, for each j∈J;(x)
h: inventory holding cost per unit of product per year;(xi)
dij: per-unit cost to ship from distribution center j to customer i, for each i∈I and for each j∈J;(xii)
β: weight factor associated with the shipment cost;(xiii)
θ: weight factor associated with the inventory cost;(xiv)
rsj: percentage of supply at distribution center j which is disrupted in scenario s, for each s∈S and for each j∈J;(xv)
πj: penalty cost for not being able to satisfy customers’ demands due to disruptions at j, per unit of demand, for each j∈J (it may be a lost-sales cost during disruptions, or the cost of filling the supply by purchasing product from a competitor on an emergency basis to compensate the disrupted supply at j);(xvi)
qs: probability that scenario s occurs, for each s∈S.Decision Variables
(i)
Xj=1, if j is selected as a DC location, and 0, otherwise, for each j∈J;(ii)
Yij=1, if customer i is assigned to a DC based at j, and 0 otherwise, for each i∈I and j∈J;(iii)
Zi=1, if customer i is not selected to be served, and 0 otherwise, for each i∈I.The total demand which is assigned to a distribution center is unknown in advance. However, the assigned demand to each distribution centerj can be obtained by the assignment decisions (Y) as follows:
(2.1)Dj=∑i∈IλiYij,
where Dj denotes the total demand that is assigned to the DC at j.


## 2.3. Working Inventory Cost

This subsection formulates the working inventory cost including costs of ordering, shipment from supplier to DCs, holding, and lost-sales. For the moment, letDj denote the unknown total demand that is assigned to the DC at j (it is obvious that Dj=∑i∈IλiYij). Also, let n be the unknown number of orders per year.Then, the expected shipment size per shipment from the supplier to DC atj is equal to Dj/n and the working inventory cost at distribution center j can be obtained by
(2.2)ajn+β(bj+cjDjn)n+∑sqs((1-rsj)Dj2nθh)+∑sqs(rsjDj2nπj).The first term of (2.2) is the fixed cost of placing n orders. The second term indicates the cost of shipping n orders of size Dj/n, assuming the shipment cost from the supplier to the distribution center j has a fixed cost bj and volume-dependent cost cj. The third term represents the cost of holding average of ∑sqs((1-rsj)Dj/2n) units. To explain this amount of average inventory, we note that if the DC at j was not subject to any disruptions, the average units to hold would be Dj/2n [29]. Also, recall that if scenario s was occurred with certainty, rsj% of supply at distribution center j would be disrupted and the average inventory would be (1-rsj)Dj/2n. Thus, knowing that each scenario s is occurred with probability of qs, the average units to hold will be ∑sqs((1-rsj)Dj/2n). Also, the average disrupted supply at distribution center j will be ∑sqs(rsjDj/2n). Therefore, the last term in (2.2) indicates the lost-sales cost due to disruptions at distribution center j, where the penalty cost for losing a unit of supply due to disruptions is denoted by πj. In order to determine the optimal number of orders, we take derivative of (2.2) with respect to n and set the derivative to zero:
(2.3)aj+βbj-∑sqs((1-rsj)Dj(2n)2θh)-∑sqs(rsjDj(2n)2πj)=0.Solving (2.3) for n, we obtain n=∑sqs(rsjπj+(1-rsj)θh)Dj/2(aj+βbj). Plugging this into (2.2), working inventory cost at distribution center j can be calculated as follows:
(2.4)2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)Dj+βcjDj.SinceDj=∑i∈IλiYij, (2.4) can be rewritten as follows:
(2.5)2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij.


## 2.4. Integrated Model

The problem is formulated as follows:(2.6)Max⁡∑i∈IRiλi(1-Zi)-(∑j∈JfjXj)-β∑j∈J∑i∈IdijλiYij-∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)
subject to
(2.7)∑j∈JYij+Zi=1∀i∈I,(2.8)Yij≤Xj∀i∈I,∀j∈J,(2.9)Xj∈{0,1}∀j∈J,(2.10)Yij∈{0,1}∀i∈I,∀j∈J.The objective function (2.6) is composed of four components. The first component indicates the sales revenue which is gained by serving the customers. The second component represents the fixed cost of locating DCs. The third component indicates the expected shipment cost from the DCs to customers. Finally, the fourth component represents the working inventory cost. Constraints (2.7) require each customer to be assigned to exactly one DC, or not to be served at all. Constraints (2.8) state that customers can only be assigned to candidate sites that are selected as DCs. Constraints (2.9) and (2.10) are binary constraints.


## 3. Solution Methods

In order to solve the model formulated in Section2, two solution methods are developed. The first solution approach is based on Lagrangian relaxation which has been largely applied in the literature of integrated supply chain design models [1, 3, 7–9, 18, 19, 25, 27]. Also, a solution method based on genetic algorithm is developed, to obtain near-optimal solutions in relatively short time. Genetic algorithm has been successfully used to solve various supply chain design problems and has proven to be a very effective heuristic procedure to solve these problems, particularly problems of large scale [21, 30–34]. It can overcome computational complexity caused by the nonlinear and stochastic objective functions to solve the model [34, 35].


### 3.1. Lagrangian Relaxation

Here, a Lagrangian relaxation is developed to solve the model due to its proven effectiveness for solving integrated supply chain design models [1, 3, 7–9, 18, 19, 25, 27]. Lagrangian relaxation approach provides both upper and lower bounds on the optimal value of the objective function. In other words, it allows the decision maker to know how far from the optimality the best found feasible solution is [36]. Figure 2 outlines the steps of the developed Lagrangian relaxation. In the following subsections, the steps of the proposed Lagrangian relaxation are explained with more details.Figure 2
Lagrangian relaxation procedure.


#### 3.1.1. Finding a Lower Bound

First, the model formulated in Section2 is converted into a standard model for which an effective Lagrangian relaxation approach exists in the literature [1, 27]. Replacing Zi in (2.6) with 1-∑j∈JYij according to (2.7), the original model formulated in Section 2 can be written as follows:(3.1)Max⁡∑j∈J∑i∈IRiλiYij-(∑j∈JfjXj)-β∑j∈J∑i∈IdijλiYij-∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)
subject to
(3.2)∑j∈JYij≤1∀i∈I,(3.3)Yij≤Xj∀i∈I,∀j∈J,(3.4)Xj∈{0,1}∀j∈J,(3.5)Yij∈{0,1}∀i∈I,∀j∈J.By changing the sign of (3.1) and rearranging the terms, the problem can be converted into a minimizing model with the following objective:(3.6)Min⁡∑j∈JfjXj+β∑j∈J∑i∈IdijλiYij+∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)-∑j∈J∑i∈IRiλiYij=∑j∈J(fjXj+∑i∈I(βdij+βcj-Ri)λiYij+2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij)=∑j∈J(fjXj+∑i∈IuiYij+∑i∈IviYij),
where(3.7)uij=(βdij+βcj-Ri)λi,vij=2(aj+βbj)λi∑sqs(rsjπj+(1-rsj)θh).Obviously the upper and lower bounds for (3.6) can easily be used as the lower and upper bounds for (3.5) by changing their signs.Relaxing constraints (3.2) with Lagrange multipliers, ωi, leads to the following Lagrangian dual problem:(3.8)Max⁡ωMin⁡X,Y∑j∈J{fjXj+∑i∈IuijYij+∑i∈IvijYij}+∑i∈Iωi(1-∑j∈JYij)=∑j∈J{fjXj+∑i∈I(uij-ωi)Yij+∑i∈IvijYij}+∑i∈Iωi
subject to(3.9)Yij≤Xj∀i∈I,∀j∈J,(3.10)Xj∈{0,1}∀j∈J,(3.11)Yij∈{0,1}∀i∈I,∀j∈J.For given values of the Lagrange multipliers,ωi, the objective is to minimize (3.8) over the decision variables Xj and Yij. This problem is decomposed by j; as a result, we need to solve the following sub-problem for each candidate location j∈J:(3.12)SPj:Ṽj=Min⁡∑i∈I(uij-ωi)Pi+∑i∈IvijPi
subject to(3.13)Pi∈{0,1}∀i∈I.In (3.12)-(3.13), the assignment variables Yij have been replaced by Pi to simplify the notation, as SPj is specific to distribution j. Note that the facility location cost, fj, is not included in Ṽj and will be added later. Subproblem SPj can be solved using the exact algorithm developed by Shen et al. [6]. Modified to our problem, their algorithm is given as follows. (1)
DefineI0={i:(uij-ωi)<0andvij=0} and I-={i:(uij-ωi)<0andvij>0}.(2)
Sort the elements ofI- in increasing order of (uij-ωi)/vij and represent the elements by 1-,2-,…,n-, respectively, where n=|I-|.(3)
Find the value ofm that minimizes(3.14)∑i∈I0(uij-ωi)Pi+∑i∈I0vijPi+∑i=1,i∈I-m(uij-ωi)Pi+∑i=1,i∈I-mvijPi.(4)
The optimal solution to subproblemSPj is obtained by Pi=1 for i∈I0, P1-=P2-=⋯=Pm-=1 for i∈I- and Pi=0 for all other i∈I.After solving subproblemSPj for each j,fj is added to the optimal objective value of Ṽj. If Ṽj+fj<0, then we select the candidate distribution center j and set Xj=1; otherwise, we set Xj=0. For each selected distribution center j (those for which Xj=1), the assignment variables Yij are the same as the optimal Pi values in subproblem SPj; for each unselected distribution center j (those for which Xj=0), Yij=0,∀i∈I.Having solved the Lagrangian problem, the optimal Lagrange multipliers are found using a standard subgradient optimization procedure [37, 38]. The optimal objective value of the Lagrangian dual problem (3.8) provides a lower bound on the optimal objective value of (3.6).


#### 3.1.2. Finding an Upper Bound

The lower bound solution obtained in each iteration of Lagrangian relaxation can violate constraints (3.2). In other words, some customers can be assigned to more than one DC in the obtained lower bound solution. Thus, in each iteration of Lagrangian relaxation, the obtained lower bound solution is converted into feasible upper bound solution using the procedure adapted from Shen et al. [6], as follows.(1)
If a customer is assigned to more than one DC, we find the DC which assigning it to the customer leads to the least objective (3.6). If the resulted objective (3.6) is less than the case that the customer is not served at all, we assign the customer to that DC. Otherwise, the customer is not assigned to any DC and is not served at all.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.If the obtained feasible solution results in a less value for (3.6) than the best known upper bound, it will be taken as the new upper bound solution. Also, it will be improved using customer reassignment algorithm.


#### 3.1.3. Customer Reassignment

Each obtained upper bound is improved using following algorithm.(1)
For each customer, we check whether the objective value (3.6) is improved if the customer is assigned, instead to another located DC, or if it is not served at all. The best improving swap is performed.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.


#### 3.1.4. Variable Fixing

At the end of Lagrangian procedure, a variable fixing technique is employed. This method uses two following rules whose detailed proofs of the validity are presented by Shu et al. [2].(1)
Let LB and UB denote the current lower and upper bounds on the solution, respectively. If no DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=0 in the optimal solution to (3.8) at some iteration), and if LB+fj+Ṽj>UB, then no DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=0.(2)
If a DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=1 in the optimal solution to (3.8) at some iteration), and if LB-(fj+Ṽj)>UB, then a DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=1.


### 3.2. Genetic Algorithm

In order to find near-optimal solutions for the model in relatively short time, here a solution approach based on genetic algorithm (GA) is developed. GA is a stochastic search and heuristic optimization technique based on the mechanism of natural genetics which has proven to be a very effective heuristic procedure to solve supply chain design problems, particularly problems of large scale [21, 30–34].GA starts with an initial set of random solution called population. Each solution in the population is called chromosome and each component of chromosome is designated by gene. The chromosomes evolve through successive iterations, called generations. During each generation, the chromosomes are evaluated, using some measures of fitness. To create next generation, new chromosomes (called offspring) are formed by crossover or mutation operators. Crossover operator combines two chromosomes from current generation, while mutation operator modifies a chromosome to form offspring.A new generation is created by selecting some of current chromosomes (called parents) and offspring based on the fitness values. Also, some chromosomes are rejected so as to keep the population size constant. Fitter chromosomes have higher probabilities of being selected. After several generations, the algorithms converge to the best chromosome, which may represent the optimum or suboptimal solution to the problem [33, 39, 40]. Figure 3 summarizes the steps of the proposed GA. In the following subsections, the developed GA is explained with more details.Figure 3
Genetic algorithm procedure.


#### 3.2.1. Chromosome Representation

In this GA-based approach, each chromosome is indicated as a single-dimensional array. Letn be the number of candidate DCs and m be the number of customers. Then, each chromosome C can be indicated by: C=(Xj,Yi)=(X1,X2,…,Xn,Y1,Y2,…,Ym), where Xj correspond to the location genes and Yi correspond to the assignment genes. Location genes indicate where DCs are located and assignment genes represent how customers are assigned to the located DCs, respectively. In particular, if Xj=1, it means that candidate site j is selected as a DC location, while if Xj=0, candidate location j is not chosen as a DC site. Also, Yi=j represents that customer i is assigned to distribution center j. If customer i is not served at all, the corresponding assignment gene takes the value of 0 (i.e., Yi=0). For example, in Figure 4, distribution centers are located at 2 and 4. It follows that customers 1 and 4 are assigned to the DC at 4 and customer 2 is allocated to the DC at 2. Also, customer 3 is not served at all.Figure 4
Chromosome structure.


#### 3.2.2. Chromosomes Fitness

The rank-based evaluation function is defined as the objective function (2.6) for the chromosomes. In fact, the objective (2.6) is calculated for each of the chromosomes. Obviously, chromosomes with better values for objective (2.6) will have the better rank.


#### 3.2.3. Crossover Operator

Crossover operator generates offspring by merging parent chromosomes. LetCk denote the chromosomes of the population for k=1,2,…,pop-size. In order to determine which of the chromosomes are selected as parents for crossover operation, the following procedure is repeated from k=1 topop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PC, where the parameter PC is the probability of crossover. Then, randomly we group the selected parents C1′,C2′,C3′,… to the pairs (C1′,C2′),(C3′,C4′),…. Without loss of generality, let us explain the crossover operator on each pair by (C1′,C2′).Crossover operator assigns each customeri in offspring chromosome either to the DC which is allocated to customer i in parent chromosome C1′, or to the DC which is assigned to customer i in parent chromosome C2′. This occurs randomly and with probability of 0.5. If a customer is allocated to an unselected candidate DC site, a DC is located in that candidate location. A sample of crossover operator is shown in Figure 5.Figure 5
Sample of crossover.


#### 3.2.4. Mutation Operator

Mutation operator modifies a chromosome to form an offspring. In order to decide which of chromosomesCk undergo mutation, the following practice is repeated from k=1 toPop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PM, where the parameter PM is the probability of mutation. A selected chromosome is modified by one of the two following types of mutation for several times. Each type of mutation is occurred with probability 0.5.The first type of mutation generates offspring by modifying the assignment genes of parent chromosome. That is, in the first type of mutation, two located DCs are selected randomly; lets and t denote them. Then, if any customer in the parent chromosome is assigned to s, that customer will be assigned to t, and if any customer is assigned to t, it will be allocated to s. The second type of mutation modifies location genes of parent chromosome to form an offspring. In other words, the second type of mutation randomly selects a location in which no DC is located; let t denotes it. Next, a DC is selected randomly from the located DCs and is named s. This type of mutation closes distribution center s and locates a DC at t instead. Then, all the customers assigned to distribution center s are allocated to distribution center t. The samples of mutation type 1 and mutation type 2 are illustrated in Figures 6 and 7, respectively.Figure 6
Sample of mutation type 1.Figure 7
Sample of mutation type 2.


## 3.1. Lagrangian Relaxation

Here, a Lagrangian relaxation is developed to solve the model due to its proven effectiveness for solving integrated supply chain design models [1, 3, 7–9, 18, 19, 25, 27]. Lagrangian relaxation approach provides both upper and lower bounds on the optimal value of the objective function. In other words, it allows the decision maker to know how far from the optimality the best found feasible solution is [36]. Figure 2 outlines the steps of the developed Lagrangian relaxation. In the following subsections, the steps of the proposed Lagrangian relaxation are explained with more details.Figure 2
Lagrangian relaxation procedure.


### 3.1.1. Finding a Lower Bound

First, the model formulated in Section2 is converted into a standard model for which an effective Lagrangian relaxation approach exists in the literature [1, 27]. Replacing Zi in (2.6) with 1-∑j∈JYij according to (2.7), the original model formulated in Section 2 can be written as follows:(3.1)Max⁡∑j∈J∑i∈IRiλiYij-(∑j∈JfjXj)-β∑j∈J∑i∈IdijλiYij-∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)
subject to
(3.2)∑j∈JYij≤1∀i∈I,(3.3)Yij≤Xj∀i∈I,∀j∈J,(3.4)Xj∈{0,1}∀j∈J,(3.5)Yij∈{0,1}∀i∈I,∀j∈J.By changing the sign of (3.1) and rearranging the terms, the problem can be converted into a minimizing model with the following objective:(3.6)Min⁡∑j∈JfjXj+β∑j∈J∑i∈IdijλiYij+∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)-∑j∈J∑i∈IRiλiYij=∑j∈J(fjXj+∑i∈I(βdij+βcj-Ri)λiYij+2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij)=∑j∈J(fjXj+∑i∈IuiYij+∑i∈IviYij),
where(3.7)uij=(βdij+βcj-Ri)λi,vij=2(aj+βbj)λi∑sqs(rsjπj+(1-rsj)θh).Obviously the upper and lower bounds for (3.6) can easily be used as the lower and upper bounds for (3.5) by changing their signs.Relaxing constraints (3.2) with Lagrange multipliers, ωi, leads to the following Lagrangian dual problem:(3.8)Max⁡ωMin⁡X,Y∑j∈J{fjXj+∑i∈IuijYij+∑i∈IvijYij}+∑i∈Iωi(1-∑j∈JYij)=∑j∈J{fjXj+∑i∈I(uij-ωi)Yij+∑i∈IvijYij}+∑i∈Iωi
subject to(3.9)Yij≤Xj∀i∈I,∀j∈J,(3.10)Xj∈{0,1}∀j∈J,(3.11)Yij∈{0,1}∀i∈I,∀j∈J.For given values of the Lagrange multipliers,ωi, the objective is to minimize (3.8) over the decision variables Xj and Yij. This problem is decomposed by j; as a result, we need to solve the following sub-problem for each candidate location j∈J:(3.12)SPj:Ṽj=Min⁡∑i∈I(uij-ωi)Pi+∑i∈IvijPi
subject to(3.13)Pi∈{0,1}∀i∈I.In (3.12)-(3.13), the assignment variables Yij have been replaced by Pi to simplify the notation, as SPj is specific to distribution j. Note that the facility location cost, fj, is not included in Ṽj and will be added later. Subproblem SPj can be solved using the exact algorithm developed by Shen et al. [6]. Modified to our problem, their algorithm is given as follows. (1)
DefineI0={i:(uij-ωi)<0andvij=0} and I-={i:(uij-ωi)<0andvij>0}.(2)
Sort the elements ofI- in increasing order of (uij-ωi)/vij and represent the elements by 1-,2-,…,n-, respectively, where n=|I-|.(3)
Find the value ofm that minimizes(3.14)∑i∈I0(uij-ωi)Pi+∑i∈I0vijPi+∑i=1,i∈I-m(uij-ωi)Pi+∑i=1,i∈I-mvijPi.(4)
The optimal solution to subproblemSPj is obtained by Pi=1 for i∈I0, P1-=P2-=⋯=Pm-=1 for i∈I- and Pi=0 for all other i∈I.After solving subproblemSPj for each j,fj is added to the optimal objective value of Ṽj. If Ṽj+fj<0, then we select the candidate distribution center j and set Xj=1; otherwise, we set Xj=0. For each selected distribution center j (those for which Xj=1), the assignment variables Yij are the same as the optimal Pi values in subproblem SPj; for each unselected distribution center j (those for which Xj=0), Yij=0,∀i∈I.Having solved the Lagrangian problem, the optimal Lagrange multipliers are found using a standard subgradient optimization procedure [37, 38]. The optimal objective value of the Lagrangian dual problem (3.8) provides a lower bound on the optimal objective value of (3.6).


### 3.1.2. Finding an Upper Bound

The lower bound solution obtained in each iteration of Lagrangian relaxation can violate constraints (3.2). In other words, some customers can be assigned to more than one DC in the obtained lower bound solution. Thus, in each iteration of Lagrangian relaxation, the obtained lower bound solution is converted into feasible upper bound solution using the procedure adapted from Shen et al. [6], as follows.(1)
If a customer is assigned to more than one DC, we find the DC which assigning it to the customer leads to the least objective (3.6). If the resulted objective (3.6) is less than the case that the customer is not served at all, we assign the customer to that DC. Otherwise, the customer is not assigned to any DC and is not served at all.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.If the obtained feasible solution results in a less value for (3.6) than the best known upper bound, it will be taken as the new upper bound solution. Also, it will be improved using customer reassignment algorithm.


### 3.1.3. Customer Reassignment

Each obtained upper bound is improved using following algorithm.(1)
For each customer, we check whether the objective value (3.6) is improved if the customer is assigned, instead to another located DC, or if it is not served at all. The best improving swap is performed.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.


### 3.1.4. Variable Fixing

At the end of Lagrangian procedure, a variable fixing technique is employed. This method uses two following rules whose detailed proofs of the validity are presented by Shu et al. [2].(1)
Let LB and UB denote the current lower and upper bounds on the solution, respectively. If no DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=0 in the optimal solution to (3.8) at some iteration), and if LB+fj+Ṽj>UB, then no DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=0.(2)
If a DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=1 in the optimal solution to (3.8) at some iteration), and if LB-(fj+Ṽj)>UB, then a DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=1.


## 3.1.1. Finding a Lower Bound

First, the model formulated in Section2 is converted into a standard model for which an effective Lagrangian relaxation approach exists in the literature [1, 27]. Replacing Zi in (2.6) with 1-∑j∈JYij according to (2.7), the original model formulated in Section 2 can be written as follows:(3.1)Max⁡∑j∈J∑i∈IRiλiYij-(∑j∈JfjXj)-β∑j∈J∑i∈IdijλiYij-∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)
subject to
(3.2)∑j∈JYij≤1∀i∈I,(3.3)Yij≤Xj∀i∈I,∀j∈J,(3.4)Xj∈{0,1}∀j∈J,(3.5)Yij∈{0,1}∀i∈I,∀j∈J.By changing the sign of (3.1) and rearranging the terms, the problem can be converted into a minimizing model with the following objective:(3.6)Min⁡∑j∈JfjXj+β∑j∈J∑i∈IdijλiYij+∑j∈J(2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij+βcj∑i∈IλiYij)-∑j∈J∑i∈IRiλiYij=∑j∈J(fjXj+∑i∈I(βdij+βcj-Ri)λiYij+2(aj+βbj)∑sqs(rsjπj+(1-rsj)θh)∑i∈IλiYij)=∑j∈J(fjXj+∑i∈IuiYij+∑i∈IviYij),
where(3.7)uij=(βdij+βcj-Ri)λi,vij=2(aj+βbj)λi∑sqs(rsjπj+(1-rsj)θh).Obviously the upper and lower bounds for (3.6) can easily be used as the lower and upper bounds for (3.5) by changing their signs.Relaxing constraints (3.2) with Lagrange multipliers, ωi, leads to the following Lagrangian dual problem:(3.8)Max⁡ωMin⁡X,Y∑j∈J{fjXj+∑i∈IuijYij+∑i∈IvijYij}+∑i∈Iωi(1-∑j∈JYij)=∑j∈J{fjXj+∑i∈I(uij-ωi)Yij+∑i∈IvijYij}+∑i∈Iωi
subject to(3.9)Yij≤Xj∀i∈I,∀j∈J,(3.10)Xj∈{0,1}∀j∈J,(3.11)Yij∈{0,1}∀i∈I,∀j∈J.For given values of the Lagrange multipliers,ωi, the objective is to minimize (3.8) over the decision variables Xj and Yij. This problem is decomposed by j; as a result, we need to solve the following sub-problem for each candidate location j∈J:(3.12)SPj:Ṽj=Min⁡∑i∈I(uij-ωi)Pi+∑i∈IvijPi
subject to(3.13)Pi∈{0,1}∀i∈I.In (3.12)-(3.13), the assignment variables Yij have been replaced by Pi to simplify the notation, as SPj is specific to distribution j. Note that the facility location cost, fj, is not included in Ṽj and will be added later. Subproblem SPj can be solved using the exact algorithm developed by Shen et al. [6]. Modified to our problem, their algorithm is given as follows. (1)
DefineI0={i:(uij-ωi)<0andvij=0} and I-={i:(uij-ωi)<0andvij>0}.(2)
Sort the elements ofI- in increasing order of (uij-ωi)/vij and represent the elements by 1-,2-,…,n-, respectively, where n=|I-|.(3)
Find the value ofm that minimizes(3.14)∑i∈I0(uij-ωi)Pi+∑i∈I0vijPi+∑i=1,i∈I-m(uij-ωi)Pi+∑i=1,i∈I-mvijPi.(4)
The optimal solution to subproblemSPj is obtained by Pi=1 for i∈I0, P1-=P2-=⋯=Pm-=1 for i∈I- and Pi=0 for all other i∈I.After solving subproblemSPj for each j,fj is added to the optimal objective value of Ṽj. If Ṽj+fj<0, then we select the candidate distribution center j and set Xj=1; otherwise, we set Xj=0. For each selected distribution center j (those for which Xj=1), the assignment variables Yij are the same as the optimal Pi values in subproblem SPj; for each unselected distribution center j (those for which Xj=0), Yij=0,∀i∈I.Having solved the Lagrangian problem, the optimal Lagrange multipliers are found using a standard subgradient optimization procedure [37, 38]. The optimal objective value of the Lagrangian dual problem (3.8) provides a lower bound on the optimal objective value of (3.6).


## 3.1.2. Finding an Upper Bound

The lower bound solution obtained in each iteration of Lagrangian relaxation can violate constraints (3.2). In other words, some customers can be assigned to more than one DC in the obtained lower bound solution. Thus, in each iteration of Lagrangian relaxation, the obtained lower bound solution is converted into feasible upper bound solution using the procedure adapted from Shen et al. [6], as follows.(1)
If a customer is assigned to more than one DC, we find the DC which assigning it to the customer leads to the least objective (3.6). If the resulted objective (3.6) is less than the case that the customer is not served at all, we assign the customer to that DC. Otherwise, the customer is not assigned to any DC and is not served at all.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.If the obtained feasible solution results in a less value for (3.6) than the best known upper bound, it will be taken as the new upper bound solution. Also, it will be improved using customer reassignment algorithm.


## 3.1.3. Customer Reassignment

Each obtained upper bound is improved using following algorithm.(1)
For each customer, we check whether the objective value (3.6) is improved if the customer is assigned, instead to another located DC, or if it is not served at all. The best improving swap is performed.(2)
The unnecessary DCs, those which no longer serve any customers after performing step 1, are closed.


## 3.1.4. Variable Fixing

At the end of Lagrangian procedure, a variable fixing technique is employed. This method uses two following rules whose detailed proofs of the validity are presented by Shu et al. [2].(1)
Let LB and UB denote the current lower and upper bounds on the solution, respectively. If no DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=0 in the optimal solution to (3.8) at some iteration), and if LB+fj+Ṽj>UB, then no DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=0.(2)
If a DC is located at candidate site locationj∈J in the Lagrangian solution (i.e., Xj=1 in the optimal solution to (3.8) at some iteration), and if LB-(fj+Ṽj)>UB, then a DC will be located at j in any optimal solution for (3.6). Thus, we fix Xj=1.


## 3.2. Genetic Algorithm

In order to find near-optimal solutions for the model in relatively short time, here a solution approach based on genetic algorithm (GA) is developed. GA is a stochastic search and heuristic optimization technique based on the mechanism of natural genetics which has proven to be a very effective heuristic procedure to solve supply chain design problems, particularly problems of large scale [21, 30–34].GA starts with an initial set of random solution called population. Each solution in the population is called chromosome and each component of chromosome is designated by gene. The chromosomes evolve through successive iterations, called generations. During each generation, the chromosomes are evaluated, using some measures of fitness. To create next generation, new chromosomes (called offspring) are formed by crossover or mutation operators. Crossover operator combines two chromosomes from current generation, while mutation operator modifies a chromosome to form offspring.A new generation is created by selecting some of current chromosomes (called parents) and offspring based on the fitness values. Also, some chromosomes are rejected so as to keep the population size constant. Fitter chromosomes have higher probabilities of being selected. After several generations, the algorithms converge to the best chromosome, which may represent the optimum or suboptimal solution to the problem [33, 39, 40]. Figure 3 summarizes the steps of the proposed GA. In the following subsections, the developed GA is explained with more details.Figure 3
Genetic algorithm procedure.


### 3.2.1. Chromosome Representation

In this GA-based approach, each chromosome is indicated as a single-dimensional array. Letn be the number of candidate DCs and m be the number of customers. Then, each chromosome C can be indicated by: C=(Xj,Yi)=(X1,X2,…,Xn,Y1,Y2,…,Ym), where Xj correspond to the location genes and Yi correspond to the assignment genes. Location genes indicate where DCs are located and assignment genes represent how customers are assigned to the located DCs, respectively. In particular, if Xj=1, it means that candidate site j is selected as a DC location, while if Xj=0, candidate location j is not chosen as a DC site. Also, Yi=j represents that customer i is assigned to distribution center j. If customer i is not served at all, the corresponding assignment gene takes the value of 0 (i.e., Yi=0). For example, in Figure 4, distribution centers are located at 2 and 4. It follows that customers 1 and 4 are assigned to the DC at 4 and customer 2 is allocated to the DC at 2. Also, customer 3 is not served at all.Figure 4
Chromosome structure.


### 3.2.2. Chromosomes Fitness

The rank-based evaluation function is defined as the objective function (2.6) for the chromosomes. In fact, the objective (2.6) is calculated for each of the chromosomes. Obviously, chromosomes with better values for objective (2.6) will have the better rank.


### 3.2.3. Crossover Operator

Crossover operator generates offspring by merging parent chromosomes. LetCk denote the chromosomes of the population for k=1,2,…,pop-size. In order to determine which of the chromosomes are selected as parents for crossover operation, the following procedure is repeated from k=1 topop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PC, where the parameter PC is the probability of crossover. Then, randomly we group the selected parents C1′,C2′,C3′,… to the pairs (C1′,C2′),(C3′,C4′),…. Without loss of generality, let us explain the crossover operator on each pair by (C1′,C2′).Crossover operator assigns each customeri in offspring chromosome either to the DC which is allocated to customer i in parent chromosome C1′, or to the DC which is assigned to customer i in parent chromosome C2′. This occurs randomly and with probability of 0.5. If a customer is allocated to an unselected candidate DC site, a DC is located in that candidate location. A sample of crossover operator is shown in Figure 5.Figure 5
Sample of crossover.


### 3.2.4. Mutation Operator

Mutation operator modifies a chromosome to form an offspring. In order to decide which of chromosomesCk undergo mutation, the following practice is repeated from k=1 toPop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PM, where the parameter PM is the probability of mutation. A selected chromosome is modified by one of the two following types of mutation for several times. Each type of mutation is occurred with probability 0.5.The first type of mutation generates offspring by modifying the assignment genes of parent chromosome. That is, in the first type of mutation, two located DCs are selected randomly; lets and t denote them. Then, if any customer in the parent chromosome is assigned to s, that customer will be assigned to t, and if any customer is assigned to t, it will be allocated to s. The second type of mutation modifies location genes of parent chromosome to form an offspring. In other words, the second type of mutation randomly selects a location in which no DC is located; let t denotes it. Next, a DC is selected randomly from the located DCs and is named s. This type of mutation closes distribution center s and locates a DC at t instead. Then, all the customers assigned to distribution center s are allocated to distribution center t. The samples of mutation type 1 and mutation type 2 are illustrated in Figures 6 and 7, respectively.Figure 6
Sample of mutation type 1.Figure 7
Sample of mutation type 2.


## 3.2.1. Chromosome Representation

In this GA-based approach, each chromosome is indicated as a single-dimensional array. Letn be the number of candidate DCs and m be the number of customers. Then, each chromosome C can be indicated by: C=(Xj,Yi)=(X1,X2,…,Xn,Y1,Y2,…,Ym), where Xj correspond to the location genes and Yi correspond to the assignment genes. Location genes indicate where DCs are located and assignment genes represent how customers are assigned to the located DCs, respectively. In particular, if Xj=1, it means that candidate site j is selected as a DC location, while if Xj=0, candidate location j is not chosen as a DC site. Also, Yi=j represents that customer i is assigned to distribution center j. If customer i is not served at all, the corresponding assignment gene takes the value of 0 (i.e., Yi=0). For example, in Figure 4, distribution centers are located at 2 and 4. It follows that customers 1 and 4 are assigned to the DC at 4 and customer 2 is allocated to the DC at 2. Also, customer 3 is not served at all.Figure 4
Chromosome structure.


## 3.2.2. Chromosomes Fitness

The rank-based evaluation function is defined as the objective function (2.6) for the chromosomes. In fact, the objective (2.6) is calculated for each of the chromosomes. Obviously, chromosomes with better values for objective (2.6) will have the better rank.


## 3.2.3. Crossover Operator

Crossover operator generates offspring by merging parent chromosomes. LetCk denote the chromosomes of the population for k=1,2,…,pop-size. In order to determine which of the chromosomes are selected as parents for crossover operation, the following procedure is repeated from k=1 topop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PC, where the parameter PC is the probability of crossover. Then, randomly we group the selected parents C1′,C2′,C3′,… to the pairs (C1′,C2′),(C3′,C4′),…. Without loss of generality, let us explain the crossover operator on each pair by (C1′,C2′).Crossover operator assigns each customeri in offspring chromosome either to the DC which is allocated to customer i in parent chromosome C1′, or to the DC which is assigned to customer i in parent chromosome C2′. This occurs randomly and with probability of 0.5. If a customer is allocated to an unselected candidate DC site, a DC is located in that candidate location. A sample of crossover operator is shown in Figure 5.Figure 5
Sample of crossover.


## 3.2.4. Mutation Operator

Mutation operator modifies a chromosome to form an offspring. In order to decide which of chromosomesCk undergo mutation, the following practice is repeated from k=1 toPop-size: generating a random number r from the interval [0,1], the chromosome Ck will be selected as a parent provided that r<PM, where the parameter PM is the probability of mutation. A selected chromosome is modified by one of the two following types of mutation for several times. Each type of mutation is occurred with probability 0.5.The first type of mutation generates offspring by modifying the assignment genes of parent chromosome. That is, in the first type of mutation, two located DCs are selected randomly; lets and t denote them. Then, if any customer in the parent chromosome is assigned to s, that customer will be assigned to t, and if any customer is assigned to t, it will be allocated to s. The second type of mutation modifies location genes of parent chromosome to form an offspring. In other words, the second type of mutation randomly selects a location in which no DC is located; let t denotes it. Next, a DC is selected randomly from the located DCs and is named s. This type of mutation closes distribution center s and locates a DC at t instead. Then, all the customers assigned to distribution center s are allocated to distribution center t. The samples of mutation type 1 and mutation type 2 are illustrated in Figures 6 and 7, respectively.Figure 6
Sample of mutation type 1.Figure 7
Sample of mutation type 2.


## 4. Computational Results

This section summarizes the computational experience with the solution approaches outlined in the previous section. Two sets of experiments were designed. The objective of the first set of experiments was to evaluate the performance of the proposed solution methods in terms of the solution quality and time. The second set of experiments was designed to study the benefits of considering facility disruptions during the supply chain design phase. The developed solution approaches were coded in Visual Basic.Net and executed on Pentium 5 computer with 1.00 GB RAM and 2.00 GHz CPU.


### 4.1. Experimental Design

The computational experiments were conducted on the 49-node, 88-node, and 150-node datasets described in Daskin [41]. These datasets have been very popular in the literature and have been used in lots of studies to validate the new solution methods [1, 3, 6–8, 24, 27, 28, 30, 34]. The 49-node dataset represents the capitals of the lower 48 United States plus Washington, DC; the 88-node data set represents the 50 largest cities in the 1990 U.S. census along with the 49-node dataset minus duplicates; the 150-node dataset contains the 150 largest cities in the 1990 U.S. census.For all three data sets, the mean of demand was obtained by dividing the population data given in Daskin [41] by 1000. Fixed costs of locating DCs (fj) were gained by dividing the fixed cost in Daskin [41] by 10 for the 49-node problem and by 100 for 88-node problem. For the 150-node problem, fixed locating costs were set to 10,000 for all the candidate DC locations. We set the per-unit cost to ship from distribution center j to customer i,dij, to the great-circle distance between these locations. Similar to Daskin et al. [1], the fixed ordering aj and shipping costs bj were set to 10 and the variable shipping cost cj was set to 5 for all DCs. Also, we set the inventory holding cost per unit of product to 1 and set the selling price at each customer to 500. The disruption scenarios were generated randomly. In particular, for each s∈S and for each j∈J,rsjrandomly was set to 0 with probability of ninety percent, or to a random number from (0,1] with probability of ten percent. The probability of occurrence associated with each scenario was drawn uniformly from (0,1] and then normalized such that the total probability of all the scenarios is equal to 1. Table 1 shows the parameters for the Lagrangian relaxation approach in the computational experiments. These parameters were set similar to Qi et al. [27]. The parameters for the genetic algorithm were set based on the optimal values suggested by Sourirajan et al. [34] and Grefenstette [42]. These parameters are given in Table 2.Table 1
Parameters for Lagrangian relaxation.
ParameterValueInitial value of the scalar used to define step size2Minimum value of the scalar10-12Maximum number of iterations before halving the scalar12Maximum number of iterations400Initial Lagrange multiplier value0Table 2
Parameters for the genetic algorithm.
ParameterValuePopulation size100Probability of crossover0.95Probability of mutation0.01Number of generations400


### 4.2. Performance of the Algorithms

Tables3, 4, and 5 present the computational results with Lagrangian relaxation on 49-node, 88-node, and 150-node problems, respectively. In each table, the first column gives the number of scenarios in the problem. The columns labeled β and θ give the values of β and θ. The columns marked LB and UB give the lower and upper bounds obtained for the profit maximizing model. The last column in each table indicates the percentage gap between the obtained upper and lower bounds and is gained by ((UB-LB)/LB)×100. In all of the cases, the gap was always less than one percent representing that bounds provided by the Lagrangian relaxation approach are very tight and can be relied upon to produce good feasible solutions.Table 3
Computational results with Lagrangian relaxation for 49-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1108537010856490.026(2)200.0050.1931632932200.30.061(3)200.0050.59222889228810.064(4)200.0051913837914104.90.029(5)200.0055736483737237.70.102(6)400.0010.1108395610842240.025(7)400.0050.1931431932001.20.061(8)400.0050.5922758923345.30.064(9)400.0051914907915496.20.064(10)400.0055741487742397.60.123(11)600.0010.1107975610800510.027(12)600.0050.1918411918992.60.063(13)600.0050.59090539096650.067(14)600.0051900592901211.60.069(15)600.0055722305723255.60.131Table 4
Computational results with Lagrangian relaxation for 88-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1221670922169070.009(2)200.0050.1218604621915040.249(3)200.0050.5218134521826320.059(4)200.0051217717021834240.286(5)200.0055210308621097390.315(6)400.0010.1221752022197000.098(7)400.0050.1218790921931680.240(8)400.0050.5218334221891980.268(9)400.0051217922721853820.282(10)400.0055210525421110110.273(11)600.0010.1221636122187390.107(12)600.0050.1218400421866800.122(13)600.0050.5217928921854440.282(14)600.0051217512421816770.300(15)600.0055209965921094050.462Table 5
Computational results with Lagrangian relaxation for 150-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1289491029018560.239(2)200.0050.1288606428884410.082(3)200.0050.5287683528790140.076(4)200.0051286898428728500.135(5)200.0055275404627667750.460(6)400.0010.1289467228960580.048(7)400.0050.1288594528883220.082(8)400.0050.5287675028800200.114(9)400.0051286884328729070.141(10)400.0055275231927650490.460(11)600.0010.1289361729353891.423(12)600.0050.1288566628880430.082(13)600.0050.5287610228794710.117(14)600.0051286790928720720.145(15)600.0055274801527750040.973Also, the computational results with genetic algorithm are presented in Tables6–8. Similarly, in each table, the first column gives the number of scenarios in the problem and the columns labeled β and θ give the values of β and θ. The column marked GA represents the objective value obtained by genetic algorithm. The column labeled UB indicates the upper bound obtained by Lagrangian relaxation for the integrated model. The percentage gap between this upper bound and the objective value obtained by genetic algorithm is given in the column marked GAP. Tables 6–8 show that the gap does not exceed one percent suggesting that the obtained solutions by genetic algorithm are close to optimal values.Table 6
Computational results with genetic algorithm for 49-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1108380910856490.170(2)200.0050.1925985.3932200.30.667(3)200.0050.5919960.89228810.316(4)200.0051913925.2914104.90.020(5)200.0055730852.1737237.70.866(6)400.0010.1108213010842240.193(7)400.0050.1931904.4932001.20.010(8)400.0050.5922998.5923345.30.038(9)400.0051915301.9915496.20.021(10)400.0055739498.9742397.60.390(11)600.0010.1107854710800510.139(12)600.0050.1918801.3918992.60.021(13)600.0050.5906391.29096650.360(14)600.0051901080.2901211.60.015(15)600.0055716398.5723255.60.948Table 7
Computational results with genetic algorithm for 88-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1220725322169070.435(2)200.0050.1217402021915040.798(3)200.0050.5218193321826320.032(4)200.0051216821621834240.697(5)200.0055209154021097390.863(6)400.0010.1221946922197000.010(7)400.0050.1217182421931680.973(8)400.0050.5218202921891980.327(9)400.0051218484921853820.024(10)400.0055210221521110110.417(11)600.0010.1219970822187390.858(12)600.0050.1216717321866800.892(13)600.0050.5217067721854440.676(14)600.0051216268021816770.871(15)600.0055209864421094050.510Table 8
Computational results with genetic algorithm for 150-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1287766229018560.834(2)200.0050.1286058728884410.964(3)200.0050.5287681928790140.076(4)200.0051286131928728500.401(5)200.0055274115027667750.926(6)400.0010.1287274528960580.805(7)400.0050.1287091728883220.603(8)400.0050.5287867928800200.047(9)400.0051285972528729070.459(10)400.0055274747527650490.636(11)600.0010.1290103529353891.170(12)600.0050.1286551528880430.780(13)600.0050.5286233828794710.595(14)600.0051285040728720720.754(15)600.0055274779027750040.981Table9 summarizes the computational results and compares the performance of Lagrangian relaxation approach with that of genetic algorithm in terms of the solution quality and time. In this table, GAP and CPU time (per second) values are averaged and reported for different values of β and θ. It follows from Table 9 that typically Lagrangian relaxation results in lower values of GAP than genetic algorithm. Thus, the performance of Lagrangian relaxation is better than genetic algorithm from the quality point of view. However, the computational time of genetic algorithm is consistently lower than that of Lagrangian relaxation. This suggests that genetic algorithm is faster than Lagrangian relaxation at solving the model.Table 9
Performance Comparison between Lagrangian relaxation and genetic algorithm.
Data setβθLagrangian relaxationgenetic algorithmAverage GAPAverage  time (sec.)Average GAPAverage time (sec.)49-node0.0010.10.0259045.80.1672931.10.0050.10.0618072.90.2326341.60.0050.50.0618073.10.237961.20.00510.054142.60.0184891.40.00550.1188193.20.7348961.388-node0.0010.10.0714512.60.4345263.40.0050.10.20375810.90.8876883.80.0050.50.20272114.30.3450743.10.00510.28950118.70.5305573.20.00550.3500598.60.5964793.9150-node0.0010.10.57010958.50.936365.30.0050.10.08233141.10.7823385.60.0050.50.10210142.90.2392875.80.00510.14034739.20.5381855.90.00550.63103352.50.8474755.2


### 4.3. The Benefits of Considering Disruptions in the Supply Chain Design Model

This section compares the performance of two different approaches for designing the supply chain networks. The first approach considers disruptions scenarios during making supply chain design decisions including location and allocation decisions, as we do in this article. That is, the first approach uses the presented model in this study to determine the supply chain design decisions. The second approach, however, makes supply chain design decisions without taking into consideration the disruptions scenarios. In order to determine the supply chain design decisions under the second method, the proposed model by Daskin et al. [1] is used. By comparing the total profits under these two approaches, the benefits of considering disruptions in the supply chain design model are implied.Tables10, 11, and 12 demonstrate the total profits of the two approaches for the same datasets in Section 4.1. The first column marked the number of scenarios,β and θ give the number of scenarios and the values of β and θ for each instance, respectively. The total profits under the first and second approaches are denoted by TP1 and TP2, respectively. The last column in each table represents the percentage difference between the total profit of the first and that of the second approach. In other words, this column indicates the amount of increase in total profit by following the first approach instead of the second approach for each instance.Table 10
The benefits of considering disruptions in the supply chain design model for 49-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.110853701017859.996.22(2)200.0050.1931632877317.8545.83(3)200.0050.5922288868887.5255.79(4)200.0051913837864855.3375.36(5)200.0055736483714093.9173.04(6)400.0010.11083956975235.21310.03(7)400.0050.1931431848347.3558.92(8)400.0050.5922758841001.6418.86(9)400.0051914907836682.4528.55(10)400.0055741487703819.465.08(11)600.0010.11079756891230.60217.46(12)600.0050.1918411778169.6415.27(13)600.0050.5909053773513.19814.91(14)600.0051900592776400.36313.79(15)600.0055722305645162.82610.68Table 11
The benefits of considering disruptions in the supply chain design model for 88-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.122167091996589.89.93(2)200.0050.121860461998483.258.58(3)200.0050.521813451998984.568.36(4)200.005121771702001254.668.08(5)200.005521030861992253.375.27(6)400.0010.122175201907732.4613.97(7)400.0050.121879091929298.1611.82(8)400.0050.521833421931384.3311.54(9)400.005121792271936025.2711.16(10)400.005521052541929465.298.35(11)600.0010.122163611715241.7822.61(12)600.0050.121840041735846.3820.52(13)600.0050.521792891722074.1720.98(14)600.005121751241730093.6320.46(15)600.005520996591737047.8917.27Table 12
The benefits of considering disruptions in the supply chain design model for 150-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.122167091939177.0312.52(2)200.0050.121860461940553.0311.23(3)200.0050.521813451945541.6110.81(4)200.005121771701946825.4110.58(5)200.005521030861914649.498.96(6)400.0010.122175201805726.5418.57(7)400.0050.121879091808963.1617.32(8)400.0050.521833421811082.1917.05(9)400.005121792271815078.1716.71(10)400.005521052541809044.7614.07(11)600.0010.122163611586471.228.42(12)600.0050.121840041589299.7127.23(13)600.0050.521792891598072.6226.67(14)600.005121751241606981.6126.12(15)600.005520996591617787.2622.95It follows from Tables10–12 that considering facility disruptions during the supply chain design phase can lead to increase in total profit. Particularly, the results show that the benefit of considering disruptions scenarios in the supply chain design model can be significant up to twenty-eight percent. There are some properties for the presented model in this study which can explain why considering disruptions in the supply chain design phase results in higher profits. First, in the optimal solution DCs are more likely to be opened at locations with lower risks of disruptions to reduce the lost-sales costs due to disruptions. Likewise, when a DC is more likely to be disrupted, fewer customers will be assigned to this DC in the optimal solution and lower lost-sales costs will be incurred. Also, when the possible disruptions are considered during the design of supply chain, fewer customers are selected to be served. That is, the optimal solution will involve more customers not served by any DC, in order to reduce the risk of incurring considerable lost-sales costs.It can be concluded that not only the facility location costs affect the location decisions, but also the rates of disruptions at facilities are important factors in determining where to locate DCs. For instance, despite low facility location cost, it is possible that a DC is not opened in a candidate site because of high rates of disruptions. Similarly, the number of customers to be served and the assignments of customers to DCs can be dependent to the rates of disruptions at the DCs. Therefore, significant increase in total profit can be realized, if the facility disruptions are considered in the supply chain design model.


## 4.1. Experimental Design

The computational experiments were conducted on the 49-node, 88-node, and 150-node datasets described in Daskin [41]. These datasets have been very popular in the literature and have been used in lots of studies to validate the new solution methods [1, 3, 6–8, 24, 27, 28, 30, 34]. The 49-node dataset represents the capitals of the lower 48 United States plus Washington, DC; the 88-node data set represents the 50 largest cities in the 1990 U.S. census along with the 49-node dataset minus duplicates; the 150-node dataset contains the 150 largest cities in the 1990 U.S. census.For all three data sets, the mean of demand was obtained by dividing the population data given in Daskin [41] by 1000. Fixed costs of locating DCs (fj) were gained by dividing the fixed cost in Daskin [41] by 10 for the 49-node problem and by 100 for 88-node problem. For the 150-node problem, fixed locating costs were set to 10,000 for all the candidate DC locations. We set the per-unit cost to ship from distribution center j to customer i,dij, to the great-circle distance between these locations. Similar to Daskin et al. [1], the fixed ordering aj and shipping costs bj were set to 10 and the variable shipping cost cj was set to 5 for all DCs. Also, we set the inventory holding cost per unit of product to 1 and set the selling price at each customer to 500. The disruption scenarios were generated randomly. In particular, for each s∈S and for each j∈J,rsjrandomly was set to 0 with probability of ninety percent, or to a random number from (0,1] with probability of ten percent. The probability of occurrence associated with each scenario was drawn uniformly from (0,1] and then normalized such that the total probability of all the scenarios is equal to 1. Table 1 shows the parameters for the Lagrangian relaxation approach in the computational experiments. These parameters were set similar to Qi et al. [27]. The parameters for the genetic algorithm were set based on the optimal values suggested by Sourirajan et al. [34] and Grefenstette [42]. These parameters are given in Table 2.Table 1
Parameters for Lagrangian relaxation.
ParameterValueInitial value of the scalar used to define step size2Minimum value of the scalar10-12Maximum number of iterations before halving the scalar12Maximum number of iterations400Initial Lagrange multiplier value0Table 2
Parameters for the genetic algorithm.
ParameterValuePopulation size100Probability of crossover0.95Probability of mutation0.01Number of generations400


## 4.2. Performance of the Algorithms

Tables3, 4, and 5 present the computational results with Lagrangian relaxation on 49-node, 88-node, and 150-node problems, respectively. In each table, the first column gives the number of scenarios in the problem. The columns labeled β and θ give the values of β and θ. The columns marked LB and UB give the lower and upper bounds obtained for the profit maximizing model. The last column in each table indicates the percentage gap between the obtained upper and lower bounds and is gained by ((UB-LB)/LB)×100. In all of the cases, the gap was always less than one percent representing that bounds provided by the Lagrangian relaxation approach are very tight and can be relied upon to produce good feasible solutions.Table 3
Computational results with Lagrangian relaxation for 49-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1108537010856490.026(2)200.0050.1931632932200.30.061(3)200.0050.59222889228810.064(4)200.0051913837914104.90.029(5)200.0055736483737237.70.102(6)400.0010.1108395610842240.025(7)400.0050.1931431932001.20.061(8)400.0050.5922758923345.30.064(9)400.0051914907915496.20.064(10)400.0055741487742397.60.123(11)600.0010.1107975610800510.027(12)600.0050.1918411918992.60.063(13)600.0050.59090539096650.067(14)600.0051900592901211.60.069(15)600.0055722305723255.60.131Table 4
Computational results with Lagrangian relaxation for 88-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1221670922169070.009(2)200.0050.1218604621915040.249(3)200.0050.5218134521826320.059(4)200.0051217717021834240.286(5)200.0055210308621097390.315(6)400.0010.1221752022197000.098(7)400.0050.1218790921931680.240(8)400.0050.5218334221891980.268(9)400.0051217922721853820.282(10)400.0055210525421110110.273(11)600.0010.1221636122187390.107(12)600.0050.1218400421866800.122(13)600.0050.5217928921854440.282(14)600.0051217512421816770.300(15)600.0055209965921094050.462Table 5
Computational results with Lagrangian relaxation for 150-node problem.
Number of scenariosβθLBUBGAP(1)200.0010.1289491029018560.239(2)200.0050.1288606428884410.082(3)200.0050.5287683528790140.076(4)200.0051286898428728500.135(5)200.0055275404627667750.460(6)400.0010.1289467228960580.048(7)400.0050.1288594528883220.082(8)400.0050.5287675028800200.114(9)400.0051286884328729070.141(10)400.0055275231927650490.460(11)600.0010.1289361729353891.423(12)600.0050.1288566628880430.082(13)600.0050.5287610228794710.117(14)600.0051286790928720720.145(15)600.0055274801527750040.973Also, the computational results with genetic algorithm are presented in Tables6–8. Similarly, in each table, the first column gives the number of scenarios in the problem and the columns labeled β and θ give the values of β and θ. The column marked GA represents the objective value obtained by genetic algorithm. The column labeled UB indicates the upper bound obtained by Lagrangian relaxation for the integrated model. The percentage gap between this upper bound and the objective value obtained by genetic algorithm is given in the column marked GAP. Tables 6–8 show that the gap does not exceed one percent suggesting that the obtained solutions by genetic algorithm are close to optimal values.Table 6
Computational results with genetic algorithm for 49-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1108380910856490.170(2)200.0050.1925985.3932200.30.667(3)200.0050.5919960.89228810.316(4)200.0051913925.2914104.90.020(5)200.0055730852.1737237.70.866(6)400.0010.1108213010842240.193(7)400.0050.1931904.4932001.20.010(8)400.0050.5922998.5923345.30.038(9)400.0051915301.9915496.20.021(10)400.0055739498.9742397.60.390(11)600.0010.1107854710800510.139(12)600.0050.1918801.3918992.60.021(13)600.0050.5906391.29096650.360(14)600.0051901080.2901211.60.015(15)600.0055716398.5723255.60.948Table 7
Computational results with genetic algorithm for 88-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1220725322169070.435(2)200.0050.1217402021915040.798(3)200.0050.5218193321826320.032(4)200.0051216821621834240.697(5)200.0055209154021097390.863(6)400.0010.1221946922197000.010(7)400.0050.1217182421931680.973(8)400.0050.5218202921891980.327(9)400.0051218484921853820.024(10)400.0055210221521110110.417(11)600.0010.1219970822187390.858(12)600.0050.1216717321866800.892(13)600.0050.5217067721854440.676(14)600.0051216268021816770.871(15)600.0055209864421094050.510Table 8
Computational results with genetic algorithm for 150-node problem.
Number of scenariosβθGAUBGAP(1)200.0010.1287766229018560.834(2)200.0050.1286058728884410.964(3)200.0050.5287681928790140.076(4)200.0051286131928728500.401(5)200.0055274115027667750.926(6)400.0010.1287274528960580.805(7)400.0050.1287091728883220.603(8)400.0050.5287867928800200.047(9)400.0051285972528729070.459(10)400.0055274747527650490.636(11)600.0010.1290103529353891.170(12)600.0050.1286551528880430.780(13)600.0050.5286233828794710.595(14)600.0051285040728720720.754(15)600.0055274779027750040.981Table9 summarizes the computational results and compares the performance of Lagrangian relaxation approach with that of genetic algorithm in terms of the solution quality and time. In this table, GAP and CPU time (per second) values are averaged and reported for different values of β and θ. It follows from Table 9 that typically Lagrangian relaxation results in lower values of GAP than genetic algorithm. Thus, the performance of Lagrangian relaxation is better than genetic algorithm from the quality point of view. However, the computational time of genetic algorithm is consistently lower than that of Lagrangian relaxation. This suggests that genetic algorithm is faster than Lagrangian relaxation at solving the model.Table 9
Performance Comparison between Lagrangian relaxation and genetic algorithm.
Data setβθLagrangian relaxationgenetic algorithmAverage GAPAverage  time (sec.)Average GAPAverage time (sec.)49-node0.0010.10.0259045.80.1672931.10.0050.10.0618072.90.2326341.60.0050.50.0618073.10.237961.20.00510.054142.60.0184891.40.00550.1188193.20.7348961.388-node0.0010.10.0714512.60.4345263.40.0050.10.20375810.90.8876883.80.0050.50.20272114.30.3450743.10.00510.28950118.70.5305573.20.00550.3500598.60.5964793.9150-node0.0010.10.57010958.50.936365.30.0050.10.08233141.10.7823385.60.0050.50.10210142.90.2392875.80.00510.14034739.20.5381855.90.00550.63103352.50.8474755.2


## 4.3. The Benefits of Considering Disruptions in the Supply Chain Design Model

This section compares the performance of two different approaches for designing the supply chain networks. The first approach considers disruptions scenarios during making supply chain design decisions including location and allocation decisions, as we do in this article. That is, the first approach uses the presented model in this study to determine the supply chain design decisions. The second approach, however, makes supply chain design decisions without taking into consideration the disruptions scenarios. In order to determine the supply chain design decisions under the second method, the proposed model by Daskin et al. [1] is used. By comparing the total profits under these two approaches, the benefits of considering disruptions in the supply chain design model are implied.Tables10, 11, and 12 demonstrate the total profits of the two approaches for the same datasets in Section 4.1. The first column marked the number of scenarios,β and θ give the number of scenarios and the values of β and θ for each instance, respectively. The total profits under the first and second approaches are denoted by TP1 and TP2, respectively. The last column in each table represents the percentage difference between the total profit of the first and that of the second approach. In other words, this column indicates the amount of increase in total profit by following the first approach instead of the second approach for each instance.Table 10
The benefits of considering disruptions in the supply chain design model for 49-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.110853701017859.996.22(2)200.0050.1931632877317.8545.83(3)200.0050.5922288868887.5255.79(4)200.0051913837864855.3375.36(5)200.0055736483714093.9173.04(6)400.0010.11083956975235.21310.03(7)400.0050.1931431848347.3558.92(8)400.0050.5922758841001.6418.86(9)400.0051914907836682.4528.55(10)400.0055741487703819.465.08(11)600.0010.11079756891230.60217.46(12)600.0050.1918411778169.6415.27(13)600.0050.5909053773513.19814.91(14)600.0051900592776400.36313.79(15)600.0055722305645162.82610.68Table 11
The benefits of considering disruptions in the supply chain design model for 88-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.122167091996589.89.93(2)200.0050.121860461998483.258.58(3)200.0050.521813451998984.568.36(4)200.005121771702001254.668.08(5)200.005521030861992253.375.27(6)400.0010.122175201907732.4613.97(7)400.0050.121879091929298.1611.82(8)400.0050.521833421931384.3311.54(9)400.005121792271936025.2711.16(10)400.005521052541929465.298.35(11)600.0010.122163611715241.7822.61(12)600.0050.121840041735846.3820.52(13)600.0050.521792891722074.1720.98(14)600.005121751241730093.6320.46(15)600.005520996591737047.8917.27Table 12
The benefits of considering disruptions in the supply chain design model for 150-node problem.
Number of scenariosβθTP1TP2Profit difference (%)(1)200.0010.122167091939177.0312.52(2)200.0050.121860461940553.0311.23(3)200.0050.521813451945541.6110.81(4)200.005121771701946825.4110.58(5)200.005521030861914649.498.96(6)400.0010.122175201805726.5418.57(7)400.0050.121879091808963.1617.32(8)400.0050.521833421811082.1917.05(9)400.005121792271815078.1716.71(10)400.005521052541809044.7614.07(11)600.0010.122163611586471.228.42(12)600.0050.121840041589299.7127.23(13)600.0050.521792891598072.6226.67(14)600.005121751241606981.6126.12(15)600.005520996591617787.2622.95It follows from Tables10–12 that considering facility disruptions during the supply chain design phase can lead to increase in total profit. Particularly, the results show that the benefit of considering disruptions scenarios in the supply chain design model can be significant up to twenty-eight percent. There are some properties for the presented model in this study which can explain why considering disruptions in the supply chain design phase results in higher profits. First, in the optimal solution DCs are more likely to be opened at locations with lower risks of disruptions to reduce the lost-sales costs due to disruptions. Likewise, when a DC is more likely to be disrupted, fewer customers will be assigned to this DC in the optimal solution and lower lost-sales costs will be incurred. Also, when the possible disruptions are considered during the design of supply chain, fewer customers are selected to be served. That is, the optimal solution will involve more customers not served by any DC, in order to reduce the risk of incurring considerable lost-sales costs.It can be concluded that not only the facility location costs affect the location decisions, but also the rates of disruptions at facilities are important factors in determining where to locate DCs. For instance, despite low facility location cost, it is possible that a DC is not opened in a candidate site because of high rates of disruptions. Similarly, the number of customers to be served and the assignments of customers to DCs can be dependent to the rates of disruptions at the DCs. Therefore, significant increase in total profit can be realized, if the facility disruptions are considered in the supply chain design model.


## 5. Conclusion

This paper has addressed a supply chain design problem where distribution centers are subject to partial and complete disruptions. The problem has been formulated as a nonlinear mixed-integer programming which maximizes the total profit for the whole supply chain. The integrated model simultaneously determines the optimal number and location of DCs, the subset of customers to serve, the assignment of customers to DCs, and the cycle-order quantities at DCs. In order to obtain near-optimal solutions with reasonable computational requirements, two heuristics based on Lagrangian relaxation and genetic algorithm have been presented. Computational results for different data sets have revealed that the proposed solution approaches are effective. Also, it has been demonstrated that the benefits of considering disruptions during supply chain design phase can be significant.In future, it would be interesting to formulate the problem when suppliers are unreliable. Also, the model can be extended to consider constraints on the maximum capacity of DCs, or on the maximum supply that can be filled by suppliers. Finally, considering routing decisions in the model makes it more useful.



---

*Source: 234324-2012-03-07.xml*