1.INTRODUCTION
In the current industrial world, transportation plays a major role in economic development of countries. Due to high frequency of transportation activities, transportation costs constitute a high percentage of logistic costs (between 30 to 60%) (Ghiani et al., 2004). Therefore, efficient transportation throughout all supply chain is highly significant. The design of supply chain network includes location of distribution centers, assigning flow among facilities and the vehicle routing problem (VRP) which could be considered simultaneously. As the number of vehicles in supply chain and FCR of vehicles increase, the destructive effects of onroad transportation on the environment will also rise. The increase of greenhouse gases has numerous risks and effects on life of people and our planet. The emission of these gases is mostly due to consumption of fuel. As a result, for preventing from the environmental pollution and reduction of destructive environmental effects of transportation, the number of ve hicles, the FCR of them and the distance a vehicle travels should be reduced.
The VRP is a complicated problem and solving it by exact methods is highly timeconsuming. In fact, by increasing the dimensions of the problem, the necessary time for computations to attain optimal solution increases exponentially. Therefore, solving large dimensions of the problem demands using heuristic and metaheuristic methods. In recent years, in addition to the heuristic methods for solving the VRP, metaheuristic methods and algorithms such as genetic algorithm, SA, VNS and etc. have been developed. Dantzig and Fulkerson (1954) used binary integer programming for solving the VRP methods. Dantzig and Ramser (1959) studied the dynamic programming solution method. Fisher (1981) offered the Lagrange Relaxation Method for solving the integer programming problems. Laport and Nobert (1987) conducted a review study and collected all of the methods applied to that date and categorized them into approximate and exact solution methods. Laport et al. (1992) extended the branch and bound methods for solving the VRP. Barmel and SimchiLevi (1997) studied branch and bound method for solving the VRPTW models. Also, Kuo and Wang, (2012), Xu et al. (2012), Salhi et al. (2013) and Xu and Jiang (2014) developed and solved different VRP problems through VNS algorithm.
Facility location problem (FLP) and VRP are the most challenging problems of supply chain management. The appearance of location routing problem (LRP) dates back to late 1970s and early 1980s. LRP has been studied by many researchers in different decades. Min et al. (1998) and Nagy and Salhi (2007) reviewed the LRP field literature and their work obtained the reference status among similar papers in the field. Recent studies in this field include the studies by Balakrishnan et al. (1987), Laporte (1988), Laporte (1989), Berman et al. (1995), Min et al. (1998), Nagy and Salhi (2007) and Prodhon and Prins (2014). MartínezSalazar et al. (2014) studied the biobjective transportation location routing problem for determining the vehicles’ optimal load and minimization of transportation costs; they used metaheuristic algorithms to solve the problem. Hamidi et al. (2012) and Hamidi et al. (2014) studied the LRP in four levels by considering multiple products. In their papers, the limitations of route length, number of vehicles, capacity of vehicles and fixed and variable costs of the problem were considered. In addition, Hamidi et al. (2014) proposed an improved and developed heuristic algorithm for a multiproduct fourechelon capacitated LRP.
Because of the NPhard nature of LRP problems, the application of heuristic and metaheuristic algorithms has drawn the attention of many researchers. Tuzun and Burke (1999) used twophase Tabu search algorithm to solve the LRP while distribution centers had indefinite capacities. During the problem solving procedure, the algorithm alternatively transferred the information between two FLP and VRP phases. Yu et al. (2010) used SA algorithm to solve CLRP. For solving the problem and validation of the proposed algorithm, they used the trial examples of LRP literature. The obtained results signified the superiority of the proposed SA algorithm compared with other algorithms introduced in the literature.
To bring the classic VRP closer to the realworld conditions and addressing the challenges frequently encountered by the managers, the present study aims to offer an integrated model in which the routing, location and supplier selection problems are simultaneously considered and pollutants emission is reduced. Due to the fact that the problem is NPhard, the metaheuristic algorithms of SA, TS, BA, and VNSSA are used in this study to solve the proposed problem. The metaheuristic algorithms employed in the present study are among the most efficient algorithms known ever in the field. In the next section, the problem statement and explanation of the proposed mathematical model are presented. In the third section, details of the proposed approaches to solve the model are explained. In the fourth section, to measure the accuracy and efficiency of the introduced mathematical model and metaheuristic algorithms a number of numerical examples based on problem structure are provided and solved by GAMS and MATLAB software. Conclusion and further suggestions are presented in the last section of the paper.
2.PROBLEM DESCRIPTION AND FORMULATION
In the present study, a review of the literature and a determination of the defects of this research field have been done to suggest an integrated and comprehensive model and as a result, to resolve some of the relevant defects. The proposed model intends to consider VRP, FLP and supplier selection problems in a threeechelon supply chain simultaneously. Moreover, by optimizing the FCR of vehicles, the model tries to reduce the destructive effects of road transportation in the chain. This supply chain includes production centers, distribution centers and clients. In this problem, each distribution center is authorized to be supported by multiple suppliers while the client is authorized to be served by one distribution center. In addition, potential locations for determining the distribution and production centers are definite. Potential routes for travelling of vehicles are defined and for each potential route, the maximum flow capacity is determined. In this model, the maximum capacity of each vehicle is fixed and definite and decision making regarding the vehicles’ capacity is associated with network design. The demand of geographically dispersed customers is presumed to be constant and defined. Each client is served by one of the existing routes and one vehicle. In the problem, limitation of timewindows is presumed to respond to the clients’ demands. Different vehicles have different FCRs as well as different transportation costs for serving the clients, meanwhile each vehicle will have certain weight and volume capacity.
Transportation increases the costs of environmental pollution, sound pollution and the costs of removing natural views. As the number of vehicles and FCR of vehicles increase, the destructive effects of road transportation on environment will also increase. The escalation of pollutants is accompanied by numerous risks and consequences for the human life and our planet. The emission of these gases is mostly due to fuel consumption. Therefore, one should optimize the sum of fuel consumption of vehicles so as to reduce the destructive effects of road transportation and the costs of fuel consumption. In this section, the suggested mathematical model is explained. In addition, the method of calculating the FCR and the linearization of the mathematical model are described.
2.1.Mathematical Notation
The mathematical notation that will be used in this paper is listed as follows:

