# Hybrid Electromagnetism-Like Algorithm for Dynamic Supply Chain Network Design under Traffic Congestion and Uncertainty

Authors: Javid Jouzdani; Mohammad Fathian
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

With the constantly increasing pressure of the competitive environment, supply chain (SC) decision makers are forced to consider several aspects of business climate. More specifically, they should take into account the endogenous features (e.g., available means of transportation, and the variety of products) and exogenous criteria (e.g., the environmental uncertainty, and transportation system conditions). In this paper, a mixed integer nonlinear programming (MINLP) model for dynamic design of a supply chain network is proposed. In this model, multiple products and multiple transportation modes, the time value of money, traffic congestion, and both supply-side and demand-side uncertainties are considered. Due to the complexity of such models, conventional solution methods are not applicable; therefore, two hybrid Electromagnetism-Like Algorithms (EMA) are designed and discussed for tackling the problem. The numerical results show the applicability of the proposed model and the capabilities of the solution approaches to the MINLP problem.


---

## Body

## 1. Introduction

Since its introduction by [1], Supply Chain Management (SCM) has been a point of interest to both researchers and practitioners. The reason might be the fact that a significant portion of budget and time in companies is spent on supply chain activities [2]. Obviously, the business units involved in these activities operate efficiently if only the whole supply chain network (SCN) is well structured [3]. Therefore, supply chain network design (SCND), as one of the most important strategic decisions, is known to be vital to survival of companies in today’s highly competitive business environment. Generally, SCND includes the determination of the location, capacity, number and the technology of the facilities, and the quantities of products transported from one element to another in order to fulfill the demand in each node of SCN [4]. However, several other factors may affect the decisions regarding the design of the SC. In this paper, the multiperiod design of a SC with multiple transportation modes and multiple products considering the time value of money, traffic congestion, and uncertainty in both demand and supply sides is investigated. In what follows the significance of these features in SCND problem is discussed.In today’s business environment, SCs are urged to be more responsible for their business activities regarding the environment and the society. The concern about these aspects of SC activities has led to the introduction of the “Sustainable Supply Chain Management (SSCM)” concepts. Despite the importance of sustainability in SCM, the relevant literature is still rare [5]. SSCM is defined as the integration of the SC’s objective considering economic, social, and environmental aspects of decisions in order to enable the long-term operation of the SC [6]. The research on sustainable supply chain network design (SSCND) under uncertainty is gaining ground in recent years. Considering traffic congestion, Jouzdani et al. [7] proposed a model for SCND under demand uncertainty and presented a dairy industry case study. Pishvaee et al. [5] proposed an accelerated Benders decomposition algorithm for sustainable supply chain network design under uncertainty and addressed a case of medical needle and syringe supply chain. They consider economic, environmental, and social objectives. Subulan et al. [8] presented a closed-loop supply chain network design for lead/acid batteries. They considered the risk of end-of-life product collection besides other aspects of SCND decision making. Sadjadi et al. [9] proposed a capacitated multiechelon, multiproduct reverse logistic network design with fuzzy demand for end-of-life products. Osmani and Zhang [10] presented a sustainable dual feedstock lignocellulosic-based bioethanol SCND model considering uncertainty in demand, supply, and prices of products. A robust closed-loop SCND (CLSCND) model under uncertainty for a case of medical device industry was investigated by [11]. Dubey and Gunasekaran [12] proposed a SSCND model and addressed a case of an Indian Company. In a recent research, Boukherroub et al. [13] proposed and integrated approach to sustainable supply chain planning. A goal programming approach was proposed by [14] for low carbon supply chain configuration for a new product. Environmental impact of supply chain and its total variable and fixed costs are considered as two objectives by Govindan et al. [15] for integrating sustainable order allocation and sustainable supply chain network strategic design with stochastic demand. They proposed a robust hybrid multiobjective metaheuristic comprised of Electromagnetism-Like Algorithm (EMA) and Variable Neighborhood Search (VNS) as a solution approach to their model. For a review of recent studies in the optimization approaches in SSCND, one may refer to a review by Eskandarpour et al. [16]. The subject of traffic congestion in SCND problem has already attracted researchers in the past few years [7, 17]. Specifically, since roads are of the most widely used transportation systems in many SCs especially for distribution of good in urban areas, the mutual effect between the transportation activities of a SC and the traffic congestion is a significant issue in SSCND. More specifically, in a SC, different vehicles have different effects on traffic congestion and, therefore, traffic congestion should also be considered in decision regarding the amounts of different products transported by different vehicles. As an example, heavy duty vehicles (e.g., trucks) are capable of transporting higher volumes of goods compared to light vehicles (e.g., pickups); however, the impact of heavy duty vehicles on traffic congestion is larger than that of light types. Therefore, considering different transportation modes, their capacities, and their impact on traffic congestion are significant issues in SCND.The effectiveness of the SCs is significantly affected by the uncertainty rooted in their complex and dynamic nature [18]. Specifically, the fluctuations and uncertainties in demand-side and supply-side in a SC are results of its inherent complex and variable characteristics. In order to capture the uncertainty in the SCs, some authors have proposed Stochastic Programming (SP) models [19–21]. However, it is not always possible to precisely determine the distribution parameters [22]. In addition, the unavailability of historical data is a major drawback of SP methods [5]. Therefore, many authors have utilized Possibilistic Programming (PP) [5, 7, 23] and Scenario-based/Robust Optimization (RO) [24–28] for modeling the epistemic uncertainty in SCND. The reader is referred to the reviews in SCND modeling uncertainty for further reading [18, 29].SCND has its roots in Facility Location Problem (FLP) and inherits its complexity [30]. Therefore, especially in large-scale instances of the problem, conventional solution methods fail to provide a solution in a reasonable amount of time. The problem becomes even more complex when several decision factors are involved. The complexity of the large-scale SCND problem has inspired researchers to propose heuristic and metaheuristic solution methods [22]. A hybrid of Particle Swarm Optimization (PSO) and Genetic Algorithm (GA) was proposed by Soleimani and Kannan [31] for tackling a CLSCND problem. Mousavi et al. [32] proposed a modified PSO algorithm for solving an integrated location and inventory control model of a two-echelon supply chain network. Roghanian and Pazhoheshfar [33] utilized a priority-based GA to solve a reverse logistics network design under uncertainty. A hybrid of Multiobjective Adaptive Memory Programming (MOAMP) and Tabu Search (TS) algorithms was proposed by Cardona-Valdés et al. [34] for a biobjective SCND stochastic problem. The authors [11] applied a hybrid of Memetic Algorithm (MA) and Variable Neighborhood Search (VNS) to solve their robust model for a closed-loop global SCND under uncertainty. A nondominated sorting genetic algorithm (NSGA-II) was utilized by Pasandideh et al. [35] to optimize a multiproduct multiperiod three-echelon supply chain problem under uncertainty. Govindan et al. [15] proposed a robust hybrid multiobjective metaheuristic comprised of Electromagnetism-Like Algorithm (EMA) and VNS as a solution approach to their SCND model.Among the metaheuristics, EMA is a relatively new population-based algorithm proposed by Birbil and Fang [36] that has been applied to many optimization problems (such as set covering [37], machine scheduling [38], flow-shop scheduling [39, 40], project scheduling [41], periodic job-shop scheduling [42] SCND [15], and travelling salesman problem [43, 44]) and has outperformed many other similar algorithms [37, 38, 41, 42, 45–47]. In addition, EMA can be easily augmented in order to further improve its performance [15, 40, 48].Simulated Annealing (SA) and Variable Neighborhood Search (VNS) are well-known, efficient, and widely used metaheuristic methods. SA was originally proposed by Kirkpatrick et al. [49] and was inspired by the annealing process of materials. SA has already been successfully applied to many problems; however, as a parameter-sensitive algorithm, it may fail to fully explore the search space; therefore, hybridization of SA is proposed to resolve this drawback [48, 50–55]. VNS, originally proposed by Mladenović and Hansen [56], utilizes more than a single neighborhood structure and switches between them in a local search process. VNS is a simple and effective algorithm; however, it may get trapped in local minima due to limited exploration; therefore, hybridization of VNS is proposed to overcome this weakness [15, 57–62].In this paper, a hybrid of EMA and SA and a combination of EMA and VNS are studied for solving the MINLP model of SCND problem under uncertainty and traffic congestion. In these hybrid algorithms the simple local search procedure in EMA is replaced by SA and VNS. The results show that the replacement leads to promising improvements. The main contributions of this research are as follows:(1)
concurrent consideration of multiple products, multiple transportation vehicles, multiple periods, time value of money, traffic congestion, and demand-side and supply-side uncertainties along with the time varying interest rate (time value of money) in SCND,(2)
proposing and studying two hybrid population-based metaheuristics (EMA-VNS and EMA-SA) for solving MINLP model of the SCND problem.In order to clarify the characteristics of our model, Table 1 provides a comparison between the present study and the recent related works. The abbreviations used in this table are defined below the table. From the literature review point of view, it can be concluded that the literature on supply-side uncertainty, traffic congestion, and multiple transportation modes in SCND is scares compared to other features. In addition, to the best of the authors knowledge, concurrent consideration of all the features mentioned in Table 1 has not been studied. Furthermore, in the majority of the previously published researches, exact algorithms (e.g., CPLEX or LINGO solvers) are utilized. These methods may fail to efficiently solve complex problems. Therefore, a lack of study on metaheuristic methods for SCND problem with the aforementioned feature may be concluded from Table 1. The literature on general SCND problem is immense; therefore, in this table, an effort is made to consider the most applicable features of SCND and most related works from the research background.Table 1
A comparison of the recent related literature and our research (chronologically sorted).
Author(s)
MP
Uncertainty
DY
CN
MM
Network flow
Solution approach
DU
SU
FF
RF
CL
Liu et al. [65]
2010
✓
✓
LR
Rentizelas and Tatsiopoulos [66]
2010
✓
✓
✓
GA + SQP
Salema et al. [67]
2010
✓
✓
✓
✓
CPLEX
Pishvaee and Torabi [23]
2010
✓
✓
✓
LINGO
Le and Lee [68]
2013
✓
✓
✓
✓
LINGO
Bai et al. [69]
2011
✓
✓
✓
LPR + LRB
Mirzapour Al-e-hashem et al. [70]
2011
✓
✓
✓
LINGO
Georgiadis et al. [19]
2011
✓
✓
✓
✓
CPLEX
Pishvaee et al. [71]
2011
✓
✓
CPLEX
Chen and Fan [72]
2012
✓
✓
✓
✓
DA
Wang et al. [73]
2011
✓
✓
CPLEX
Kostin et al. [74]
2012
✓
✓
✓
✓
SAA
Almansoori and Shah [75]
2012
✓
✓
✓
✓
CPLEX
Bashiri et al. [76]
2012
✓
✓
✓
CPLEX
Wang and Watada [77]
2012
✓
✓
PSO
Sadjady and Davoudpour [78]
2012
✓
✓
✓
LR
Vahdani et al. [79]
2013
✓
✓
GAMS
Baghalian et al. [80]
2013
✓
✓
✓
LINGO
Cardoso et al. [21]
2013
✓
✓
✓
✓
CPLEX
Ramezani et al. [25]
2013
✓
✓
✓
✓
CPLEX
Hatefi and Jolai [81]
2014
✓
✓
✓
CPLEX
Jindal and Sangwan [82]
2013
✓
✓
✓
✓
LINGO
Jouzdani et al. [7]
2013
✓
✓
✓
✓
LINGO
Ramezani et al. [25]
2013
✓
✓
✓
SRA
Soleimani et al. [83]
2013
✓
✓
✓
GA
Yang and Liu [22]
2015
✓
✓
✓
MA
Hasani et al. [11]
2014
✓
✓
✓
✓
✓
MA + VNS
Jabbarzadeh et al. [84]
2014
✓
✓
✓
✓
LINGO
Mousavi et al. [32]
2014
✓
✓
✓
Modified PSO
Osmani and Zhang [10]
2014
✓
✓
✓
✓
SAA
Pishvaee et al. [5]
2014
✓
✓
✓
Benders Dec.
Boukherroub et al. [13]
2015
✓
✓
✓
GA
Sadjadi et al. [9]
2014
✓
✓
✓
MA
Shabani et al. [85]
2014
✓
✓
✓
✓
CPLEX
Brandenburg [14]
2015
✓
✓
✓
GA + PSO
Dubey and Gunasekaran [12]
2014
✓
✓
CPLEX
Brandenburg [14]
2015
✓
✓
✓
✓
✓
CPLEX
Boukherroub et al. [13]
2015
✓
✓
✓
✓
CPLEX
Govindan et al. [15]
2015
✓
✓
EMA + VNS
Pasandideh et al. [35]
2015
✓
✓
✓
✓
✓
NSGA-II
Current Research
✓
✓
✓
✓
✓
✓
✓
Hybrid EMAs
MP: multiple products, DU: demand uncertainty, SU: supply uncertainly, DY: dynamic model (multiple periods), CN: congestion, MM: multiple transportation modes, FF: forward network flow, RF: reverse network flow, and CL: closed-loop network.The rest of this paper is organized as follows. The next section provides a brief discussion on the concept of superiority and inferiority for fuzzy triangular numbers. In Section3, the problem is presented and Section 4 discusses the solution methods. The results of numerical experiments are provided in Section 5 and, finally, Section 6 concludes the paper and presents guidelines for future research.


