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 supplyside and demandside uncertainties are considered. Due to the complexity of such models, conventional solution methods are not applicable; therefore, two hybrid ElectromagnetismLike 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.
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 longterm 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 closedloop supply chain network design for lead/acid batteries. They considered the risk of endoflife 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 endoflife products. Osmani and Zhang [10] presented a sustainable dual feedstock lignocellulosicbased bioethanol SCND model considering uncertainty in demand, supply, and prices of products. A robust closedloop 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 ElectromagnetismLike 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 demandside and supplyside 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 Scenariobased/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 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 largescale 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 twoechelon 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 CardonaValdé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 closedloop global SCND under uncertainty. A nondominated sorting genetic algorithm (NSGAII) was utilized by Pasandideh et al. [35] to optimize a multiproduct multiperiod threeechelon supply chain problem under uncertainty. Govindan et al. [15] proposed a robust hybrid multiobjective metaheuristic comprised of ElectromagnetismLike Algorithm (EMA) and VNS as a solution approach to their SCND model.
Among the metaheuristics, EMA is a relatively new populationbased algorithm proposed by Birbil and Fang [36] that has been applied to many optimization problems (such as set covering [37], machine scheduling [38], flowshop scheduling [39, 40], project scheduling [41], periodic jobshop 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 wellknown, 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 parametersensitive 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 demandside and supplyside uncertainties along with the time varying interest rate (time value of money) in SCND,(2)proposing and studying two hybrid populationbased metaheuristics (EMAVNS and EMASA) 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 supplyside 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.
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.
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 , then the superiority of a fuzzy triangular number over a fuzzy triangular number is defined by the following equation:In the above equation, is the center and and are the left and right spreads of , respectively, and is the center and and are the left and right spreads of , respectively. Similarly, the inferiority of to 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.
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 twoechelon 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
: 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. Parameters
: interest rate in time period , : monetary value of time in time period , : parameter in the link performance function, : parameter in the link performance function, : capacity of vehicle , : traffic congestion coefficient of vehicle , : the annual maintenance cost of a facility in candidate location in time period , : the opening cost of a facility in candidate location in time period , : the closing 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 , : the left spread of the triangular fuzzy demand for product , in demand node , in time period , : 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 , : the left spread of the triangular fuzzy supply for product , in facility node , in time period , : the right spread of the triangular fuzzy supply for product , in facility node , in time period , : free flow travel time from node to node in time period , : traffic capacity of the link from node to node in time period , : 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 , : 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, : 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.
3.3.4. Decision Variables
: 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 .
3.3.5. Auxiliary Variable
: the inferiority of the demand in a demand node for product in time period to the amount of that product transported to that node, : 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 , : 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 , : total traffic flow from node to node , : the total amount of product transported to a node in time period , : the total amount of product transported from a node in time period .
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: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) asEquation (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 ElectromagnetismLike 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 ElectromagnetismLike 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] aswhere 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 aswhere 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], flowshop scheduling [88], jobshop 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 parametersensitive 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 5dimensional 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, , which is equal to 1 if a facility is operating in node in time period and equals 0 otherwise. If , then no facility is opened in node in time period . Therefore, we have , , , , and apparently according to the definition of , we set , , , .
4.5. The EMASA 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 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 EMASA 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 aswhere 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 . 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 aswhere 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.
4.6. The EMAVNS 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 solution is , 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 . Similar to what was mentioned above, the infeasible generated neighbors are discarded. Therefore, the general scheme of the EMAVNS 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.

5. Numerical Results
5.1. A Small Test Problem
In order to investigate the proposed model, a small test problem of size , 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.
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.
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 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 EMASA and EMAVNS, the local search procedure in EMA is replaced with TS in EMATS hybrid. In each problem category, 50 test problems are randomly generated and solved 30 times by EMASA, EMAVNS, and EMATS hybrids and the original EMA. The hybrid of EMA and VNS has been shown to outperform other wellknown 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, EMAVNS 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 Figures 4, 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.
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 equation, and 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. The RPI is given byin which and are the best and the worst known objective function values for problem instance , respectively, and is the best objective function value for problem obtained by using algorithm .
In order to compare the algorithm by using the MIC criterion, a MannWhitney test (twosample Wilcoxon test) is utilized. The MannWhitney 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 MannWhitney 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 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 EMASA with that of EMA is rejected with 99.95% confidence. Considering the alternative hypothesis in (27), it can be concluded that EMASA is superior to EMA according to MIC criterion.
From Table 8, it can be concluded that the difference between the MICs of EMASA and EMAVNS is not statistically significant. However, according to this table, EMAVNS and EMASA are both improvements of the EMA in terms of MIC criterion. Furthermore, the MICs of EMASA and EMAVNS hybrids are higher than that of the hybrid of EMA and TS (EMATS) and the difference is statistically significant (at least with 90% confidence). In addition, although the EMATS 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 EMAVNS hybrid is preferable in our case; however, if the target would be a realtime application, then the EMASA 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 supplyside and demandside 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 ElectromagnetismLike 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 realworld 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.
Appendix
The data for the small test problem are presented here. The traffic capacity, transportation costs, link flows, and so forth, from a city to itself, can be interpreted as the local transportation properties such as total capacity of urban roads within the city. Similar discussion can be presented regarding the transportation costs, basic link flow, and free flow travel times.
It should be noted that the aforementioned data are generated randomly. In addition, the monetary value of time and the interest rate are assumed to be constant for both periods and equal 0.01 and 0.1, respectively. and are set to 0.15 and 4, respectively. Furthermore, all violation cost coefficients for demand and supply are set to 1 and the spreads for all fuzzy numbers are set to 0.1. The data are summarized in Tables 9, 10, 11, 12, 13, 14, 15, and 16.
Conflict of Interests
The authors declare that there is no conflict of interests.