(1) Sets and indices

I : set of potential producers (i∈I) ; i index of potential producers

J : set of potential distribution centers (j∈J) ; j index of potential distribution centers

K : set of customers (k∈K) ; k index of customers

V : set of vehicles belonging to distribution centers (v∈V) ; v index of vehicles belonging to distribution centers

E, F : set of producers and distribution centers (E, F = J ∪ K) ; e, f index of producers and distribution centers


(2) Parameters

D_{j} : total demands for distribution center j, for all j ∈ J

D_{k} : total demands for customer k, for all k ∈ K

d_{ij} : distance from producer i to distribution center j (Kilometer), for all i∈ I, j ∈ J

d_{ef} : distance from point e to point f (Kilometer), for all j ∈ J, j ∈ J

CapV^{M} : capacity of vehicles belonging to producers

$Cap{V}_{v}^{D}$ : capacity of vehicle v belonging to distribution centers, for all v∈V

C a p_{i} : production capacity at producer i, for all i∈I

C a p_{j} : supply capacity at distribution center j, for all j ∈ J

t_{ij} : mean duration of moving from producer i to distribution center j for all i∈I, j∈J

t_{ef} : mean duration of moving from point e to point f, for all e∈E, f∈F

${t}_{j}^{*}$ : mean time of unloading in distribution center j for all j∈J.

τ_{j}: maximum time allowed for procurement demand of distribution center j, for all j∈J

S_{k} : duration of serving the customer k, for all k ∈ K

$\left[{L}_{k},\hspace{0.17em}{U}_{k}\right]$ : range of time window to serve the customer k, for all k ∈ K

TC_{ij} : transportation cost of each unit of product, from producer i to distribution center j for all i∈I, j∈J

TC_{efv} : transportation cost of each unit of product, from point e to point f by vehicle v for all e∈E, f ∈ F

FCR_{ij} : fuel consumption rate form producer i to distribution center j (Lit/km) for all i∈I, j∈J

FCR_{e f v} : fuel consumption rate from point e to point f by vehicle v (Lit/km) for all e∈E, f ∈ F, v∈V

CG : the cost of losing reputation due to the violation of time window constraint per unit time

CW : the waiting costs of vehicles until starting service

C_{0} : the cost of each unit of vehicle fuel (liter)

H_{o j} : the cost of inventory holding for distribution center j, for all j∈J

Λ_{j} : the cost of commodity shortcoming for distribution center j, for all j∈J