## 2. Inferiority and Superiority of Fuzzy Triangular Numbers

Fuzzy numbers may be utilized for modeling the uncertainty when the data is imprecise. The ease of application and the wide range of real problems in which they can be expressed through fuzzy triangular numbers are a good reason for their popularity among both researchers and practitioners.In 2007, van Hop proposed an approach for comparing the fuzzy triangular numbers [63]. According to this theorem, if P~≤Q~, then the superiority of a fuzzy triangular number P~=u,a,b over a fuzzy triangular number Q~=v,c,d is defined by the following equation:(1)SQ~,P~=v-u+d-b2.In the above equation, u is the center and a and b are the left and right spreads of P~=u,a,b, respectively, and v is the center and c and d are the left and right spreads of Q~=v,c,d, respectively. Similarly, the inferiority of P~ to Q~ is calculated by the following equation:(2)IP~,Q~=v-u+a-c2.The concepts of superiority and inferiority can be suitably utilized for modeling the constraints of an optimization problem where a fuzzy triangular number is involved [64].In this paper, the products demand and supply capacities are considered as triangular fuzzy numbers and the superiority and inferiority concepts are utilized for modeling “the inferiority of the demand for a product in a demand node in a time period to the amount of that product transported to that node” and “the superiority of the capacity for a product in a supply node in a time period over the amount of that product transported from that node.” These are discussed in more details in Sections3.3 and 3.4.


## 3. The Problem

In this section, a definition of the problem is presented, the assumptions based on which the mathematical model is built are provided, and the mathematical model is studied.


### 3.1. Problem Definition

In this paper, our focus is on facility location and transportation planning in SCND. Specifically, we assume a two-echelon supply chain which utilizes multiple vehicle types to transport multiple products to its customers in multiple planning periods. Here, for the sake of simplicity, two echelons are considered; however, the number of echelons may be extended to create a more realistic and complex model with some minor modifications. The problem is the location and relocation of facilities and determining the amounts of products transported from each facility to each customer in each planning period such that the total cost is minimized. The total cost consists of the facility location and relocation cost, the transportation cost, the traffic congestion cost, and the uncertainty cost. Here, the traffic congestion cost is calculated through the formula provided by the US Bureau of Public Roads (BPR) (formulated in (16)) and is expressed through the monetary value of time and the uncertainty cost is determined through the concepts of inferiority/superiority of the amounts of products transported from the facilities to the demand nodes to/over the uncertain demand and supply amounts as triangular fuzzy numbers. The formulations of these concepts are provided later in this section.


### 3.2. Assumptions

Assumptions determine the scope of the problem and define its capabilities and limitations. Current research is based on the following assumptions.(1)
The total number of nodes in the supply network, the total number of transportation vehicle types, the total number of products, the total number of periods in the planning horizon are known and fixed.(2)
The products are produced and shipped in continuous quantities.(3)
Each facility has a limited capacity for each product in each period. The capacity is subject to uncertainty and can be expressed through fuzzy triangular numbers.(4)
Each demand node has a limited demand for each product in each period. The demand is subject to uncertainty and can be expressed through fuzzy triangular numbers.(5)
All parameters except the capacity of vehicles and their traffic congestion coefficients are time varying.


### 3.3. Nomenclature

In order to facilitate the understanding of the mathematical model, the sets, parameters, and the decision variables are introduced in this section.


#### 3.3.1. Sets

N
: the set of network nodes (including candidate facility locations and demand points),
T: the set of time periods in the planning horizon,
P: the set of products,
M: the set of transportation modes.


#### 3.3.2. Subscripts

i: subscript for candidate facility location (i∈N),
j: subscript for demand point (j∈N),
t: subscript for time period (t∈T),
p: subscript for product (p∈P),
m: subscript for transportation vehicle (m∈M).


#### 3.3.3. Parameters

IR
t: interest rate in time period t∈T,
MVT
t: monetary value of time in time period t∈T,
α: parameter in the link performance function,
β: parameter in the link performance function,
CAP
m: capacity of vehicle m∈M,
TCC
m: traffic congestion coefficient of vehicle m∈M,
MC
i
t: the annual maintenance cost of a facility in candidate location i∈N in time period t∈T,
OC
i
t: the opening cost of a facility in candidate location i∈N in time period t∈T,
CC
i
t: the closing cost of a facility in candidate location i∈N in time period t∈T,
D
j
p
t
c: the center of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
ls: the left spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
rs: the right spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
S
i
p
t
c: the center of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
ls: the left spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
rs: the right spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
FFTT
i
j
t: free flow travel time from node i∈N to node j∈N in time period t∈T,
TCAP
i
j
t: traffic capacity of the link from node i∈N to node j∈N in time period t∈T,
BF
i
j
t: basic traffic flow of the link from node i∈N to node j∈N in time period t∈T,
C
i
j
m
t: the transportation cost for vehicle m∈M for link from node i∈N to node j∈N in time period t∈T,
C
I
D
,
x
j
p
t
in: the unit cost of inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
C
S
S
,
x
i
p
t
out: the unit cost of superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node.


#### 3.3.4. Decision Variables

y
i
t: a binary variable which is equal to 1 if a facility is operating in i∈N in time period t∈T and equals 0 otherwise,
x
i
j
m
p
t: the amount of product p∈P transported from node i∈N to node j∈N by means of transportation vehicle m∈M in time period t∈T.


#### 3.3.5. Auxiliary Variable

λ
I
D
,
x
j
p
t
in: the inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
λ
S
S
,
x
i
p
t
out: the superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node,
z
i
t
o: being 1 if a facility is opened in candidate supply location i∈N in time period t∈T,
z
i
t
c: being 1 if a facility is closed in candidate supply location i∈N in time period t∈T,
n
i
j
m
p
t: the number of vehicles of type m∈M transporting product p∈P from node i∈N to node j∈N in time period t∈T minus 1,
z
i
j
m
t: the total number of vehicles of type m∈M transporting goods from node i∈N to node j∈N in time period t∈T,
z
i
j
m
p
t
l: the auxiliary variable used form linearization, determining the amount of product p∈P transported by vehicle of type m∈M from node i∈N to node j∈N in time period t∈T if a facility is operating in node i∈N,
FLW
i
j
t: total traffic flow from node i∈N to node j∈N,
x
j
p
t
in: the total amount of product p∈P transported to a node j∈N in time period t∈T,
x
i
p
t
out: the total amount of product p∈P transported from a node i∈N in time period t∈T.


