# Designing a Fuzzy Strategic Integrated Multiechelon Agile Supply Chain Network

Authors: Morteza Abbasi; Reza Hosnavi; Mehrdad Mohammadi
Journal: Journal of Optimization
 (2013)

---

## Abstract

This paper integrates production, distribution and logistics activities at the strategic decision making level, where the objective is to design a multiechelon supply chain network considering agility as a key design criterion. A network with five echelons of supply chains including suppliers, plants, distribution centers, cross-docks, and customer zones is addressed in this paper. The problem has been mathematically formulated as a biobjective optimization model that aims to minimize the cost (fixed and variable) and maximize the plant flexibility and volume flexibility. A novel multiobjective parallel simulating annealing algorithm (MOPSA) is proposed to obtain the Pareto-optimal solutions of the problem. The performance of the proposed solution algorithm is compared with two well-known metaheuristics, namely, nondominated sorting genetic algorithm (NSGA-II) and Pareto archive evolution strategy (PAES). Computational results show that MOPSA outperforms the other metaheuristics.


---

## Body

## 1. Introduction

The ever-growing global market and impact of competitive environment caused to receive more attention to supply chain for producing and delivering products at low cost, high quality, and short lead time. A supply chain can be defined as a network of suppliers, plants, warehouses, distribution centers, and channels between them which are organized to acquire/achieve raw materials, convert them to finished products, and distribute the finished products to customers; on the other hand, supply chain is a chain that links each element of the manufacturing and supply process from raw materials through to the customer. In most researches, production, distribution and logistics are considered exclusive while considering an integrated approach has more benefit for network. Integrated approach is based on integrating decisions of different functions into a single optimization model. For example, in the supply chain management, demands are satisfied by keeping the service level at a deterministic level then minimizing costs, however; in an integrated approach the costs are minimized and simultaneously the service level is improved by redesigning its supply chain network.The supply chain network design decisions are classified into three categories based on their importance and the length of the planning horizon. First, strategic (long-range) planning that lasts more than one year; decisions about company selection and facility location are generally placed in this category. Second, tactical (medium) planning that is revised every few months includes decisions about production, inventory, and logistics. Finally, operational or short-range planning involves within-days or shift decisions; raw material, semifinished and finished product flows in the network are operational decisions that are easily classified in this class.In competitive environment, companies cope with several challenges such as speed of designing, manufacturing, and distributing, while they have to consider simultaneously higher product efficiency and lower operational cost. In this situation, the concept of agile supply chain was proposed to response to these challenges. Agility is used for showing the flexibility and speed in coping with innovative products and unpredictable demands. An agile company should be able to operate in an uncertain and challenging market and use the opportunities that are continuously appearing in competitive environment. Besides, responsiveness is a major factor of an agile supply chain that is as important as a high level of efficiency in cost, quality, and smooth operations flow. In an agile manufacturing when an opportunity appears, the virtual organization can be formed to attain the objective. The virtual organization, according to Goldman et al. [1], is defined as an organizational tool for agile competitors who are integrated by sharing core competencies and resources to produce a particular product that could not be made by each of competitors. Because of short-live durations of virtual organization, strategic and tactical planning should be considered together in the agile manufacturing paradigm.Uncertainties involved in the supply chain optimization can be categorized into two main groups: (1) environmental and (2) system [2]. Environmental uncertainties are attributed to the performance of each member of supply chain consisting of suppliers, manufacturers, and so forth. System uncertainties are attributed to routine processes in the supply chain such as production and distribution. The above-mentioned uncertainties adversely contribute to the quality of decisions made in strategic, tactical, and operational levels of supply chain. With respect to the alluded points and necessity of considering uncertainties in adopting various decisions in the supply chain networks, researchers, in tune with designing their supply chain networks, have strived to include foregoing two aspects of uncertainties to better represent the practical features of real-world problems. Thus, for this purpose, advanced optimization methodologies are desired. Hence, many researchers have recently made attempts to model the uncertainty in supply chains [3–5].In the real world, after designing a supply chain network, the costs, demands, distances, times, and other parameters of the problem may change due to uncertain environment. For example, with variability in the transportation time from an origin to a destination, there is a possibility that a flow may not be delivered on time. A failure in on-time delivery may result in a huge and nonmeasureable cost, such as lost-opportunity cost and lost-sale cost due to customer churning, to compensate an unsatisfied customer.In this paper, we consider a multiechelon supply chain network design problem consisting of five echelon and agility as a design criterion with an integrated approach. The objective of this model is to determine the best suppliers and plants in upstream level and determine the number of distribution and cross-dock that should be opened in downstream level to satisfy demands in lower cost and higher efficiency. Usage of cross-dock caused to decrease the construction and setup costs; this concept is used in this paper at the intermediate level between distribution and customers zones. Also, in order to cope with uncertain environment, some parameters of the model such as demand and cost are considered as fuzzy parameters.The rest of the paper is organized as follows. Section2 presents related literature. In Section 3, the fuzzy multiobjective mathematical model for agile supply chain management (ASCM) is presented. Defuzzification of the fuzzy model is proposed in Section 4. Section 5 explains the proposed multiobjective metaheuristic algorithm followed by experimental results in Section 6 followed by conclusion is Section 7.


## 2. Literature Review

Many studies have been done in traditional supply chains dealing with supplier selection and allocation problem of facilities for satisfying demand. Samadhi and Hoang [6] proposed a three-stage hierarchical procedure for ranking suppliers based on strategic, manufacturing, and logistics factors. De Boer et al. [7] suggested an outranking method for evaluating the supplier by use of the multicriteria factors. Jayaraman [8] considered a classical capacitated warehouse problem. In this paper, the number and location of warehouses are determined and the demand of the customers is satisfied by allocating the customers to warehouses at the minimum cost. Tsiakis et al. [9] addressed a multiproduct, multi-echelon supply chain network under uncertainty. It is assumed that the locations of manufacturing sites are fixed although the number, location, and capacity of warehouses and distribution centers should be determined from a set of potential locations. The objective of this model is to minimize the total cost of the network. This model is solved by using branch-and-bound techniques. Shankar et al. [10] presented a multiobjective optimization of single-product for a multi-echelon supply chain network including suppliers, plants, distribution centers (DCs), and customer zones (CZs). The goal of this model is to find the number and location of plants, the flow of raw materials from suppliers to plants, and the quantity of products to be shipped from plants to DCs, from DCs to CZs with optimizing the objectives simultaneously. The objectives minimize location and shipment costs while the maximum customer demands should be met. They proposed multiobjective hybrid particle swarm optimization (MOHPSO) algorithm for solving their model.Integrated analysis in supply chain network is an approach based on the integrating the decisions of different functions such as supply, processing, distribution, inventory management, production planning, and facilities location on a single optimization model for a production-distribution system. Sarmiento and Nagi [11] presented a review on integrated model of production-distribution system with considering transportation system. This paper consists of two sections: analyses of production, distribution and inventory planning and the inventory/routing problem. This paper considers decisions at the strategic and tactical levels. Alonso-Ayuso et al. [12] presented a two-stage stochastic supply chain planning. This model determines the size of plant, product selection, assigning product to plants, and vendor selection. The strategic decisions are made in the first stage and tactical decisions are made in second stage. The objective of this model is to maximize the expected profit of the supply chain. Teo and Shu [13] proposed a distribution network design problem with integrating transportation and inventory function. The goal of this model is to determine the number and location of warehouses, assign the retailers to warehouses, and determine the optimal inventory policy. In this model, there is a trade-off between inventory cost, direct shipment cost, and facility location cost, and the objective is to minimize the multi-echelon inventory, transportation, and facility location costs.Elhedhli and Coffin [14] presented an integrated production-distribution problem involving three echelons: plants, distributions, and customer zones. The objective of this model is to minimize production, distribution, handling, transportation costs, and the fixed costs of opening distribution centers. The solution method for this model is based on two-level Lagrangian relaxation approach and branch-and-price algorithm. Amiri [15] presented a distribution network design problem and considered the multiple levels of capacities for plants and warehouses in supply chain system. In this model, the best strategies for distributing the product from the plants to the warehouses and from the warehouses to the customers are selected. The objective is to minimize total costs of the distribution network by determining the optimum numbers, locations and capacities of plants and warehouses. Van Hoek et al. [16], considered the concept of agility in supply chain. In their model, the agile supply chain consists of customer sensitivity, process integration, and information integration. Dotoli et al. [17] proposed a single and multiobjective optimization model for agile supply chain network. For demonstrating the ability of reconfiguring a flexible system economically, quickly, and environmentally responsible, an Internet-based distributed manufacturing system is considered in this model. Sha and Che [18] offered a way to make manufacturing systems more agile and competitive. They developed a partner selection and production-distribution planning model based on the analytic hierarchy process (AHP), multiattribute utility theory (MAUT), and integer programming (IP), for the virtual integration. They used branch and bound for solving this model. Wu and Barnes [19] presented a review on agile supply chain partner decision making. They concentrated on particular attention given to those methods that are especially relevant for use in agile supply chains.Pan and Nagi [20] proposed a supply chain design problem with uncertain demand in an agile manufacturing system. They considered the integrated optimization of logistics and production costs. The main objective is to select one company in each echelon. A scenario-based approach is used to handle the uncertainty of demand. The formulation is a robust optimization model with considering the expected total costs, variable cost of uncertain demand, and the penalty of/for unmet demands. The heuristic approach is used for solving this problem. Bachlaus et al. [21] proposed an integrated multi-echelon supply chain network with considering agility. This paper presented a five-echelon supply chain model consisting of suppliers, plants, distribution centers, cross-docks, and customer zones and integrated production, distribution and logistics activities. This model is formulated as a multiobjective optimization model that minimizes the fixed and variable costs and maximizes the plant flexibility and volume flexibility. They used a weighted-sum approach for solving their multiobjective problem while this method always caused suboptimality. They proposed HTPSO algorithm for solving the model.Cochran and Marquez Uribe [22] presented a modified set covering formulation for capacity planning and equipment selection in supply chains. In this model, use/need of multi-functional equipment to deal with changing demand is considered. Usage of this equipment is costly, but it is necessary to compete with widespread technological developments among competitors and have flexibility in meeting varying supply chain demands. Costantino et al. [23] considered the manufacturing system (MSC) composing of different stages in supply chain planning. They proposed a technique for the strategic management of the supply chain planning and allowed to have an ability of reconfiguration of the chain for the improvement of the MSC agility. They considered the supplier capacity constraints and the level of demands and maximum production capacities with single/multiple sourcing.Pan and Nagi [24] presented a supply chain network design problem in an agile manufacturing system. They considered a multiple echelons and multiple period model under the situation of multiple customers that have heavy demands. The strategic and tactical decisions are integrated. The main goal of this model is the selection of one or more companies in each echelon, production, inventory, and transportation with minimizing the total operational costs including fixed costs, production, raw material holding, finished products holding, and transportation costs under production and transportation capacity limits. They used Lagrangian heuristic for solving this paper.The agile supply chain that is presented in this paper consists of cross-dock site at the location between distribution and customer zones. This model integrated a production, distribution, and logistics and formulated as a multiobjective optimization model. Due to the complexity of the proposed integrated model, a well-known multiobjective metaheuristic algorithm, namely, multiobjective parallel simulated annealing (MOPSA) is developed for obtaining Pareto solutions of the model. Two other multiobjective metaheuristic algorithms, called NSGA-II and PAES, are considered in order to show the higher performance of the proposed MOPSA.This paper has been prepared with the following innovations:(i)
designing a fuzzy integrated agile supply chain at the strategic level by integrating production, distribution, and logistics across the whole chain,(ii)
presenting a multiobjective optimization model with considering objectives of fixed and variable costs and volume flexibility,(iii)
developing a new multiobjective metaheuristic algorithm, namely, MOPSA,(iv)
extracting nondominated Pareto front solutions of the multiobjective model.