N: a large number


(3) Decision variables

γ_{ij} : amount of flow from producer i to distribution center j, for all i∈I, j∈J

γ_{efv} : amount of flow from point e to point f by vehicle v, for all e∈E, f∈F, k∈K

In_{j} : inventory of distribution center j at the end of the course, for all j∈J

Br_{j} : shortage of distribution center j at the end of the course, for all j∈J

B_{k ,v} : arrival time of vehicle v to customer k, for all k∈K, v∈V

ZW_{k ,v} : waiting time of vehicle v for serve to the customer k, for all k∈K, v∈V

ZG_{k,v} : the amount of violation from time window of vehicle v to serve to the customer k, for all k∈K, v∈V

Car_{jv} : 1: if vehicle v allocated to distribution center j, for all j∈J, v∈V
0: otherwise;

R_{e, fv} : 1: if point e will be before point f upon vehicle v, for all e∈E, f ∈F, v∈V
0: otherwise;

Y_{ij} : 1: if producer i allocated to distribution center j, for all i∈I, j∈J
0: otherwise;

Y_{jk} : 1: if customer k allocated to distribution center j, for all j∈J, k ∈ K
0: otherwise;

2.2.The Mathematical Model Formulation
The proposed model is formulated as the following mixed integer nonlinear programming (MINLP) model:
Objective function
Equation (1) signifies the objective function with nine terms. The first expression refers to transportation costs of each unit of product from producer i to distribution center j. The second expression includes transportation costs of each unit of product from distribution center j to the client k by vehicle v. The third expression refers to transportation costs of each unit of product from client e to the client f by the vehicle v. It should be noted that the transportation costs of each unit of product do not cover fuel consumption costs and they only include driver’s fee, depreciation or maintenance of the vehicle. The fuel consumption costs of the vehicles are considered in fourth to sixth expressions. The seventh expression includes the inventory cost and costs of shortage of a product in distribution center j. Finally, the eighth and ninth expressions include penalties of deviation from timewindows (waiting costs or delay costs).
The equations (2)(33) show the model limitations. Equation (2) guarantees that the extent of passing flow from producer i to distribution center j is less than capacity of producer i and capacity of vehicles. Equation (3) signifies that there is a passing flow from producer i to distribution center j just when producer i is assigned to distribution center j. The equation (4) guarantees that the total flow sent from producer i to all assigned distribution centers is less than the capacity of producer i. The equation (5) guarantees that total extent of flow assigned to producer i is less than the capacity of producer i. The equation (6) manifests that total flow sent to distribution center j should be equal with its demand. The equation (7) guarantees that the extent of passing flow from distribution center j to all clients is less than or equal to the distribution center j capacity. The equations (8) and (9) guarantee that each client is served only by one distribution center. The equation (10) guarantees that the extent of flow of vehicle v from distribution center j to client k is less than or equal to the distribution center j and vehicle v capacities respectively. The equation (11) denotes that if the k^{th} client is not assigned to the distribution center j, there will be no flow between distribution center j and k^{th} client. The equation (12) guarantees that the total passing flow from distribution center j is less than the capacity of distribution center j. The equation (13) shows that the total passing flow from the distribution centers to k^{th} client is equal to the demand by the k^{th} client. The equations (14) and (15) show that an individual route cannot be used by multiple distribution centers (depot) or vehicles. In fact, they guarantee that each distribution center and each vehicle use only one route to respond to the orders. The equations (16)(19) show that each node of the client will have one input flow and one output flow at most and the client should be serviced from only one vehicle. Actually, they guarantee that each client is served by one vehicle and one distribution center only. The equation (20) shows that if one flow (route) leads to k^{th} client, one flow (route) should leave the client too. Additionally, it guarantees that each vehicle moves from one distribution center and returns to the same distribution center after service and delivery of orders to clients. The equation (21) refers to the fact that k^{th} client is assigned to distribution center j if and only if one route exists between the k^{th} client and the distribution center j. The equation (22) suggests that all demands responded by vehicle v should be less than or equal to its capacity. The equation (23) suggests that the level of flow sent to k^{th} client should be precisely equal to the demand by the k^{th} client. The equation (24) guarantees that the whole flow between distribution centers and the clients responded by a vehicle does not exceed the capacity of the vehicle. The equation (25) shows that the vehicle v is assigned to distribution center j. If the vehicle v is assigned to distribution center j, the equation (26) indicates that the difference between the two demands by the distribution center j and total flow sent from distribution center j is equal to the difference between the inventory and shortage of product distribution center j. The equation (27) shows that the duration in which the i^{th} producer is authorized to deliver the order of distribution center j should be less than maximum permissible duration. The equations (28) and (29) determine the time that the ve hicle v arrives at the k^{th} client. The durations ef t and ij t are calculated based on the average speed limits of the routes ef and ij. They are considered as parameters of the proposed model. In addition, for estimation of${t}_{j}^{*}$the mean duration of unloading in distribution center j is determined and${t}_{j}^{*}$is regarded as the parameter of the proposed model. The equations (30) and (31) determine the deviation from time windows for serving the clients. The equations (32) and (33) show that R_{e, f v}, R_{j, k v}, Y_{ij}, Y_{jk} are binary decision variables. Besides, they show that the decision variables ${B}_{k,\hspace{0.17em}v},\hspace{0.17em}{\gamma}_{ij},\hspace{0.17em}{\gamma}_{jk},\hspace{0.17em}{\gamma}_{e,\hspace{0.17em}f\hspace{0.17em}v}\hspace{0.17em}\text{and}\hspace{0.17em}Z{W}_{k,v},\hspace{0.17em}Z{G}_{k,v}$ are greater than or equal to zero.
2.3.FCR Computation
The FCR computation of vehicles is highly complex and dependent to numerous factors. Barth et al. (2005), Barth and Boriboonsomsin (2009) and Bektas and Laporte (2011) surveyed the FCR computation methods for the vehicles. In the present paper, the simplified form of their computations is used for determining FCR. The FRC calculation method is briefly discussed below.(34)(35)
In above, k indicates engine friction factor (kilo joule /rev/liter), N indicates engine speed (Rev/second), V indicates engine displacement, P indicates engine power output (kilowatt), η indicates amount of diesel motor proficiently, τ indicates fuel to air mass ratio, π indicates heating value of one kind of diesel fuel (kilo joule/gram), ε indicates vehicle drive train efficiency, P_{a} indicates necessary power of vehicle engine for discharging energy caused by accessory parts of vehicles (like air condition and cooler) and P_{TF} indicates total force necessary for pulling a vehicle (kilowatt) and can be calculated as follows:(36)(37)(38)
In above, M indicates total mass of vehicle accompanied by load (kilogram), v indicates speed of vehicle (meter/second), a indicates acceleration of vehicle (meter/ second2), g indicates gravitational constant (9.81 meter/ second2), θ indicates road angle, A indicates frontal surface area of the vehicle (meter^{2}), ρ indicates air density (kilogram/meter^{3}), C_{d} indicates coefficient of the aero dynamic drag, C_{r} indicates coefficient of the rolling resistance, α indicates special fixed for each vehicle route, β indicates special fixed vehicle. Consequently P_{TF} becomes as follows:(39)
It is assumed that all the parameters for a vehicle traveling along a path or arc remain fixed only the speed and load of the vehicle change along a route. In other word given vehicle with average speed v = v_{ij} during distance d_{ij} and gradient of the route θ = θ_{ij} carries the total mass M. M is equal to total weight of empty vehicle and load mass which the vehicle carries along the route (i, j). Where λ =1/(1,000.εη) and ω = τ / (π.ξ) are constants and ξ is the conversion factor from unit (gram per second) to unit (liter per second). Therefore the fuel consumption rate formula for each arc (i, j) is as follows:(40)
Then, final FCR function in unit of liter per second is as follows:(41)
Function of FCR per distance traveled is Ushaped; hence, the calculation of vehicle’s optimal speed with lowest FCR would be possible. Description of vehicle’s function behavior is shown in Figure 1. This function consists of two components, the first one is the phrase kNVω / v which is represented by the dashed line, and the second one is the phrase $\left(M\alpha v\lambda +\beta {v}^{3}\lambda +{P}_{a}/n\right)\omega /v$ which is represented by the dotted line in Figure 1. As illustrated in Figure 1, by increasing the vehicle speed, the amount of kNVω / v will decrease and $\left(M\alpha v\lambda +\beta {v}^{3}\lambda +{P}_{a}/n\right)\omega /v$ will increase. The illustration shown in FCR’s function Figure depends on different elements and factors such as the amount of vehicle load, road gradient, etc. It can be observed that the optimal vehicle speed is about 45 kilometers per hour. So for maximum fuel efficiency of the vehicles and also reducing the amount of pollutants and harmful gases emissions, a limit of 45 kilometers per hour can be recommended for the vehicles speed. In Figure 1 a chart of fuel consumption rate is shown based on the speed function (Bektas and Laporte, 2011). In this chart the speed of vehicle varies from 20 kilometer/hour to 120 kilometer/ hour. The FCR function depicted in Figure 1 is the sum of kNVω / v and $\left(M\alpha v\lambda +\beta {v}^{3}\lambda +{P}_{a}/\eta \right)\omega /v$ components.
As a result, the final formula of $FCR\left(v\right)=kNV\omega /v+\left(M\alpha v\lambda +\beta {v}^{3}\lambda +{P}_{a}/\eta \right)\omega /V$ through which the parameters FCR_{ij} and FCR_{efv} are calculated for suggested model in this paper. The unit of the formula is liter per kilometer.
2.4.Linearization of the Suggested Model
The proposed model is nonlinear and because nonlinear models are harder to solve than linear ones, in the present section we try to linearize the proposed model through linearizing the nonlinear expressions. There are two types of nonlinear expressions in the model. In the first type, two binary variables are multiplied together and in the second type a binary variable is multiplied by continuous variable.