### 3.4. Mathematical Model

In what follows, the dynamic multiproduct multimode SCND problem under transportation cost uncertainty is formulated according to the notations introduced above:(3a)
m
i
n
∑
i
∈
N
∑
t
∈
T
y
i
t
MC
i
t
∏
n
=
1
t
1
+
IR
n
+
∑
i
∈
N
∑
t
∈
T
z
i
t
o
OC
i
t
∏
n
=
1
t
1
+
IR
n
+
∑
i
∈
N
∑
t
∈
T
z
i
t
c
CC
i
t
∏
n
=
1
t
1
+
IR
n
(3b)
+
∑
i
∈
N
∑
j
∈
N
∑
m
∈
M
∑
t
∈
T
C
i
j
m
t
z
i
j
m
t
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(3c)
+
∑
i
∈
N
∑
j
∈
N
∑
t
∈
T
MVT
t
×
FLW
i
j
t
×
t
FLW
i
j
t
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
(3d)
+
∑
j
∈
N
∑
p
∈
P
∑
t
∈
T
C
I
D
,
x
j
p
t
in
×
λ
I
D
,
x
j
p
t
in
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(3e)
+
∑
i
∈
N
∑
p
∈
P
∑
t
∈
T
C
S
S
,
x
i
p
t
out
×
λ
S
S
,
x
i
p
t
out
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
(4)
subject to
n
i
j
m
p
t
=
x
i
j
m
p
t
CAP
m
+
1
∀
i
,
j
∈
N
,
∀
p
∈
P
,
∀
m
∈
M
,
∀
t
∈
T
f
f
f
f
f
f
f
f
(5)
z
i
j
m
t
=
∑
p
∈
P
n
i
j
m
p
t
∀
i
,
j
∈
N
,
∀
m
∈
M
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
(6)
FLW
i
j
t
=
BF
i
j
t
+
∑
m
∈
M
TCC
m
×
z
i
j
m
t
∀
i
,
j
∈
N
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
i
i
f
f
f
(7)
∑
j
∈
I
∑
m
∈
M
x
i
j
m
p
t
-
∑
i
∈
I
∑
m
∈
M
x
j
i
m
p
t
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
=
y
i
t
×
x
i
p
t
out
-
x
i
p
t
in
∀
i
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(8)
λ
I
D
,
x
j
p
t
in
=
x
j
p
t
in
-
D
j
p
t
c
+
D
j
p
t
ls
2
∀
j
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
i
f
f
f
f
f
f
f
f
f
f
f
(9)
λ
S
S
,
x
i
p
t
out
=
-
x
i
p
t
out
+
S
i
p
t
c
+
S
i
p
t
rs
2
∀
i
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
(10)
z
i
t
o
=
y
i
t
×
1
-
y
i
t
-
1
∀
i
∈
N
,
∀
t
∈
T
,
t
≥
2
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
f
f
f
f
f
f
(11)
z
i
t
o
=
y
i
t
∀
i
∈
N
,
t
=
1
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(12)
z
i
t
+
1
c
=
1
-
y
i
t
+
1
×
y
i
t
∀
i
∈
N
,
∀
t
∈
T
,
t
≤
T
-
1
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(13)
z
i
t
c
=
y
i
t
∀
i
∈
N
,
t
=
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(14)
y
i
t
,
z
i
t
c
,
z
i
t
o
∈
0,1
∀
i
∈
N
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(15)
x
i
j
m
p
t
≥
0
∀
i
,
j
∈
N
,
∀
p
∈
P
,
∀
m
∈
M
,
∀
t
∈
T
.
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
fIn the above model, (3a) calculates the total facility cost which consists of annual maintenance cost, opening cost, and closing cost of facilities in the SC. Equation (3b) presents the total transportation cost. Total traffic congestion cost is calculated in (3c) in which the following formula is provided by the US Bureau of Public Roads (BPR) as(16)tFLWijt=FFTTijt1+αFLWijtTCAPijtβ,∀i,j∈N,∀t∈T.Equation (3d) calculates the cost of inferiority of the demand in demand nodes to the amounts of products transported to those nodes. Similarly, (3e) is for obtaining the cost of superiority of the capacity in supply nodes over the amounts of products transported from those nodes.Constraint (4) and Constraint (5) determine the total number of each vehicle type used for transporting goods from each node to another in each time period. Constraint (6) is for obtaining the flow of each link on each time period based on its basic flow and the congestion caused by the fleet of the SC. Constraint (7) presents the flow balance in each node. Constraint (8) calculates the inferiority of the demand in each demand node for each product in each time period to the amount of that product transported to that node. Similarly, Constraint (9) gives the superiority of the capacity in each supply node for each product in each time period over the amount of that product transported from that node. In Constraint (10), binary variables are defined to determine if a facility is opened in each node in each time period and Constraint (11) initializes the variables for the first time period in the planning horizon. Similarly, Constraint (12) calculates binary variables that determine if a facility is closed in each node in each time period and Constraint (13) determines the values for the final time period in the planning horizon. Finally, Constraint (14) and Constraint (15) determine the types of variables.


## 3.1. Problem Definition

In this paper, our focus is on facility location and transportation planning in SCND. Specifically, we assume a two-echelon supply chain which utilizes multiple vehicle types to transport multiple products to its customers in multiple planning periods. Here, for the sake of simplicity, two echelons are considered; however, the number of echelons may be extended to create a more realistic and complex model with some minor modifications. The problem is the location and relocation of facilities and determining the amounts of products transported from each facility to each customer in each planning period such that the total cost is minimized. The total cost consists of the facility location and relocation cost, the transportation cost, the traffic congestion cost, and the uncertainty cost. Here, the traffic congestion cost is calculated through the formula provided by the US Bureau of Public Roads (BPR) (formulated in (16)) and is expressed through the monetary value of time and the uncertainty cost is determined through the concepts of inferiority/superiority of the amounts of products transported from the facilities to the demand nodes to/over the uncertain demand and supply amounts as triangular fuzzy numbers. The formulations of these concepts are provided later in this section.


## 3.2. Assumptions

Assumptions determine the scope of the problem and define its capabilities and limitations. Current research is based on the following assumptions.(1)
The total number of nodes in the supply network, the total number of transportation vehicle types, the total number of products, the total number of periods in the planning horizon are known and fixed.(2)
The products are produced and shipped in continuous quantities.(3)
Each facility has a limited capacity for each product in each period. The capacity is subject to uncertainty and can be expressed through fuzzy triangular numbers.(4)
Each demand node has a limited demand for each product in each period. The demand is subject to uncertainty and can be expressed through fuzzy triangular numbers.(5)
All parameters except the capacity of vehicles and their traffic congestion coefficients are time varying.


## 3.3. Nomenclature

In order to facilitate the understanding of the mathematical model, the sets, parameters, and the decision variables are introduced in this section.


### 3.3.1. Sets

N
: the set of network nodes (including candidate facility locations and demand points),
T: the set of time periods in the planning horizon,
P: the set of products,
M: the set of transportation modes.


### 3.3.2. Subscripts

i: subscript for candidate facility location (i∈N),
j: subscript for demand point (j∈N),
t: subscript for time period (t∈T),
p: subscript for product (p∈P),
m: subscript for transportation vehicle (m∈M).


### 3.3.3. Parameters

IR
t: interest rate in time period t∈T,
MVT
t: monetary value of time in time period t∈T,
α: parameter in the link performance function,
β: parameter in the link performance function,
CAP
m: capacity of vehicle m∈M,
TCC
m: traffic congestion coefficient of vehicle m∈M,
MC
i
t: the annual maintenance cost of a facility in candidate location i∈N in time period t∈T,
OC
i
t: the opening cost of a facility in candidate location i∈N in time period t∈T,
CC
i
t: the closing cost of a facility in candidate location i∈N in time period t∈T,
D
j
p
t
c: the center of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
ls: the left spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
rs: the right spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
S
i
p
t
c: the center of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
ls: the left spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
rs: the right spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
FFTT
i
j
t: free flow travel time from node i∈N to node j∈N in time period t∈T,
TCAP
i
j
t: traffic capacity of the link from node i∈N to node j∈N in time period t∈T,
BF
i
j
t: basic traffic flow of the link from node i∈N to node j∈N in time period t∈T,
C
i
j
m
t: the transportation cost for vehicle m∈M for link from node i∈N to node j∈N in time period t∈T,
C
I
D
,
x
j
p
t
in: the unit cost of inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
C
S
S
,
x
i
p
t
out: the unit cost of superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node.


### 3.3.4. Decision Variables

y
i
t: a binary variable which is equal to 1 if a facility is operating in i∈N in time period t∈T and equals 0 otherwise,
x
i
j
m
p
t: the amount of product p∈P transported from node i∈N to node j∈N by means of transportation vehicle m∈M in time period t∈T.


### 3.3.5. Auxiliary Variable

λ
I
D
,
x
j
p
t
in: the inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
λ
S
S
,
x
i
p
t
out: the superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node,
z
i
t
o: being 1 if a facility is opened in candidate supply location i∈N in time period t∈T,
z
i
t
c: being 1 if a facility is closed in candidate supply location i∈N in time period t∈T,
n
i
j
m
p
t: the number of vehicles of type m∈M transporting product p∈P from node i∈N to node j∈N in time period t∈T minus 1,
z
i
j
m
t: the total number of vehicles of type m∈M transporting goods from node i∈N to node j∈N in time period t∈T,
z
i
j
m
p
t
l: the auxiliary variable used form linearization, determining the amount of product p∈P transported by vehicle of type m∈M from node i∈N to node j∈N in time period t∈T if a facility is operating in node i∈N,
FLW
i
j
t: total traffic flow from node i∈N to node j∈N,
x
j
p
t
in: the total amount of product p∈P transported to a node j∈N in time period t∈T,
x
i
p
t
out: the total amount of product p∈P transported from a node i∈N in time period t∈T.


