FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 18, N o 1, 2020, pp. 107 - 120 https://doi.org/10.22190/FUME190318021M © 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper  THE VEHICLE ROUTING PROBLEM WITH STOCHASTIC DEMANDS IN AN URBAN AREA – A CASE STUDY Danijel Marković, Goran Petrović, Žarko Ćojbašić, Aleksandar Stanković Faculty of Mechanical Engineering, University of Niš, Serbia Abstract. The vehicle routing problem with stochastic demands (VRPSD) is a combinatorial optimization problem. The VRPSD looks for vehicle routes to connect all customers with a depot, so that the total distance is minimized, each customer visited once by one vehicle, every route starts and ends at a depot, and the travelled distance and capacity of each vehicle are less than or equal to the given maximum value. Contrary to the classical VRP, in the VRPSD the demand in a node is known only after a vehicle arrives at the very node. This means that the vehicle routes are designed in uncertain conditions. This paper presents a heuristic and meta-heuristic approach for solving the VRPSD and discusses the real problem of municipal waste collection in the City of Niš. Key Words: Vehicle Routing Problem, Stochastic Demands, Municipal Waste, Heuristic and Meta-heuristic 1. INTRODUCTION Waste and mail collection, delivery of newspapers, milk, bread and postal packages, distribution of buses and planes in their networks of lines, all of these represent problems that researchers and traffic experts face daily in practice. Inadequate collection and transport, as functions of municipal waste management, result in enormous economic and environmental losses, and they serve as a major incentive to many researchers in their quest for appropriate systemic solutions. The set of transport means supporting the municipal waste collection and transport process most often comprises the majority of the vehicle fleet of public utility companies, which usually amounts to 50–70% of all transport units. This process is also dominant in those business systems that integrate a number of public utilities. In such systems, 15–40% of transport units support the municipal waste collection and transport process. Optimizing and applying a combination of heuristic and meta- Received March 18, 2019 / Accepted May 28, 2019 Corresponding author: Danijel Marković Faculty of Mechanical Engineering, University of Niš, Aleksandra Medvedeva 14, 18000 Niš, Serbia E-mail: danijel.markovic@masfak.ni.ac.rs 108 D. MARKOVIĆ, G. PETROVIĆ, Ţ. ĆOJBAŠIĆ, A. STANKOVIĆ heuristic methods in only one such system can lead to the reduction in mechanization fuels of 10 ÷ 25% [1]. Collecting municipal waste in urban areas can be observed as solving a vehicle routing problem. Such a problem is known in the literature as the waste collection vehicle routing problem. The vehicle routing problem (VRP) is one of the most challenging problems of combinatorial optimization. It was first mentioned in 1959 by Dantzig and Ramser [2]. Since then, the VRP has been increasingly applied in the solution of various problems and it bears a great economic importance for the reduction in operating costs of distribution systems [3, 4]. With the aim of satisfying real problems for the VRP solution, several constraints are commonly introduced in the solution of the problem, such as a larger number of depots, different types of vehicles (homogeneous and heterogeneous), different types of customer demands (deterministic and stochastic), infrastructural limitations (one-way streets, prohibited roads), types of service performance (pickup, delivery and mixed), etc. If all these constraints are taken into consideration, the VRP becomes much more complicated to solve and it falls under the NP difficult problems category. Because of this, different variants of the VRP problem have been introduced in the literature. The basic model of the vehicle routing problem is the capacitated vehicle routing problem – CVRP [5]. In the CVRP customer demands are deterministic, known in advance and cannot be separated. Vehicles are identical and they have a common starting point, while the only constraint is the vehicle capacity. The goal function expresses the demand to minimize total costs. Various vehicle routing problems are derived from the CVRP such as: Fig. 1 Basic vehicle routing problems and their relations For the sake of a simplified presentation, Fig. 1 uses abbreviations. The detailed explanation of Fig. 1 is as follows: distance–constrained capacitated vehicle routing problem – DCCVPR [6]; vehicle routing problem with backhauls – VRPB [7]; vehicle routing problem with time window – VRPTW [8], vehicle routing problem with pickup and delivery – VRPPD [9]; vehicle routing problem with backhauls and time window – VRPBTW [10]; vehicle routing problem with pickup and delivery and time window – VRPPDTW [11] In the above vehicle routing problems, the demands of the given transport network are deterministic, i.e. known in advance. However, in certain cases these demands can be random variables, that is, they become a stochastic quantity, which results in the standard The Vehicle Routing Problem with Stochastic Demands in an Urban Area – a Case Study 109 CVRP being expanded into a capacitated vehicle routing problem with stochastic demands – CVRPSD [12]. In the past few years, numerous researchers have applied heuristic and meta-heuristic methods to solve the CVRPSD [13, 14, 15, 16]. Examples of stochastic demands can be found in various transport activities. One such example is the municipal waste collection considered in this paper. Namely, here we discuss a capacitated vehicle routing problem with stochastic demands – CVRPSD, for municipal waste collection in urban areas. 2. MATHEMATICAL FORMULATION OF THE CVRPSD In practice, municipal waste collection in urban areas can be observed as a capacitated vehicle routing problem with stochastic demands – CVRPSD. This means that the amount of waste in nodes for the given transport network is a randomly variable quantity. The amount of waste may vary depending on the season and it can be known only after a vehicle arrives at a certain node to be served. For this reason, it is very hard to design the routes in the classical manner. Solving the CVRPSD for municipal waste collection (CVRPSD-MWC) thus encompasses finding their routes for the given transport network with minimal costs while meeting the following conditions:  there is only one depot and each route begins and ends in that depot,  the locations of the depot and the nodes are known,  the amount of waste in each node is a stochastic variable with normal distribution,  the capacity of the waste collection vehicle is known,  the sum of the amounts of municipal waste in a single route must not be greater than the vehicle capacity,  each node must be visited only once. 2.1 The chance-constrained model of the CVRPSD The Chance-Constrained (Ch-C) method is one of the main methods for stochastic optimization under various conditions of uncertainty. The Ch-C method allows for the probability of meeting a certain constraint to be above a necessary level. In other words, this method limits the allowable region thus reaching a high level of solution reliability. In recent times, numerous researchers have applied the Ch-C method to solving certain variants of the VRP [17, 18, 19, 20]. In this paper the Ch-C model for the CVRPSD observes one level of constraint, and that is the amount of waste. Namely, it is expected that the amount of waste per node is smaller than the vehicle capacity, with the level of reliability α. Using the Ch-C method to solve the CVRPSD, the goal function can be defined as follows:      n i n j ijij xcF 0 0 min (1) with the constraints:    n i ij njx 0 ,...,2,1,1 (2) 110 D. MARKOVIĆ, G. PETROVIĆ, Ţ. ĆOJBAŠIĆ, A. STANKOVIĆ      n i n i jiij njxx 0 0 ,...,1,0,0 (3)      n j n j jj zxx 1 0 1 00 2 (4)      Vi Vj iji QxqP )( (5) 0 ≤ qij ≤ Q, i = 0,1,..., n; j = 0,1,..., n (6) zi{0, 1}, i = 0,1,..., n (7) xij{0, 1}, i= 0,1,..., n; j = 0,1,..., n (8) The minimization function for the CVRPSD is shown in Eq. (1), where cij is the transport costs of the vehicle between node i and node j; i, jV it is assumed that cij=dij The constraint given in Eq. (2) shows that each node must be visited only once, while the constraint given in Eq. (3) presents the continuation of the vehicle flow, i.e. the fact that after serving node j the vehicle must leave that same node. The constraint given in Eq. (4) shows that each vehicle must start its route in the depot and end it there. Using the Ch-C conditions one can assure that the amount of collected waste on the route is smaller than the vehicle capacity with known probability (P) as shown in constraint (5). The capacity constraint (6) shows that the vehicle load never exceeds the vehicle capacity. Furthermore, Q is the maximum capacity of the vehicle, while qij is the capacity of the vehicle after visiting node i, and before visiting node j. Constraints (7) and (8) define the intervals of variables zi and xij. It follows that: 0, Vji therwiseo0, j node the visits and inode the after vehicle the If1, xij     mkVi therwiseo inodethevisitsvehicletheIf zi     , ,0 ,1 0 Constraint (5) can be solved by applying the Ch-C conditions. It is assumed that the amount of waste per each node is a random variable with normal distribution, which can be presented as: qi ~ N(μi , σi 2 ) (9) where μi is the total expected amount of waste for the i-th node, σi 2 is the standard deviation (variance) from the amount of waste for the i-th node. Parameters μi and σi 2 can be written using Eqs. (10) and (11):      n i iji n j i xqE 0 0 )]([ (10)      n i iji n j i xqVar 0 0 2 )]([ (11) where E(qi xij) is the mathematical expectation of normal distribution, while Var (qi xij) is the variance, i.e. the normal distribution scaling parameter. The Vehicle Routing Problem with Stochastic Demands in an Urban Area – a Case Study 111 If the expected customer demand is presented in the following way: QxqExqE n i iji n j iji     0 0 )]([)( (12) and the standard deviation as:      n i iji n j iji xqVarxqVar 0 0 )]([)( (13) using Eqs. (12) and (13) one can rework the Ch-C condition with constraint (5) into Eq. (14) [21].                             n i iji n j n i iji n j xqVar QxqE P 0 0 0 0 )]([ )]([ (14) It is important to emphasize that Eq. (14) holds if and only if Eq. (15) holds as well: 0 01 0 0 [ ( )] ( ) [ ( )] n n i ij i j n n i ij i j E q x Q Var q x            (15) Eq. (15) can be written as a deterministic equivalent: 1 0 0 0 0 ( ) [ ( )] [ ( )] n n n n i ij i ij i j i j Var q x E q x Q          (16) where Φ is the standard function of normal distribution, while Φ -1 is the inverse function of function Φ . As the amount of waste is assumed to be a random variable with normal distribution, this means that the waste collection vehicle routes should be designed under the conditions of uncertainty regarding the amount of waste, i.e. the demand value in nodes. The next section of the paper presents the results of the routing optimization for the observed problem. Parameter α takes the value of 0.8 for the optimization of the case study routes [22]. 3. DEFINING THE MODEL AND METHOD FOR THE CVRPSD The CVRPSD model considered in this paper is defined by a transport network that comprises one depot and 29 nodes (Fig. 2). The transport network presents “area” 103 according to the division of the territory of the City of Niš by the PUC “Mediana-Niš”. The nodes in the transport network present the locations of the containers as defined by the coordinates, i.e. the latitude and longitude (Tab. 1). Apart from the coordinates, Tab. 2 provides the number of containers per each node for the observed transport network. In the transport network, the first and the final node (the depot) is marked with “1”. The other nodes of the transport network are numbered from 2 to 30. 112 D. MARKOVIĆ, G. PETROVIĆ, Ţ. ĆOJBAŠIĆ, A. STANKOVIĆ Fig. 2 The transport network of container nodes for municipal waste collection Table 1 Coordinates and number of containers for the observed transport network Location number Latitude Longitude Number of containers per location Depot 43.319256 21.919682 2 43.322794 21.913082 4 3 43.322464 21.914317 4 4 43.324196 21.914412 2 5 43.324696 21.916001 2 6 43.323830 21.916709 3 7 43.323338 21.917632 5 8 43.322615 21.918220 1 9 43.323829 21.921069 2 10 43.322712 21.920759 2 11 43.322109 21.922545 1 12 43.321878 21.921097 1 13 43.321168 21.917535 4 14 43.322144 21.918179 2 15 43.322471 21.915755 3 16 43.319104 21.920206 1 17 43.319854 21.920548 2 18 43.320413 21.921150 2 19 43.320076 21.922198 3 20 43.320751 21.921739 2 21 43.320552 21.923959 2 22 43.322529 21.910364 1 23 43.322181 21.911683 1 24 43.321311 21.913388 1 25 43.321119 21.913767 1 26 43.320082 21.916993 1 27 43.319579 21.917587 1 28 43.319241 21.918091 1 29 43.319097 21.918627 1 30 43.318779 21.919317 1 The Vehicle Routing Problem with Stochastic Demands in an Urban Area – a Case Study 113 To solve the model, waste containers for municipal waste collection, so-called semi- underground containers, with the capacity of 3m 3 are installed in the transport network. This model does not consider the optimal locations and the number of containers. The number and location of containers are selected on the basis of the previous positions of waste containers determined by the PUC “Mediana-Niš”. On the locations where there are two or more containers, their positions are defined by a single node, i.e. a single coordinate. The waste collection vehicle departs the depot and it is always assumed empty. Waste collection should be performed with only one vehicle. This vehicle possesses a superstructure with a telescopic crane adapted to semi-underground waste containers. It is predicted that the waste collection in “area” 103 takes place three times a week (on Mondays, Wednesdays and Fridays). The matrix of shortest distances is symmetrical, and its elements represent the shortest possible real distance between the node pairs for the transport network. The amount of municipal waste in each node is stochastic, i.e. randomly variable. For the optimization purposes, the municipal waste collection route was monitored at specific time intervals. The monitoring was performed in ten instances for each node, during different seasons. The assessment of how full the containers were was made at each transport network node and this was recorded in a Tab. 2. This table was filled on the basis of the routing card. Table 2 Intervals of monitoring the amount of waste assessment per transport network node Node Interval of the amount of waste assessment 1 2 3 4 5 6 7 8 9 10 1 12.8 13.6 12.5 13.8 13.0 12.6 14.2 13.9 14.4 13.3 2 13.3 14.4 13.9 12.6 14.2 13.0 13.8 12.5 13.6 12.8 3 6.9 6.2 6.8 6.4 6.6 7.2 7.0 7.1 6.3 6.5 4 6.4 6.8 6.2 6.9 6.5 6.3 7.1 7.0 7.2 6.6 5 10.3 9.4 10.2 9.6 10.0 10.8 10.4 10.7 9.5 9.7 6 16.6 18.0 17.4 15.8 17.8 16.2 17.2 15.6 17.0 16.0 7 3.2 3.4 3.1 3.4 3.2 3.2 3.6 3.5 3.6 3.3 8 6.9 6.2 6.8 6.4 6.6 7.2 7.0 7.1 6.3 6.5 9 6.6 7.2 7.0 6.3 7.1 6.5 6.9 6.2 6.8 6.4 10 3.1 3.4 3.6 3.7 3.1 3.0 3.4 3.4 3.5 3.2 11 3.4 3.1 3.4 3.2 3.3 3.6 3.5 3.6 3.2 3.2 12 13.3 14.4 13.9 12.6 14.2 13.0 13.8 12.5 13.6 12.8 13 6.2 6.7 7.2 7.4 6.2 6.0 6.8 6.9 7.0 6.3 14 9.6 10.2 9.4 10.3 9.7 9.5 10.7 10.4 10.8 10.0 15 3.3 3.6 3.5 3.2 3.6 3.2 3.4 3.1 3.4 3.2 16 6.9 6.2 6.8 6.4 6.6 7.2 7.0 7.1 6.3 6.5 17 6.4 6.8 6.2 6.9 6.5 6.3 7.1 7.0 7.2 6.6 18 10.3 9.4 10.2 9.6 10.0 10.8 10.4 10.7 9.5 9.7 19 6.6 7.2 7.0 6.3 7.1 6.5 6.9 6.2 6.8 6.4 20 6.4 6.8 6.2 6.9 6.5 6.3 7.1 7.0 7.2 6.6 21 3.3 3.6 3.5 3.2 3.6 3.2 3.4 3.1 3.4 3.2 22 3.1 3.4 3.6 3.7 3.1 3.0 3.4 3.4 3.5 3.2 23 3.4 3.4 3.1 3.4 3.2 3.2 3.6 3.5 3.6 3.3 24 3.3 3.6 3.5 3.2 3.6 3.2 3.4 3.1 3.4 3.2 25 3.4 3.1 3.4 3.2 3.3 3.6 3.5 3.6 3.2 3.2 26 3.2 3.4 3.1 3.4 3.2 3.2 3.6 3.5 3.6 3.3 27 3.4 3.1 3.4 3.2 3.3 3.6 3.5 3.6 3.2 3.2 28 3.3 3.6 3.5 3.2 3.6 3.2 3.4 3.1 3.4 3.2 29 3.4 3.1 3.5 3.4 3.2 3.2 3.6 3.5 3.6 3.3 30 12.8 13.6 12.5 13.8 13.0 12.6 14.2 13.9 14.4 13.3 114 D. MARKOVIĆ, G. PETROVIĆ, Ţ. ĆOJBAŠIĆ, A. STANKOVIĆ The routing card contains the following information: the numerical marker of the container node; the name of the place where the containers are located; the time of arrival, pickup and departure; the number of containers per node and the assessment of how full the container is; the accessibility of the container; the note. 3.1. Stochastic optimization problems In stochastic optimization, optimization parameters have a random character described by the methods from the theory of probability and statistics. When one observes the constraint functions whose parameters have a random character, it is not certain that a selection of control variables will ensure that the function is satisfied. Therefore, a new optimization task is defined to demand that the probability of a constraint being met should be greater than a predetermined value. If the parameters in constraint functions have a random character (random quantities), then they can be described using, for example, normal distribution, Poisson distribution or Gamma distribution. 3.2. Heuristics and meta-heuristics for the CVRPSD When solving a vehicle routing problem, heuristic methods are used to construct routes, with the construction and improvement of routes being performed iteratively in relation to the goal function. Bearing in mind that this paper deals with the determination of the most favorable (optimal) routes of municipal waste collection vehicles, as well as considering all the observed constraints, it can be safely assumed that the sufficiently good solutions were determined by applying the global optimization methods. On the basis of research, the Clarke and Wright’s savings algorithm is considered in this paper as the representative of the approaches of the constructive heuristic methods [23]. Meta-heuristics is conceived as a means of solution of complex optimization problems where other optimization methods cannot provide an efficient and economical solution to the optimization problem at hand. One of the basic representatives of meta-heuristics is the genetic algorithm [24]. Today these methods are regarded as belonging to the most practical approaches to solving various complex problems [25], which is particularly related to the solution of numerous real problems that are combinatorial in nature, such as the vehicle routing problem itself. It can generally be said that meta-heuristics is a higher level of heuristics. Out of the group of meta-heuristic methods used for the solution of municipal waste collection vehicle routing problems, this paper employs the 2-OPT local search [26] and Simulated Annealing – SA [27]. 3.3. Stochastic simulation for computing the expected value and probability check The first step in the optimization of the CVRPSD is to compute the expected value of the amount of municipal waste (μi) and check probability (β). This step is necessary due to the stochastic character of the amount of municipal waste per transport network node. Based on the input data on the assessed amount of waste per node, one can compute the expected values of the amount of municipal waste (μi) for each node of the transport network since the distribution is known, i.e. normal distribution. After the expected value of the amount of municipal waste is computed, variance (σi 2 ) gets computed as well. To compute the mathematical expectancy and variance Procedure 1 was used, and its pseudo code is shown in Algorithm 1. The Vehicle Routing Problem with Stochastic Demands in an Urban Area – a Case Study 115 Algorithm 1: Procedure 1 Start Define the assessed amount of waste per transport network node; Define Q; for n = 1; n ≤ nuk; n = n + 1; sumn = 0; for i =1; i ≤ 10; i = i + 1; sumn = sumn + qni; pi = qni / sumn: end for end for for n = 1; n ≤ nuk; n = n + 1; compute En; compute σn; compute Φ(α); if   n n QE     probability condition = TRUE; else probability condition = FALSE; end if end for end Here, sum is the amount of waste per transport network node, qni is the assessed amount of waste in the node. When these two parameters are computed, then probability (β) is checked. Algorithm 1 presents the procedures to check the probability. The last step in this procedure is the checking of the Ch-C condition, and if this condition is met, the procedure continues (TRUE). In the opposite case the procedure is stopped (FALSE). 3.4. Initial solution The next step in the CVRPSD optimization is the formation of the initial solution. Bearing in mind that this is a stochastic problem, in line with the previous explanation, the problem is reduced to the solution of the deterministic problem by applying Eqs. (15) and (16). The C-W savings algorithm was used to obtain the initial solution. In the application of the C-W savings algorithm parameter qi was substituted with parameter μi. Procedure 2 presents the pseudo code for the C-W savings algorithm (Algorithm 2). Algorithm 2: Procedure 2 start Define distance matrix; Define Q; Call Procedure 1; Compute s; Sort s' in a non-increasing sequence; Form a partial route; Expected demand = µi; for all savings from sequence if (probability condition == TRUE) if met operative constraints 116 D. MARKOVIĆ, G. PETROVIĆ, Ţ. ĆOJBAŠIĆ, A. STANKOVIĆ if Expected demand + µi ≤ Q Expected demand = µi + Expected demand Form route end if end if end if end for Vehicle fullness = Expected demand; Print routes; Print vehicle fullness; end 3.5. 2-OPT search and SA meta-heuristic for the CVRPSD The first improvement of the initial CVRPSD solution was performed by applying the 2-OPT local search. During the improvement of the initial solution, the number of iterations was varied (1e3 and 1e6). The initial solution was improved by applying the 2- OPT local search algorithm. The pseudo code of the 2-OPT algorithm for the improvement of the initial solution is presented in Algorithm 3. Algorithm 3: 2-OPT algorithm for the improvement of the initial solution start Load initial route; U0 = initial route length; for (i = 1; i ≤ n - 2; i = i + 1) for (j = i + 2; j ≤ n; j = j + 1) U' = d(i, j) + d(i+1, j+1) - d(i, i+1)- d(j, j+1); if (U'< U) U' = U; end if end for end for end The next algorithm used to optimize the CVRPSD was the SA algorithm. The solution obtained by applying the SA algorithm largely depends on adjusting the parameters of the algorithm itself. However, this paper does not consider the selection of optimal parameters for the given problem but uses the recommended parameters. The parameters of the SA algorithm used to solve the CVRPSD are [27]: initial temperature T0 = 100, temperature reduction factor α = 0.8. The pseudo code of the SA algorithm for the solution of the CVRPSD is presented by Algorithm 4. Algorithm 4: SA algorithm for solving CVRPSD start Define model; Load initial solution U obtained by C-W savings algorithm; Define SA algorithm parameters; Best solution = U; T = T0; The Vehicle Routing Problem with Stochastic Demands in an Urban Area – a Case Study 117 for it1 = 1; it1