Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.16 No.4 pp.574-589
DOI : https://doi.org/10.7232/iems.2017.16.4.574

Incorporating Vehicle Routing, Location and Supplier Selection Problems for Reducing Pollutants Emission

Ebrahim Teimoury, Seyed Omid Hashemi Amiri*, Farbod Ketabchi
Department of Industrial Engineering, Iran University of Science and Technology, Tehran, Iran
Management of Technology (MOT), Iran University of Science and Technology, Tehran, Iran
Corresponding Author, Omid_hashemi@alumni.iust.ac.ir
20170926 20171001

ABSTRACT

Vehicle routing, location and supplier selection problems are the most practical and challenging problems of the supply chain management. The particular attention to these problems is due to their high practicality in real world and the difficulty of solving. Therefore, the present paper attempts to offer a combined model of vehicle routing, location and supplier selection problems with more compatibility with the real-world problems. Moreover, the proposed model will attempt to reduce the fuel consumption rate (FCR) of vehicles and pollutant emissions simultaneously. Because the problem is NP-hard, meta-heuristic algorithms including simulated annealing (SA), Tabu search algorithm (TS), bat algorithm (BA) and variable neighborhood search algorithm based on simulated annealing (VNS-SA) are used for solving the large dimensions of the problem. For smaller dimensions of the problem, several numerical examples are generated and solved by GAMS Software and meta-heuristic algorithms of SA, TS, BA, and VNS-SA. Then, the results are compared for verification of the solution methods’ efficiency. In the rest of the paper, a real example is offered as an example of real-world problems and solved by metaheuristic algorithms. More clearly, the core problem consists of assigning producers to distribution centers, vehicle routing and service of distribution centers to clients in a supply chain of natural honey.