## 3.3.1. Sets

N
: the set of network nodes (including candidate facility locations and demand points),
T: the set of time periods in the planning horizon,
P: the set of products,
M: the set of transportation modes.


## 3.3.2. Subscripts

i: subscript for candidate facility location (i∈N),
j: subscript for demand point (j∈N),
t: subscript for time period (t∈T),
p: subscript for product (p∈P),
m: subscript for transportation vehicle (m∈M).


## 3.3.3. Parameters

IR
t: interest rate in time period t∈T,
MVT
t: monetary value of time in time period t∈T,
α: parameter in the link performance function,
β: parameter in the link performance function,
CAP
m: capacity of vehicle m∈M,
TCC
m: traffic congestion coefficient of vehicle m∈M,
MC
i
t: the annual maintenance cost of a facility in candidate location i∈N in time period t∈T,
OC
i
t: the opening cost of a facility in candidate location i∈N in time period t∈T,
CC
i
t: the closing cost of a facility in candidate location i∈N in time period t∈T,
D
j
p
t
c: the center of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
ls: the left spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
D
j
p
t
rs: the right spread of the triangular fuzzy demand for product p∈P, in demand node j∈N, in time period t∈T,
S
i
p
t
c: the center of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
ls: the left spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
S
i
p
t
rs: the right spread of the triangular fuzzy supply for product p∈P, in facility node i∈N, in time period t∈T,
FFTT
i
j
t: free flow travel time from node i∈N to node j∈N in time period t∈T,
TCAP
i
j
t: traffic capacity of the link from node i∈N to node j∈N in time period t∈T,
BF
i
j
t: basic traffic flow of the link from node i∈N to node j∈N in time period t∈T,
C
i
j
m
t: the transportation cost for vehicle m∈M for link from node i∈N to node j∈N in time period t∈T,
C
I
D
,
x
j
p
t
in: the unit cost of inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
C
S
S
,
x
i
p
t
out: the unit cost of superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node.


## 3.3.4. Decision Variables

y
i
t: a binary variable which is equal to 1 if a facility is operating in i∈N in time period t∈T and equals 0 otherwise,
x
i
j
m
p
t: the amount of product p∈P transported from node i∈N to node j∈N by means of transportation vehicle m∈M in time period t∈T.


## 3.3.5. Auxiliary Variable

λ
I
D
,
x
j
p
t
in: the inferiority of the demand in a demand node j∈N for product p∈P in time period t∈T to the amount of that product transported to that node,
λ
S
S
,
x
i
p
t
out: the superiority of the capacity in a supply node i∈N for product p∈P in time period t∈T over the amount of that product transported from that node,
z
i
t
o: being 1 if a facility is opened in candidate supply location i∈N in time period t∈T,
z
i
t
c: being 1 if a facility is closed in candidate supply location i∈N in time period t∈T,
n
i
j
m
p
t: the number of vehicles of type m∈M transporting product p∈P from node i∈N to node j∈N in time period t∈T minus 1,
z
i
j
m
t: the total number of vehicles of type m∈M transporting goods from node i∈N to node j∈N in time period t∈T,
z
i
j
m
p
t
l: the auxiliary variable used form linearization, determining the amount of product p∈P transported by vehicle of type m∈M from node i∈N to node j∈N in time period t∈T if a facility is operating in node i∈N,
FLW
i
j
t: total traffic flow from node i∈N to node j∈N,
x
j
p
t
in: the total amount of product p∈P transported to a node j∈N in time period t∈T,
x
i
p
t
out: the total amount of product p∈P transported from a node i∈N in time period t∈T.


## 3.4. Mathematical Model

In what follows, the dynamic multiproduct multimode SCND problem under transportation cost uncertainty is formulated according to the notations introduced above:(3a)
m
i
n
∑
i
∈
N
∑
t
∈
T
y
i
t
MC
i
t
∏
n
=
1
t
1
+
IR
n
+
∑
i
∈
N
∑
t
∈
T
z
i
t
o
OC
i
t
∏
n
=
1
t
1
+
IR
n
+
∑
i
∈
N
∑
t
∈
T
z
i
t
c
CC
i
t
∏
n
=
1
t
1
+
IR
n
(3b)
+
∑
i
∈
N
∑
j
∈
N
∑
m
∈
M
∑
t
∈
T
C
i
j
m
t
z
i
j
m
t
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(3c)
+
∑
i
∈
N
∑
j
∈
N
∑
t
∈
T
MVT
t
×
FLW
i
j
t
×
t
FLW
i
j
t
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
(3d)
+
∑
j
∈
N
∑
p
∈
P
∑
t
∈
T
C
I
D
,
x
j
p
t
in
×
λ
I
D
,
x
j
p
t
in
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(3e)
+
∑
i
∈
N
∑
p
∈
P
∑
t
∈
T
C
S
S
,
x
i
p
t
out
×
λ
S
S
,
x
i
p
t
out
∏
n
=
1
t
1
+
IR
n
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
f
f
f
f
f
f
f
(4)
subject to
n
i
j
m
p
t
=
x
i
j
m
p
t
CAP
m
+
1
∀
i
,
j
∈
N
,
∀
p
∈
P
,
∀
m
∈
M
,
∀
t
∈
T
f
f
f
f
f
f
f
f
(5)
z
i
j
m
t
=
∑
p
∈
P
n
i
j
m
p
t
∀
i
,
j
∈
N
,
∀
m
∈
M
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
(6)
FLW
i
j
t
=
BF
i
j
t
+
∑
m
∈
M
TCC
m
×
z
i
j
m
t
∀
i
,
j
∈
N
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
i
i
f
f
f
(7)
∑
j
∈
I
∑
m
∈
M
x
i
j
m
p
t
-
∑
i
∈
I
∑
m
∈
M
x
j
i
m
p
t
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
=
y
i
t
×
x
i
p
t
out
-
x
i
p
t
in
∀
i
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(8)
λ
I
D
,
x
j
p
t
in
=
x
j
p
t
in
-
D
j
p
t
c
+
D
j
p
t
ls
2
∀
j
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
i
f
f
f
f
f
f
f
f
f
f
f
(9)
λ
S
S
,
x
i
p
t
out
=
-
x
i
p
t
out
+
S
i
p
t
c
+
S
i
p
t
rs
2
∀
i
∈
N
,
∀
p
∈
P
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
(10)
z
i
t
o
=
y
i
t
×
1
-
y
i
t
-
1
∀
i
∈
N
,
∀
t
∈
T
,
t
≥
2
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
f
f
f
f
f
f
(11)
z
i
t
o
=
y
i
t
∀
i
∈
N
,
t
=
1
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(12)
z
i
t
+
1
c
=
1
-
y
i
t
+
1
×
y
i
t
∀
i
∈
N
,
∀
t
∈
T
,
t
≤
T
-
1
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(13)
z
i
t
c
=
y
i
t
∀
i
∈
N
,
t
=
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(14)
y
i
t
,
z
i
t
c
,
z
i
t
o
∈
0,1
∀
i
∈
N
,
∀
t
∈
T
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
(15)
x
i
j
m
p
t
≥
0
∀
i
,
j
∈
N
,
∀
p
∈
P
,
∀
m
∈
M
,
∀
t
∈
T
.
f
f
f
f
f
f
f
f
f
f
f
f
f
i
f
f
f
f
f
fIn the above model, (3a) calculates the total facility cost which consists of annual maintenance cost, opening cost, and closing cost of facilities in the SC. Equation (3b) presents the total transportation cost. Total traffic congestion cost is calculated in (3c) in which the following formula is provided by the US Bureau of Public Roads (BPR) as(16)tFLWijt=FFTTijt1+αFLWijtTCAPijtβ,∀i,j∈N,∀t∈T.Equation (3d) calculates the cost of inferiority of the demand in demand nodes to the amounts of products transported to those nodes. Similarly, (3e) is for obtaining the cost of superiority of the capacity in supply nodes over the amounts of products transported from those nodes.Constraint (4) and Constraint (5) determine the total number of each vehicle type used for transporting goods from each node to another in each time period. Constraint (6) is for obtaining the flow of each link on each time period based on its basic flow and the congestion caused by the fleet of the SC. Constraint (7) presents the flow balance in each node. Constraint (8) calculates the inferiority of the demand in each demand node for each product in each time period to the amount of that product transported to that node. Similarly, Constraint (9) gives the superiority of the capacity in each supply node for each product in each time period over the amount of that product transported from that node. In Constraint (10), binary variables are defined to determine if a facility is opened in each node in each time period and Constraint (11) initializes the variables for the first time period in the planning horizon. Similarly, Constraint (12) calculates binary variables that determine if a facility is closed in each node in each time period and Constraint (13) determines the values for the final time period in the planning horizon. Finally, Constraint (14) and Constraint (15) determine the types of variables.


## 4. The Solution Algorithms

The proposed algorithms are hybridization of the Electromagnetism-Like Algorithm (EMA) with Simulated Annealing (SA) algorithm and the Variable Neighborhood Search (VNS), respectively. The hybrid algorithms and their aforementioned elements (i.e., EMA, SA, and VNS) are discussed in what follows.


### 4.1. The Electromagnetism-Like Algorithm

