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

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., availablemeans of transportation, and the variety of products) and exogenous criteria (e.g., the environmental uncertainty, and transportation system conditions). In this paper, amixed integer nonlinear programming (MINLP)model for dynamic design of a supply chain network is proposed. In thismodel,multiple products andmultiple transportationmodes, the time value ofmoney, 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.


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 endof-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][20][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 Scenariobased/Robust Optimization (RO) [24][25][26][27][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 largescale 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 prioritybased 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.
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][51][52][53][54][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][58][59][60][61][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, 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 Section 3, 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.

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 = (, , ) over a fuzzy triangular number Q = (V, , ) is defined by the following equation: In the above equation,  is the center and  and  are the left and right spreads of P = (, , ), respectively, and V is the center and  and  are the left and right spreads of Q = (V, , ), respectively.Similarly, the inferiority of P to Q is calculated by the following equation: 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 Sections 3.3 and 3.4.

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.

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.

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.
(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. : the total amount of product  ∈  transported from a node  ∈  in time period  ∈ .
In 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 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.

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.

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 Algorithm 1  is the number of solution points (particles),  is the maximum number of the algorithm iterations,  is the maximum number of local search steps, and  and  are local counter variables.
The original scheme of EMA is presented in Algorithm 1 in which the forces are calculated according to Coulomb law [86] as where   is the force particle  and particle  exert on each other,   is the distance between those particles,  is the number of particles, and   and   are the electrical charges of particle  and particle , respectively.In EMA, the charges of the particles are calculated as where () is the objective function value for solution  and  is the dimension of the solution space.The force exerted on each particle  is then obtained according to the following equation: 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 ‖  ‖ is the norm of the force vector: 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.

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][54][55].

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][58][59][60][61].

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, , of size || × || × || × || × ||, where || represents the number of elements in the set , is considered.In , each element [, , , , ] has the same interpretation as   .It should be noted that  also contains the information about the decision variable   .For example, consider the binary decision variable,  23 , which is equal to 1 if a facility is operating in node 2 ∈  in time period 3 ∈  and equals 0 otherwise.If  23 = 0, then no facility is opened in node 2 ∈  in time period 3 ∈ .Therefore, we have  23 = 0, ∀ ∈ , ∀ ∈ , ∀ ∈ , and apparently according to the definition of , we set [2, , , , 3] = 0, ∀ ∈ , ∀ ∈ , ∀ ∈ .

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, , in SA is a solution   =  1 () that is generated by adding a matrix , with randomly generated element, of the same size as  to .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.
In Algorithm 2, the force particle  and particle  exert on each other,   , is calculated as where  is obtained through the following equation and other elements are the same as defined in (17): It should be noted that, in EMA, we have  = 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  is calculated as where  worst 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   and the movement of each particle are the same as in EMA.
The neighborhood structure in VNS is simple: a neighbor to a solution  is   =  2 (), which is generated by adding a matrix  of the same size as  to .The elements of  are normally distributed random variables with a mean of their corresponding elements in  and a variance of  which is an increasing function of, , the VNS iteration number, denoted by V().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.

A Small Test Problem.
In order to investigate the proposed model, a small test problem of size 4 × 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 2.95 units of product 1 1.95 units of product 1 1.59 units of product 1 0.95 units of product 1 and 1.95 units of product 2 and 1.95 units of product 2 and 0.81 units of product 2 and 2.95 units of product 2 using 41 vehicles of type 2 using 40 vehicles of type 1 using 20 vehicles of type 2 using 33 vehicles of type 2 0.36 units of product 1 and 0.14 units of product 2 using 5 vehicles of type 2 3.31 units of product 1 and 1.09 units of product 2 using 44 vehicles of type 1 Time period 1 Time period 2 Figure 1: The illustration of the solution to the small problem.
the Appendix.The LINGO code is available upon requesting the corresponding author.
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 Table 5, 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.
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 Table 5 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.

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 Table 6.
In 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 Table 7 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.
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  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: In the above   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.The  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  value for the following test: where   is the median of the MIC of the algorithm in the column and   is that of the one mentioned in the row of the table.For example, the significance for the test  0 :  EMA-SA =  EMA  1 :  EMA-SA >  EMA (27) for 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  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.
From Table 8, 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 Table 8 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.

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 Hybrid EMAsMP: 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.

Algorithm 2 :
General procedure of the EMA-SA hybrid.

5 Figure 2 :Figure 3 :
Figure 2: The evolution of the objective function of the small problem instance for the algorithms.

Figure 4 :
Figure 4: The evolution of the objective function of the large problem instance for the algorithms.

Table 1 :
A comparison of the recent related literature and our research (chronologically sorted).
rs : the right spread of the triangular fuzzy demand for product  ∈ , in demand node  ∈ , in time period  ∈ ,    : the center of the triangular fuzzy supply for product  ∈ , in facility node  ∈ , in time period  ∈ ,  ls  : the left spread of the triangular fuzzy supply for product  ∈ , in facility node  ∈ , in time period  ∈ ,  rs  : the right spread of the triangular fuzzy supply for product  ∈ , in facility node  ∈ , in time period  ∈ , FFTT  : free flow travel time from node  ∈  to node  ∈  in time period  ∈ , TCAP  : traffic capacity of the link from node  ∈  to node  ∈  in time period  ∈ , BF  : basic traffic flow of the link from node  ∈  to node  ∈  in time period  ∈ ,   : the transportation cost for vehicle  ∈  for link from node  ∈  to node  ∈  in time period  ∈ ,  (, in  ) : the unit cost of inferiority of the demand in a demand node  ∈  for product  ∈  in time period  ∈  to the amount of that product transported to that node,  (, out  ) : the unit cost of superiority of the capacity in a supply node  ∈  for product  ∈  in time period  ∈  over the amount of that product transported from that node.  : a binary variable which is equal to 1 if a facility is operating in  ∈  in time period  ∈  and equals 0 otherwise,   : the amount of product  ∈  transported from node  ∈  to node  ∈  by means of transportation vehicle  ∈  in time period  ∈ . (, in  ) : the inferiority of the demand in a demand node  ∈  for product  ∈  in time period  ∈  to the amount of that product transported to that node,  (, out     : being 1 if a facility is closed in candidate supply location  ∈  in time period  ∈ ,   : the number of vehicles of type  ∈  transporting product  ∈  from node  ∈  to node  ∈  in time period  ∈  minus 1,   : the total number of vehicles of type  ∈  transporting goods from node  ∈  to node  ∈  in time period  ∈ ,    : the auxiliary variable used form linearization, determining the amount of product  ∈  transported by vehicle of type  ∈  from node  ∈  to node  ∈  in time period  ∈  if a facility is operating in node  ∈ , FLW  : total traffic flow from node  ∈  to node  ∈ ,  in  : the total amount of product  ∈  transported to a node  ∈  in time period  ∈ ,  out : the set of network nodes (including candidate facility locations and demand points), : the set of time periods in the planning horizon, : the set of products, : the set of transportation modes.3.3.2.Subscripts: subscript for candidate facility location ( ∈ ), : subscript for demand point ( ∈ ), : subscript for time period ( ∈ ), : subscript for product ( ∈ ), : subscript for transportation vehicle ( ∈ ).3.3.3.ParametersIR  : interest rate in time period  ∈ , MVT  : monetary value of time in time period  ∈ , : parameter in the link performance function, : parameter in the link performance function, CAP  : capacity of vehicle  ∈ , TCC  : traffic congestion coefficient of vehicle  ∈ , MC  : the annual maintenance cost of a facility in candidate location  ∈  in time period  ∈ ,  : the center of the triangular fuzzy demand for product  ∈ , in demand node  ∈ , in time period  ∈ ,  ls  : the left spread of the triangular fuzzy demand for product  ∈ , in demand node  ∈ , in time period  ∈ , ) : the superiority of the capacity in a supply node  ∈  for product  ∈  in time period  ∈  over the amount of that product transported from that node,    : being 1 if a facility is opened in candidate supply location  ∈  in time period  ∈ ,

Table 2 :
The obtained results for the global optimal solution to the small test problem.

Table 3 :
Facility open/close status for planning periods.

Table 4 :
Product flows for optimal solution.

Table 5 :
The impact of traffic congestion in SCND problem.

Table 6 :
The size characteristics of small, medium, and large test problems.

Table 7 :
The comparison of the algorithms regarding the objective function and computational time.

Table 8 :
The significance for Mann-Whitney test for MIC of the algorithms.
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.
equation, RPI  and CPU  are the Relative Percent Improvement (RPI) and the CPU for algorithm  in a problem instance , respectively. and  are the set of algorithms and the set of problems, respectively.

Table 10 :
The data for links traffic capacities.

Table 11 :
The data for links basic flows.