초록


    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 on-road 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 time-consuming. 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 meta-heuristic methods. In recent years, in addition to the heuristic methods for solving the VRP, meta-heuristic 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 Simchi-Levi (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ínez-Salazar et al. (2014) studied the bi-objective transportation location routing problem for determining the vehicles’ optimal load and minimization of transportation costs; they used meta-heuristic 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 four-echelon capacitated LRP.

    Because of the NP-hard nature of LRP problems, the application of heuristic and meta-heuristic algorithms has drawn the attention of many researchers. Tuzun and Burke (1999) used two-phase 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 real-world 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 NP-hard, the meta-heuristic algorithms of SA, TS, BA, and VNS-SA are used in this study to solve the proposed problem. The meta-heuristic 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 meta-heuristic 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 three-echelon 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 time-windows 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 (iI) ; i index of potential producers

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

      • K : set of customers (kK) ; k index of customers

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

      • E, F : set of producers and distribution centers (E, F = JK) ; e, f index of producers and distribution centers

    • (2) Parameters

      • Dj : total demands for distribution center j, for all jJ

      • Dk : total demands for customer k, for all kK

      • dij : distance from producer i to distribution center j (Kilometer), for all i∈ I, jJ

      • def : distance from point e to point f (Kilometer), for all jJ, jJ

      • CapVM : capacity of vehicles belonging to producers

      • C a p V v D : capacity of vehicle v belonging to distribution centers, for all vV

      • C a pi : production capacity at producer i, for all iI

      • C a pj : supply capacity at distribution center j, for all jJ

      • tij : mean duration of moving from producer i to distribution center j for all iI, jJ

      • tef : mean duration of moving from point e to point f, for all eE, fF

      • t j * : mean time of unloading in distribution center j for all jJ.

      • τj: maximum time allowed for procurement demand of distribution center j, for all jJ

      • Sk : duration of serving the customer k, for all kK

      • [ L k , U k ] : range of time window to serve the customer k, for all kK

      • TCij : transportation cost of each unit of product, from producer i to distribution center j for all iI, jJ

      • TCefv : transportation cost of each unit of product, from point e to point f by vehicle v for all eE, fF

      • FCRij : fuel consumption rate form producer i to distribution center j (Lit/km) for all iI, jJ

      • FCRe f v : fuel consumption rate from point e to point f by vehicle v (Lit/km) for all eE, fF, vV

      • 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

      • C0 : the cost of each unit of vehicle fuel (liter)

      • Ho j : the cost of inventory holding for distribution center j, for all jJ

      • Λj : the cost of commodity shortcoming for distribution center j, for all jJ

      • N: a large number

    • (3) Decision variables

      • γij : amount of flow from producer i to distribution center j, for all iI, jJ

      • γefv : amount of flow from point e to point f by vehicle v, for all eE, fF, kK

      • Inj : inventory of distribution center j at the end of the course, for all jJ

      • Brj : shortage of distribution center j at the end of the course, for all jJ

      • Bk ,v : arrival time of vehicle v to customer k, for all kK, vV

      • ZWk ,v : waiting time of vehicle v for serve to the customer k, for all kK, vV

        Z W k , v = { ( L k B k , v ) L k B k , v 0 L k B k , v k K , v V D

      • ZGk,v : the amount of violation from time window of vehicle v to serve to the customer k, for all kK, vV

        Z G k , v = { ( B k , v + S k U k ) U k B k , v + S k 0 U k B k , v + S k k K , v V D

      • Carjv : 1: if vehicle v allocated to distribution center j, for all jJ, vV

        0: otherwise;

      • Re, fv : 1: if point e will be before point f upon vehicle v, for all eE, fF, vV

        0: otherwise;

      • Yij : 1: if producer i allocated to distribution center j, for all iI, jJ

        0: otherwise;

      • Yjk : 1: if customer k allocated to distribution center j, for all jJ, kK

        0: otherwise;

    2.2.The Mathematical Model Formulation

    The proposed model is formulated as the following mixed integer nonlinear programming (MINLP) model:

    Objective function

    M i n i m i z e Z = i I j J T C i j γ i j + j J k K v V D T C e f v γ e f v + e K f K v V D T C e f v γ e f v + i I j J C 0 . F C R i j . d i j . Y i j + j J k K v V D C 0 . F R C j k v . d j k . R j , k v + e K f K v V D C 0 . F R C e f v . d e f . R e , f v + j J ( I n j . H o j + B r j . Λ j ) + k K v V D Z W k , v . C W + k K v V D Z G k , v . C G
    (1)

    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 time-windows (waiting costs or delay costs).

    Subject to(17)(18)

    γ i j min { C a p i , C a p V M } i I , j J
    (2)

    γ i j Y i j N i I , j J
    (3)

    j J γ i j C a p i i I
    (4)

    j J D j Y i j C a p i i I
    (5)

    i I γ i j = D j j J
    (6)

    k K D k Y j k C a p j j J
    (7)

    k K Y j k = 1 j J
    (8)

    j J Y j k 1 k K
    (9)

    γ j k v min { C a p j , C a p V v D } j J , k K , v V D
    (10)

    γ j k v Y j k N j J , k K , v V D
    (11)

    k K v V D γ j k v C a p j j J
    (12)

    j J v V D γ j k v = j J v V D D k k K
    (13)

    j J v V D R j , k v 1 k K
    (14)

    j J k K R j , k v 1 v V D
    (15)

    k K v V D R k , e v 1 e E
    (16)

    e E v V D R k , e v = 1 k K
    (17)

    k K v V D R e , k v 1 e E
    (18)

    e E v V D R e , k v = 1 k K
    (19)

    f F R e , f v = f F R f , e v e E , v V D
    (20)

    e K R f , e v + e K R e , k v Y j k 1 j J , k K , v V D
    (21)

    e E k K R e , k v D k C a p V v D v V D
    (22)

    e E v V D γ e , f v = e K v V D D f R e , f v f K
    (23)

    j J k K γ j k v + e K f K γ e , f , v C a p V v D v V D
    (24)

    C a r j v = k K R j , k v γ j J , v V
    (25)

    D j k K v V D γ j , k v e K f K v V D γ e , f v j J = I n j B r j
    (26)

    ( t i j + t j * ) . Y i j τ j i I , j J
    (27)

    j J t j k Y j k = j J v V D B k , v R j , k v k K
    (28)

    ( B f , v B e , v S e t e , f ) R e , f v = 0 e K , f K , v V D
    (29)

    Z W k , v = M a x ( ( L k B k , v ) , ) v V D , k K
    (30)

    Z G k , v = M a x ( ( B k , v + S k U k ) , 0 ) v V D , k K
    (31)

    R e , f v , R j , k v , Y i j , Y j k { 0 , 1 } i I , j J k K , v V D
    (32)

    γ i j , γ j k v γ e , f v B k , v Z W k , v , Z G k , v 0 i I , j J k K , v V D
    (33)

    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 kth client is not assigned to the distribution center j, there will be no flow between distribution center j and kth 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 kth client is equal to the demand by the kth 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 kth 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 kth client is assigned to distribution center j if and only if one route exists between the kth 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 kth client should be precisely equal to the demand by the kth 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 ith 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 kth 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 Re, f v, Rj, k v, Yij, Yjk are binary decision variables. Besides, they show that the decision variables B k , v , γ i j , γ j k , γ e , f v and Z W k , v , 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)

    F C R ( k N V + P / η ) τ / π 1 )
    (34)

    P = ( P T F / ε + P a )
    (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, Pa indicates necessary power of vehicle engine for discharging energy caused by accessory parts of vehicles (like air condition and cooler) and PTF indicates total force necessary for pulling a vehicle (kilowatt) and can be calculated as follows:(36)(37)(38)

    P T F = ( M a + M g sin θ + 0.5 C d A ρ v 2 + M g C r cos θ ) . v / 1 , 000
    (36)

    α = ( a + g sin θ + g C r cos θ )
    (37)

    β = ( 0.5 C d A ρ )
    (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 (meter2), ρ indicates air density (kilogram/meter3), Cd indicates coefficient of the aero- dynamic drag, Cr indicates coefficient of the rolling resistance, α indicates special fixed for each vehicle route, β indicates special fixed vehicle. Consequently PTF becomes as follows:(39)

    P T F = ( α M v + β v 3 ) / 1 , 000
    (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 = vij during distance dij 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)

    F C R ( v ) = ω ( k N V + M α v λ + β v 3 λ + P a / η ) 2 )
    (40)

    Then, final FCR function in unit of liter per second is as follows:(41)

    F C R ( v ) = ω ( k N V + M α v λ + β v 3 λ + P a / η ) / v = k n V ω / v + ( M α v λ + β v 3 λ + P a / η ) ω / v
    (41)

    Function of FCR per distance traveled is U-shaped; 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 ( M α v λ + β v 3 λ + P a / n ) ω / 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 ( M α v λ + β v 3 λ + P a / n ) ω / 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 ( M α v λ + β v 3 λ + P a / η ) ω / v components.

    As a result, the final formula of F C R ( v ) = k N V ω / v + ( M α v λ + β v 3 λ + P a / η ) ω / V through which the parameters FCRij and FCRefv 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 non-linear 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 non-linear expressions. There are two types of non-linear 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 = x1 × x2 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 non-linear limitations into a set of linear limitations.(42)(43)(44)

      Z X 1
      (42)

      Z X 2
      (43)

      Z X 1 + X 2 1
      (44)

    • (2) Linearization of the product of a binary and a continuous variable

      Let’s suppose Z = x1× x2 is multiplication of binary variable x1 by continuous variable x2 . 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 non-linear limitations into a set of linear limitations.(45)(46)(47)

      Z X 2
      (45)

      Z N × X 1
      (46)

      Z X 2 N ( 1 X 1
      (47)

      In the suggested model, objective function and limitations (28) and (29) include non-linear expressions. In the following, linearization of non-linear expressions is done.

    2.4.1.Linearization of the Objective Function

    The objective function has non-linear expressions. In these non-linear 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 non-linear sections of the expressions are placed in brackets. For linearization of non-linear expressions of objective function, three continuous variables F ij , F jkv , and F 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)

    (54),(57),(59)

    i I j J C 0 . d i j . { F C R i j . Y i j }
    (48)

    i I j J C 0 . d i j . Γ i j
    (49)

    Γ i j F C R i j i I , j J
    (50)

    Γ i j N Y i j i I , j J
    (51)

    Γ i j F C R i j N ( 1 Y i j ) i I , j J
    (52)

    j J k K v V D C 0 . d j k . { F C R j k v . R j , k v }
    (53)

    j J k K v V D C 0 . d j k . Γ j , k v
    (54)

    Γ j k v F C R j k v j J , k K , v V D
    (55)

    Γ j k v N . R j k v j J , k K , v V D
    (56)

    Γ j k v F C R j k v N . ( 1 R j , k v ) j J , k K , v V D
    (57)

    e K f K v V D C 0 . d e f . { F C R e f v . R e , f v }
    (58)

    e K f K v V D C 0 . d e f . Γ e f v
    (59)

    Γ e f v F C R e f v e K , f K , v V D
    (60)

    Γ e f v N . R e f v e K , f K , v V D Γ e f v F C R e f v N . ( 1 R e , f v )
    (61)

    e K , f K , v V D
    (62)

    2.4.2.Linearization of the non-linear constraints

    In the suggested model, the limitations (28) and (29) have non-linear expressions. In these non-linear expressions, a binary variable is multiplied by a continuous variable. The rewritten limitations are represented in equations (63) and (68) in which the non-linear parts are represented in brackets. To linearize the non-linear limitation (28), the continuous variable β1ef is defined while for non-linear limitation (29), two continuous variables β2ef and β3ef 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 Rj,kv could be 1 at most, we could consider j J R j , k v as a binary variable.(65)(66)(70)(71)(72)(73)(74)

    j J t j k Y j k = v V D { B k , v j J R j , k v } k K
    (63)

    j J t j k Y j k = v V D β 1 k , v k K
    (64)

    β 1 k v B k , v k K , v V D
    (65)

    β 1 k v N . j J R j , k v k K , v V D
    (66)

    β 1 k v B k , v N ( 1 j J R j , k v ) k K , v V D
    (67)

    { B f , v R e , f v } { B e , v R e , f v } = { S e + t e f } R e f , v e K , f K , v V D
    (68)

    β 2 f v β 3 e v = ( S e + t e , f ) R e , f v e K , f K , v V D
    (69)

    β 2 f v B f , v f K , v V D
    (70)

    β 2 f v N R e , f v e K , f K , v V D
    (71)

    β 2 f v B f , v N ( 1 R e , f v ) e K , f K , v V D
    (72)

    β 3 e v B e v e K , v V D
    (73)

    β 3 e v N R e , f v e K , f K , v V D
    (74)

    β 3 e v B e , v N ( 1 R e , f v ) e K , f K , v V D
    (75)

    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 NP-hard 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 meta-heuristic algorithms SA, TS, BA and VNS-SA 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 (VNS-SA)

    The main idea of the VNS-SA 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 S0 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 S0 to the neighbor solution S1 based on kth 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 S1 to local optimal solution SLocal Optimal . Now, if the optimal solution SLocal Optimal is better than the current solution S0, 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 VNS-SA algorithm, the neighborhood search is executed by using of SA algorithm. In a certain temperature, the cycle of VNSSA algorithm continues until we have kkmax. Then, the algorithm temperature should be updated through T = α×T. The VNS-SA 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 VNS-SA 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 ith 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 meta-heuristic 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 Vi in a random manner. The location of bat i is represented by xi. During search for hunt, the bat changes its loudness (Ai) and pulse propagation rate (ri). 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 (A0) to a smaller one (Amin). The bats fly randomly throughout the location xi, with the speed Vi, and fixed frequency of fmin, with variable wavelength of λ and loudness A0 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 [ f min , f max ] . 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 xi, speed Vi, loudness Ai and pulse propagation rate ri.

    • 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 ri, then:

        • 2.2.1. A solution is selected among the best solutions.

        • 2.2.2. Then, a solution is generated in neighborhood of the best selected solution (random step).

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

        • 2.3.1. A new solution is selected. Then, pulse propagation rate ri w increases and loudness Ai reduces.

      • 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 meta-heuristic 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 meta-heuristic methods in attaining optimal solutions will be investigated through a comparison of the results of the meta-heuristic methods with those of the exact method. In addition, in this section, the performance of the meta-heuristic 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 “(Answer-Best 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 NP-Hard, 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 meta-heuristic algorithms (SA, BA, and VNS-SA). 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 meta-heuristic algorithms (SA, BA and VNS-SA). It means that TS algorithm has more incremental trend than the other methods. The incremental trend for solving SA, BA and VNS-SA 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 VNS-SA. In comparison with TS and GAMS methods, the VNS-SA 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 VNS-SA algorithm becomes more suitable than GAMS method; though, the TS algorithm still offers better results. The SA algorithm has less processing time than VNS-SA but as shown in Figure 6, the VNS-SA 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, VNS-SA 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 real-world problems, the time-windows 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 time-windows and maximization of the clients’ satisfaction. Due to the fact that we face a mixed integer non-linear programming model, by increasing in the binary variables, the processing time of solving the problem will increase exponentially and hence solving such problems in real-world scale through exact methods is not justifiable. Because the proposed model is NP-hard, computing time of exact methods for solving the large dimensions of the problem is not feasible and as a result SA, TS, BA and VNS-SA 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 VNS-SA algorithm is the one which offers better quality of results in reasonable processing time. In fact, VNS-SA 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 real-world problems. Some of these ideas are mentioned as follows. It is suggested to develop a model with multiple products with better compatibility with the real-world. 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.

    Figure

    IEMS-16-574_F1.gif

    FCR based on function of speed.

    IEMS-16-574_F2.gif

    Pseudo code of VNS-SA Algorithm.

    IEMS-16-574_F3.gif

    Pseudo code of BA.

    IEMS-16-574_F4.gif

    Objective function for small dimensions of problem.

    IEMS-16-574_F5.gif

    Objective function for medium and large dimensions of problem.

    IEMS-16-574_F6.gif

    Solution time.

    IEMS-16-574_F7.gif

    Improvement trend for SA algorithm.

    IEMS-16-574_F8.gif

    Improvement trend for TS algorithm.

    IEMS-16-574_F9.gif

    Improvement trend for BA algorithm.

    IEMS-16-574_F10.gif

    Improvement trend for VNS-SA algorithm.

    Table

    Dimensions of numerical examples

    Lower and upper limits of numerical examples parameters

    Computational results for small, medium and large dimensions

    Computational results for small, medium and large dimensions

    REFERENCES

    1. BalakrishnanA. WardJ.E. WongR.T. (1987) Integrated facility location and vehicle routing models: Recent work and future prospects. , Am. J. Math. Manage. Sci., Vol.7 (1-2) ; pp.35-61
    2. BarmelJ. Simchi-LeviD. (1997) On the effectiveness of set covering formulation for the vehicle routing problem with time windows. , Oper. Res., Vol.45 (2) ; pp.295-301
    3. BarthM. BoriboonsomsinK. (2009) Energy and emissions impacts of a freeway-based dynamic ecodriving system. , Transp. Res. Part D Transp. Environ., Vol.14 (6) ; pp.400-410
    4. BarthM. YoungloveT. ScoraG. (2005) Development of a Heavy-Duty Diesel Modal Emissions and Fuel Consumption Model, California PATH Program, Institute of Transportation Studies, University of California at Berkeley,
    5. BektasT. LaporteG. (2011) The pollution-routing problem. , Transp. Res., Part B: Methodol., Vol.45 (8) ; pp.1232-1250
    6. BermanO. JailletP. Simchi-LeviD. , DreznerZ. (1995) Facility Location: A Survey of Applications and Methods., Springer, ; pp.427-452
    7. DantzigG. RamserJ. (1959) The truck dispatching problem. , Manage. Sci., Vol.6 (1) ; pp.80-91
    8. DantzigG.B. FulkersonD.R. (1954) Minimizing the number of tankers to meet a fixed schedule. , Naval Research Logistics Quarterly, Vol.1 (3) ; pp.217-222
    9. FisherM. (1981) The Lagrange relaxation method for solving integer programming problems. , Manage. Sci., Vol.27 (1) ; pp.1-18
    10. GhianiG. LaporteG. MusmannoR. (2004) Introduction to Logistics System Planning and Control., John Wiley & Sons Inc.,
    11. GloverF. WoolseyL. (1974) Converting the 0-1 polynomial programming problem to a 0-1 linear program. , Operation Research, Vol.22 (1) ; pp.180-182
    12. HamidiM. FarahmandK. SajjadiS. (2012) Modeling a four-layer location routing problem. , International Journal of Industrial Engineering Computations, Vol.3 (1) ; pp.43-52
    13. HamidiM. FarahmandK. SajjadiS. NygardK. (2014) A heuristic algorithm for a multi-product four-layer capacitated location-routing problem. , International Journal of Industrial Engineering Computations, Vol.5 (1) ; pp.87-100
    14. KuoY. WangC-C. (2012) A variable neighborhood search for the multi-depot vehicle routing problem with loading cost. , Expert Syst. Appl., Vol.39 (8) ; pp.6949-6954
    15. LaportG. NobertY. (1987) Exact algorithm for the vehicle routing problem. , surveys in Combination oriel Optimization, Amsterdam, Vol.132 ; pp.147-184
    16. LaporteG. (1989) A survey of algorithms for locationrouting problems. , InvestigaciA3n Operativa, Vol.1 ; pp.93-123
    17. LaportG. MercureH. NobertY. (1992) A Branch and Bound Algorithm for a Class of Asymmetrical Vehicle Routeing Problems. , J. Oper. Res. Soc., Vol.43 (5) ; pp.469-481
    18. LenstraJ.K. Rinnooy KanA.H. (1981) Complexity of vehicle routing and scheduling problems. , Networks, Vol.11 (2) ; pp.221-227
    19. MartA-nez-SalazarI.A. MolinamJ. A?ngel-BelloF. GA3mezT. CaballeroR. (2014) Solving a bi-objective transportation location routing problem by metaheuristic algorithms. , Eur. J. Oper. Res., Vol.234 (1) ; pp.25-36
    20. MinH. JayaramanV. SrivastavaR. (1998) Combined location-routing problems: A synthesis and future research directions. , Eur. J. Oper. Res., Vol.108 (1) ; pp.1-15
    21. NagyG. SalhiS. (2007) Location-routing: Issues, models and methods. , Eur. J. Oper. Res., Vol.177 (2) ; pp.649-672
    22. NorouziN. Tavakkoli-MoghaddamR. GhazanfariM. AlinaghianM. SalamatbakhshA. (2012) A new multiobjective competitive open vehicle routing problem solved by particle swarm optimization. , Netw. Spat. Econ., Vol.12 (4) ; pp.609-633
    23. ProdhonC. PrinsC. (2014) A survey of recent research on location-routing problems. , Eur. J. Oper. Res., Vol.238 (1) ; pp.1-17
    24. SalhiS. ImranA. WassanN. (2013) The multidepot vehicle routing problem with hetero-geneous vehicle fleet: Formulation and a variable neighborhoodsearch implementation. , Computers & Operations Research, 52(Part B), ; pp.315-325
    25. TsaiP.W. PanJ.S. LiaoB.Y. TsaiM.J. IstandaV. (2012) Bat algorithm inspired algorithm for solvIncorporating ing numerical optimization problems. , Appl. Mech. Mater., Vol.148-149 ; pp.134-137
    26. TuzunD. BurkeL.I. (1999) A two-phase tabu search approach to the location routing problem. , Eur. J. Oper. Res., Vol.116 (1) ; pp.87-99
    27. XuY. JiangW. (2014) An improved variable neighborhood search algorithm for multi-depot heterogeneous vehicle routing problem based on hybrid operators. , International Journal of Control and Automation, Vol.7 (3) ; pp.299-316
    28. XuY. WangL. YangY. (2012) A new variable neighborhood search algorithm for the multi-depot heterogeneous vehicle routing problem with time windows. , Electron. Notes Discrete Math., Vol.39 (1) ; pp.289-296
    29. YangX. (2008) Nature-inspired metaheuristic algorithms., Luniver Press,
    30. YangX. (2010) A new metaheuristic bat-inspired algorithm. , Nature Inspired Cooperative Strategies for Optimization, Vol.284 ; pp.65-74
    31. YangX. GandomiA. (2012) Bat algorithm: A novel approach for global engineering optimization. , Eng. Comput., Vol.29 (5) ; pp.464-483
    32. YuV.F. LinS-W. LeeW. TingC-J. (2010) A simulated annealing heuristic for the capacitated location- routing problem. , Comput. Ind. Eng., Vol.58 (2) ; pp.288-299
    오늘하루 팝업창 안보기 닫기