The EMA was originally proposed by Birbil and Fang [36]. The main idea in EMA was originated from the behavior of electromagnetic particles and the forces they exert on each other based on their distance and charge. The analogy here is that the solutions in the search space are considered as particles and their fitness values are regarded as their charges.In Algorithm1  popSize is the number of solution points (particles), maxIter is the maximum number of the algorithm iterations, LSIter is the maximum number of local search steps, and i and j are local counter variables.Algorithm 1:General procedure of the original EMA.
(1)Initialize the particles.
(2)i←1
(3)WHILE  i<maxIter  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)        Perform a local search.
(7)j←j+1
(8)END WHILE
(9)Calculate forces.
(10)     Move particles.
(11)i←i+1
(12)END WHILEThe original scheme of EMA is presented in Algorithm1 in which the forces are calculated according to Coulomb law [86] as(17)Fij=qiqjrij2,∀i,j∈1,2,…,m,where Fij is the force particle i and particle j exert on each other, rij is the distance between those particles, m is the number of particles, and qi and qj are the electrical charges of particle i and particle j, respectively. In EMA, the charges of the particles are calculated as(18)qi=e-nfxi-fxbest/∑k=1mfxk-fxbest,∀i∈1,2,…,m,where fx is the objective function value for solution x and n is the dimension of the solution space. The force exerted on each particle i is then obtained according to the following equation:(19)Fi=∑j=1j≠imxj-xiFij,ifqi<qj∑j=1j≠imxi-xjFij,ifqi≥qj,∀i∈1,2,…,m.The forces are used to move the particles in the solution space such that the updated particle is in the feasible solution space. The particles are moved according to the total force exerted on them through the following equation in which λ randomizes the particle movement and Fi is the norm of the force vector:(20)xi⟵xi+λFiFi,∀i∈1,2,…,m.Originally, EMA is designed for solving continuous optimization problems with bounded variables; however, many discrete problems such as machine scheduling [38, 87], cell formation [45], project scheduling [41], flow-shop scheduling [88], job-shop scheduling [42], assembly sequences planning [28], set covering [37], and travelling salesman problem [43, 44] are solved using the EMA.


### 4.2. The Simulated Annealing Algorithm

Simulated Annealing was introduced by Kirkpatrick et al. [49]. The method is inspired by metal annealing process in which the material is heated in order to prepare the material for atomic restructure and then gradually cooled down to reach the desired atomic configuration. Here, the analogy is between the temperature in the real annealing process and the temperature in SA. The algorithm is capable of escaping the local optima by allowing the acceptance of worse solutions in each iteration based on its temperature. More specifically, the higher the temperature, the more the chance of the algorithm moving to a worse solution in hope of finding the optimum point in the search space. SA is capable of efficiently and effectively solving combinatorial optimization problem and is easy to implement [89]. However, as a parameter-sensitive trajectory method, SA may fail to fully explore the search space. Therefore, it is combined with other metaheuristics to compensate its lack of exploration [48, 53–55].


### 4.3. The Variable Neighborhood Search Algorithm

Trajectory methods deal with a single solution in each step; however, they have shown their capability in exploiting the promising subsets of search space in order to reach high quality solutions. VNS is a simple and yet effective trajectory search algorithm originally proposed by Mladenović and Hansen [56]. VNS utilizes more than a single neighborhood structure and switches between the structures in a local search process. Contrary to many other metaheuristics, VNS needs few or sometimes no parameters. Despite its capabilities, VNS may be trapped in local minima due to limited exploration. Hence, its hybridization with other metaheuristics is proposed in the literature and is recently applied to solve problems [15, 57–61].


### 4.4. Solution Representation

One of the major issues in any search algorithm is the solution representation. In order to encode solutions into a particle, a 5-dimensional matrix,X, of size N×N×M×P×T, where A represents the number of elements in the set A, is considered. In X, each element Xi,j,m,p,t has the same interpretation as xijmpt. It should be noted that X also contains the information about the decision variable yit. For example, consider the binary decision variable, y23, which is equal to 1 if a facility is operating in node 2∈N in time period 3∈T and equals 0 otherwise. If y23=0, then no facility is opened in node 2∈N in time period 3∈T. Therefore, we have x2jmp3=0, ∀j∈N, ∀m∈M, ∀p∈P, and apparently according to the definition of X, we set X2,j,m,p,3=0, ∀j∈N, ∀m∈M, ∀p∈P.


### 4.5. The EMA-SA Hybrid

EMA results in satisfying solutions when combined with SA and has been successfully applied to discrete problems [40] and continuous search spaces [48]. In this section, an algorithm which benefits from both the EMA and SA is discussed. More specifically, the local search procedure in the original EMA is replaced by the SA. In addition, in order to provide the algorithm with a more guided search, the concept of temperature is incorporated into the movement of particles such that, in the early iterations of the algorithm, the particles move faster in order to increase the exploration of the search space, while in the final iterations, the temperature is dropped and the particles move slower in order to provide the intensification of the search procedure.A neighborhood to a solution,S, in SA is a solution S′=N1S that is generated by adding a matrix D, with randomly generated element, of the same size as S to S. By applying this type of neighborhood structure, some of the generated solutions may be infeasible. In such cases, the solutions are simply discarded. The general scheme of the EMA-SA hybrid for a minimization problem can be expressed as shown in Algorithm 2.Algorithm 2:General procedure of the EMA-SA hybrid.
(1)Initialize the parameters T0, Tf, maxRep, LSIter and the initial population
(2)T←T0
(3)WHILE  T≥Tf  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)i←1
(1)WHILE  i<maxRep  DO
(2)FOReach particle S in the population DO
(3)Generate a solution S′=N1S in the neighborhood S
(4)Δ←fS′-fS
(5)IF  Δ≤0  THEN
(6)S′←S
(7)ELSE
(8)Generate a random number r∈U0,1.
(9)IF  r≤e-Δ/T  THEN
(10)S′←S
(11)END IF
(12)END IF
(13)END FOR
(14)i←i+1
(15)END WHILE
(16)j←j+1
(17)END WHILE
(18)     Calculate forces.
(19)     Move particles.
(20)   Reduce T.
(21)END WHILEIn Algorithm2, the force particle i and particle j exert on each other, Fij, is calculated as(21)Fij=kqiqjrij2,∀i,j∈1,2,…,m,where k is obtained through the following equation and other elements are the same as defined in (17):(22)k=1+T-TfT0-Tf.It should be noted that, in EMA, we have k=1. The reason for using the coefficient in (22) is to incorporate the concept of temperature into the movement of particles. More specifically, as the temperature drops the force that the particles exert on each other decreases resulting in a more intense search in comparison with the early iterations of the algorithm. In addition, in a minimization problem, the charge for a particle i is calculated as(23)qi=nfxworst-fxifxworst-fxbest,∀i∈1,2,…,m,where xworst is the particle with the worst (largest) objective function and other elements are the same as defined in (18). The above equation not only considers the closeness to the best solution but also takes into account the distance from the worst particle. The calculation of Fi and the movement of each particle are the same as in EMA.


### 4.6. The EMA-VNS Hybrid

EMA and VNS have already been successfully combined for tackling problems [15]. Similar to what was mentioned in the previous section, we combine the EMA and VNS by replacing the simple local search in EMA by the VNS.The neighborhood structure in VNS is simple: a neighbor to a solutionS is S′=N2S, which is generated by adding a matrix D of the same size as S to S. The elements of D are normally distributed random variables with a mean of their corresponding elements in S and a variance of σ which is an increasing function of, i, the VNS iteration number, denoted by vi. Similar to what was mentioned above, the infeasible generated neighbors are discarded. Therefore, the general scheme of the EMA-VNS hybrid for a minimization problem can be expressed in Algorithm 3 in which other calculations are as in what was mentioned regarding the original EMA.Algorithm 3:General procedure of the EMA-VNS hybrid.
(1)Initialize the parameters maxIter, maxRep, LSIter and the initial population
(2)t←1
(3)WHILE  t<maxIter  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)i←1
(7)WHILE  i<maxRep  DO
(8)FOReach particle S in the population DO
(9)Generate a solution S′=N2S in the neighborhood S
(10)Δ←fS′-fS
(11)IF  Δ≤0  THEN
(12)S′←S
(13)ELSE
(14)σ=vi
(15)END IF
(16)END FOR
(17)i←i+1
(18)END WHILE
(19)j←j+1
(20)END WHILE
(21)      Calculate forces.
(22)     Move particles.
(23)t←t+1
(24)END WHILE


## 4.1. The Electromagnetism-Like Algorithm

The EMA was originally proposed by Birbil and Fang [36]. The main idea in EMA was originated from the behavior of electromagnetic particles and the forces they exert on each other based on their distance and charge. The analogy here is that the solutions in the search space are considered as particles and their fitness values are regarded as their charges.In Algorithm1  popSize is the number of solution points (particles), maxIter is the maximum number of the algorithm iterations, LSIter is the maximum number of local search steps, and i and j are local counter variables.Algorithm 1:General procedure of the original EMA.
(1)Initialize the particles.
(2)i←1
(3)WHILE  i<maxIter  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)        Perform a local search.
(7)j←j+1
(8)END WHILE
(9)Calculate forces.
(10)     Move particles.
(11)i←i+1
(12)END WHILEThe original scheme of EMA is presented in Algorithm1 in which the forces are calculated according to Coulomb law [86] as(17)Fij=qiqjrij2,∀i,j∈1,2,…,m,where Fij is the force particle i and particle j exert on each other, rij is the distance between those particles, m is the number of particles, and qi and qj are the electrical charges of particle i and particle j, respectively. In EMA, the charges of the particles are calculated as(18)qi=e-nfxi-fxbest/∑k=1mfxk-fxbest,∀i∈1,2,…,m,where fx is the objective function value for solution x and n is the dimension of the solution space. The force exerted on each particle i is then obtained according to the following equation:(19)Fi=∑j=1j≠imxj-xiFij,ifqi<qj∑j=1j≠imxi-xjFij,ifqi≥qj,∀i∈1,2,…,m.The forces are used to move the particles in the solution space such that the updated particle is in the feasible solution space. The particles are moved according to the total force exerted on them through the following equation in which λ randomizes the particle movement and Fi is the norm of the force vector:(20)xi⟵xi+λFiFi,∀i∈1,2,…,m.Originally, EMA is designed for solving continuous optimization problems with bounded variables; however, many discrete problems such as machine scheduling [38, 87], cell formation [45], project scheduling [41], flow-shop scheduling [88], job-shop scheduling [42], assembly sequences planning [28], set covering [37], and travelling salesman problem [43, 44] are solved using the EMA.