(1) Linearization of the product of two binary variables
Let’s suppose that Z = x_{1} × x_{2} and multiplication of two binary variables is 1. In other words, Z is equal to 1 if both variables are assigned the value 1; otherwise, it will be zero. Based on the study by Norouzi et al. (2012), we could use the following auxiliary limitations to turn nonlinear limitations into a set of linear limitations.(42)(43)(44)

(2) Linearization of the product of a binary and a continuous variable
Let’s suppose Z = x_{1}× x_{2} is multiplication of binary variable x_{1} by continuous variable x_{2} . In this case, if the binary variable is 1, the variable Z is equal with value of continuous variable. Otherwise, it will be zero. Based on the study by Glover and Woolsey (1974), we could use the following auxiliary limitations to turn nonlinear limitations into a set of linear limitations.(45)(46)(47)
In the suggested model, objective function and limitations (28) and (29) include nonlinear expressions. In the following, linearization of nonlinear expressions is done.
2.4.1.Linearization of the Objective Function
The objective function has nonlinear expressions. In these nonlinear expressions, a binary variable is multiplied by a continuous variable. For better clarity, the nonlinear expressions of objective function are rewritten in equations (48), (53) and (58) and nonlinear sections of the expressions are placed in brackets. For linearization of nonlinear expressions of objective function, three continuous variables ${\text{F}}_{\text{ij}}$ , ${\text{F}}_{\text{jkv}}$ , and ${\text{F}}_{\text{efv}}$ are defined for the model and nine limitations of (50)(52), (55)(57) and (60)(62) are added to limitations of the model.(49)(51)
2.4.2.Linearization of the nonlinear constraints
In the suggested model, the limitations (28) and (29) have nonlinear expressions. In these nonlinear expressions, a binary variable is multiplied by a continuous variable. The rewritten limitations are represented in equations (63) and (68) in which the nonlinear parts are represented in brackets. To linearize the nonlinear limitation (28), the continuous variable β1_{ef} is defined while for nonlinear limitation (29), two continuous variables β2_{ef} and β3_{ef} are defined in the model. In addition, the equations (64)(67) and (69)(75) are added to model limitations. Considering the fact that in limitation 36, the sum of discrete variables R_{j,kv} could be 1 at most, we could consider $\sum}_{j\notin J}{R}_{j,k\hspace{0.17em}v$ as a binary variable.(65)(66)(70)(71)(72)(73)(74)
3.SOLUTION METHODOLOGY
The VRP is among problems with hard complexity (Lenstra and Rinnooy Kan, 1981). So, the proposed model in the present paper is among the NPhard integer programming and due to the complexity and largeness of dimensions, such problem is solved through using the approximate methods. This leads to proper approximate solutions in reasonable computation duration. To solve the suggested model, the metaheuristic algorithms SA, TS, BA and VNSSA are used. In the following, the metaheuristic algorithms used are discussed in further details.
3.1.Hybrid Variable Neighborhood Search Algorithm Based on Simulated Annealing (VNSSA)
The main idea of the VNSSA algorithm is the replacement of neighborhoods structure for inhibition of entrapment in local optimal during algorithm searching. This algorithm systematically searches the solution space by changing the neighbors of the current solution. The algorithm starts by obtaining the initial solution S_{0} and initial solution is presumed to be the current solution of the system. The algorithm is comprised of two phases of shaking and local search. The algorithm starts with first neighborhood structure k = 1. First, the algorithm moves from current solution S_{0} to the neighbor solution S_{1} based on k^{th} neighborhood structure in shaking phase and will randomly select a neighbor among the closest neighbors. Then, in local search phase, the algorithm moves from the solution S_{1} to local optimal solution S_{Local Optimal} . Now, if the optimal solution S_{Local Optimal} is better than the current solution S_{0}, the solution will be substituted and we get back to the first neighborhood structure k = 1. Otherwise, we move to the next neighborhood structure k = k+1 and the steps are iterated. In fact, searching for better optimal solution is continued. In VNSSA algorithm, the neighborhood search is executed by using of SA algorithm. In a certain temperature, the cycle of VNSSA algorithm continues until we have k ≤ k_{max}. Then, the algorithm temperature should be updated through T = α×T. The VNSSA algorithm continues until reaching the last neighborhood structure and final temperature. In the end, when termination criterion is satisfied the best solution will be represented.
To solve the problem by VNSSA algorithm, five neighborhood structures are used. The structures that are used in this paper are among the most common neighborhood structures for solving the VRPTW problems. These five neighborhood structures include the following operators:Figure 2

The first operator is insertion. At first, the two elements i and j are randomly selected. Then, the operator removes the i^{th} element and places it in position of element j+1 (one cell after element j).

The second operator is reversion. In this operator, two elements i and j are randomly selected. Then, the arrangement of the elements i, j and all the elements in between are reversed.

The third operator is swap. In this operator, the two elements i and j are randomly chosen. Then, the elements i and j replace their positions.

The fourth operator is route swap in which instead of two elements i and j, two rows of elements are randomly selected. Then, the two rows exchange their positions.

Finally, the fifth operator is twice a route swap. In this operator, four rows of elements are randomly selected instead of two rows. Then, first and fourth elements and second and third elements will exchange their position.
3.2.Bat Algorithm (BA)
The BA algorithm is a metaheuristic algorithm which was initially introduced by Yang (2008) and was extended later in 2010 (Yang, 2010). The studies by Yang and Gandomi (2012) and Tsai et al. (2012) on applications of BA algorithm showed that the algorithm is suitable for optimization problems and can obtain proper results.
The bats use sonic locating their hunt. The bats produce significantly high sonic pulse and listen to its reflection from objects. The delay in reflection of sound and variation in returned sound volume enable bats to detect the distance, direction and even speed of their hunt. The logic of BA algorithm is as described below:
Each virtual bat i flies with speed V_{i} in a random manner. The location of bat i is represented by x_{i}. During search for hunt, the bat changes its loudness (A_{i}) and pulse propagation rate (r_{i}). The search for finding hunt is randomly done and selection of the best hunt continues until one of the termination conditions is satisfied. In this algorithm, it is presumed that all bats use the ability of sonic locating to determine distance and to distinguish the difference between hunt and barriers. The loudness changes from a large, fixed and positive value (A_{0}) to a smaller one (A_{min}). The bats fly randomly throughout the location x_{i}, with the speed V_{i}, and fixed frequency of f_{min}, with variable wavelength of λ and loudness A_{0} to find their hunt.
In addition to these assumptions, other approximations are used in the design of BA algorithm. For example, the frequency f generally is in the range $\left[{f}_{\mathrm{min}},\hspace{0.17em}{f}_{\mathrm{max}}\right]$. The range of wavelength λ should be selected in compatibility with dimensions of problem and then, it should be changed into a smaller size.Figure 3
The steps of the BA algorithm are:

1. Initialization of location x_{i}, speed V_{i}, loudness A_{i} and pulse propagation rate r_{i}.

2. Iteration of the following steps until satisfaction of termination condition:

2.1. generation new solution by adjusting loudness and updating speed and location

2.2. If a random number between 0 to 1 is larger than pulse rate r_{i}, then:

2.3. If a between 0 and 1 is less than loudness (A_{i}) and objective function of bat i is less than the best bat, then:

2.4. The bats are ranked and the best solutions are selected.

4.COMPUTATIONAL EXPERIMENTS
4.1.Numerical Example
In this section, the accuracy of the suggested mathematical model and the performance of the proposed metaheuristic methods will be evaluated. To do this, a set of examples, for small, medium and large dimensions of the problem is designed. In the following section, the ability of metaheuristic methods in attaining optimal solutions will be investigated through a comparison of the results of the metaheuristic methods with those of the exact method. In addition, in this section, the performance of the metaheuristic methods in solving problems will be compared with each other. To do the mentioned tests, we use the GAMS and MATLAB software. For running the software, a computer with Core i5 processor (2.53 GHz) and 4GB RAM will be used.
4.1.1.Generated Problem Instances
In this section, 20 numerical examples of small, medium and large scale are implemented as shown in Table 1. In Table 1, the first column represents the number of problem. The second to fifth columns represent the number of producers, number of distribution centers, number of clients and number of vehicles respectively. To generate the numerical examples, the uniform probability distribution is used. To do this, a lower and an upper bound are considered for uniform distribution of each parameter. The values of parameters are randomly determined. The limits of the uniform distribution that each parameter follows are shown in Table 2.
4.1.2.Computational Results and Findings
The results of solving the generated examples are presented in Table 3 and Table 4. To solve each numerical example, each individual algorithm and the solution method is run for ten times. The analysis of results and evaluation of performance of algorithms as well as solution methods are performed based on three performance measures. The measures included mean values of objective function, mean calculation duration and percentage of algorithm error. The error value of each algorithm is determined through the formula “(AnswerBest Answer)/ Best Answer).” In this expression, “Answer” is the value reported by the algorithm and “Best Answer” is the best value for the problem. Due to being NPHard, the duration of solving mathematical model increases exponentially as dimensions of the problem rises. For solving the problems, if duration of problem solution exceeds 144000 seconds, running of model will stop and the best obtained solution is considered as the solution for the problem. This procedure is done for the numerical examples 17 to 20 and GAMS Software is incapable of obtaining optimal solution in the due time (144,000 seconds).
As shown in Table 3 and Table 4, in small dimen sions of the problem (examples 1 to 10) GAMS software offers better results than other solution methods and mean error percentage from optimal solution is inconsiderable. By the increase in dimensions of the problem (examples 10 to 20), TS algorithm offers better results than other solution methods and its mean error percentage from optimal solution will be 1.99 percent which signifies high efficiency of the algorithm. Figure 4 and Figure 5 show the values of objective function of numerical examples.
The results of numerical examples and Figure 4 and Figure 5 show that with the increase in dimensions of the problem, the objective function of Tabu search algorithm offers better results compared with GAMS Software and other metaheuristic algorithms (SA, BA, and VNSSA). In Figure 6, the processing times of solving numerical examples are represented.
By increasing the dimensions of numerical examples, the space of searching for optimal solution increases for all solution methods. This leads to a situation in which by increasing the dimension of the problems, the processing time of solving numerical examples will increase. The analysis of processing time of numerical examples in Table 3 and Table 4 as well as the Figure 4 shows that with increasing in problem dimensions, the processing time by GAMS Software will increase expo nentially, so that GAMS Software cannot obtain the optimal solution for numerical examples 17 to 20 in maximum permissible duration (144,000 seconds). That is why the best solutions obtained for the problems 17 to 20 up to 144,000 seconds are defined as the optimal solutions. The TS algorithm, as shown in Figure 4, has smaller slope and lower incremental trend. However, the slope of TS algorithm is more than other metaheuristic algorithms (SA, BA and VNSSA). It means that TS algorithm has more incremental trend than the other methods. The incremental trend for solving SA, BA and VNSSA algorithms declines respectively and in this regard, SA algorithm has the shortest processing time compared with other solution methods.
Figure 7 to Figure 10 show that SA, BA and VNSSA algorithms are improving throughout different iterations of the algorithm to get the optimal solution. In the Figures, the horizontal axis of the diagram represents the iterations of algorithm and vertical axis of diagram represents the best solution up to a certain iteration. The review of these diagrams let one to infer that metaheuristic algorithms which have been used in the present study are improving in consecutive iterations of the algorithm to get the optimal solution. Figure 7 to Figure 10 show the improvement trend of the algorithms for numerical example 20.Figure 8Figure 9
Based on the obtained results of solving the numerical examples, we can conclude that for solving the small dimensions of the problems, GAMS Software offers better optimal solutions in reasonable processing time but with the increase in the problem dimensions, its processing time increases exponentially and quality of its results will drop gradually. Compared with other metaheuristic algorithms, TS algorithm also offers more suitable results and maintains the quality of obtained results as the dimensions of the problem increase. But the increasing in dimensions of the problems adds to its processing time by incremental slope. In comparison with other meta heuristic algorithms, SA algorithm offers less desirable results but with the increase in dimensions of the problem, the processing time will be having less incremental slope. In fact, the algorithm offers more suitable processing time for solving larger dimensions of the problem compared with other solution methods. The BA algorithm does not offer desirable results for small dimensions of the problems but with the increase in problem dimensions, the quality of its results will be more suitable than the SA and GAMS methods. Moreover, with the increase in problem dimensions, the processing time of BA algorithm has less incremental rate compared with the TS and GAMS, but will show higher incremental rate compared with SA and VNSSA. In comparison with TS and GAMS methods, the VNSSA algorithm offers less suitable results in small dimensions of the problem but with the increase in the dimensions of the problem the quality of results offered by the VNSSA algorithm becomes more suitable than GAMS method; though, the TS algorithm still offers better results. The SA algorithm has less processing time than VNSSA but as shown in Figure 6, the VNSSA has less incremental slope of processing time compared with GAMS, TS and BA methods.
We observed that with the increase in the dimensions of a problem, some methods such as TS offer better results than other methods meanwhile the processing time of TS has higher incremental rate than other methods. In addition, we observed that with the increase in dimensions of the problem, SA algorithm has more desirable processing time compared with other methods with more undesirable quality of results than other methods. As a result, we sought for a method that offers more acceptable quality of results with the increase in the problem dimensions meanwhile having shorter processing time in comparison with methods producing similar quality of results. By reviewing all findings, one could conclude that VNSSA algorithm is the one which offers suitable results in more reasonable processing time. In fact, VNSSA algorithm establishes a balance between quality of results and processing time.
5.CONCLUSION AND FURTHER SUGGESTIONS
In the present study, the proposed model intends to consider the problems of routing, location and supplier selection in an integrated way so that the vehicles᾽ fuel consumption rate (FCR) and their pollutants emission are optimized simultaneously. To apply the proposed model in realworld problems, the timewindows limitation was imposed for serving the clients. The objectives of the models are minimization of total costs of the three echelon supply chain (suppliers, distribution centers, and clients), optimization of pollutants emission and penalties of deviation from timewindows and maximization of the clients’ satisfaction. Due to the fact that we face a mixed integer nonlinear programming model, by increasing in the binary variables, the processing time of solving the problem will increase exponentially and hence solving such problems in realworld scale through exact methods is not justifiable. Because the proposed model is NPhard, computing time of exact methods for solving the large dimensions of the problem is not feasible and as a result SA, TS, BA and VNSSA algorithms were used for solving the larger dimensions of the problem. For the performance evaluation of the algorithms, some numerical examples in different dimensions (small, medium and large) were generated and results of algorithms were compared with each other. The objective of comparing the results is to determine a method which offers acceptable results as the dimensions of problem get larger meanwhile taking less time to solve the problem compared with methods with similar quality of results. Finally, the survey of all methods led to the conclusion that VNSSA algorithm is the one which offers better quality of results in reasonable processing time. In fact, VNSSA algorithm establishes a balance between quality of the results and the required processing time. In real world, the present study is highly practical and applicable for many organizations. Reducing the costs of supply chain and fuel consumption rate of vehicles leads to reducing total costs of transportation and damage to the environment. Meanwhile, decreasing the travelling distance by a vehicle, the number of accidents will decrease. In future, potential innovations and new ideas will complement this research field to make it more compatible with the realworld problems. Some of these ideas are mentioned as follows. It is suggested to develop a model with multiple products with better compatibility with the realworld. In addition, for making the problem more real and getting closer to optimal solution, the problem parameters could be presumed to be uncertain and the model should be considered in uncertainty space. In addition, VRPs have numerous economic, social, and environmental aspects. Therefore, modeling the VRPs by considering different aspects simultaneously will contribute in getting problems closer to real world.