## 3. Fuzzy Multiobjective Mathematical Model

In this section, an integrated agile supply chain network model has been mathematically formulated. The model represents the concept of agility in the supply chain network that decides its ability for rapid responding to the customer demands. The proposed model consists of 5 echelons: suppliers, plants, distribution centers, cross-docks, and customer zones, and these echelons are linked together with a forward flow of material and backward flow of information. The aim of this paper is to integrate the production, distribution and logistics activities at minimum cost, maximum plant flexibility, and volume flexibility. This model optimizes the material flow and determines the optimal number of suppliers, plants, and distribution centers and assignments of cross-docks to customer zones.The transportation costs and capacities are modeled as fuzzy parameters in the proposed formulation. The following notation is used in the formulation of the agile supply chain network design (ASCND) model.Indices
s
=
1
,
2
,
…
,
S
:Set of suppliers
i
=
1
,
2
,
…
,
I
: Set of plants
j
=
1
,
2
,
…
,
J
: Set of possible distribution centers
k
=
1
,
2
,
…
,
K
: Set of possible cross-docks
d
=
1
,
2
,
…
,
D
: Set of customer zones
p
=
1
,
2
,
…
,
P
: Product type
r
=
1
,
2
,
…
,
R
: Raw material type.Notations
Q
p
i: Quantity of product p produced at plant i
q
p
i
j: Quantity of product p transported from plant i to distribution center j
q
r
s
i: Quantity of raw material r transported from supplier s to plant i
A
L
k: Agility level of cross-dock k
A
L
j: Agility level of distribution center  j
A
L
i: Agility level of plant  i
A
L
s: Agility level of supplier  s
Cap
k: Capacity of cross-dock k to handle product families
Ca
p
j: Capacity of distribution center j to handle product families
CDF: Cross-dock flexibility
C
m
k: Cost to supply product p from cross-dock k which would be used by the customer zone m
DCF: Distribution center flexibility
D
m
p: Demand from customer zone m for product p
DCVF
m: Distribution volume flexibility of customer zone m
C
k: Fixed operating cost to open cross-dock k
F
i: Fixed cost for plant i
F
j: Fixed cost for distribution center j
L
B
h: Lower bound of high agility performance index (0.8)
U
B
h: Upper bound of high agility performance index (1.0)
P: Total number of products to be manufactured
K: Total number of cross-docks
J: Total number of distribution centers
L
B
l: Lower bound of low agility performance index (0.4)
U
B
l: Upper bound of low agility performance index (0.6)
M: Total number of customer zones
L
B
m: Lower bound of medium agility performance index (0.6)
UB
m: Upper bound of medium agility performance index (0.8)
l
p
i: Minimum production volume for product p at plant i
l
j: Minimum throughput at distribution center j
CD: Maximum number of cross-docks to be opened
DC: Maximum number of distribution centers to be opened
U
i
p: Maximum production volume for product p at plant i
U
j: Maximum throughput at distribution center j
norm(·): Represents normalized objective (·)
Z: Final integrated objective function
I: Total number of plants
Cap
i: Production capacity for each plant i
Cap
r
s: Production capacity of supplier s for raw material r
PF: Plant flexibility
PVF
m: Plant volume flexibility
STp
i: Standard units at plant i per unit of product p
R: Total number of raw material to be supplied
S: Total number of suppliers
SF: Supplier flexibility
ST
j
p: Standard units at distribution center j per unit of product p
TC
m: Total cost of m customer zone
C
r
s: Unit cost of raw material r for supplier s
C
p
j: Unit cost of throughput (handling and inventory) for product p at distribution center j
C
p
i: Unit production cost for product pat plant i
U
R
r
p: Utilization rate for each raw material r per unit of product p
C
j
k
p: Unit transportation cost from distribution center j to cross-dock k for product p
C
i
j
p: Unit transportation cost from plant i to distribution center j for product p
C
r
s
i: Unit transportation cost from supplier s to plant i for raw material r.Decision Variables
(1)
X
i
=
{
1
if
plant
i
is
opened
0
otherwise
,
Y
j
=
{
1
if
distribution
j
is
opened
0
otherwise
,
V
k
=
{
1
if
cross-dock
k
is
opened
0
otherwise
,
R
k
j
p
=
{
1
if
cross-dock
k
is
assigned
to
distribution
j
for
product
p
0
otherwise
,
A
m
k
p
=
{
1
if
customer
zone
m
is
assigned
to
cross-dock
k
for
product
p
0
otherwise
,
N
s
=
{
1
if
supplier
s
is
selected
0
otherwise
,
(2)
T
C
m
=
[
∑
i
=
1
I
∑
s
=
1
S
∑
r
=
1
R
(
C
r
s
i
~
+
C
r
s
~
)
·
q
r
s
p
]
+
[
∑
i
I
F
i
·
x
i
+
∑
p
P
∑
i
I
C
p
i
~
·
q
p
i
+
∑
p
P
∑
i
I
∑
j
J
C
p
i
j
~
·
q
p
i
j
]
+
[
∑
k
K
F
k
·
y
k
+
∑
p
P
∑
j
J
∑
k
K
C
p
j
k
~
·
R
k
j
p
+
∑
m
M
∑
p
P
∑
k
K
∑
j
J
C
p
j
~
·
D
m
p
·
R
k
j
p
·
A
m
k
p
]
+
[
∑
j
J
F
k
·
v
k
+
∑
m
M
∑
p
P
∑
k
K
C
m
k
~
·
A
m
k
p
]
.
The objective function (2) minimizes the total cost including fixed and variable costs, and are included four components, that is, (a) purchasing cost of raw material and transportation cost from suppliers to plants, (b) the fixed and variable costs of plant operations and transportation cost, (c) the variable cost of handling and inventory of products at DCs, and transportation cost to ship from plant to DCs, and (d) the transportation cost to ship products from DCs to cross-docks and distribution costs to transfer product families from cross-docks to customer zones:
(3)SF=∑p=1P∑r=1R∑s=1S(Cap~rs·Ns-D~mp)∀m,(4)PF=∑pP∑iI(Cap~i·xi-D~mp)∀m,(5)DCF=∑pP∑jJ(Cap~j·yj-D~mp)∀m,(6)CDF=∑pP∑kK(Cap~k·Vk-D~mp)∀m.
The flexibility objective is divided into two types: distribution volume flexibility and delivery flexibility. Specifically, distribution volume flexibility can be defined as the minimum of the difference between the capacity of suppliers, plants, distribution centers, and cross-docks with total demand of products. In (3), supplier flexibility is equal to the difference between supplier’s capacities with total demand of products ordered by customer at zone m. In (4), plant flexibility is the difference between plant’s capacities with total demand of products ordered by customer at zone m, and distribution center flexibility and cross-dock flexibility are defined by the same definitions:
(7)DVFm=min(SF,PF,DCF,CDF)∀m,(8)PVFm=∑iI(xi·Cap~i-∑pPqpi·stpi).
In (7), distribution volume flexibility is defined as the minimum of supplier flexibility, plant flexibility, distribution center flexibility, and cross-dock flexibility. Equation (8) describes plant volume flexibility as the difference between utilization plant capacities with available plant capacity. Finally, the three objectives of the model are as follows:
(9)minZ1=∑m=1MTCm,maxZ2=∑m=1MDVFm,maxZ3=∑m=1MPVFm.Because of the different dominance of cost and flexibility, the normalized approach is considered to integrate both objectives together. The general objective function consists of normalized cost, distribution volume flexibility, and plant flexibility:(10)∑iIqrsi≤Cap~rs∀r,s,(11)∑sSURrp·qrsi≤∑sSCap~rs∀r,p,i,(12)∑pPstpi·qpi≤cap~i·xi∀i,(13)lpi·xi≤qpi≤upi·xi∀i,p,(14)lj·yj≤∑k∑pstpj·D~mp·Rkjp·Amkp≤uj·yj∀j,(15)∑kKAmkp=1∀m,p,(16)qpi=∑jJqpij∀p,i,(17)∑iI∑jJqpij=∑mMD~mp∀p,(18)∑iIqpij=∑mMD~mp·Rkjp·Amkp∀p,j,(19)qpij,qpi,qrsi≥0,∀r,p,i,j,k,(20)∑jJRkjp=Vk∀k,p,(21)∑mM∑pPd~mp·Amkp≤cap~k∀k,(22)∑mM∑pPRkjp·cap~k≤cap~j∀j,(23)Amkp≤Vk∀m,p,k,(24)Rkjp≤yj∀p,k,j,(25)∑kKVk≤DC,(26)∑jJyj≤CD,(27)LBl≤ALs+ALi+ALj+ALk≤UBl,(28)LBm≤ALs+ALi+ALj+ALk≤UBm,(29)LBu≤ALs+ALi+ALj+ALk≤UBu.Constraint (10) ensures that the total number of raw material shipped to plants by each supplier should not exceed the supplier capacity. Constraint (11) checks the raw material supplied by the supplier meets the production requirements. Constraint (12) checks whether the plant capacity, and the range of the quantity of products to be produced at a plant is considered in Constraint (13). Constraint (14) ensures the minimum and maximum capacities for distribution centers and the customer zones that are assigned to the open DCs in this limited area. Constraint (15) enforces each customer zone to be assigned to exactly one cross-dock. Constraint (16) ensures that the amount of products transfer from plant to distribution centers is equal to the quantity of products that produced at that plant. Constraint (17) ensures that all demands are satisfied and constraint (18) checks the demand requirements at each distribution center. Constraint (19) enforces the nonnegativity for all variables. Constraint (20) ensures that cross-docks should be assigned to only open distribution centers. Constraints (21) and (22) ensure the capacity restriction for distribution centers’ cross-docks. Constraint (23) imposes that open cross-dock can satisfy customer demand for products. Constraint (24) ensures that there will be a flow from cross-dock to open distribution centers if the distribution centers are opened. Because of limited capital availability, there is an upper bound for the number of DCs and CDs: a maximum number of cross-docks and distribution centers to be opened. Constrains (25) and (26) ensure that the number of opening cross-docks and distribution centers should not exceed the maximum number of them. Constraints (27), (28), and (29) are modeled to decide the feasibility of agile supply chain. The agility is categorized into three classes: high, medium, and low. Constraint (27) ensures that the agility of selected supplier, plant, distribution centers, and cross-docks should be lesser than the upper bound and greater than its lower bound of low-level agility. Constraints (28) and (29) are similarly defined for medium- and high-level agility performance.


## 4. Defuzzification Procedure

Jiménez et al. [25] proposed a method to convert the possibilistic mixed integer programming model into an equivalent auxiliary crisp model. Jiménez et al.’s [25] method has some advantages as follows.(i)
This is an efficient method for solving the fuzzy linear problems because it keeps linearity and the number of the objective functions does not increase.(ii)
This is a general ranking method which can be used in ranking method which can be applied to different kinds of membership functions such as triangular, trapezoidal, and nonlinear ones in both symmetric and asymmetric forms.(iii)
This method has strong mathematical concepts such as expected interval and expected value of fuzzy numbersIn the following, some concepts that are used in the Jimenez method are initially defined, and the defuzzification procedure is clearly delineated.


### 4.1. Concepts Definition

A fuzzy number is a fuzzy seta~ on the real line R whose membership function μa~ is semicontinuous, such that
(30)r=μa~(x)={0∀x∈(-∞,a1)fa(x)increasingon[a1,a2]1∀x∈[a2,a3]ga(x)decreasingon[a3,a4]0∀x∈(a4,∞).Anr-cut of a fuzzy number a~ is defined by ar={x∈S∣μa~(x)≥r} where S is a universe. Since μa~ is continuous semi-continuous, the r-cuts are closed and bounded intervals, and we represent them by ar=[fa-1(r),ga-1(r)]; the expected interval of a fuzzy number a~, noted EI (a~), can be defined by following equation:
(31)EI(a~)=[E1a,E2a]=[∫01fa-1(r),∫01ga-1(r)].
The expected value of a fuzzy number a~ noted EV(a~) is defined as the half-point of its expected interval:
(32)EV(a~)=E1a+E2a2.Consideringa~, b~ as a fuzzy numbers and γ,δ as nonnegative real numbers, the principle to aggregate fuzzy numbers can be computed by the following equation:
(33)EI(δa~+γb~)=δEI(a~)+γEI(b~),EV(δa~+γb~)=δEV(a~)+γEV(b~).Definition 1.
For any pair of fuzzy numbersa~,b~, the degree in which considering a~ is bigger thanb~ is the following:
(34)μM(a~,b~)={0ifE2a-E1b<0E2a-E1bE2a-E1b-(E1a-E2b)if0∈[E1a-E2b]1ifE2a-E1b>0.
When μM(a~,b~)≥∝, it is said that a~is bigger than or equal to b~ at least in a degree ∝ and is represented by a~≥b~∝.For fuzzy linear constraintai~x≥bi~, a decision vector x∈R is said to be feasible if
(35)min{μM(ai~,bi~)}=∝,
where ai~=(ai1~,ai2~,…,ain~).That is to say(36)ai~x≥b~i∝,i=1,…,m.
For fuzzy linear constraint ai~x=bi~, a decision vector x∈R is said to be feasible if
(37)a~≥b~∝/2,a~≤b~∝/2.
These equations can be rewritten as follows:
(38)∝2≤μM(a~,b~)≤1-∝2.


### 4.2. Defuzzified Model

In this method, each individual fuzzy parameter that exists in each objective function is replaced by itsexpected value, and for defuzzifying constraints, the concept of ∝ feasibility is employed:
(39)minZ1=∑m=1M{[∑jJFk.vk+∑mM∑pP∑kKEV(Cmk)·Amkp][∑i=1I∑s=1S∑r=1R(EV(Crsi)+EV(Crs))·qrsp]+[∑pP∑iI∑jJEV(Cpij)·qpij∑iIFi·xi+∑pP∑iIEV(Cpi)·qpi+∑pP∑iI∑jJEV(Cpij)·qpij]+[∑pP∑kK∑jJEV(Cpj)·EV(Dmp)·Rkjp·Amkp∑kKFk·yk+∑pP∑jJ∑kKEV(Cpjk)·Rkjp+∑mM∑pP∑kK∑jJEV(Cpj)·EV(Dmp)·Rkjp·Amkp∑mM∑pP∑kK∑jJEV(Cpj)]+[∑jJFk.vk+∑mM∑pP∑kKEV(Cmk)·Amkp]}maxZ2=∑m=1MDVFmmaxZ1=∑m=1MPVFms.t.:SF≤∑p=1P∑r=1R∑s=1S([α2×E2D~mp+(1-α2)×E1D~mp]Caprs·Ns-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,SF≥∑p=1P∑r=1R∑s=1S([(1-α2)×E2D~mp+α2×E1D~mp]Caprs·Ns-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,PF≤∑pP∑iI((1-α2)×E1D~mpCapi·xi-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,PF≥∑pP∑iI((1-α2)Capi·xi-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,DCF≤∑pP∑jJ((1-α2)Capj·yj-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,DCF≥∑pP∑jJ([(1-α2)×E2D~mp+α2×E1D~mp]Capj·yj-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,CDF≤∑pP∑kK((1-α2)Capk·Vk-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,CDF≥∑pP∑kK((1-α2)Capk·Vk-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,DVFm=min(SF,PF,DCF,CDF)∀m,PVFm≤∑iI(∑pPqpi·stpixi·[(1-α2)×E2Capi+α2×E1Capi]-∑pPqpi·stpi),PVFm≥∑iI(∑pPqpi·stpixi·[α2×E2Capi+(1-α2)×E1Capi]-∑pP(qpi·stpi)),∑iIqrsi≤[(1-α)×E2Caprs+α×E1Caprs]∀r,s,∑sSURrp·qrsi≤∑sS[(1-α)×E2Caprs+α×E1Caprs]∀r,p,i,∑pPstpi·qpi≤[(1-α)×E2Caprs+α×E1Caprs]·xi∀i,lpi·xi≤qpi≤upi·xi∀i,p,lj·yj≤∑kK∑pPstpj·[(1-α)×E2Dm+α×E1Dm]·Rkjp·Amkp∀j,∑kK∑pPstpj·[α×E2Dm+(1-α)×E1Dm]·Rkjp·Amkp≤uj·yj∀j,∑kKAmkp=1∀m,p,qpi=∑jJqpij∀p,i,∑iI∑jJqpij≤∑mM[(1-α2)×E2Dm+α2×E1Dm]∀p,∑iI∑jJqpij≥∑mM[α2×E2Dm+(1-α2)×E1Dm]∀p,∑iIqpij≤∑mM[(1-α2)×E2Dm+α2×E1Dm]·Rkjp·Amkp∀p,j,∑iIqpij≥∑mM[α2×E2Dm+(1-α2)×E1Dm]·Rkjp·Amkp∀p,j,qpij,qpi,qrsi≥0,∀r,p,i,j,k,∑jJRkjp=Vk∀k,p,∑mM∑pPdmp·Amkp≤[(1-α)×E2Capk+α×E1Capk]∀k,∑mM∑pPRkjp·[α×E2Caprs+(1-α)×E1Caprs]≤[(1-α)×E2Capj+α×E1Capj]∀j,Amkp≤Vk∀m,p,k,Rkjp≤yj∀p,k,j,∑kKVk≤DC,∑jJyj≤CD,LBl≤ALs+ALi+ALj+ALk≤UBl,LBm≤ALs+ALi+ALj+ALk≤UBm,LBu≤ALs+ALi+ALj+ALk≤UBu.


## 4.1. Concepts Definition

A fuzzy number is a fuzzy seta~ on the real line R whose membership function μa~ is semicontinuous, such that
(30)r=μa~(x)={0∀x∈(-∞,a1)fa(x)increasingon[a1,a2]1∀x∈[a2,a3]ga(x)decreasingon[a3,a4]0∀x∈(a4,∞).Anr-cut of a fuzzy number a~ is defined by ar={x∈S∣μa~(x)≥r} where S is a universe. Since μa~ is continuous semi-continuous, the r-cuts are closed and bounded intervals, and we represent them by ar=[fa-1(r),ga-1(r)]; the expected interval of a fuzzy number a~, noted EI (a~), can be defined by following equation:
(31)EI(a~)=[E1a,E2a]=[∫01fa-1(r),∫01ga-1(r)].
The expected value of a fuzzy number a~ noted EV(a~) is defined as the half-point of its expected interval:
(32)EV(a~)=E1a+E2a2.Consideringa~, b~ as a fuzzy numbers and γ,δ as nonnegative real numbers, the principle to aggregate fuzzy numbers can be computed by the following equation:
(33)EI(δa~+γb~)=δEI(a~)+γEI(b~),EV(δa~+γb~)=δEV(a~)+γEV(b~).Definition 1.
For any pair of fuzzy numbersa~,b~, the degree in which considering a~ is bigger thanb~ is the following:
(34)μM(a~,b~)={0ifE2a-E1b<0E2a-E1bE2a-E1b-(E1a-E2b)if0∈[E1a-E2b]1ifE2a-E1b>0.
When μM(a~,b~)≥∝, it is said that a~is bigger than or equal to b~ at least in a degree ∝ and is represented by a~≥b~∝.For fuzzy linear constraintai~x≥bi~, a decision vector x∈R is said to be feasible if
(35)min{μM(ai~,bi~)}=∝,
where ai~=(ai1~,ai2~,…,ain~).That is to say(36)ai~x≥b~i∝,i=1,…,m.
For fuzzy linear constraint ai~x=bi~, a decision vector x∈R is said to be feasible if
(37)a~≥b~∝/2,a~≤b~∝/2.
These equations can be rewritten as follows:
(38)∝2≤μM(a~,b~)≤1-∝2.


## 4.2. Defuzzified Model

In this method, each individual fuzzy parameter that exists in each objective function is replaced by itsexpected value, and for defuzzifying constraints, the concept of ∝ feasibility is employed:
(39)minZ1=∑m=1M{[∑jJFk.vk+∑mM∑pP∑kKEV(Cmk)·Amkp][∑i=1I∑s=1S∑r=1R(EV(Crsi)+EV(Crs))·qrsp]+[∑pP∑iI∑jJEV(Cpij)·qpij∑iIFi·xi+∑pP∑iIEV(Cpi)·qpi+∑pP∑iI∑jJEV(Cpij)·qpij]+[∑pP∑kK∑jJEV(Cpj)·EV(Dmp)·Rkjp·Amkp∑kKFk·yk+∑pP∑jJ∑kKEV(Cpjk)·Rkjp+∑mM∑pP∑kK∑jJEV(Cpj)·EV(Dmp)·Rkjp·Amkp∑mM∑pP∑kK∑jJEV(Cpj)]+[∑jJFk.vk+∑mM∑pP∑kKEV(Cmk)·Amkp]}maxZ2=∑m=1MDVFmmaxZ1=∑m=1MPVFms.t.:SF≤∑p=1P∑r=1R∑s=1S([α2×E2D~mp+(1-α2)×E1D~mp]Caprs·Ns-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,SF≥∑p=1P∑r=1R∑s=1S([(1-α2)×E2D~mp+α2×E1D~mp]Caprs·Ns-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,PF≤∑pP∑iI((1-α2)×E1D~mpCapi·xi-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,PF≥∑pP∑iI((1-α2)Capi·xi-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,DCF≤∑pP∑jJ((1-α2)Capj·yj-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,DCF≥∑pP∑jJ([(1-α2)×E2D~mp+α2×E1D~mp]Capj·yj-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,CDF≤∑pP∑kK((1-α2)Capk·Vk-[α2×E2D~mp+(1-α2)×E1D~mp])∀m,CDF≥∑pP∑kK((1-α2)Capk·Vk-[(1-α2)×E2D~mp+α2×E1D~mp])∀m,DVFm=min(SF,PF,DCF,CDF)∀m,PVFm≤∑iI(∑pPqpi·stpixi·[(1-α2)×E2Capi+α2×E1Capi]-∑pPqpi·stpi),PVFm≥∑iI(∑pPqpi·stpixi·[α2×E2Capi+(1-α2)×E1Capi]-∑pP(qpi·stpi)),∑iIqrsi≤[(1-α)×E2Caprs+α×E1Caprs]∀r,s,∑sSURrp·qrsi≤∑sS[(1-α)×E2Caprs+α×E1Caprs]∀r,p,i,∑pPstpi·qpi≤[(1-α)×E2Caprs+α×E1Caprs]·xi∀i,lpi·xi≤qpi≤upi·xi∀i,p,lj·yj≤∑kK∑pPstpj·[(1-α)×E2Dm+α×E1Dm]·Rkjp·Amkp∀j,∑kK∑pPstpj·[α×E2Dm+(1-α)×E1Dm]·Rkjp·Amkp≤uj·yj∀j,∑kKAmkp=1∀m,p,qpi=∑jJqpij∀p,i,∑iI∑jJqpij≤∑mM[(1-α2)×E2Dm+α2×E1Dm]∀p,∑iI∑jJqpij≥∑mM[α2×E2Dm+(1-α2)×E1Dm]∀p,∑iIqpij≤∑mM[(1-α2)×E2Dm+α2×E1Dm]·Rkjp·Amkp∀p,j,∑iIqpij≥∑mM[α2×E2Dm+(1-α2)×E1Dm]·Rkjp·Amkp∀p,j,qpij,qpi,qrsi≥0,∀r,p,i,j,k,∑jJRkjp=Vk∀k,p,∑mM∑pPdmp·Amkp≤[(1-α)×E2Capk+α×E1Capk]∀k,∑mM∑pPRkjp·[α×E2Caprs+(1-α)×E1Caprs]≤[(1-α)×E2Capj+α×E1Capj]∀j,Amkp≤Vk∀m,p,k,Rkjp≤yj∀p,k,j,∑kKVk≤DC,∑jJyj≤CD,LBl≤ALs+ALi+ALj+ALk≤UBl,LBm≤ALs+ALi+ALj+ALk≤UBm,LBu≤ALs+ALi+ALj+ALk≤UBu.


## 5. Proposed Multiobjective Algorithm

In addition to validating the higher performance and justifing the use of MOPSA in under-hand HLP, the proposed MOPSA was compared with several well-known metaheuristic algorithms existing in the literature (such as NSGA-II, NRGA, MOTS, MOSS, MOIWO, and PAES) while the proposed MOPSA had better performance. However, in this paper, the comparison with the best two of them (i.e., NSGA-II and PAES) is reported. For more details about NSGA-II and PAES, the interested readers are referred to [26, 27].


### 5.1. SA in General

Simulated annealing is a local search-based heuristic that is capable of escaping from being trapped into a local optimum by accepting, with small probability, worse solutions during its iterations. It has been applied successfully to a wide variety of highly complicated combinatorial optimization problems as well as various real-world problems. SA was introduced by Metropolis et al. [28] and used by Kirkpatrick et al. [29] in optimization problems as a first time. The concept of the method is adopted from the “annealing” process used in the metallurgical industry. Annealing is the process by which slow cooling is applied to metals to produce better aligned, low-energy-state crystallization. The optimization procedure of SA searches for a (near) global minimum mimicking the slow cooling procedure in the physical annealing process. It starts from a random initial solution. At each iteration, a new solution is taken from the predefined neighborhood of the current solution.The objective function value of this new solution is then compared with that of the current best solution in order to determine if an improvement has been achieved. If the objective function value of the new solution is better, that is, being smaller in the case of minimization, the new solution becomes the current solution from which the search continues by proceeding with a new iteration. A new solution with a degraded (larger) objective function value may also be accepted as the new current solution, with a small probability determined by the Boltzmann function,exp(-Δ/kT), where Δ is the difference of objective function values between the current solution and the new solution, k is a predetermined constant, and T is the current temperature. The basic idea is not to restrict the search to those solutions that decrease the objective function value but also allow moves that increase the objective function value. This mechanism may avoid the procedure being trapped prematurely in a local minimum.


### 5.2. Proposed MOPSA

In this paper, as a first time in agile supply chain management, we propose a new concept of simulated annealing algorithm adopted in recent studies used in optimization problems, in which, unlike classical SA, more than one solution are used for searching the solution space to obtain Pareto-optimal solutions. So, this kind of SA is called parallel SA (PSA). This PSA is able to search the solution space widely and may also obtain nondominated Pareto solutions in less computational CPU time. Now, the proposed MOPSA will be described in detail.


### 5.3. Solution Representation

One important decision in designing a metaheuristic method is to decide how to represent and relate solutions in an efficient way to the searching space. Representation should be easy to decode to reduce the cost of the algorithm. The integer coding procedure is adopted to encode the solution, where each solution in the algorithm represents the number of decision variables used in the proposed model. For example, if there are 4 suppliers, 5 plants, 6 DCs, 6 CDs, and 12 CZs for 5 raw materials and 5 products and the objective is to designn agile supply chains, then the length of the solution will be the sum of the decision variables used in the n supply chains to be designed. Each cell of the solution represents the integer value corresponding to the supplier, plant, distribution center, cross-dock, and customer zones for specified amount of raw material and products.


### 5.4. Initial Solution

The proposed PSA algorithm, unlike classical SA, is started with more than one initial for searching the solution space. The number of initial solutions, namely,nPop, is an input parameter of the PSA algorithm. The initial solutions are created randomly according to Section 5.3.


### 5.5. Neighborhood Search

We use three consecutive steps for searching the initial solution neighborhood for obtaining the better solution at each iteration. The proposed three steps include mutation, assimilation and crossover that will be described next in detail.


#### 5.5.1. Mutation Step

The first step of the proposed MOPSA is the mutation of the initial solutions to a predefined numbernMutate. For mutating the initial solution, three different procedures are applied including (a) swap, (b) reversion, and (c) inversion mutation as shown in Figure 1. In the swap mutation, the places of two random selected bits are exchanged. In the reversion mutation, a random part of the solution is selected and its permutation is reversed. Finally, in the inversion mutation, one bit is chosen randomly and its value is replaced with new random value.Figure 1
Three steps of mutation.


#### 5.5.2. Assimilation Step

At the end of the mutation step, we havenPop group of mutated solution with a number of nMutate. For example, if nPop and nMutate are equal to 5 and 10, respectively, then we have totally 50 new mutated solutions which are included in 5 groups with 10 solutions in each group. Now, the mutated solutions, called colonies, start moving toward the best solution in each group, called imperialist, by the assimilation operator of imperialist competitive algorithm (ICA) [24] as illustrated in Figure 2, in which d is the distance between colony and the imperialist and x is the amount of moving of colony toward the imperialist which is uniformly distributed between 0 and β×d, in which β is a number greater than 1. For moving colony toward the imperialist in other direction, an θ angle is applied with the uniform distribution between –γ and γ. It has been estimated that π/4 is the most proper value of γ. The best solution in each population group is chosen according to nondominance strategy and crowding distance metric.Figure 2
Moving colonies toward the imperialist with a random angleθ.


#### 5.5.3. Crossover Step

In this step, all the new created solutions in the two previous steps are merged for being combined according to the crossover operator of the genetic algorithm. For applying the crossover on the solutions, three different kinds of crossover are adopted on a pair of solutions, including: (a) one-point crossover, (b) two-point crossover and (c) uniform crossover as they are shown in Figure3. It must be mentioned that operators a and b are applied for both continous and discrete parts, but operator c is applied for both parts differently. In this step, a binary tournament selection is applied for choosing parents.Figure 3
Three steps of crossover.


### 5.6. Next Population

At each iteration, the proposed MOPSA needs initial solutions to apply three main neighborhood search strategies again for obtaining nondominated Pareto solutions. In this step, the nondominance strategy and crowding distance metric are applied for choosing thenPop better solutions among the whole new created solutions till now to compare them with previous nPop initial solutions correspondingly for choosing the next iteration’s initial solutions. Three different cases may occur when two corresponding solutions are be compared in order to being chosen for next iteration initial solutions as follows.(1)
The new solution will be chosen for next iteration if dominates the old one.(2)
If two solutions are not dominated with each other, a solution with higher crowding distance metric will be chosen for next iteration.(3)
If the old solution dominates the new one, the probable acceptance function of simulated annealing will be adopted according to the next section.


### 5.7. Probability Acceptance Function

A new solution with a degraded (larger) objective function value may also be accepted as the new current solution, with a small probability determined by the Boltzmann function,exp(-Δ/kT), where Δ is the difference of objective function values between the current solution and the new solution, k is a predetermined constant, and T is the current temperature. The Boltzmann function and Δ in proposed MOPSA are as follows:
(40)P=e-Δf/T,(41)Δf=|(f1(x)-f1(y))f1(x)+(f2(x)-f2(y))f2(x)|,
where x and y are the new and old solutions, respectively. Also, f1 and f2 are the first and the second objective functions, respectively. Variable T, in the Boltzmann function, is the current temperature that is calculated as follows:
(42)Tn+1=α·Tn,
where α is a constant factor between 0 and 1 for decreasing the temperature at each next iteration. According to above description, a worse solution will be accepted for next iteration when a created random value r is less than probability value P calculated according to (40).


### 5.8. Stopping Criteria

In the proposed MOPSA, the number of function Calls-NFCs is the stopping criteria.The pseudocode of proposed MOPSA is presented in Pseudocode1.Pseudocode 1
NFC← 0
set the parameters of PSA (nPop, nMutate, pCrossover, α, β, T0)
create initial solution←  nPop
terminate←  false
while  (terminate  =  false)  do
mutate each initial solution←nMutate
find the best solution (imperialist)← non-dominated sorting & crowding distance
update NFC
assimilate mutated solutions (colonies) toward imperialist←β
merge whole new created solution
apply crossover←  pCrossover
update NFC
find better new solutions←nPop
update NFC
if  (new solution dominates old solution)  then
accept the new solution
else if  (no one dominate the other one) then
calculate crowding distance-CD of each solution
if  (CD of new solution > CD of old solution)  then
accept the new solution
end if
else
apply probable acceptance function←P=e-Δf/T←Δf & T
create random valuer
if  (r<P)  then
accept new solution
else
accept old solution
end if
end if
update NFC
T←α×T
if  (NFC = predefined value)  then
terminate  =true
end if
end while


## 5.1. SA in General

Simulated annealing is a local search-based heuristic that is capable of escaping from being trapped into a local optimum by accepting, with small probability, worse solutions during its iterations. It has been applied successfully to a wide variety of highly complicated combinatorial optimization problems as well as various real-world problems. SA was introduced by Metropolis et al. [28] and used by Kirkpatrick et al. [29] in optimization problems as a first time. The concept of the method is adopted from the “annealing” process used in the metallurgical industry. Annealing is the process by which slow cooling is applied to metals to produce better aligned, low-energy-state crystallization. The optimization procedure of SA searches for a (near) global minimum mimicking the slow cooling procedure in the physical annealing process. It starts from a random initial solution. At each iteration, a new solution is taken from the predefined neighborhood of the current solution.The objective function value of this new solution is then compared with that of the current best solution in order to determine if an improvement has been achieved. If the objective function value of the new solution is better, that is, being smaller in the case of minimization, the new solution becomes the current solution from which the search continues by proceeding with a new iteration. A new solution with a degraded (larger) objective function value may also be accepted as the new current solution, with a small probability determined by the Boltzmann function,exp(-Δ/kT), where Δ is the difference of objective function values between the current solution and the new solution, k is a predetermined constant, and T is the current temperature. The basic idea is not to restrict the search to those solutions that decrease the objective function value but also allow moves that increase the objective function value. This mechanism may avoid the procedure being trapped prematurely in a local minimum.


## 5.2. Proposed MOPSA

In this paper, as a first time in agile supply chain management, we propose a new concept of simulated annealing algorithm adopted in recent studies used in optimization problems, in which, unlike classical SA, more than one solution are used for searching the solution space to obtain Pareto-optimal solutions. So, this kind of SA is called parallel SA (PSA). This PSA is able to search the solution space widely and may also obtain nondominated Pareto solutions in less computational CPU time. Now, the proposed MOPSA will be described in detail.


## 5.3. Solution Representation

One important decision in designing a metaheuristic method is to decide how to represent and relate solutions in an efficient way to the searching space. Representation should be easy to decode to reduce the cost of the algorithm. The integer coding procedure is adopted to encode the solution, where each solution in the algorithm represents the number of decision variables used in the proposed model. For example, if there are 4 suppliers, 5 plants, 6 DCs, 6 CDs, and 12 CZs for 5 raw materials and 5 products and the objective is to designn agile supply chains, then the length of the solution will be the sum of the decision variables used in the n supply chains to be designed. Each cell of the solution represents the integer value corresponding to the supplier, plant, distribution center, cross-dock, and customer zones for specified amount of raw material and products.


## 5.4. Initial Solution

The proposed PSA algorithm, unlike classical SA, is started with more than one initial for searching the solution space. The number of initial solutions, namely,nPop, is an input parameter of the PSA algorithm. The initial solutions are created randomly according to Section 5.3.


## 5.5. Neighborhood Search

We use three consecutive steps for searching the initial solution neighborhood for obtaining the better solution at each iteration. The proposed three steps include mutation, assimilation and crossover that will be described next in detail.


### 5.5.1. Mutation Step

The first step of the proposed MOPSA is the mutation of the initial solutions to a predefined numbernMutate. For mutating the initial solution, three different procedures are applied including (a) swap, (b) reversion, and (c) inversion mutation as shown in Figure 1. In the swap mutation, the places of two random selected bits are exchanged. In the reversion mutation, a random part of the solution is selected and its permutation is reversed. Finally, in the inversion mutation, one bit is chosen randomly and its value is replaced with new random value.Figure 1
Three steps of mutation.


### 5.5.2. Assimilation Step

At the end of the mutation step, we havenPop group of mutated solution with a number of nMutate. For example, if nPop and nMutate are equal to 5 and 10, respectively, then we have totally 50 new mutated solutions which are included in 5 groups with 10 solutions in each group. Now, the mutated solutions, called colonies, start moving toward the best solution in each group, called imperialist, by the assimilation operator of imperialist competitive algorithm (ICA) [24] as illustrated in Figure 2, in which d is the distance between colony and the imperialist and x is the amount of moving of colony toward the imperialist which is uniformly distributed between 0 and β×d, in which β is a number greater than 1. For moving colony toward the imperialist in other direction, an θ angle is applied with the uniform distribution between –γ and γ. It has been estimated that π/4 is the most proper value of γ. The best solution in each population group is chosen according to nondominance strategy and crowding distance metric.Figure 2
Moving colonies toward the imperialist with a random angleθ.


### 5.5.3. Crossover Step

In this step, all the new created solutions in the two previous steps are merged for being combined according to the crossover operator of the genetic algorithm. For applying the crossover on the solutions, three different kinds of crossover are adopted on a pair of solutions, including: (a) one-point crossover, (b) two-point crossover and (c) uniform crossover as they are shown in Figure3. It must be mentioned that operators a and b are applied for both continous and discrete parts, but operator c is applied for both parts differently. In this step, a binary tournament selection is applied for choosing parents.Figure 3
Three steps of crossover.


## 5.5.1. Mutation Step

The first step of the proposed MOPSA is the mutation of the initial solutions to a predefined numbernMutate. For mutating the initial solution, three different procedures are applied including (a) swap, (b) reversion, and (c) inversion mutation as shown in Figure 1. In the swap mutation, the places of two random selected bits are exchanged. In the reversion mutation, a random part of the solution is selected and its permutation is reversed. Finally, in the inversion mutation, one bit is chosen randomly and its value is replaced with new random value.Figure 1
Three steps of mutation.


## 5.5.2. Assimilation Step

At the end of the mutation step, we havenPop group of mutated solution with a number of nMutate. For example, if nPop and nMutate are equal to 5 and 10, respectively, then we have totally 50 new mutated solutions which are included in 5 groups with 10 solutions in each group. Now, the mutated solutions, called colonies, start moving toward the best solution in each group, called imperialist, by the assimilation operator of imperialist competitive algorithm (ICA) [24] as illustrated in Figure 2, in which d is the distance between colony and the imperialist and x is the amount of moving of colony toward the imperialist which is uniformly distributed between 0 and β×d, in which β is a number greater than 1. For moving colony toward the imperialist in other direction, an θ angle is applied with the uniform distribution between –γ and γ. It has been estimated that π/4 is the most proper value of γ. The best solution in each population group is chosen according to nondominance strategy and crowding distance metric.Figure 2
Moving colonies toward the imperialist with a random angleθ.


## 5.5.3. Crossover Step

In this step, all the new created solutions in the two previous steps are merged for being combined according to the crossover operator of the genetic algorithm. For applying the crossover on the solutions, three different kinds of crossover are adopted on a pair of solutions, including: (a) one-point crossover, (b) two-point crossover and (c) uniform crossover as they are shown in Figure3. It must be mentioned that operators a and b are applied for both continous and discrete parts, but operator c is applied for both parts differently. In this step, a binary tournament selection is applied for choosing parents.Figure 3
Three steps of crossover.


## 5.6. Next Population

At each iteration, the proposed MOPSA needs initial solutions to apply three main neighborhood search strategies again for obtaining nondominated Pareto solutions. In this step, the nondominance strategy and crowding distance metric are applied for choosing thenPop better solutions among the whole new created solutions till now to compare them with previous nPop initial solutions correspondingly for choosing the next iteration’s initial solutions. Three different cases may occur when two corresponding solutions are be compared in order to being chosen for next iteration initial solutions as follows.(1)
The new solution will be chosen for next iteration if dominates the old one.(2)
If two solutions are not dominated with each other, a solution with higher crowding distance metric will be chosen for next iteration.(3)
If the old solution dominates the new one, the probable acceptance function of simulated annealing will be adopted according to the next section.


## 5.7. Probability Acceptance Function

A new solution with a degraded (larger) objective function value may also be accepted as the new current solution, with a small probability determined by the Boltzmann function,exp(-Δ/kT), where Δ is the difference of objective function values between the current solution and the new solution, k is a predetermined constant, and T is the current temperature. The Boltzmann function and Δ in proposed MOPSA are as follows:
(40)P=e-Δf/T,(41)Δf=|(f1(x)-f1(y))f1(x)+(f2(x)-f2(y))f2(x)|,
where x and y are the new and old solutions, respectively. Also, f1 and f2 are the first and the second objective functions, respectively. Variable T, in the Boltzmann function, is the current temperature that is calculated as follows:
(42)Tn+1=α·Tn,
where α is a constant factor between 0 and 1 for decreasing the temperature at each next iteration. According to above description, a worse solution will be accepted for next iteration when a created random value r is less than probability value P calculated according to (40).


## 5.8. Stopping Criteria

In the proposed MOPSA, the number of function Calls-NFCs is the stopping criteria.The pseudocode of proposed MOPSA is presented in Pseudocode1.Pseudocode 1
NFC← 0
set the parameters of PSA (nPop, nMutate, pCrossover, α, β, T0)
create initial solution←  nPop
terminate←  false
while  (terminate  =  false)  do
mutate each initial solution←nMutate
find the best solution (imperialist)← non-dominated sorting & crowding distance
update NFC
assimilate mutated solutions (colonies) toward imperialist←β
merge whole new created solution
apply crossover←  pCrossover
update NFC
find better new solutions←nPop
update NFC
if  (new solution dominates old solution)  then
accept the new solution
else if  (no one dominate the other one) then
calculate crowding distance-CD of each solution
if  (CD of new solution > CD of old solution)  then
accept the new solution
end if
else
apply probable acceptance function←P=e-Δf/T←Δf & T
create random valuer
if  (r<P)  then
accept new solution
else
accept old solution
end if
end if
update NFC
T←α×T
if  (NFC = predefined value)  then
terminate  =true
end if
end while


## 6. Experimental Results

This section evaluates the performance of the proposed MOPSA. The performance of the proposed MOPSA is compared with well-known multiobjective evolutionary algorithms (MOEAs), namely, NSGA-II and PAES. At first, a brief discussion about the implementation of NSGA-II and PAES is presented.


### 6.1. Parameters Setting

It is well known that the quality of an algorithm is significantly influenced by the values of its parameters. In this section, for optimizing the behavior of the proposed algorithms, appropriate tuning of their parameters has been carried out. For this purpose, response surface methodology (RSM) is employed. RSM is defined as a collection of mathematical and statistical method-based experiential, which can be used to optimize processes. Regression equation analysis is used to evaluate the response surface model. First of all, parameters of each algorithm that statistically have significant impact on the algorithms are recognized. To select the values that result in solutions with high quality, we consider problems in two different sizes including Small-S and Large-L sizes. To identify significant parameters, two levels for each parameter are considered. Each factor is measured at two levels, which can be coded as −1 when the factor is at its low level (L) and +1 when the factor is at its high level (H). The coded variable can be defined as follows:
(43)Xi=ri-((h+l)/2)((h-l)/2),
where Xi and ri are coded variable and natural variable, respectively. h and l represent high level and low level of factor. Factors and their levels are shown in Table 1.Table 1
Parameters and their levels for small and large sizes.
Factor
T
0
α
n
M
u
t
a
t
e
n
P
o
p
p
C
r
o
s
s
o
v
e
r
β
Problem size
S
L
S
L
S
L
S
L
S
L
S
L
Lower limit
5
5
0.75
0.75
5
8
3
5
0.4
0.4
1
1
Upper limit
15
15
0.95
0.95
15
20
8
12
0.8
0.8
3
3After developing regression models for each problem size separately, tuned parameters of proposed MOPSA have been shown in Table2.Table 2
Tuned parameters of proposed MOPSA.
Factor
T
0
α
n
M
u
t
a
t
e
n
P
o
p
p
C
r
o
s
s
o
v
e
r
β
Problem size
S
L
S
L
S
L
S
L
S
L
S
L
Tuned value
10
13
0.84
0.91
10
16
5
6
0.5
0.7
1.8
2For the proposed MOPSA, the NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.(i)
NSGA-II assumptions and parameters’ value.
(a)
The initial population is randomly generated.(b)
One of these operators, namely, one-point crossover, two-point crossover, and uniform crossover, is selected randomly as a crossover operator.(c)
One of these operators, namely, inversion, swap, and reversion, is selected randomly as mutation operator.(d)
The crossover and mutation ratios are set to 0.8 and 0.2, respectively, using RSM method.(e)
The number of the initial population is set to 200 and 300 for small- and large-size problems, respectively.(f)
The NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.(ii)
PAES assumptions and parameters’ value.
(a)
The size of archive is 150.(b)
One of these operators, namely, inversion, swap, and reversion, is selected randomly as revolution.(c)
The NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.


### 6.2. Comparison Metrics

To validate the reliability of the proposed MOPSA, the following four comparison metrics are taken into account.(1)
Quality metrics: this metric is simply measured by putting together the nondominated solutions found by algorithms, and the ratios between nondominated solutions are reported.(2)
Mean ideal distance (MID): the closeness between Pareto solutions and ideal point (f1best,f2best) is determined by using MID. The equation of MID is defined by
(44)MID=(∑i=1n(f1i-f1bestf1,totalmax-f1,totalmin)2+(f2i-f2bestf2,totalmax-f2,totalmin)2)×(n)-1,wheren is the number of nondominated solutions and fi,totalmax and fi,totalmin are the maximum and minimum values of each fitness functions among the all nondominated solutions obtained by the algorithms, respectively. Regarding this definition, an algorithm with a lower value of MID has a better performance.(3)
Diversification metric (DM): this metric measures the spread of the Pareto solutions set and is calculated by(45)DM=(maxf1i-minf1if1.totalmax-f1.totalmin)2+(maxf2i-minf2if2.totalmax-f2.totalmin)2.Regarding this definition, an algorithm with a higher value of DM has a better performance.(4)
Spacing metric (SM): this metric measures the uniformity of the spread of the nondominated set solutions. This metric is defined by(46)SM=∑i=1n-1|d--di|(n-1)d-,where di is the Euclidean distance between consecutive solutions in the obtained nondominated set of solutions and d- is the average of these distances. Regarding this definition, an algorithm with a lower value of SM has a better performance.


### 6.3. Data Generation

In order to evaluate the performance of the proposed MOPSA, several test problems are generated in this subsection. The number of potential locations for all five echelons is distributed randomly between 5 and 100 nodes. Other parameters are presented as in Table3. It should be mentioned that these ranges of data are obtained from real steel industry, and some of the ranges are also generated randomly because of lack of information.Table 3
Data generation and their distribution.
Parameters
AL
k
AL
j
AL
i
AL
s
Cap
k
Cap
j
C
m
k
D
m
p
C
k
F
i
F
j
Level
~(0.05, 0.3)
~(100, 500)
~(100, 400)
~(10, 20)
~(20, 50)
~(103, 104)
~(105, 106)
~(103, 104)
Parameters
l
p
i
l
j
U
i
p
U
j
Cap
i
Cap
r
s
ST
p
i
ST
j
p
C
r
s
Level
~(100, 200)
~(100, 300)
~(300, 500)
~(250, 450)
~(300, 500)
~(200, 500)
~(5, 10)
~(5, 10)
~(10, 20)
Parameters
C
p
j
C
p
i
UR
r
p
C
j
k
p
C
i
j
p
C
r
s
i
Level
~(10, 20)
~(20, 40)
~(0.6, 1)
~(10, 20)
~(10, 20)
~(10, 20)


### 6.4. Computational Results

The proposed MOPSA is applied to a number of test problems, and its performance is compared with NSGA-II and PAES. Table4 shows the comparison between the proposed MOPSA algorithm and NSGA-II and PAES regarding QM and SM comparison metrics. Besides, Table 5 shows the results of DM and MID comparison metrics. Both Tables 4 to 5 show that the proposed MOPSA is superior to NSGA-II and PAES in each test problem.(i)
The proposed MOPSA can achieve a greater number of Pareto-optimal solutions with higher qualities than both NSGA-II and PAES.(ii)
The proposed MOPSA provides nondominated solutions that have less average values of the spacing metric.(iii)
These data reveal that nondominated solutions obtained by the proposed MOPSA are more uniformly distributed in comparison with NSGA-II.(iv)
The average values of the diversification metric in our proposed PSO are considerably greater than those of NSGA-II (i.e., PSO finds nondominated solutions that diverse more).(v)
In most of the test problems, the values of MID in the proposed MOPSA are smaller than those of NSGA-II and PAES.(vi)
The computational time of the proposed MOPSA was less than 10 minutes even for the largest test problem.Table 4
Quality and spacing metrics.
Problem number
Quality metric (QM)
Spacing metric (SM)
PAES
NSGA-II
MOPSA
NSGA-II
PAES
MOPSA
1
0
0.08
0.92
0.558
0.879
0.444
2
0.09
0.18
0.73
0.521
0.937
0.315
3
0
0.20
0.80
0.409
0.810
0.385
4
0
0
1
0.395
0.762
0.373
5
0
0.10
0.90
0.608
0.726
0.359
6
0
0.00
1
0.359
0.932
0.282
7
0
0.15
0.85
0.608
0.761
0.274
8
0
0
1
0.595
0.816
0.335
9
0
0.25
0.75
0.870
0.866
0.303
10
0
0.10
0.90
0.672
0.768
0.243
11
0
0
1
0.574
0.893
0.586
12
0
0
1
0.862
0.845
0.303
13
0
0
1
0.674
0.745
0.394
14
0.08
0.17
0.75
0.624
0.935
0.476
15
0.06
0.2
0.74
0.870
0.730
0.314
16
0
0
1
0.517
0.788
0.521
17
0
0
1
0.657
0.771
0.442
18
0
0.12
0.88
0.601
0.859
0.276
19
0
0
1
0.552
0.727
0.486
20
0
0
1
0.562
0.895
0.610
21
0
0
1
0.772
0.819
0.636
22
0
0
1
0.703
0.812
0.481
23
0
0.04
0.96
0.641
0.835
0.423
24
0
0.10
0.90
0.774
0.686
0.359
25
0
0
1
0.531
0.808
0.224
26
0
0
1
0.724
0.767
0.512
27
0
0.05
0.95
0.721
0.866
0.335
28
0
0
1
0.669
0.794
0.572
29
0
0
1
0.555
0.527
0.303
30
0
0
1
0.621
0.862
0.394Table 5
Diversification and mean ideal distance metrics.
Problem number
Diversity metric (DM)
Mean ideal distance (MID)
PAES
NSGA-II
MOPSA
PAES
NSGA-II
MOPSA
1
0.929
0.872
1.379
0.862
0.569
0.255
2
0.606
1.065
1.143
0.698
0.673
0.258
3
0.617
0.772
1.300
0.795
0.743
0.302
4
0.667
0.805
0.971
0.893
0.621
0.480
5
0.860
0.758
1.111
0.622
0.797
0.317
6
0.893
0.754
1.358
0.776
0.527
0.282
7
0.859
1.048
1.296
0.724
0.596
0.245
8
0.780
0.932
1.380
0.793
0.653
0.519
9
0.819
0.920
1.228
0.679
0.518
0.312
10
0.718
0.758
0.917
0.828
0.718
0.239
11
0.898
1.042
1.325
0.899
0.667
0.430
12
0.675
0.949
1.367
0.656
0.659
0.227
13
0.875
0.840
1.240
0.835
0.749
0.385
14
0.673
0.905
1.279
0.658
0.758
0.203
15
0.771
1.033
1.212
0.898
0.737
0.372
16
0.793
1.022
1.240
0.673
0.606
0.447
17
0.648
0.724
1.098
0.895
0.655
0.282
18
0.836
0.860
1.084
0.638
0.636
0.245
19
0.690
0.911
1.394
0.669
0.481
0.319
20
0.754
0.867
0.918
0.607
0.468
0.312
21
0.833
0.963
1.343
0.782
0.496
0.239
22
0.700
0.951
1.357
0.633
0.649
0.430
23
0.716
0.817
1.298
0.722
0.647
0.227
24
0.847
0.873
0.949
0.866
0.571
0.385
25
0.706
0.706
1.031
0.764
0.572
0.203
26
0.930
1.094
1.068
0.711
0.486
0.372
27
0.994
0.767
1.240
0.662
0.610
0.437
28
0.892
0.742
0.968
0.732
0.639
0.270
29
0.737
0.849
1.261
0.887
0.532
0.334
30
0.834
0.779
0.953
0.637
0.534
0.471


## 6.1. Parameters Setting

It is well known that the quality of an algorithm is significantly influenced by the values of its parameters. In this section, for optimizing the behavior of the proposed algorithms, appropriate tuning of their parameters has been carried out. For this purpose, response surface methodology (RSM) is employed. RSM is defined as a collection of mathematical and statistical method-based experiential, which can be used to optimize processes. Regression equation analysis is used to evaluate the response surface model. First of all, parameters of each algorithm that statistically have significant impact on the algorithms are recognized. To select the values that result in solutions with high quality, we consider problems in two different sizes including Small-S and Large-L sizes. To identify significant parameters, two levels for each parameter are considered. Each factor is measured at two levels, which can be coded as −1 when the factor is at its low level (L) and +1 when the factor is at its high level (H). The coded variable can be defined as follows:
(43)Xi=ri-((h+l)/2)((h-l)/2),
where Xi and ri are coded variable and natural variable, respectively. h and l represent high level and low level of factor. Factors and their levels are shown in Table 1.Table 1
Parameters and their levels for small and large sizes.
Factor
T
0
α
n
M
u
t
a
t
e
n
P
o
p
p
C
r
o
s
s
o
v
e
r
β
Problem size
S
L
S
L
S
L
S
L
S
L
S
L
Lower limit
5
5
0.75
0.75
5
8
3
5
0.4
0.4
1
1
Upper limit
15
15
0.95
0.95
15
20
8
12
0.8
0.8
3
3After developing regression models for each problem size separately, tuned parameters of proposed MOPSA have been shown in Table2.Table 2
Tuned parameters of proposed MOPSA.
Factor
T
0
α
n
M
u
t
a
t
e
n
P
o
p
p
C
r
o
s
s
o
v
e
r
β
Problem size
S
L
S
L
S
L
S
L
S
L
S
L
Tuned value
10
13
0.84
0.91
10
16
5
6
0.5
0.7
1.8
2For the proposed MOPSA, the NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.(i)
NSGA-II assumptions and parameters’ value.
(a)
The initial population is randomly generated.(b)
One of these operators, namely, one-point crossover, two-point crossover, and uniform crossover, is selected randomly as a crossover operator.(c)
One of these operators, namely, inversion, swap, and reversion, is selected randomly as mutation operator.(d)
The crossover and mutation ratios are set to 0.8 and 0.2, respectively, using RSM method.(e)
The number of the initial population is set to 200 and 300 for small- and large-size problems, respectively.(f)
The NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.(ii)
PAES assumptions and parameters’ value.
(a)
The size of archive is 150.(b)
One of these operators, namely, inversion, swap, and reversion, is selected randomly as revolution.(c)
The NFCs stopping criteria were set on 30000 and 100000 for small- and large-size problems, respectively.


## 6.2. Comparison Metrics

To validate the reliability of the proposed MOPSA, the following four comparison metrics are taken into account.(1)
Quality metrics: this metric is simply measured by putting together the nondominated solutions found by algorithms, and the ratios between nondominated solutions are reported.(2)
Mean ideal distance (MID): the closeness between Pareto solutions and ideal point (f1best,f2best) is determined by using MID. The equation of MID is defined by
(44)MID=(∑i=1n(f1i-f1bestf1,totalmax-f1,totalmin)2+(f2i-f2bestf2,totalmax-f2,totalmin)2)×(n)-1,wheren is the number of nondominated solutions and fi,totalmax and fi,totalmin are the maximum and minimum values of each fitness functions among the all nondominated solutions obtained by the algorithms, respectively. Regarding this definition, an algorithm with a lower value of MID has a better performance.(3)
Diversification metric (DM): this metric measures the spread of the Pareto solutions set and is calculated by(45)DM=(maxf1i-minf1if1.totalmax-f1.totalmin)2+(maxf2i-minf2if2.totalmax-f2.totalmin)2.Regarding this definition, an algorithm with a higher value of DM has a better performance.(4)
Spacing metric (SM): this metric measures the uniformity of the spread of the nondominated set solutions. This metric is defined by(46)SM=∑i=1n-1|d--di|(n-1)d-,where di is the Euclidean distance between consecutive solutions in the obtained nondominated set of solutions and d- is the average of these distances. Regarding this definition, an algorithm with a lower value of SM has a better performance.


## 6.3. Data Generation

In order to evaluate the performance of the proposed MOPSA, several test problems are generated in this subsection. The number of potential locations for all five echelons is distributed randomly between 5 and 100 nodes. Other parameters are presented as in Table3. It should be mentioned that these ranges of data are obtained from real steel industry, and some of the ranges are also generated randomly because of lack of information.Table 3
Data generation and their distribution.
Parameters
AL
k
AL
j
AL
i
AL
s
Cap
k
Cap
j
C
m
k
D
m
p
C
k
F
i
F
j
Level
~(0.05, 0.3)
~(100, 500)
~(100, 400)
~(10, 20)
~(20, 50)
~(103, 104)
~(105, 106)
~(103, 104)
Parameters
l
p
i
l
j
U
i
p
U
j
Cap
i
Cap
r
s
ST
p
i
ST
j
p
C
r
s
Level
~(100, 200)
~(100, 300)
~(300, 500)
~(250, 450)
~(300, 500)
~(200, 500)
~(5, 10)
~(5, 10)
~(10, 20)
Parameters
C
p
j
C
p
i
UR
r
p
C
j
k
p
C
i
j
p
C
r
s
i
Level
~(10, 20)
~(20, 40)
~(0.6, 1)
~(10, 20)
~(10, 20)
~(10, 20)


## 6.4. Computational Results

The proposed MOPSA is applied to a number of test problems, and its performance is compared with NSGA-II and PAES. Table4 shows the comparison between the proposed MOPSA algorithm and NSGA-II and PAES regarding QM and SM comparison metrics. Besides, Table 5 shows the results of DM and MID comparison metrics. Both Tables 4 to 5 show that the proposed MOPSA is superior to NSGA-II and PAES in each test problem.(i)
The proposed MOPSA can achieve a greater number of Pareto-optimal solutions with higher qualities than both NSGA-II and PAES.(ii)
The proposed MOPSA provides nondominated solutions that have less average values of the spacing metric.(iii)
These data reveal that nondominated solutions obtained by the proposed MOPSA are more uniformly distributed in comparison with NSGA-II.(iv)
The average values of the diversification metric in our proposed PSO are considerably greater than those of NSGA-II (i.e., PSO finds nondominated solutions that diverse more).(v)
In most of the test problems, the values of MID in the proposed MOPSA are smaller than those of NSGA-II and PAES.(vi)
The computational time of the proposed MOPSA was less than 10 minutes even for the largest test problem.Table 4
Quality and spacing metrics.
Problem number
Quality metric (QM)
Spacing metric (SM)
PAES
NSGA-II
MOPSA
NSGA-II
PAES
MOPSA
1
0
0.08
0.92
0.558
0.879
0.444
2
0.09
0.18
0.73
0.521
0.937
0.315
3
0
0.20
0.80
0.409
0.810
0.385
4
0
0
1
0.395
0.762
0.373
5
0
0.10
0.90
0.608
0.726
0.359
6
0
0.00
1
0.359
0.932
0.282
7
0
0.15
0.85
0.608
0.761
0.274
8
0
0
1
0.595
0.816
0.335
9
0
0.25
0.75
0.870
0.866
0.303
10
0
0.10
0.90
0.672
0.768
0.243
11
0
0
1
0.574
0.893
0.586
12
0
0
1
0.862
0.845
0.303
13
0
0
1
0.674
0.745
0.394
14
0.08
0.17
0.75
0.624
0.935
0.476
15
0.06
0.2
0.74
0.870
0.730
0.314
16
0
0
1
0.517
0.788
0.521
17
0
0
1
0.657
0.771
0.442
18
0
0.12
0.88
0.601
0.859
0.276
19
0
0
1
0.552
0.727
0.486
20
0
0
1
0.562
0.895
0.610
21
0
0
1
0.772
0.819
0.636
22
0
0
1
0.703
0.812
0.481
23
0
0.04
0.96
0.641
0.835
0.423
24
0
0.10
0.90
0.774
0.686
0.359
25
0
0
1
0.531
0.808
0.224
26
0
0
1
0.724
0.767
0.512
27
0
0.05
0.95
0.721
0.866
0.335
28
0
0
1
0.669
0.794
0.572
29
0
0
1
0.555
0.527
0.303
30
0
0
1
0.621
0.862
0.394Table 5
Diversification and mean ideal distance metrics.
Problem number
Diversity metric (DM)
Mean ideal distance (MID)
PAES
NSGA-II
MOPSA
PAES
NSGA-II
MOPSA
1
0.929
0.872
1.379
0.862
0.569
0.255
2
0.606
1.065
1.143
0.698
0.673
0.258
3
0.617
0.772
1.300
0.795
0.743
0.302
4
0.667
0.805
0.971
0.893
0.621
0.480
5
0.860
0.758
1.111
0.622
0.797
0.317
6
0.893
0.754
1.358
0.776
0.527
0.282
7
0.859
1.048
1.296
0.724
0.596
0.245
8
0.780
0.932
1.380
0.793
0.653
0.519
9
0.819
0.920
1.228
0.679
0.518
0.312
10
0.718
0.758
0.917
0.828
0.718
0.239
11
0.898
1.042
1.325
0.899
0.667
0.430
12
0.675
0.949
1.367
0.656
0.659
0.227
13
0.875
0.840
1.240
0.835
0.749
0.385
14
0.673
0.905
1.279
0.658
0.758
0.203
15
0.771
1.033
1.212
0.898
0.737
0.372
16
0.793
1.022
1.240
0.673
0.606
0.447
17
0.648
0.724
1.098
0.895
0.655
0.282
18
0.836
0.860
1.084
0.638
0.636
0.245
19
0.690
0.911
1.394
0.669
0.481
0.319
20
0.754
0.867
0.918
0.607
0.468
0.312
21
0.833
0.963
1.343
0.782
0.496
0.239
22
0.700
0.951
1.357
0.633
0.649
0.430
23
0.716
0.817
1.298
0.722
0.647
0.227
24
0.847
0.873
0.949
0.866
0.571
0.385
25
0.706
0.706
1.031
0.764
0.572
0.203
26
0.930
1.094
1.068
0.711
0.486
0.372
27
0.994
0.767
1.240
0.662
0.610
0.437
28
0.892
0.742
0.968
0.732
0.639
0.270
29
0.737
0.849
1.261
0.887
0.532
0.334
30
0.834
0.779
0.953
0.637
0.534
0.471


## 7. Conclusion

This paper integrates production, distribution and logistics activities at the strategic decision making level, where the objective is to design a multi-echelon supply chain network considering agility as a key design criterion. Agility performance index has been conceptualized as a key design criterion which is decided by the manufacturers. A network with five echelons of supply chains including suppliers, plants, distribution centers, cross-docks, and customer zones is addressed in this paper. The problem has been mathematically formulated as a multiobjective optimization model that aims to minimize the cost (fixed and variable) and maximizes the plant flexibility and volume flexibility. Also, in order to cope with uncertain environment, some parameter of the model such as demand and cost are considered as fuzzy parameters. To validate our proposed MOPSA, various test problems were designed to evaluate its performance and reliability in comparison with two conventional multiobjective evolutionary algorithms, known as NSGA-II and PAES. Experimental results showed the higher performance of the proposed MOPSA in all of four comparison metrics. The proposed MOPSA has higher quality in the whole test problems which makes more near optimal solution for the decision makers.



---

*Source: 543604-2013-06-19.xml*