## 4.2. The Simulated Annealing Algorithm

Simulated Annealing was introduced by Kirkpatrick et al. [49]. The method is inspired by metal annealing process in which the material is heated in order to prepare the material for atomic restructure and then gradually cooled down to reach the desired atomic configuration. Here, the analogy is between the temperature in the real annealing process and the temperature in SA. The algorithm is capable of escaping the local optima by allowing the acceptance of worse solutions in each iteration based on its temperature. More specifically, the higher the temperature, the more the chance of the algorithm moving to a worse solution in hope of finding the optimum point in the search space. SA is capable of efficiently and effectively solving combinatorial optimization problem and is easy to implement [89]. However, as a parameter-sensitive trajectory method, SA may fail to fully explore the search space. Therefore, it is combined with other metaheuristics to compensate its lack of exploration [48, 53–55].


## 4.3. The Variable Neighborhood Search Algorithm

Trajectory methods deal with a single solution in each step; however, they have shown their capability in exploiting the promising subsets of search space in order to reach high quality solutions. VNS is a simple and yet effective trajectory search algorithm originally proposed by Mladenović and Hansen [56]. VNS utilizes more than a single neighborhood structure and switches between the structures in a local search process. Contrary to many other metaheuristics, VNS needs few or sometimes no parameters. Despite its capabilities, VNS may be trapped in local minima due to limited exploration. Hence, its hybridization with other metaheuristics is proposed in the literature and is recently applied to solve problems [15, 57–61].


## 4.4. Solution Representation

One of the major issues in any search algorithm is the solution representation. In order to encode solutions into a particle, a 5-dimensional matrix,X, of size N×N×M×P×T, where A represents the number of elements in the set A, is considered. In X, each element Xi,j,m,p,t has the same interpretation as xijmpt. It should be noted that X also contains the information about the decision variable yit. For example, consider the binary decision variable, y23, which is equal to 1 if a facility is operating in node 2∈N in time period 3∈T and equals 0 otherwise. If y23=0, then no facility is opened in node 2∈N in time period 3∈T. Therefore, we have x2jmp3=0, ∀j∈N, ∀m∈M, ∀p∈P, and apparently according to the definition of X, we set X2,j,m,p,3=0, ∀j∈N, ∀m∈M, ∀p∈P.


## 4.5. The EMA-SA Hybrid

EMA results in satisfying solutions when combined with SA and has been successfully applied to discrete problems [40] and continuous search spaces [48]. In this section, an algorithm which benefits from both the EMA and SA is discussed. More specifically, the local search procedure in the original EMA is replaced by the SA. In addition, in order to provide the algorithm with a more guided search, the concept of temperature is incorporated into the movement of particles such that, in the early iterations of the algorithm, the particles move faster in order to increase the exploration of the search space, while in the final iterations, the temperature is dropped and the particles move slower in order to provide the intensification of the search procedure.A neighborhood to a solution,S, in SA is a solution S′=N1S that is generated by adding a matrix D, with randomly generated element, of the same size as S to S. By applying this type of neighborhood structure, some of the generated solutions may be infeasible. In such cases, the solutions are simply discarded. The general scheme of the EMA-SA hybrid for a minimization problem can be expressed as shown in Algorithm 2.Algorithm 2:General procedure of the EMA-SA hybrid.
(1)Initialize the parameters T0, Tf, maxRep, LSIter and the initial population
(2)T←T0
(3)WHILE  T≥Tf  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)i←1
(1)WHILE  i<maxRep  DO
(2)FOReach particle S in the population DO
(3)Generate a solution S′=N1S in the neighborhood S
(4)Δ←fS′-fS
(5)IF  Δ≤0  THEN
(6)S′←S
(7)ELSE
(8)Generate a random number r∈U0,1.
(9)IF  r≤e-Δ/T  THEN
(10)S′←S
(11)END IF
(12)END IF
(13)END FOR
(14)i←i+1
(15)END WHILE
(16)j←j+1
(17)END WHILE
(18)     Calculate forces.
(19)     Move particles.
(20)   Reduce T.
(21)END WHILEIn Algorithm2, the force particle i and particle j exert on each other, Fij, is calculated as(21)Fij=kqiqjrij2,∀i,j∈1,2,…,m,where k is obtained through the following equation and other elements are the same as defined in (17):(22)k=1+T-TfT0-Tf.It should be noted that, in EMA, we have k=1. The reason for using the coefficient in (22) is to incorporate the concept of temperature into the movement of particles. More specifically, as the temperature drops the force that the particles exert on each other decreases resulting in a more intense search in comparison with the early iterations of the algorithm. In addition, in a minimization problem, the charge for a particle i is calculated as(23)qi=nfxworst-fxifxworst-fxbest,∀i∈1,2,…,m,where xworst is the particle with the worst (largest) objective function and other elements are the same as defined in (18). The above equation not only considers the closeness to the best solution but also takes into account the distance from the worst particle. The calculation of Fi and the movement of each particle are the same as in EMA.


## 4.6. The EMA-VNS Hybrid

EMA and VNS have already been successfully combined for tackling problems [15]. Similar to what was mentioned in the previous section, we combine the EMA and VNS by replacing the simple local search in EMA by the VNS.The neighborhood structure in VNS is simple: a neighbor to a solutionS is S′=N2S, which is generated by adding a matrix D of the same size as S to S. The elements of D are normally distributed random variables with a mean of their corresponding elements in S and a variance of σ which is an increasing function of, i, the VNS iteration number, denoted by vi. Similar to what was mentioned above, the infeasible generated neighbors are discarded. Therefore, the general scheme of the EMA-VNS hybrid for a minimization problem can be expressed in Algorithm 3 in which other calculations are as in what was mentioned regarding the original EMA.Algorithm 3:General procedure of the EMA-VNS hybrid.
(1)Initialize the parameters maxIter, maxRep, LSIter and the initial population
(2)t←1
(3)WHILE  t<maxIter  DO
(4)j←1
(5)WHILE  j<LSIter  DO
(6)i←1
(7)WHILE  i<maxRep  DO
(8)FOReach particle S in the population DO
(9)Generate a solution S′=N2S in the neighborhood S
(10)Δ←fS′-fS
(11)IF  Δ≤0  THEN
(12)S′←S
(13)ELSE
(14)σ=vi
(15)END IF
(16)END FOR
(17)i←i+1
(18)END WHILE
(19)j←j+1
(20)END WHILE
(21)      Calculate forces.
(22)     Move particles.
(23)t←t+1
(24)END WHILE


## 5. Numerical Results

### 5.1. A Small Test Problem

In order to investigate the proposed model, a small test problem of size4×4×2×2×2, that is, 4 network nodes, 2 transportation modes, 2 products, and 2 periods, is studied. The problem is coded into LINGO 9.0 and solved by using a PC equipped with an Intel Atom CPU N455 @1.66 GHz and 2.00 GB of RAM running Windows 7 Starter operating system. The obtained cost components are presented in Table 2 and the optimal solution is depicted in Tables 3 and 4 and illustrated in Figure 1. It should be noted that the demand for both products in period 1 in City 3 is satisfied by the facility opened in that node. It is interesting to see that the configuration of the network is changed in the second time period in a way that the demand in City 1 is partially satisfied by products initially transported from City 3 to City 4 and then to City 1. The obtained solution may be justified by the model parameter data presented in the Appendix. The LINGO code is available upon requesting the corresponding author.Table 2
The obtained results for the global optimal solution to the small test problem.
Cost component
Value
Total system cost
184.15
Total facility maintenance cost
3.39
Total facility opening cost
101.02
Total facility closing cost
6.01
Total demand violation cost
—
Total supply violation cost
11.79
Total transportation cost
37.22
Total congestion cost
24.73Table 3
Facility open/close status for planning periods.
Facilities
Period 1
Period 2
City 1
✗
✗
City 2
✓
✗
City 3
✓
✓
City 4
✗
✗Table 4
Product flows for optimal solution.
Product flows
Period 1
Period 2
Mode 1
Mode 2
Mode 1
Mode 2
From
To
Product 1
Product 2
Product 1
Product 2
Product 1
Product 2
Product 1
Product 2
City 1
City 1
—
—
—
—
—
—
—
—
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
—
—
—
—
City 2
City 1
—
—
2.95
1.95
—
—
—
—
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
1.95
1.95
—
—
—
—
—
—
City 3
City 1
—
—
—
—
—
—
1.59
0.81
City 2
—
—
—
—
—
—
0.95
2.95
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
3.31
1.09
—
—
City 4
City 1
—
—
—
—
—
—
0.36
0.14
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
—
—
—
—Figure 1
The illustration of the solution to the small problem.In order to shed light on the impact of congestion in SCND, the model is solved not considering the traffic congestion cost. The results, presented in Table5, indicate that neglecting the congestion costs leads to an increase in the total system cost. Moreover, by not considering the congestion cost, total transportation cost decreases while the total congestion cost increases. This translates to the significant fact that slightly cheaper but more congested routes are selected for transportation of products, and/or cheaper transportation modes that have higher congestion costs are preferred.Table 5
The impact of traffic congestion in SCND problem.
Cost component
Value
Value (neglecting congestion)
Relative percent change
Total system cost
184.15
189.21
2.75%
Total facility maintenance cost
3.39
3.39
0.00%
Total facility opening cost
101.02
101.02
0.00%
Total facility closing cost
6.01
6.01
0.00%
Total demand violation cost
—
—
0.00%
Total supply violation cost
11.79
11.79
0.00%
Total transportation cost
37.22
35.78
−4.02%
Total congestion cost
24.73
31.22
20.79%Here, the monetary value of time plays a significant role. To show its effect, in another experiment the monetary value of time is increased from 0.01 to 0.05 for the same test problem. In this case, the result of neglecting the congestion cost is a 3.66% decrease in total transportation cost and 13.87% and 34.00% increase in the total system cost and total congestion cost, respectively. It is notable that the decrease in total transportation cost is smaller and the increase in the total system cost and the congestion cost is greater in comparison to what is presented in Table5 where monetary value of time is 0.01. Therefore, the higher the monetary value of time is, the more significant it is to consider the congestion cost in SCND.


