Route Optimization for Garbage Trucks in City Environments
submitted in partial fulfillment for the degree of master of science Martijn P.M. de Jong
10774807
master information studies data science
faculty of science university of amsterdam
2019-07-04
Internal Supervisor External Supervisor Title, Name Dr. Frank M. Nack Dr. Jeroen Silvis
Affiliation UvA, FNWI, IvI Municpality of Amsterdam Email nack@uva.nl j.silvis@amsterdam.nl
Route Optimization for Garbage Trucks in City Environments
Martijn P.M. de Jong
martijn.dejong@mail.com University of Amsterdam
ABSTRACT
Collecting waste of citizens is a costly and yet often an inefficient process. An increase in efficiency can lead to a great decrease in cost. In this thesis, it is studied to which extent automatic planning dynamic routes for collecting waste can increase efficiency based on historical weight data. The filling level for each site is estimated on a certain point in time based on the historical data. This estimate together with the driving time between sites formed the input for calculating routes for each working day. A linear- and quadratic programming model are implemented. The first model consists of two steps where first the containers to be emptied are determined based on filling level and second the routes are calculated by min-imizing time. The second model combines these two steps in a single objective function, where containers are determined based on filling level and travel time, which maximizes profit. Both models are applied in a case study for the city of Amsterdam. An offline simulation for the first quarter of 2019 is conducted and compared to the current planning. The second model performs better on both kilograms collected per minute and per kilometer driven than the first model. Both models perform better on kilograms per minute than the current planning.
KEYWORDS
Vehicle Routing Problem, Waste Collection, Optimization, Simplex
1
INTRODUCTION
Collecting waste of citizens is an important task for each govern-ment and thereby municipalities. In 1900 worldwide fewer than 300.000 tonnes of waste per day was created by urban residents. By 2000 this was more than three million tonnes and by 2025 the prediction is that it will be doubled. Enough to create a line of garbage trucks of 5.000 kilometers every day [7]. Different fractions of waste such as glass and paper are also separated for collection which results in an even more complex logistical operation. Govern-ments spend a significant part of their budget on collecting waste. According to a report of the European Commission, collecting resid-ual waste costs between 42 and 126 euro per tonne1. One way to reduce cost, and nuisance, is to generate less waste as a society. As mentioned, this is not in line with the predictions. A small reduc-tion of the cost per tonnage through optimizareduc-tion will lead to a significant decrease in overall costs. This is not only beneficial for the government as an institute, but also for the citizens since this concerns money of tax-payers.
Besides costs, garbage collection has also other effects on ur-ban environments, such as CO2 emissions, unpleasant smells of garbage trucks, holding up traffic in narrow streets, and the possi-bility of creating unsafe situations in traffic. Planning of picking up waste consists of multiple challenges which include selecting 1http://ec.europa.eu/environment/waste/studies/pdf/eucostwaste.pdf
which containers have to be picked up and what routes should be chosen. These challenges do also interact with each other, selecting different routes can lead to being able to visit different containers. The plan of the Municipality of Amsterdam to improve the cur-rent waste collection process in the city forms the basis for the project outlined in this paper. Every year on average 381 kilo waste per citizen is collected in this city[15], which leads to notable costs. Efficiency improvements have been made in Amsterdam over the years by, for example, adding containers into the process where citizens dispose their waste which led to a reduction in the number of stops a garbage truck has during a tour. Routes are adjusted based on the input of the truck drivers. The planning of these routes can however still be improved. Each container in the city has a 7 or 14-day pickup cycle. The filling speed of these containers, however, differs from site to site and do not exactly match these standard cycles. At each pickup, the net weight of the container is registered. This data is stored, but currently not used for planning purposes. This historical weight data gives the opportunity to set the next step in the improvement of collecting garbage and planning routes. Finding routes between multiple connected points is also known as the Vehicle Routing Problem [18]. The estimated filling level of a container, based on the historical weight data, can serve as input for a routing model with the goal to prevent visiting containers that have a low filling level. This could result in using resources more efficiently and a reduction in costs. According to a study by Ramos et al. [14], adding distance between points, to determine which sites to visit will lead to more optimal routes. The models used in this study are however tested in a provincial environment and not a city environment, where travel time is considered to be a better estimate of ’distance’ between sites.
The current situation in Amsterdam, available data and study of Ramos et al. lead to the following research question and sub-questions:
Main Research question:
To what extent can the automatic planning of dynamic routes for waste collection based on historical weight data increase efficiency? Sub questions:
(1) Can historical weight data be used to make a stable estimate of the total collected waste?
(2) Does adding travel time to determine the priority of sites lead to a more optimal solution?
(3) Can the method for waste collection of Ramos et al. be applied to city environments given specific constraints?
2
RELATED WORK
In this section, related literature is reviewed, started with an overall explanation of the Vehicle Routing Problem followed by filling level calculation. Finally, several studies on waste collection are discussed.
2.1
Vehicle Routing Problem
The problem of finding a set of optimal routes given a set of points and roads is known as the Vehicle Routing Problem (VRP) [18] which is part of the research field of Network Optimization. In this case the network consists of sites as nodes and roads as the edges between the sites. The ’depot’, the place where the trucks are emptied, is the starting- and endpoint of each route. The VRP has been first described by Dantzig and Ramser [4] as in the simplest form as a generalization of the Travelling Salesman problem, which consists of only one route as output.
As the aim is to find the routes with the lowest costs or the high-est profit, finding the solution to this problem is NP-hard [11]. This results in the fact that the explored problem space is limited and thus the absolute optimum cannot be guaranteed. This also leads to applying all kinds of algorithms and heuristics to this problem, such as branch-and-bound and centre tree as described in the research of Laporte. [10]. Laporte defines the vehicle routing problem as follow: The Vehicle Routing Problem (VRP) can be described as the problem of designing optimal delivery or collection routes from one or several depots to a number of geographically scattered cities or customers, subject to side constraints.
Laporte, page 1 [10]
The most common side constraints are capacity constraints and time window constraints, which introduces respectively the Capac-ity Vehicle Routing Problem (CVRP) and Vehicle Routing Problem with Time Window (VRPTW) [6]. The capacity constraint is de-pendent on the maximum capacity of a vehicle and the cost/weight of the transported goods, in this case the weight of waste. This introduces weights on the nodes of the graph as well. The time win-dows are applicable if delivery or pickup points are only available in specified time slots, e.g. if the truck has to be at a certain location during breaks.
2.2
Filling Level Calculation
In previous studies on garbage collection two types of input sources for the current fill level of a container are used: (1) real-time sensor data, as in the study of Ramos et al. [14], and (2) prediction based on historical data, e.g. Nuortio et al. [13]. In this study the second option is explored because sensor data is not available for each container in Amsterdam.
In essence this is trend prediction based on historical data which is used in a wide range of research fields, such as climate variability [8] and stock prediction [9]. Most of the inventory models assume that the demand distribution and the values of the parameters are known as described by Akcay et al. [1]. However, in practice this not always applies. Thereby most inventory models are based on fixed periods e.g. a week and not to estimate the inventory on a
certain day. This will limit the dynamic aspect of the desired model where the historical data can influence the cycles of single sites.
At the Municipality of Amsterdam previous research on the estimation of the filling level of containers showed that a regression model, based on a model of Winston [22], could be used where the mean amount of waste per day for a site is multiplied by the number of days the site is not emptied. On top of this estimate a safety margin is added. When the safety margin reaches the threshold, the truck should visit the site [20].
2.3
Waste Collection
There are multiple studies that are conducted on waste collection which is a combination of filling levels calculation and a vehicle routing problem. Not only the method differs in each study, but also the problem definition because the process of collecting waste is different in each city.
A good overview of research in this field is provided by Nuortio et al. [13]. They also conducted a case study in Eastern Finland. They first created feasible routes for each day and then applied several techniques such as Guided Variable Neighborhood Thresholding as meta-heuristic to improve these routes. This approach is however based on fixed pick-up cycles for each container. Eisenstein and Iyer [5] used a Markov decision algorithm based on the probability of collecting all waste in a sequence of blocks in Chicago depending on the number of dump visits. This method is really specific for a city which is divided into clear blocks of buildings. Another approach is to first create one big single route and then divide this route in smaller feasible routes based on, for example, capacity constraints. This route-first, cluster-second method is among others, applied by Maurão and Almeida [12].
Another, widely used, method is linear optimization. Chang et al. applied this method in Taiwan [3], Tung in Vietnam, (Hanoi) [19] and Ramos et al. in Portugal [14]. Linear optimization (LO), also called linear programming, is a method to find an optimum in a mathematical defined feasible area. An objective function is defined which e.g. minimizes cost or maximizes profit. This objective is subjected to one or more constraints. Together these define the feasible area which is enclosed by the x- and y-axis and linear lines given by the constraints. This method is very flexible due to the fact that single constraints can be added and removed easily.
The study of Ramos et al. applies three different LO-approaches to a case study. The first approach has as input sensor data about the filling level of the containers and a symmetric distance matrix. Based on the fill levels the containers that have to be picked up on a certain day are determined. These containers form a fully connected graph where the edges represent the roads between pickup points with the distance as weight. With the help of linear optimization the best set of routes is found. The objective is to minimize the costs, subject to constraints such as the maximum capacity of the truck. This process repeats each day. The second approach uses the same input data, but instead of determining the containers that have to be picked up, the objective function is changed to maximizing the profit which is based on the collected waste (profit) and distance travelled (cost). In this solution, constraints are added concerning the service level and allowed percentage of the total number of bins to overflow. This model is more complex and is better described as a 2
Mixed Integer Linear Programming (MILP) solution. This algorithm is executed each day as well. In the third approach, a heuristic is added to determine on which day the algorithm of approach 2 should be run to maximize the overall profit. The results are more kilogram waste per kilometer is collected, more profit and less distance travelled. The third approach scores the best on kilogram per km ratio, profit and distance travelled, but higher on vehicle usage rate in comparison to the first and second approach.
3
CASE STUDY
In this section background information and requirements are dis-cussed. The goal is to increase efficiency and thereby reduce the overall cost of garbage collection. The cost is expressed in time spent on collecting, distance driven by the trucks and the amount of waste collected.
3.1
Background
The Municipality of Amsterdam has the task of picking up the waste of citizens. Until recently each of the districts in the city man-aged the process by itself. The new situation aims for a centralized managed collecting process.
In each district separate (underground) containers are placed where citizens can dispose general waste, glass, paper, plastic and textile. An exception is the city centre where waste is picked up at the door. Each container is visited every seven or fourteen days. This results in a fixed schedule. The containers are placed in a well, a concrete pit in the ground. A group of wells that are placed next to each other is called a site.
When a container is being emptied, the container is picked up from the well with an arm on the garbage truck. Before emptying in the back of the truck the container is weighted while it hangs on the arm. After emptying the container is weighted a second time to determine the net weight. Each fraction of waste (e.g. general waste and glass) is picked up individually. An exception are the containers in the Northern part of the city where containers are swapped with a new container when picked up.
For this case study it has been chosen to focus on glass waste since glass has the most constant conversion ratio of kilograms to cubic meters. This conversion is necessary to calculate the fill level which has as input the weight data in kilograms and the maximum capacity/volume in cubic meters. General waste, paper and plastic have much more variety in this conversion. Geographically, the case study area limited to district West of the city based on the available historical data as will be discussed in the approach section.
3.2
Requirements
For a feasible planning multiple requirements have to be met. These requirements are either formed by physical limitations or by poli-cies. First of all, it is undesirable that a container overflows. This means that a container has to be picked up before it is totally full. An overflowing container will lead to additional placement of waste around the container which has to be picked up by another team of the municipality. This results in more overall costs, complaints of people living near the container and unnecessary garbage on the street.
Second, a truck has a maximum capacity of 8000 kilograms of glass waste, which has a safety margin included. If this limit is reached, the truck has to empty the back of the truck at the dump in the Western port area of the city.
There are also requirements regarding the working hours of the truck drivers. A truck driver is available for 8.5 hours each day at most. This is including 6.5 hours of driving time and picking up garbage, 30 minutes coffee break (09:00 - 09:30), 1-hour lunch break (12:00 - 13:00) and 15 minutes for personal hygiene (15:15 - 15:30). The drivers start their working day at 7:00 A.M., all other times as mentioned are fixed as well. The drivers do all start and end at the same time and have breaks together at the same location.
For weighing and emptying the container the driver has 1.5 minute per container plus an additional 1.5 minutes arrival time per site.
Each truck used for collecting glass waste in West is parked at the break location (Nieuw-Zeelandweg). This is the starting and end point of each day. All other routes start and end at the dump which is close to the break location (Aziëhavenweg).
In the approach section several more requirements regarding the routes are introduced in the form of constraints.
4
APPROACH
This section will describe the approach used for the case study on the glass garbage pickup in the city of Amsterdam.
4.1
Available Data
The municipality has multiple data sources combined in a single API which is publicly available via the data API portal2. This API is split up in two parts. The first part covers the data from the municipality itself and contains a description of each site, which is the geographical place where multiple containers can be placed. Each location of a site is described with a coordinate and the address of the nearest building. These sites contain one or multiple wells that can each hold one container of a specific dimension. For each site is also described how many containers are present and for which fraction of waste.
The second part covers data provided by third parties which are Welvaarts3, Enevo4and Sidcon5. In this study, the data from the supplier Welvaarts is used. This supplier provides weight data from each moment a container is picked up. This is stored with additional information, such as a time stamp, fraction and site id. Weight data is linked to the nearest site with a container of the fraction that is picked up based on the GPS location of the truck. The quality of this data is however not perfect. If the container does not hang free or still enough during weighing, the result is inaccurate. Thereby drivers can also be in a rush or weight moments are skipped when the equipment is broken. This all leads to inaccurate or missing data points in this weight data and thus has an influence on the filling level calculation.
For comparison to the current situation, the API data is not sufficient. Due to missing data points, the total time and distance 2https://api.data.amsterdam.nl/afval/
3Welvaarts: https://kilogram.nl
4Enevo: https://totalwastesystems.nl/enevo-one/ 5Sidcon: https://sidcon.nl
of the current situation cannot be calculated based on this data. The routes that are driven by the trucks are also tracked with GPS, but this data is also incomplete due to failure in the signal and/or other technical defects and thus not suitable to determine time and distance. There is also no access to the historical mileage of the trucks and there is no tachograph installed in the trucks.
However, additional data is available on the current planning. This is not publicly available via the API but as internal documents. The current planning is available as lists of addresses of the sites for each working day of the week. These addresses slightly differ from the addresses in the data from the API since a different source is used to create these routes. The breaks and visits of the dump are not included in these routes. Also, the net weights of the trucks at the dump sites are available with time stamp via Maltha6. The trucks are weighted before and after emptying which results in the net weight of waste that is picked up on a single route and with help of the time stamps also the number of visits to the dump on a day.
4.2
Input Data Processing
4.2.1 Filling Level and Weight.
The filling levels are required to determine which containers have to be emptied on a specific day. Therefore the current estimated weight of the container and the volume is needed. The volume is available via the API, the current weight is estimated with the help of the historical weight data. Due to the fact that until recently each district in the city managed the waste collection by themselves, there are large differences in the number of historical weight mo-ments per site between districts. Also, the frequency of collecting waste at a certain site has an influence on the number of weight points in the data set since weighting only takes place if a container is emptied. At the start of this research in March 2019, district West has on average the most historical weight points for glass per site (47.4) followed by East (31) and New-West (18.5). This is since the start of weighting which differs for each district, with the oldest record from begin 2017. Based on these numbers, district West will be used for this case study.
For each site in West where a container for glass waste is placed, general data is collected from the API, such as site id, coordinates, number of glass containers and volume of the containers. In total 161 sites are registered in the API which is less than the 208 sites on the current planning. This has an effect on the comparison of the performance of the models and the current situation and therefore a correction should be applied. This will be discussed in section 5.2.2. Next, for each site all historical weight moments are requested at the API. This is only required once for this study and all data is gathered in just three minutes. If this model will be put in production, new historical data should be loaded every day.
For each site, the number of days between a single weight mo-ment and the previous weight momo-ment is calculated. The first record cannot be used because the preceding interval cannot be calculated and is therefore filtered out. Next, for each site the outliers are filtered out of the historical net weights. A value is considered an outlier if it deviates two times the standard deviation from the 6http://www.maltha.nl/
mean. Examples of these outliers are net weights of (almost) 0 kilo-grams which appear in the dataset if the container is not weighted after emptying or extremely high values which can be the result of a measurement if the container did not hang still on the arm. These values are filtered out because they will not contribute to a realistic sight on the behaviour of the specific site and will lead to unnecessary calculations in the next step. In total 4.8% of the weight moments were filtered out.
According to the planning, each container has to be picked up every seven days in district West. In practice this could differ a bit to six or eight days if the actual tour e.g. falls on a holiday. In the data however are gaps of much more days, up to over 100 days. If the container was not picked up in these long periods, the net weight of the next pick up moment should be proportional high. This is however not the case since the extreme high net weights are already removed as an outlier. The values of net weight are similar to the periods of 6/7/8 days. These large time periods will lead to a too low estimate of weight per day if the registered periods are used. The assumption is made that in these long periods the waste is picked up, but not weighted. The registered net weight is thus based on a similar number of days as the other registered net weights and therefore considered to be seven days as planned. With the help of the number of days between each pair of successive weight moments an average weight per day can be calculated for each period. On this set of values the same outlier filter is applied as on the total net weights to filter out the slipped through outliers, which filters out another 4.0% of the values.
After the outlier filtering process, a set is established of the average amount of waste per day disposed for each site. With this data, the estimated weight on a day in the future can be calculated. A random sample of the size of the number of days a site is not emptied is taken from the set of average weights per day. The sum of this sample is the estimated weight. See also formula 1 where d is the number of days the site is not emptied. This method incorporates the fluctuation of filling speed by taking a different amount of kilograms for each day of the period, and therefore this method is chosen over a simple linear function. To ensure the same outcome of the simulation, the day number is set as a seed for the random sampling (25 for 25th of January).
sum([random kд/day 1, random kд/day 2, ...,d]) (1) This estimated amount of waste on a site is used for deciding which sites to visit and if a truck can empty the container based on the maximum capacity of the truck. Because glass containers in the city do not all have the same volume, a certain weight means for one container almost full and for another just half full. To decide which sites to visit, the filling level is determined with the help of the estimated weight. For the conversion of weight to volume, 300 kg/m3is used as given by Stimular7. With this number the filling level is calculated as follows:
f illinд level = (estimated weiдht/300) / volume (2) If there are two or more glass containers present on a single site, the containers are seen as one large container for the amount of
7http://www.duurzamebedrijfsvoeringoverheden.nl/themas/afval/hoeveelheden.html
waste and maximum capacity. This is the current working method and it has been chosen to apply the same rule in this study. 4.2.2 Distance and Time Matrix.
To lower the needed computational power to solve the problem of finding routes, the routes from each site to each site and between each site and the break- and dump location are pre-computed. This way the algorithm does not have to compute the best routes over the whole road network of the city. The matrices are asymmetric since a truck cannot always drive the same route from A to B and from B to A because of for example one-way streets. For retrieving the matrices the Bing Maps Routes API8is used. This API provides in the Enterprise version an option for truck routing which takes properties of the vehicle such as dimensions into account. During this study this version could not be used. Yet is chosen for this API over alternatives such as Google Maps with an eye on the future when the enterprise version might be available. Thereby, the algorithm is designed to use retrieved matrices from any source.
The API has a limitation of 2500 result combinations per request, which means per request only 50 * 50 geolocations can be retrieved, regardless of the license. Therefore the requests have to be made in batches. Both the distance and time matrix are represented and saved as key/value pairs with a tuple of site ids as key and the distance/time as value:
{(siteid1, siteid2) : distance, ...}
The travel time and distance are based on normal traffic con-ditions without traffic jams. Historical traffic information is not available for this study. In the models that will be introduced later, the travel time is taken as input for optimization. To provide a more realistic travel time, an additional 25% is added to the travel time to and from each node in the network. This number is based on the TomTom Traffic Index for Amsterdam [17] which indicates that in the last year the average tour in Amsterdam took 25% longer than usual due to congestion. If the model will be used for real planning, the distance and travel time matrix can be retrieved real time right before the routes are calculated.
4.3
Route Optimization
4.3.1 Time Windows.
The fixed times of coffee and lunch break divide the working day into three parts for the truck drivers. In the two models introduced in the next sections, this partition will be handled in different ways. In the first model the routes will be calculated for the whole day. The driver can pause the routes any time to take a break. It is taken into account that the driver has to drive to the break location at most three times a day (two breaks and end of the day) and from the break location to a point in the network at most twice a day (after both breaks). For this, the mean driving time from a node in the network of that specific day to and from the break location is used in the route calculation.
In the second model the day consists of three time windows, each with limited time. The duration of the route of the first time window, from 7 A.M. till 9 A.M., is 115 minutes at a maximum which is 120 minutes minus five minutes to drive from the dump 8
https://docs.microsoft.com/en-us/bingmaps/rest-services/routes/calculate-a-distance-matrix
to the break location. The second time window has a limit of 140 (150-2*5) minutes and the third 125 (135-2*5) minutes. In this model routes are calculated for each separate time window in a day. 4.3.2 Model 1: Minimizing Time.
The first model consists of two steps. First, the sites that need to be visited are determined. Containers are not supposed to overflow and therefore it is chosen to pick up containers that exceed an estimated filling level of 75% which gives a 25% margin to prevent overflow. This is based on a goal set by the Municipality and a good target which will be a huge improvement to the current situation of 17.12%. Sites that reach a fill level of 75% or higher on the next workday are selected. On a Friday the filling level of containers on Monday are calculated because the waste in district West is not picked up during the weekend. Sites with little historical weight data are also added to the selection if the last visit has been seven days ago. This is the current pick up cycle. These sites will gather more historical records and will be planned in more dynamically if there are enough historical data points. A site is considered to have too little available data if it has up to and including 10 historical weight moments. This number is chosen in order to guarantee enough possible variation for calculating the estimated weight for periods around the seven days of the current cycle. For example, for a period of eight days, eight different amounts of weights can be summed to simulate the variation of disposed waste. There are in total five sites in the data set which have little available data.
Each model is simulated for two scenarios. In the first scenario the set of sites is determined as described above, in the second scenario sites are added that are not visited in the last 14 days. This scenario is also taken into account because the drivers do not only pick up the waste, but also check if the containers are in a good condition and can notify other teams if there is bulky waste disposed right next to the containers.
With this set of selected sites, the calculated amount of waste and the time matrix, a directed graph can be created with weighted nodes and arcs. Because of the limited available calculation power and increasing mathematical complexity for each extra node in the network, a maximum of 35 containers is set for the graph. First all sites that have to be picked up due to exceeding the filling level threshold are added to the graph, second the sites that have limited historical data and if this total amount of sites does not exceed 35, the sites will be supplemented with sites that are exceeding the 14-day period.
The second step is to calculate the routes to pick up these con-tainers with the help of linear optimization principles. A schematic overview of the process of model 1 can be found in Figure 1.
Figure 1: Steps Model 1.
The model for finding routes requires the following input: a list of sites (nodes) without the depot (N), a list of nodes including the depot (V), maximum capacity of the vehicle (Q), a set of key/value 5
pairs with site id as key and estimated weight as value (q), a list of all possible roads in the form of tuples (A) and a set of key/value pairs of costs in time per road (c) where the key is a tuple from A and the value the travel time. The list of sites is determined as described earlier in this section.
Besides these variables, two more model variables are needed. The first one is a binary variable which keeps track of active and inactive arcs (x). For each possible arc in A, x can be 1 if the arc is active, thus a truck drives this arc from one site to the other, or 0 if this arc is not used. A second variable (continuous) is used to keep track of the picked up weight by the truck (u). This variable has an upper bound which is equal to the maximum load of the truck defined as Q.
The objective of this model is to minimize travel time. Travel time is chosen over travel distance since it could be beneficial to drive more kilometers, e.g. over the highway, which costs less time. Time savings can also lead to better deployment of the workforce since there is time over to spend on other jobs. The objective is programmed as follow:
minimize(sum(c[i, j] ∗ x[i, j] + 1.5 ∗ x[i, j] + (3) 1.5 ∗nr_containers[j] ∗ x[i, j] f or i, j in A) +
4 ∗break_travel_time)
In this formula, c[i,j] is the cost in time of travelling from node i to j. This is multiplied by x[i,j] which is 0 if arc from site i to j is not active and 1 if the arc is active. Furthermore, 1.5 minutes are added for each visited site and an additional 1.5 minutes for each container on these sites. Lastly, four times the average driving time from all nodes in the network of the specific day from and to the break location is added. If all containers can be picked up before one of the two breaks, too much break travel time is added in this formula. This is added in the first place to make sure that a full schedule will fit in a full day, but will be corrected after route calculation. For the results of the performance of this model, the break travel time that is added and was not needed is subtracted from the total.
This objective function is subjected to multiple constraints. First of all, each node has to be visited and only once. This requirement is met by the following two constraints. The first constraint makes sure that each site in N should only have one active incoming arc (sum must be equal to 1). The second constraint makes sure each site has one active outgoing arc, see:
sum(x[i, j] f or j in V i f j! = i) == 1 f or i in N (4) sum(x[i, j] f or i in V i f i! = j) == 1 f or j in N (5) Since the variable ’u’ has an upper bound equal to the maximum capacity of the truck, the capacity constraint is already met. How-ever, the check has to be made that a route does not exceed this maximum. This is done by making use of an indicator constraint which only applies if x[i,j] is 1, if the arc i to j is active. The con-straint itself makes sure that the new total amount of waste at the end node of the arc is the sum of the total waste at the begin node of the arc plus the amount of waste at the end arc. This applies to each possible arc in the set A where the begin or end node is not the depot (0). Thereby, a container can only be emptied totally and
the truck can not dump waste from the truck in a container. This constraint is met by making sure that the cumulative amount of picked up waste at node i is at least as much as the amount of waste that has to be picked up at i for each node in N and thus no waste is dumped from the truck in the container.
indicator_constraint(x[i, j],u[i] + q[j] == u[j]) (6) f or i, j in A i f i! = 0 and j! = 0
u[i] >= q[i] f or i in N (7) The routes also have to fit in the total available time on a day. This is the sum of the time before the coffee break, between coffee and lunch break and after lunch. Which is in total 405 minutes. The maximum travel time is implemented with the following constraint which makes sure that the objective is less than the total available time:
sum(c[i, j] ∗ x[i, j] + 1.5 ∗ x[i, j] + 1.5 ∗ nr_containers[j] ∗ x[i, j] (8) f or i, j in A) + 4 ∗ break_travel_time <= 405) 4.3.3 Model 2: Maximize Profit.
The second model has a lot in common with model 1 but does not consist of the two steps, see Figure 2. The input network consists of all sites. The objective function is based on maximizing the profit instead of minimizing the cost, inspired by the study of Ramos et al. [14]. This profit can be described by income minus cost and makes it possible to decide which site to visit based on multiple conditions, in this case filling level and travel time. Ramos et al. used the amount of waste as income since revenues are generated from selling the waste and the cost is based on travel distance. This results in the fact that a site with more waste becomes more interesting to visit than another site with less waste with the same distance. In this case study the sites have a different number of containers with different capacities which results in different maximums. It could thus be the case that a site which contains less waste has a higher priority to pick up than a site which more kilograms of waste. Therefore in this model the filling level is used as income. For the cost, it is chosen to use the time as in the first model. Both the income and time are values between 0 and 1. Where income is a percentage and cost is the normalized time where the site which has the furthest travel time has value 1.
Figure 2: Steps Model 2.
This model can also pick up waste at sites that are not above a certain threshold as in the first model, but still worth to pick up considering filling level and travel time. The sites with filling levels above the threshold of 75% should be picked up regardless of the travel time to this site and thus should be profitable to visit. The same holds for sites that have less than 10 historical weight records and in the second scenario for sites that exceed the maximum period of 14 days. In order to make sites with higher filling levels way more profitable than the sites with low filling levels, the filling level 6
is squared. Sites that have to be visited due to one of the above-mentioned rules must have a good profit and therefore do not have costs in the form of time. Sites that are not mandatory to visit get penalty by subtracting the travel time as a cost. A container that has a relative distance of 0.3 in the network and a filling level of 0.5 has thus a ’profit’ of -0.05 and is unattractive to pick up. If the filling level is 0.7, the profit is 0.19 and thus will be considered to pick up, because it will raise the profit. If this container had a distance of 0.6 relative in the network the profit was -0.11 and thus not attractive. Now, if the container has distance 0.3 and exceeded the threshold of 0.75 with a filling level of 0.8, the profit will be 0.64 and thus compared to the other examples really profitable to pick up.
The filling level of a site with little available data and/or exceed-ing the 14-day rule is adjusted to 75% in order to make it more attractive. The amount of waste is however still based on the little historical data that is available since this is the best estimator for these containers. The objective function is as shown below where fadj is the adjusted filling level, cnorm the normalized travelling time and p the penalty of either 0 or 1.
maximize(sum(f adj[j] ∗ f adj[j] ∗ x[i, j] − (9) cnorm[i, j] ∗ x[i, j] ∗ p[j] f or i, j in A)
In the first model all nodes in the network had to be picked up, but in this model each node can be visited at most once or can be skipped. Constraint 1 and 2 have to be relaxed where each node has zero or one active incoming arc and zero or one active outgoing arc:
sum(x[i, j] f or j in V i f j! = i) <= 1 f or i in N (10) sum(x[i, j] f or i in V i f i! = j) <= 1 f or j in N (11) This can lead to a solution of fragmented active arcs which do not form a closed route. The sum of active incoming arcs must be equal to the sum of active outgoing arcs for each node in the network. The following constraint has to be added to make sure closed routes are found:
sum(x[j, i] f or j in V i f j! = i) == (12) sum(x[i, j] f or j in V i f j! = i) f or i in V
As of the last constraint, the routes do have to fit in the current time window. In contrast to the first model the routes are now specifically calculated for a certain time window and not the whole day. This is handled by limiting the sum of travel time, 1.5 minute for each site and 1.5 minute for each container. This sum must be less or equal to the adjusted, as described earlier, available time in the current time window.
sum(c[i, j] ∗ x[i, j] + 1.5 ∗ x[i, j] + (13) 1.5 ∗nr_containers[j] ∗ x[i, j] f or i, j in A) <=
timewindow_dict[current_timewindow])
The remaining containers that do not fit in the route due to the maximum capacity of a truck (Q) or the time limit will be picked up in the next time window.
Because the network contains all available sites, solving this mathematical problem is considered to be much more intensive than at the first model. For the simulation is chosen to only select
sites with a filling level of 55 percent or higher for the network. This way the size of the network decreases and therefore the total problem space. This means that sites with a filling level of less than 55% are not visited, sites with a filling level between 55% and 75% are considered to be visited and above 75% are most likely emptied.
4.4
Implementation
The processing of the data and implementation of the models is done in Python 3 [21]. For the implementation of the routing models, the DoCplex library version 2.9.141 of IBM [16] is used. This library is a Python layer for accessing the Cplex Solver of IBM Optimization Studio9. For using this solver on large mathematical problems a license is needed. In this study, the academic license of version 12.9.0.0 is used which is free to request for students and researchers. The implementation in Python of model 1, and with that of the second model, is based on an implementation presented by Hernán Cáceres [2]. All simulations are run locally on a 2.6 GHz Intel Core i7 quad-core processor.
5
SIMULATION
To evaluate the models and compare these models with each other and the current planning, an offline simulation is used. The first quarter of 2019 is chosen as simulation period. This way all histori-cal data from 23-02-2017, the first weight moment, till 01-01-2019 could be used for the filling level calculation, which are 8520 records in total. For each day the filling levels are calculated, together with the estimated weight. For the first model for each day and for the second model for each time window in a day, routes are calculated. Between the time windows on a single day the filling levels are not recalculated, only the sites that are emptied are adjusted to be empty.
The simulation differs on an important point from a real-life implementation. In the real-life implementation, the number of historical weight moments will increase each time the site is visited. Therefore the containers with little historical data will be planned dynamically if the number of weight moments exceeds 10. During the offline simulation no new historical data is added if a site is visited, because for these specific sites the estimation of filling level is still based on the few data points available. Adding the estimations as new ’historical’ data in the simulation will lead to an unreliable estimation for the dynamic route planning. This results in the fact that sites with less than 10 historical weight moments are picked up in a seven-day cycle for the whole simulation period. The simulation is run by a Python script. For the first day, the initial estimates are calculated based on the available historical data. This is stored in a data frame and saved as a Pickle file. After calculating the routes, the estimates of the filling level and weight of the emptied containers are adjusted and saved for the next day/time window. Also, the solution output of the solver is saved, which contains the objective value, active arcs and the load weight in the truck at each visited site.
5.1
Solving Time
The solver runs until it has explored the whole problem space. The required time for this can not be determined in advance. A time 9https://www.ibm.com/products/ilog-cplex-optimization-studio
limit of 60 seconds is applied to the solver for each set of routes that have to be calculated. This number is chosen based on a sample of situations that is tested for 15, 30, 60 and 120 seconds. Between 120 and 60 seconds the improvement was almost nothing. This means for the first model each day takes at most 60 seconds to calculate, which comes to a total of 67 minutes for three months. This can be shorter if a day no routes have to be calculated or an optimal solution is found earlier by the solver. For the second model this time limit of 60 seconds is set for each time window. The calculation for one day can thus take up to three minutes, which results in 3.5 hours at most in total. Again this is the most extreme situation.
5.2
Evaluation
The two created models in two scenarios and the current planning are compared on several factors. In both models the time spent picking up waste plays an important role. This is one of the factors to compare with. In addition, there is also comparison on amount of waste collected in kilograms, the distance driven by the trucks, the mean filling level of the sites that are visited, the total number of visits of the dump location and number of days that the truck does not have to drive. Finally, the models are compared with the current planning on kilograms waste collected per minute and per kilometer driven because this is the best estimator of efficiency and cost reduction.
5.2.1 Comparing Models.
The two models can be compared in a uniform way since both models have the exact same input, namely the starting point of all the historical data up to 2019, and weight estimation method. There is randomness involved in the simulation because the sum of a random sample of a set of ’kilograms per day’ is used as an estimated weight. By using a seed as described in the approach section and by taking the average of the above metrics over a period of three months, the influence of this randomness is limited. For the metrics based on average, the days on which the trucks do not have to drive in the simulation are also included in the calculation. 5.2.2 Compare to Current Planning.
Comparing the models with the help of the simulation to the current situation is a bit more complicated due to the incomplete data and difference in API-data and the planning as described earlier. Chosen is to determine the spend time and distance based on the current planning. This way the planning of the two models is compared with the current planning. With the help of the planned route list for each day of the week and Bing Maps, the routes are calculated which results in distance and travel time for each day of the week. An additional 25% is added to the travel time, the same as with the calculated travel time for the simulation. Also, the travel time to and from the break location is added as described earlier in the approach section. Since the limited Bing Maps edition without truck routes support is used, the same bias is introduced for both the routes of the models as for the current situation. Because the planning consists of 208 sites and the API of 161, the final averages of sites, minutes and kilometers per day are corrected by multiplying by the factor of 161/208.
Evaluation of the calculation of the amount of waste at a site is also a challenge with the incomplete data in the API. Since a
majority of the sites miss weight moments, there is chosen to use the data of Maltha of total net weight of the glass that was dumped because this gives the most accurate insight in the amount of waste that had to be picked up in this period. In the first quarter of 2019 also only 161 sites have a weight registered in the API, assumed is that in the real situation only the 161 sites are visited instead of the planned 208. There is thus no need to apply a correction factor. The total amount of waste that is collected in the simulation should be around the total amount that is collected in the real situation.
5.3
Results
In this section the results of the filling level calculation and the sim-ulation of the two models are discussed on the metrics as introduced in the approach section.
5.3.1 Simulation Models.
Figure 3 shows, among others, the results of the first model, which minimizes time, for the two different scenarios. Fewer sites are visited on average each day in the first scenario than in scenario 2, nine sites in comparison to 19 sites per day. The average time spent on collecting the waste also increases with 47% from 158 minutes to 233 minutes per day and driven kilometers increases with 32% from 78.23 to 103.37 kilometers per day. As a result of picking up more containers which are expected to not have reached the threshold of 75% due to the 14-day rule, the mean filling level of containers that are emptied decreases from 75.56% to 48.81%, which is a decrease of 36%. The mean filling level for optimized containers has a higher value since this only contains the containers that are emptied because of the exceeded 75% threshold, 80.47% in scenario 1 and 85.51% in scenario 2. The second scenario results in eight more rides to the dump location (from 116 to 124) than in the first scenario. The second scenario has two days in the simulation period where no containers have to be picked up, compared to 0 in the first scenario.
The second model is simulated for the same two scenarios as model 1. The same behaviour can be seen between the two scenarios. The mean number of sites per days increases from almost 11 in scenario 1 to 20 in scenario 2. The time spent on collecting increases with 45% from 108.82 to 157.03 minutes per day and amount of kilometers with 20% from 53.94 to 63.66 on average per day. The mean filling level is in scenario 1 67.95 and scenario 2 43.68, a decrease of 36% comparable to the decrease seen in model 1. The filling levels of the optimized containers are 71.05% and 66.05% for scenario 1 and 2 respectively. The dump site is visited 12 more times in the second scenario and in the first scenario, 108 compared to 96. Six days no containers have to be emptied in the first scenario comparison to seven days in scenario 2.
In the simulation of the second model, scenario 2, two sites were not emptied on time on the first day of the simulation period. In this simulation also nine sites (seven on the first day) that were not visited for 14 days were visited one day later. In all other simulations, all containers were picked up on time.
5.3.2 Statistics Current Planning.
Figure 3 also shows the statistics of the current planning in the first quarter of 2019. The mean sites per day is based on the route lists and corrected for the number of sites in the API, which results in 8
Figure 3: Results of Simulations of Model 1 and 2 in both scenarios compared to current situation. 33,45 sites per day. Corrected mean time and distance per day are
based on Bing Maps and are respectively 287,64 minutes and 60,66 kilometers. The mean filling level is based on the available data via the API and results in 17.12%. In total 99 visits to the dump are registered by Maltha in the simulation period. At least one weight moment on each working day was registered and thus every day a truck had to pick up glass waste.
5.3.3 Amount of Waste.
The amount of waste collected in each simulation is different. Figure 4 describes the total amounts in 1000 kilograms and the mean kilogram collected per minute and per kilometer for each model and scenario. In the first model, second scenario, the most waste is collected compared to the other three simulations and the real situation. Per minute both models perform better in both scenarios (4.73, 3.49, 6.48, 4.90 respectively) than the current situation of 2.59k kg per minute. Per kilometer only the second model in the first scenario performs better with 13.32k kg/km compared to the current situation of 12.29k kg/km. Model 2 collects slightly less in scenario two with 12.09k kg/km and the first model collects 9.56k kg/km and 7.88k kg/km respectively for both scenarios.
Figure 4: Kilograms waste collected per minute and per kilo-meter per model and scenario.
Figure 5: Mean distance and time per day of the week from Monday till Friday for Model 1 and 2, scenario 1.
5.3.4 Week Planning.
In both models a difference of the mean amount of kilometers driven and time spend on collecting is seen for each day of the week. Notable is the high means for both metrics on Friday, see Figure 5. The planned routes in the simulations are more dynamic than the current routes, e.g. not every Monday the same route is planned. The spread of sites that has to be visited is also larger than the current static routes as can be seen in Figure 6.
Figure 6: Planning on Mondays for Model 2, scenario 1: Cur-rent Planning (L), 21-01-19 (M), 04-03-19 (R). Same color be-longs to the same route.
6
DISCUSSION
The second model has a lower mean in distance and time spent on collecting in comparison with the first model. This seems counter-intuitive since this model can visit other sites on top of the sites that are visited in the first model. Visiting more sites would lead to more distance travelled and longer time spending. However, on a larger period this really pays off by days that a truck does not have to drive at all. Picking up waste from a site which did not reach the hard threshold of 75% yet, does have an advantage if the site is close to sites that have to be emptied.
The higher means on distance and time on Fridays can be ex-plained by the fact that on Friday also all the sites have to be visited that are exceeding the thresholds of filling level or number of days not emptied on Saturday and Sunday. Further research could inves-tigate the possibilities of spreading this ’weekend load’ over the week for a more balanced work week.
The mean filling level is slightly lower for the second model compared to the first model despite the fact of outperforming the first model on other factors. This also has a benefit because the chance of overflowing will be less.
The first scenario performs better for both models. This can be explained by the fact that in the second scenario containers with a cycle longer than 14 days has to be picked up before the threshold of the filling level is reached.
The better performance on the amount of waste collected per minute in comparison with the current planning was expected. The current planning is not dynamic and historical weight moments were not taken into account. Fewer containers have to be picked up which results in the decrease of travel time. The increase in driven kilometers in three of the four simulations is notable and can be explained by the fact that the spread of sites that has to be visited is larger as can be seen in Figure 6.
6.1
Limitations Data
The historical data that is used is limited and incomplete. As seen in the results, the mean filling level of containers is lower for all picked up containers than for the optimized containers that have enough historical data available and thus can be picked up based on the estimated filling level. This could indicate that with more, and more complete, historical data the mean filling level will increase. The limited available data about the current situation and plan-ning in the simulation period and the difference in the number of sites between the API and current planning makes the comparison of the models to the current situation hard and therefore determin-ing the improvement. The amount of waste collected is compared to the real collected amount of waste in this period with the as-sumption that this only covers the containers registered in the API. The distance and time are compared to the current planning cor-rected for the difference in the number of sites between API and the planning. Ideally, the comparison would be made based on the same data set and the same calculation methods. There is, therefore, room for improvement of the comparison of the performance of the models to the current situation.
In further studies of the implementation of this model in real-life, the route network has to be adjusted for trucks where the dimensions and weight of trucks are taken into account. This can easily be implemented by creating net matrices as input for the algorithm.
6.2
Filling Level Calculation
The used method for calculation and evaluation of the estimated filling level of a site on a certain point in time is not advanced. In the simulations for both models, the same method is used and therefore this has no influence on the comparison. The variation of the total amount of waste collected in the simulations could be explained by the randomness that is involved in the calculation of the estimated amount of waste at a site. Since the total amount is not consistently lower or higher than in the real situation in the four simulations, a constant bias could not be concluded. The better performing model, model 2, has in both scenarios less waste collected than model 1 which could reflect in the final comparison. It would be interesting to further study this result in a following study.
Based on the fact that the total amount of waste collected in the simulations is around the amount registered at Maltha, the used method seems to be a reasonable estimation. Interesting would be to compare and validate this method on individual site level if there is more (complete) data available. The filling level calculation could also be extended with factors such as seasonality and holidays in further research.
6.3
Other Limitations
Limited computing power of one CPU is used for running the simulations. With more computer power, or setting a greater time limit, a larger part of the problem space could be explored. This could lead to more efficient routes and an overall more optimal solution.
For further research it could be interesting to investigate the influence of changing the fixed parameters that are used in this study, such as the 75% threshold for filling level, the minimum of 10 historical weight values, maximum graph size of 35 nodes in the first model, the seven day cycle for sites with low amount of historical records and the 14-day rule. These parameters are chosen based on the explained arguments, but further research could possibly lead to better results.
6.4
Scalability & Generizability
The case study is performed on a small district compared to the whole city and only for one fraction of waste. However, the model could be used for other districts or cities as well as long as the context is similar. The first part of the process, the filling level calculation, would work the best on glass fraction since it has the most constant conversion ratio from weight to volume. The second part, the route calculation, could be used in combination with any other filling level calculation method, or with the input of sensors which e.g. are present in a lot of plastic containers. The results in increasing efficiency compared to the current situation depends of course on the efficiency of the current situation and will, therefore, be different in each case.
Increasing the network size will lead to an increase in possible solutions and thus a bigger mathematical problem. For model 1 a maximum network size of 35 is chosen. It would be interesting to conduct further research on the performance of these models on a larger scale with more computational power.
7
CONCLUSIONS
The main research question, to what extent can the automatic plan-ning of dynamic routes for waste collection based on historical weight data increase efficiency, can be answered with an increase up to 149% in the amount of waste collected per minute. This also results in a 7,7% increase in the amount of waste collected per kilo-meter. This is the result of a simulation in the city of Amsterdam with a model that determines which sites to visit based on filling level and travel time and maximizes profit. With this result must be mentioned that there could be made improvements in the com-parison between the model performance and the current situation. Besides the model which maximizes profit, another model which determines the sites only based on filling level and minimizes time is also tested in this study. The maximizing profit model performs better on time spent and kilometers driven, and the amount of waste collected per minute and per kilometer compared to the minimizing time model. This holds for two different scenarios, wherein the second scenario each site should be visited at least once in every 14 days as an additional constraint. It can be concluded that the maximizing profit model creates a better planning for cost reduction than the minimizing time model. Compared to the current planning in district West in Amsterdam, both models led to collecting a 10
greater amount of waste per minute. Only the maximizing time model also led to collecting more waste per kilometer in the first scenario.
The used method to estimate the amount of waste in a container based on historical data results in a similar total amount of picked-up waste in the simulations as is collected in the real situation in the same period. Therefore the answer to sub-question 1 is yes. Both models determine the sites to visit in different ways. The second model, which also uses travel time, leads to better results and therefore the second sub-question can be answered positively as well. The results show that the models of Ramos et al. can, with some adjustments, be applied in city environments as well, which answers the third sub-question.
The results indicate that applying these models can lead to a decrease in several factors influencing the cost of waste collection in city environments. A few suggestions for further research based on obtained insights are made as well, such improving the weight estimation method, implement truck-specific routes and compare the performance of the models in the simulation to the real situation based on the same source data. This makes this study valuable for both the research field and the Municipality of Amsterdam.
REFERENCES
[1] Alp Akcay, Bahar Biller, and Sridhar Tayur. 2011. Improved inventory targets in the presence of limited historical demand data. Manufacturing & Service Operations Management 13, 3 (2011), 297–309.
[2] Hernán Cáceres. 2018. CPLEX Python. Capacitated vehcile routing problem. https://youtu.be/-hGL39jdtQE Accessed: 2019-04-17.
[3] Ni-Bin Chang, HY Lu, and YL Wei. 1997. GIS technology for vehicle routing and scheduling in solid waste collection systems. Journal of environmental engineering 123, 9 (1997), 901–910.
[4] George B Dantzig and John H Ramser. 1959. The truck dispatching problem. Management science 6, 1 (1959), 80–91.
[5] Donald D Eisenstein and Ananth V Iyer. 1997. Garbage collection in Chicago: a dynamic scheduling model. Management Science 43, 7 (1997), 922–933. [6] Nasser A El-Sherbeny. 2010. Vehicle routing with time windows: An overview
of exact, heuristic and metaheuristic methods. Journal of King Saud University-Science 22, 3 (2010), 123–131.
[7] Daniel Hoornweg, Perinaz Bhada-Tata, and Chris Kennedy. 2013. Environment: Waste production must peak this century. Nature News 502, 7473 (2013), 615. [8] Walter Immerzeel. 2008. Historical trends and future predictions of climate
variability in the Brahmaputra basin. International Journal of Climatology 28, 2 (2008), 243–254.
[9] K Senthamarai Kannan, P Sailapathi Sekar, M Mohamed Sathik, and P Arumugam. 2010. Financial stock market forecast using data mining techniques. In Proceedings of the International Multiconference of Engineers and computer scientists, Vol. 1. 4. [10] Gilbert Laporte. 1992. The vehicle routing problem: An overview of exact and approximate algorithms. European journal of operational research 59, 3 (1992), 345–358.
[11] Jan Karel Lenstra and AHG Rinnooy Kan. 1981. Complexity of vehicle routing and scheduling problems. Networks 11, 2 (1981), 221–227.
[12] M Cândida Mourão and M Teresa Almeida. 2000. Lower-bounding and heuristic methods for a refuse collection vehicle routing problem. European Journal of operational research 121, 2 (2000), 420–434.
[13] Teemu Nuortio, Jari Kytöjoki, Harri Niska, and Olli Bräysy. 2006. Improved route planning and scheduling of waste collection and transport. Expert systems with applications 30, 2 (2006), 223–232.
[14] Tânia Rodrigues Pereira Ramos, Carolina Soares de Morais, and Ana Paula Barbosa-Povoa. 2018. The smart waste collection routing problem: Alterna-tive operational management approaches. Expert Systems with Applications 103 (2018), 146–158.
[15] Ministerie van Infrastructuur en Waterstaat Rijkswaterstaat. 2019. Afvalmoni-tor 2017. https://afvalmoniAfvalmoni-tor.databank.nl/Jive/Jive?cat_open=Gemeentelijk% 20niveau Accessed: 2019-06-24.
[16] The IBM Decision Optimization on Cloud team. 2019. DoCplex. https://pypi.org/ project/docplex/
[17] TomTom. 2019. TomTom Traffic Index Amsterdam. https://www.tomtom.com/en_ gb/traffic-index/amsterdam-traffic#statistics Accessed: 2019-04-15.
[18] Paolo Toth and Daniele Vigo. 2002. The vehicle routing problem. SIAM.
[19] Dang Vu Tung and Anulark Pinnoi. 2000. Vehicle routing–scheduling for waste collection in Hanoi. European Journal of Operational Research 125, 3 (2000), 449–468.
[20] Bijl B. van der Habbou I. Duijvestijn D. Uijl J. den, Oskam J. 2017. Container-netwerken Amsterdam: Voorspellen vulgraden containers. Internal Report of Study Conducted by Students of HvA.
[21] Guido Van Rossum and Fred L Drake Jr. 1995. Python tutorial. Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands.
[22] Wayne L Winston and Jeffrey B Goldberg. 2004. Operations research: applications and algorithms. Vol. 3. Thomson/Brooks/Cole Belmont Calif.