### 5.2. Comprehensive Experiments

In order to justify the proposed model and the performance of the solution algorithm, several test problems are solved and the results are presented in this section. The test problems are categorized into three groups based on their sizes: small, medium, and large. The sizes of these test problems are shown in Table6.Table 6
The size characteristics of small, medium, and large test problems.
Category
N
M
P
T
Min
Max
Min
Max
Min
Max
Min
Max
Small
10
14
3
5
5
9
3
4
Medium
15
20
6
9
10
15
5
6
Large
11
25
10
13
16
20
7
10In order to further study different hybrids of EMA, a third algorithm comprised of EMA and Tabu Search (TS) is added to the comparison. Following the same structure explained for EMA-SA and EMA-VNS, the local search procedure in EMA is replaced with TS in EMA-TS hybrid. In each problem category, 50 test problems are randomly generated and solved 30 times by EMA-SA, EMA-VNS, and EMA-TS hybrids and the original EMA. The hybrid of EMA and VNS has been shown to outperform other well-known algorithms in the context of SCND [15]. Therefore, to demonstrate the performance of the proposed algorithms, they are used in the comparison. All algorithms are coded by MATLAB R2013a. In each problem category, the parameters of the algorithms are tuned for the largest problem in the category by using a factorial Design of Experiments (DOE) and the Response Optimizer in MINITAB 16.The overall performances of algorithms are compared in Table7 which provides the results for the experiments for small, medium, and large problems and the corresponding mean objective function values, the mean CPU times, and the gap from the best results in each criterion. The statistical significance of the gaps is tested by means of Wilcoxon and the results indicate that the minimum significance is 5% for the gaps. In addition, in this table, the best value in each criterion for each problem group is indicated by boldface numbers.Table 7
The comparison of the algorithms regarding the objective function and computational time.
Problem
Algorithm
OBJ
TIME
OBJ GAP
CPU GAP
Small
EMA-SA
30203.0
83.15
0.000%
18.055%
EMA
39772.9
68.13
24.061%
0.000%
EMA-VNS
30276.3
108.82
0.242%
37.388%
EMA-TS
36967.8
75.91
18.299%
10.238%
Medium
EMA-SA
112212.6
1499.71
0.242%
12.167%
EMA
123348.4
1317.24
9.248%
0.000%
EMA-VNS
111941.0
1541.92
0.000%
14.571%
EMA-TS
121293.9
1493.37
7.711%
11.794%
Large
EMA-SA
237645.1
9048.42
0.004%
15.623%
EMA
250421.9
7634.80
5.106%
0.000%
EMA-VNS
237636.1
9825.95
0.000%
22.300%
EMA-TS
247392.1
10963.4
3.944%
30.362%
All
EMA-SA
124831.5
3477.21
0.053%
15.156%
EMA
135961.7
2950.21
8.235%
0.000%
EMA-VNS
124765.1
3754.09
0.000%
21.413%
EMA-TS
133328.5
4090.02
6.423%
27.868%Generally, considering the solution quality, EMA-VNS outperforms other algorithms (at the cost of a larger computational time), while with less computational effort (at the cost of a lower solution quality), EMA performs better than other algorithms.In order to demonstrate the evolution of the objective function, a problem instance from each problem category is solved by using all four algorithms. The algorithms are started with the same initial population of solutions. The evolution of the objective function value of the small, medium, and large problem instance in each algorithm is depicted in Figures4, 3, and 2, respectively. These figures show that the replacement of the local search procedure in the EMA with SA and VNS results in promising improvements in the quality of the obtained solutions.Figure 2
The evolution of the objective function of the small problem instance for the algorithms.Figure 3
The evolution of the objective function of the medium problem instance for the algorithms.Figure 4
The evolution of the objective function of the large problem instance for the algorithms.In order to further investigate the performance of the four algorithms, a slightly modified version of the Marginal Improvement per CPU (MIC) criterion, originally proposed by Osman [90], is utilized. The MIC for an algorithm α in a problem instance ρ is calculated by the following equation:(24)MICαρ=RPIαρCPUαρ,∀α∈A,∀ρ∈P.In the above equation, RPIαρ and CPUαρ are the Relative Percent Improvement (RPI) and the CPU for algorithm α in a problem instance ρ, respectively. A and P are the set of algorithms and the set of problems, respectively. The RPI is given by(25)RPIαρ=cWρ-cBαρcWρ-cBρ×100,∀α∈A,∀ρ∈Pin which cBρ and cWρ are the best and the worst known objective function values for problem instance ρ, respectively, and cBαρ is the best objective function value for problem ρ obtained by using algorithm α.In order to compare the algorithm by using the MIC criterion, a Mann-Whitney test (two-sample Wilcoxon test) is utilized. The Mann-Whitney test is a nonparametric hypothesis test which is utilized to determine whether two populations have the same population median. This test does not require the data to be sampled from normally distributed populations. TheP values for Mann-Whitney test for the medians of MICs for each pair of algorithms in each problem category are given in Table 8. In this table, each value represents the P value for the following test:(26)H0:ηc=ηrH1:ηc>ηr,where ηc is the median of the MIC of the algorithm in the column and ηr is that of the one mentioned in the row of the table. For example, the significance for the test (27)H0:ηEMA-SA=ηEMAH1:ηEMA-SA>ηEMAfor the small problem category is equal to 0.0005 and can be found in the fifth row and third column of Table 8. From this P value, it can be concluded that the null hypothesis of equality of the MIC median for EMA-SA with that of EMA is rejected with 99.95% confidence. Considering the alternative hypothesis in (27), it can be concluded that EMA-SA is superior to EMA according to MIC criterion.Table 8
The significance for Mann-Whitney test for MIC of the algorithms.
Problem
Algorithm
EMA-SA
EMA-VNS
EMA-TS
Small
EMA-SA
EMA-VNS
0.0625
EMA-TS
0.0050
0.0129
EMA
0.0005
0.0085
0.4001
Medium
EMA-SA
EMA-VNS
0.4725
EMA-TS
0.0087
0.0105
EMA
0.0001
0.0001
0.1352
Large
EMA-SA
EMA-VNS
0.4047
EMA-TS
0.0000
0.0000
EMA
0.0002
0.0002
0.1086
All
EMA-SA
EMA-VNS
0.2568
EMA-TS
0.0000
0.0000
EMA
0.0000
0.0000
0.1947From Table8, it can be concluded that the difference between the MICs of EMA-SA and EMA-VNS is not statistically significant. However, according to this table, EMA-VNS and EMA-SA are both improvements of the EMA in terms of MIC criterion. Furthermore, the MICs of EMA-SA and EMA-VNS hybrids are higher than that of the hybrid of EMA and TS (EMA-TS) and the difference is statistically significant (at least with 90% confidence). In addition, although the EMA-TS performs better than EMA considering the objective function value, due to its relatively high computational effort, the MIC of this hybrid is not higher than that of EMA.The results presented in Table8 show that the hybrid of EMA and VNS outperforms other algorithms in medium and large instances (the most interesting) in terms of solution quality. Given the timeframe of the SCND problem, the execution time is not a concern; therefore, the EMA-VNS hybrid is preferable in our case; however, if the target would be a real-time application, then the EMA-SA hybrid may be considered as the superior algorithm.


## 5.1. A Small Test Problem

In order to investigate the proposed model, a small test problem of size4×4×2×2×2, that is, 4 network nodes, 2 transportation modes, 2 products, and 2 periods, is studied. The problem is coded into LINGO 9.0 and solved by using a PC equipped with an Intel Atom CPU N455 @1.66 GHz and 2.00 GB of RAM running Windows 7 Starter operating system. The obtained cost components are presented in Table 2 and the optimal solution is depicted in Tables 3 and 4 and illustrated in Figure 1. It should be noted that the demand for both products in period 1 in City 3 is satisfied by the facility opened in that node. It is interesting to see that the configuration of the network is changed in the second time period in a way that the demand in City 1 is partially satisfied by products initially transported from City 3 to City 4 and then to City 1. The obtained solution may be justified by the model parameter data presented in the Appendix. The LINGO code is available upon requesting the corresponding author.Table 2
The obtained results for the global optimal solution to the small test problem.
Cost component
Value
Total system cost
184.15
Total facility maintenance cost
3.39
Total facility opening cost
101.02
Total facility closing cost
6.01
Total demand violation cost
—
Total supply violation cost
11.79
Total transportation cost
37.22
Total congestion cost
24.73Table 3
Facility open/close status for planning periods.
Facilities
Period 1
Period 2
City 1
✗
✗
City 2
✓
✗
City 3
✓
✓
City 4
✗
✗Table 4
Product flows for optimal solution.
Product flows
Period 1
Period 2
Mode 1
Mode 2
Mode 1
Mode 2
From
To
Product 1
Product 2
Product 1
Product 2
Product 1
Product 2
Product 1
Product 2
City 1
City 1
—
—
—
—
—
—
—
—
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
—
—
—
—
City 2
City 1
—
—
2.95
1.95
—
—
—
—
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
1.95
1.95
—
—
—
—
—
—
City 3
City 1
—
—
—
—
—
—
1.59
0.81
City 2
—
—
—
—
—
—
0.95
2.95
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
3.31
1.09
—
—
City 4
City 1
—
—
—
—
—
—
0.36
0.14
City 2
—
—
—
—
—
—
—
—
City 3
—
—
—
—
—
—
—
—
City 4
—
—
—
—
—
—
—
—Figure 1
The illustration of the solution to the small problem.In order to shed light on the impact of congestion in SCND, the model is solved not considering the traffic congestion cost. The results, presented in Table5, indicate that neglecting the congestion costs leads to an increase in the total system cost. Moreover, by not considering the congestion cost, total transportation cost decreases while the total congestion cost increases. This translates to the significant fact that slightly cheaper but more congested routes are selected for transportation of products, and/or cheaper transportation modes that have higher congestion costs are preferred.Table 5
The impact of traffic congestion in SCND problem.
Cost component
Value
Value (neglecting congestion)
Relative percent change
Total system cost
184.15
189.21
2.75%
Total facility maintenance cost
3.39
3.39
0.00%
Total facility opening cost
101.02
101.02
0.00%
Total facility closing cost
6.01
6.01
0.00%
Total demand violation cost
—
—
0.00%
Total supply violation cost
11.79
11.79
0.00%
Total transportation cost
37.22
35.78
−4.02%
Total congestion cost
24.73
31.22
20.79%Here, the monetary value of time plays a significant role. To show its effect, in another experiment the monetary value of time is increased from 0.01 to 0.05 for the same test problem. In this case, the result of neglecting the congestion cost is a 3.66% decrease in total transportation cost and 13.87% and 34.00% increase in the total system cost and total congestion cost, respectively. It is notable that the decrease in total transportation cost is smaller and the increase in the total system cost and the congestion cost is greater in comparison to what is presented in Table5 where monetary value of time is 0.01. Therefore, the higher the monetary value of time is, the more significant it is to consider the congestion cost in SCND.


## 5.2. Comprehensive Experiments

In order to justify the proposed model and the performance of the solution algorithm, several test problems are solved and the results are presented in this section. The test problems are categorized into three groups based on their sizes: small, medium, and large. The sizes of these test problems are shown in Table6.Table 6
The size characteristics of small, medium, and large test problems.
Category
N
M
P
T
Min
Max
Min
Max
Min
Max
Min
Max
Small
10
14
3
5
5
9
3
4
Medium
15
20
6
9
10
15
5
6
Large
11
25
10
13
16
20
7
10In order to further study different hybrids of EMA, a third algorithm comprised of EMA and Tabu Search (TS) is added to the comparison. Following the same structure explained for EMA-SA and EMA-VNS, the local search procedure in EMA is replaced with TS in EMA-TS hybrid. In each problem category, 50 test problems are randomly generated and solved 30 times by EMA-SA, EMA-VNS, and EMA-TS hybrids and the original EMA. The hybrid of EMA and VNS has been shown to outperform other well-known algorithms in the context of SCND [15]. Therefore, to demonstrate the performance of the proposed algorithms, they are used in the comparison. All algorithms are coded by MATLAB R2013a. In each problem category, the parameters of the algorithms are tuned for the largest problem in the category by using a factorial Design of Experiments (DOE) and the Response Optimizer in MINITAB 16.The overall performances of algorithms are compared in Table7 which provides the results for the experiments for small, medium, and large problems and the corresponding mean objective function values, the mean CPU times, and the gap from the best results in each criterion. The statistical significance of the gaps is tested by means of Wilcoxon and the results indicate that the minimum significance is 5% for the gaps. In addition, in this table, the best value in each criterion for each problem group is indicated by boldface numbers.Table 7
The comparison of the algorithms regarding the objective function and computational time.
Problem
Algorithm
OBJ
TIME
OBJ GAP
CPU GAP
Small
EMA-SA
30203.0
83.15
0.000%
18.055%
EMA
39772.9
68.13
24.061%
0.000%
EMA-VNS
30276.3
108.82
0.242%
37.388%
EMA-TS
36967.8
75.91
18.299%
10.238%
Medium
EMA-SA
112212.6
1499.71
0.242%
12.167%
EMA
123348.4
1317.24
9.248%
0.000%
EMA-VNS
111941.0
1541.92
0.000%
14.571%
EMA-TS
121293.9
1493.37
7.711%
11.794%
Large
EMA-SA
237645.1
9048.42
0.004%
15.623%
EMA
250421.9
7634.80
5.106%
0.000%
EMA-VNS
237636.1
9825.95
0.000%
22.300%
EMA-TS
247392.1
10963.4
3.944%
30.362%
All
EMA-SA
124831.5
3477.21
0.053%
15.156%
EMA
135961.7
2950.21
8.235%
0.000%
EMA-VNS
124765.1
3754.09
0.000%
21.413%
EMA-TS
133328.5
4090.02
6.423%
27.868%Generally, considering the solution quality, EMA-VNS outperforms other algorithms (at the cost of a larger computational time), while with less computational effort (at the cost of a lower solution quality), EMA performs better than other algorithms.In order to demonstrate the evolution of the objective function, a problem instance from each problem category is solved by using all four algorithms. The algorithms are started with the same initial population of solutions. The evolution of the objective function value of the small, medium, and large problem instance in each algorithm is depicted in Figures4, 3, and 2, respectively. These figures show that the replacement of the local search procedure in the EMA with SA and VNS results in promising improvements in the quality of the obtained solutions.Figure 2
The evolution of the objective function of the small problem instance for the algorithms.Figure 3
The evolution of the objective function of the medium problem instance for the algorithms.Figure 4
The evolution of the objective function of the large problem instance for the algorithms.In order to further investigate the performance of the four algorithms, a slightly modified version of the Marginal Improvement per CPU (MIC) criterion, originally proposed by Osman [90], is utilized. The MIC for an algorithm α in a problem instance ρ is calculated by the following equation:(24)MICαρ=RPIαρCPUαρ,∀α∈A,∀ρ∈P.In the above equation, RPIαρ and CPUαρ are the Relative Percent Improvement (RPI) and the CPU for algorithm α in a problem instance ρ, respectively. A and P are the set of algorithms and the set of problems, respectively. The RPI is given by(25)RPIαρ=cWρ-cBαρcWρ-cBρ×100,∀α∈A,∀ρ∈Pin which cBρ and cWρ are the best and the worst known objective function values for problem instance ρ, respectively, and cBαρ is the best objective function value for problem ρ obtained by using algorithm α.In order to compare the algorithm by using the MIC criterion, a Mann-Whitney test (two-sample Wilcoxon test) is utilized. The Mann-Whitney test is a nonparametric hypothesis test which is utilized to determine whether two populations have the same population median. This test does not require the data to be sampled from normally distributed populations. TheP values for Mann-Whitney test for the medians of MICs for each pair of algorithms in each problem category are given in Table 8. In this table, each value represents the P value for the following test:(26)H0:ηc=ηrH1:ηc>ηr,where ηc is the median of the MIC of the algorithm in the column and ηr is that of the one mentioned in the row of the table. For example, the significance for the test (27)H0:ηEMA-SA=ηEMAH1:ηEMA-SA>ηEMAfor the small problem category is equal to 0.0005 and can be found in the fifth row and third column of Table 8. From this P value, it can be concluded that the null hypothesis of equality of the MIC median for EMA-SA with that of EMA is rejected with 99.95% confidence. Considering the alternative hypothesis in (27), it can be concluded that EMA-SA is superior to EMA according to MIC criterion.Table 8
The significance for Mann-Whitney test for MIC of the algorithms.
Problem
Algorithm
EMA-SA
EMA-VNS
EMA-TS
Small
EMA-SA
EMA-VNS
0.0625
EMA-TS
0.0050
0.0129
EMA
0.0005
0.0085
0.4001
Medium
EMA-SA
EMA-VNS
0.4725
EMA-TS
0.0087
0.0105
EMA
0.0001
0.0001
0.1352
Large
EMA-SA
EMA-VNS
0.4047
EMA-TS
0.0000
0.0000
EMA
0.0002
0.0002
0.1086
All
EMA-SA
EMA-VNS
0.2568
EMA-TS
0.0000
0.0000
EMA
0.0000
0.0000
0.1947From Table8, it can be concluded that the difference between the MICs of EMA-SA and EMA-VNS is not statistically significant. However, according to this table, EMA-VNS and EMA-SA are both improvements of the EMA in terms of MIC criterion. Furthermore, the MICs of EMA-SA and EMA-VNS hybrids are higher than that of the hybrid of EMA and TS (EMA-TS) and the difference is statistically significant (at least with 90% confidence). In addition, although the EMA-TS performs better than EMA considering the objective function value, due to its relatively high computational effort, the MIC of this hybrid is not higher than that of EMA.The results presented in Table8 show that the hybrid of EMA and VNS outperforms other algorithms in medium and large instances (the most interesting) in terms of solution quality. Given the timeframe of the SCND problem, the execution time is not a concern; therefore, the EMA-VNS hybrid is preferable in our case; however, if the target would be a real-time application, then the EMA-SA hybrid may be considered as the superior algorithm.


## 6. Conclusion and Future Works

In this paper, the designing of a multiproduct multimode multiperiod supply chain network considering traffic congestion and both supply-side and demand-side uncertainties is investigated and an optimization model is proposed. In order to justify the proposed model, a small numerical problem is solved by LINGO. The experiments with the numerical example show that neglecting the congestion cost in SCND results in preferring the cheaper but more congested routes for transportation of products, and/or cheaper transportation vehicles that have higher congestion costs are given preference. This fact shows the significance of considering traffic congestion in SCND.In addition, to tackle the complexity of the larger problems, two hybrid Electromagnetism-Like Algorithms are presented and several test problems are solved. The numerical experiments results show that proposed algorithms are promising.The outcomes of the experiment with the small test problem yield the applicability of the proposed model in many real-world situations. Therefore, one of future studies may be comprehensive real case study and sensitivity analyses of an application of the proposed model. Uncertainty may exist in other parameters such as transportation costs. Therefore, a possible future research may be the study of uncertainty in parameters other than those considered in this research.



---

*Source: 2498616-2016-01-27.xml*