Help with school bus routing on a map with known bus stop coordinates

[AsharFatmi]

I want to generate all the Routes for a school with the following given data:-

• Capacity of the busses
• Location [lat,lng] of the school
• Location [lat,lng] of all the stops.

I need to generate the most optimal routes based on the stops, so was thinking of applying heuristics for the solution (Gillet & Miller Sweep heuristic).
All [Lat, Lng] data are available.

• Blue Marker in the map is the school Location
• Yellow Markers are the Stops

?

How can I use the VeRyPy for this ?

[yorak]

Dear Ashar,

definitely you can. What you need to do, is to calculate a distance matrix based on location [lat,lng] information of the school and stops. You can refer to the example in the VeRyPy repostory: https://github.com/yorak/VeRyPy/blob/master/examples/bing_maps_example.py

The example script takes a tab separated file with a list of addresses (one for each row) as an input, but you could modify it to accept coordinates. I'd propose adding a new command line flag to input such file. If you need any help implementing the feature to the bing_maps_example.py file, please let me know. Also, note that you need to get a Bing Maps API key to calculate the distances in the road network: https://www.bingmapsportal.com

Please note that the developer key only allows 50 locations (https://docs.microsoft.com/en-us/bingmaps/rest-services/routes/calculate-a-distance-matrix). From the map it seems you have more? If you have the time, this can be done in pieces - or you can use another map server that can calculate larger distance matrices for you.

Hope this allows you to continue. Please let me know if you hit any roadblocks.

[AsharFatmi]

Hi Yorak,

Thank you for the prompt response, I really appreciate it.

Just a few queries before I dive into this:-

• The distance matrix has to be the onRoad travel distance right? if yes then I can create my own distance matrix using distance from the normal route and directions APIs increasing the run time of course [ Services like OSRM etc. ]?
• Can you explain the format of the .vrp file as if I create my own distance matrix I would have to create the .vrp file and I can't understand the structure of data in that file
• The distance matrix calculation has to be from school to all the individual stops only or also between all the stops? one to many or all to all?

Thanks in advance

[yorak]

Hi,

no problem, I'm just happy if I could be of assistance. Regarding your questions:

The distance matrix has to be the onRoad travel distance right? if yes then I can create my own distance matrix using distance from the normal route and directions APIs increasing the run time of course [ Services like OSRM etc. ]?
The distance matrix calculation has to be from school to all the individual stops only or also between all the stops? one to many or all to all?

Yes, the distance matrix is 1+n x 1+n matrix D, where n is the number of bus stops and the 1 being the school. The values are the distances or travel times it takes to move from any stop to another (or straight from/to the school). Hence, it contains precalculated distances/travel times for all possible combinations of locations. When computing such a matrix, be sure to use a API that is capable of effectively computing such distances as the matrix size is quadratic to the number of stops and can grow quite large. Here, I'd recommend against rolling out your own solution for the on road distance calculation, because it is another, quite extensive, can of worms (but, if you are interested, please check bidirectional Dijkstra, R-trees, and contraction hierarchies). Hence, if you already use a MAP API, please check if it provides a distance matrix computation service.

edit: In OSRM this seems to be called table service, docs https://github.com/Project-OSRM/osrm-backend/blob/master/docs/http.md#table-service

Can you explain the format of the .vrp file as if I create my own distance matrix I would have to create the .vrp file and I can't understand the structure of data in that file.

The VRP file is TSPLIB95 format and how the matrix is given is determined by the EDGE_WEIGHT_TYPE : EXPLICIT and EDGE_WEIGHT_FORMAT attributes. If you set the EDGE_WEIGHT_FORMAT : FULL_MATRIX then you can give the full matrix separated by white space, row by row. Note that then you can input asymmetric distances where it might take shorter time to get from A to B than from B to A. However, note that some VeRyPy CVRP algorithms might not fully support such asymmetric distances.

Another option for using the TSPLIB .vrp files is to import the algorithms as a Python module as shown below:

from classic_heuristics.gillet_miller_sweep import gillet_miller_init
from util import sol2routes

# Get this using e.g. OSRM distance matrix / table API. Be sure to have the school at index 0
full_distance_matrix = ...  
# Again, first is the school, and the rest are indexed as the distance matrix and show how 
# many pupils are expected to be picked up at each stop (only 4 stops in this example)
pupils_on_each_bus_stop = [0,12,5,3,1] 
capacity_of_the_bus = 42 # or whatever
minimize_number_of_buses = False # default optimizes total travel time / length

solution = gillet_miller_init(
    D=full_distance_matrix, 
    d=pupils_on_each_bus_stop, 
    C=capacity_of_the_bus,
    minimize_K=minimize_number_of_buses )

for route_idx, route in enumerate(sol2routes(solution)):
    print("Route #%d : %s"%(route_idx+1, route))

[AsharFatmi]

Hi,

Thanks for the prompt and detailed info, I will implement this and get back to you.

[yorak]

Hi,

one more thing that occurred to me when thinking about this. I spend too much time in my youth in a school bus. Although they should not, they were often overfull. The issue in modeling this is that in CVRP the C constraint is hard, i.e., it cannot be violated, even a little bit. Also, some classical algorithms allow minimizing K (the number of vehicles), but not all. Hence, after you have managed to compute the distance matrix and produce an initial solution with Gillett-Miller Sweep algorithm, I'd recommend looking into adapting the algorithm SG82-LR3OPT, where the Lagrangian relaxation makes the constraints soft and 3-opt* (inter-route 3-opt) finds the moves. Just be sure to add another relaxed constraint calculation for K here. If your problem instances are fairly large (lets say over 100 bus stops), you might also need to split the problem or improve the performance of LR3OPT.

I just wanted to share the insight.

-Jussi

[AsharFatmi]

Hi,

Thanks for the insight it is very helpful.

Although Bus overloading is what I am trying to remove so I don't think this is applicable to my case. If a bus can transport 30 students I only want 30 students assigned to that bus.

Having said that I do want variable C (capacity) for the busses as different sizes of busses can have different capacities of students.

For reference let us say we have 3 busses with [20, 30, 50] students capacity for the same school.

Any idea as to how to handle this case?

Thanks in Advance

[yorak]

Oh, what you are describing is the Heterogeneous Fleet VRP (HFVRP). By default, VeRyPy expects all of the vehicles to be the same, but modifying the code to support multiple different vehicle types is possible (although, unfortunately, not trivial). However, please create a separate branch for this. Modifying the code is algorithm specific. I guess the most straightforward approach to make the necessary changes would be to use the API as in my previous example but give the capacity C as a list of capacities, try to run the algorithm, and fix the places where exceptions are thrown.

If you are mostly interested in Gillet & Miller, then, looking at the code you need to:

1. add a vehicle_idx field to routedata.py
2. assign an increasing vehicle_idx for each new current_route in sweep.py (i.e., here and here)
3. If the C is a list (or similar, check with hasattr(C, '__getitem__')), instead of straight comparison, look up the vehicle capacity using the vehicle_idx in the sweep.py capacity constraint check.
4. Finally, you need to check all permutations of capacity constraints adding yet another for loop layer to the algorithm. You can use the built-in from itertools import permutations. It is usual that many vehicles have the same capacity, so at the beginning of the function it is advisable to prune the duplicates HFC_permutations = list(set(permutations(C)).

There are some open questions such as what to do if the algorithm runs out of vehicles during the sweep? After you get the basic Sweep to work as implemented insweep.py, you can look into the improvement procedure as implemented in the file gillet_miller_sweep.py. However, some of the algorithms might be easier to convert to HFVRP such as the LR3OPT I mentioned earlier.

[AsharFatmi]

I am not fixated on using the Gillet & Miller I would be happy to try the LR3OPT if it gives better results although running out of Vehicles is a serious limitation for HFVRP and kind of counter-intuitive as it would require back and forth human intervention in the system.

Before I implement the HFVRP I am trying to do the CVRP. Unfortunately, I was not able to build the LKH in windows so now I am working on Linux.

from classic_heuristics.gillet_miller_sweep import gillet_miller_init
from util import sol2routes

# Get this using e.g. OSRM distance matrix / table API. Be sure to have the school at index 0
full_distance_matrix = ...  
# Again, first is the school, and the rest are indexed as the distance matrix and show how 
# many pupils are expected to be picked up at each stop (only 4 stops in this example)
pupils_on_each_bus_stop = [0,12,5,3,1] 
capacity_of_the_bus = 42 # or whatever
minimize_number_of_buses = False # default optimizes total travel time / length

solution = gillet_miller_init(
    D=full_distance_matrix, 
    d=pupils_on_each_bus_stop, 
    C=capacity_of_the_bus,
    minimize_K=minimize_number_of_buses )

for route_idx, route in enumerate(sol2routes(solution)):
    print("Route #%d : %s"%(route_idx+1, route))

In the above code, the parameter "points" is missing. I added it as follows:-

points = [[24.419954, 54.4416255], [24.440368, 54.59438], [24.423042, 54.525564]]
solution = gillet_miller_init(
    points=points,
    D=full_distance_matrix,
    d=pupils_on_each_bus_stop,
    C=capacity_of_the_bus,
    minimize_K=minimize_number_of_buses)

I am getting the following error:-

Traceback (most recent call last):
  File "test.py", line 35, in <module>
    minimize_K=minimize_number_of_buses)
  File "/home/ubuntu/VeRyPy/classic_heuristics/gillet_miller_sweep.py", line 336, in gillet_miller_init
    intra_route_improvement=_improvement_callback)
  File "/home/ubuntu/VeRyPy/classic_heuristics/sweep.py", line 485, in sweep_init
    callback_data = pcds_callback(D,d,C,L,sweep)
  File "/home/ubuntu/VeRyPy/classic_heuristics/gillet_miller_sweep.py", line 88, in _pack_datastructures_callback
    return (len(D), D,produce_nn_list(D),d,C,L,node_to_pos,pos_to_node, avr)
  File "/home/ubuntu/VeRyPy/util.py", line 38, in produce_nn_list
    NN_D[i] = sorted(enumerate(D[i,:]), key=itemgetter(1))
TypeError: list indices must be integers or slices, not tuple```

How do I have to pass the points array?

[yorak]

Hi, the error suggests that the D (full_distance_matrix) is not a 2D Numpy array of size (1+n) x (1+n). Could you verify this?

During testing you could use "as the crow flies" distances. To compute these, you can use the haversine formula that comes with VeRyPy. An usage example below:

import numpy as np
from cvrp_io import _haversine 

points = [[24.419954, 54.4416255], [24.440368, 54.59438], [24.423042, 54.525564]]

n_incl_depot = len(points)
geodesic_D = np.zeros( (n_incl_depot, n_incl_depot) )

for i, p1 in enumerate(points):
    for j, p2 in enumerate(points):
        geodesic_D[i,j] = 0.0 if p1==p2 else _haversine(p1,p2)
    
print(geodesic_D)

[AsharFatmi]

Hi,

I passed the exact response that came from the OSRM.

full_distance_matrix = [
    [
        0,
        7562.1,
        23971.5
    ],
    [
        7545.1,
        0,
        18324.4
    ],
    [
        12008.1,
        17728.9,
        0
    ]
]

I converted into NumPy array and it worked thanks.

[AsharFatmi]

Hi,

So I ran the Gillet & Miller on a dummy (test) data, basically, one school with 2 stops having 4 and 5 students respectively.

School Loc = [24.419954, 54.4416255]
stops = [[24.459459, 54.321149], [24.440368, 54.59438]]

I ran the following code:-

from classic_heuristics.gillet_miller_sweep import gillet_miller_init
from util import sol2routes
import numpy

points = [[24.419954, 54.4416255], [
    24.459459, 54.321149], [24.440368, 54.59438]]
# Get this using e.g. OSRM distance matrix / table API. Be sure to have the school at index 0
full_distance_matrix = [
    [
        0,
        14997,
        17728.9
    ],
    [
        15737.5,
        0,
        32034.8
    ],
    [
        18324.4,
        33151.2,
        0
    ]
]
full_distance_matrix = numpy.array(full_distance_matrix)

# Again, first is the school, and the rest are indexed as the distance matrix and show how
# many pupils are expected to be picked up at each stop (only 4 stops in this example)
pupils_on_each_bus_stop = [0, 4, 5]
capacity_of_the_bus = 10  # or whatever
minimize_number_of_buses = False  # default optimizes total travel time / length

solution = gillet_miller_init(
    points=points,
    D=full_distance_matrix,
    d=pupils_on_each_bus_stop,
    C=capacity_of_the_bus,
    minimize_K=minimize_number_of_buses)

for route_idx, route in enumerate(sol2routes(solution)):
    print("Route #%d : %s" % (route_idx+1, route))

As you can see the bus capacity I entered 10. I got the following result:-
Route #1 : [0, 1, 2, 0]

I the above image you can see the blue marker is the school and the two stops labeled as 1 and 2.

Result:-

There was one route created containing both the stops as the stop1 + stop2 = 9 student < Bus Capacity (10).

Issue:-

Realistically the distance between the two stops is a lot. one bus should not be assigned to both the stops. I don't think any school in the world will have a bus go to both stops it's not practical. The total distance of the route being 68 KM.

Is there any way to handle this scenario?
I have not fixated on Gillet & Miller any other algorithm that can handle this scenario is also welcomed.

Thanks in advance

[yorak]

I'm very happy to hear that you have been able to resolve the initial issues and get the pipeline from coordinates to an optimized routes works. It is easy to build on that foundation. While it has already been used in several other research and R&D projects, I'm aware there is room for improvement in getting started with VeRyPy. Unfortunately, I've been, and currently am, working on entirely different topics, and have not been able to find larger blocks of time to work on VeRyPy. I've applied a research grant for that purpose so hopefully I can make some improvements in the near future.

Regarding the issue: that is typical. The way to solve it with VeRyPy is to use the maximum route duration/length constraint L in conjunction with the capacity constraint C. This allows you to avoid such "nonsensical" solutions.

A more granular way to model this would be to per delivery travel time constraints (that is, the travel time is measured for each pickup/delivery) or multicriteria optimization where the pupil travel time is taken into account in the cost function. In fact, this is usually what is desired in bus routing: Without such modified cost function, it may be that a pupil is picked up from a bus stop near the school, only to be transported around the city when other pupils are picked up to before delivering them to the school. In effect the first pupils would've been better off just walking to school.

If we consider VeRyPy, it might be that the simplest way to implement this would be to use the construction heuristics as is and offload the travel duration minimization to improvement heuristics (local search). The delta calculations ( e.g., for 2-opt* ) need to be adapted to also account for the transportation times of the students. To make this performant, you probably need to add some auxiliary data to the RouteData structure. Again, this is not entirely trivial modification which also means it is not CVRP anymore. Use of an specialized git branch for such extension is advisable.

Finally, please note that I pushed some changes to master make sure nobody uses the inter-route 3-opt by accident. You might consider pulling from my upstream repository.

[AsharFatmi]

I will try using the maximum route duration/length constraint L and check if it works for my case.

With respect to modeling based per delivery travel time constraints (that is, the travel time is measured for each pickup/delivery) this is basically Duration Matrix (1+n)*(1+n) where n is the number of stops and 1 is the school. [Correct me if I am wrong]

so we will result with a 2D NumPy array similar to the distance matrix for the cost function.

What should be the next step?

NOTE- I am not completely familiar with your code, so I am not sure where travel duration minimization is happening and how to factor in the new cost function.

[yorak]

Yes, I think that by using the L constraint you can already force the routes to different buses and, thus, to be more balanced. Also, each experiment you do increases your understanding of the problem and VeRyPy code base which you can leverage in choosing the next step.

Be warned that my intuition says that measuring the total transportation time cost of the pupils becomes a quite bit more complex as one has to consider all the next stops the bus route will take. Hence, the cost of traveling an edge is dependent on who the bus is transporting. Such context dependent travel costs are tricky to compute effectively and the travel time cost is a weighted sum over the selected edges (values in D) before offloading at the school. Drawing this out with pen and paper might help.

Another way I can come up for solving the problem of excessive transportation times is to allow the buses to begin their routes far away from the school. The drawback would be that you need to specify a number of starting locations for the bus and implement a specialized version of your chosen algorithm among the ones VeRyPy implements for such Open VRP variant. If you think about it, making driving when bus is empty next to free as opposed to making the travel costs expensive with pupils on board will create similar (but not the same!) routings to the OVRP alternative.

I see that the major issue for in adapting VeRyPy for the bus routing use is that the computation of cost and constraint checks are not done in an unified way between the different algorithms. By avoiding of using such abstraction, I was able to reach the main goal in writing the algorithms as close to the original descriptions (my eye was on replications, not extendability). Hence, I cannot point a single point where the modifications need to be made, sorry.

Perhaps some more specialized solver would be a better fit to you? If you would like to continue with VeRyPy, my proposal is that you first look what you can do with vanilla VeRyPy. Then, if the solutions given by, e.g., Gillet and Miller, Paessens savings and perhaps Stewart and Golden (LR3OPT) are usable to you (however, be sure to set both the capacity C and maximum route duration L!), you might be able to make an informed decision if investing the time to modify best of these algorithms to heterogeneous fleet and load dependent travel times.

[AsharFatmi]

Regarding using the L constraint, I didn't find any documentation hinting at values that L can take on.
Can you help me with the format of how to pass the values?

• How to pass duration/length as L.
• Unit of measure in which it is expecting the value.

Thanks in advance

[yorak]

Adding L was a bit of an afterthought. Some classical CVRPs had them and I decided to add the support for these, which required implementing the maximum route duration/length constraint. Hence, all construction algorithms (that is, the *_init(..) functions) have the optional parameter L. It uses the same unit as the distance matrix D and specifies the maximum duration/length for each route. Therefore, this should do it:

solution = gillet_miller_init(
    D=full_travel_time_matrix, 
    d=pupils_on_each_bus_stop, 
    C=capacity_of_the_bus,
    L=maximum_route_duration)

As you can see, I'd recommend optimizing using travel time instead of distance. By doing so, you can use stopping times. VeRyPy currently supports only one and same stopping time for all stops. However, adding support for a stopping time that is dependent on the number of pupils picked up from a bus stop is trivial. Just add half of the stopping time to outgoing and half to incoming edges of each bus stop. In practice you need to add those to each row/column of the distance/travel time matrix. Then, given that D is in time units and u is a Numpy array of stopping times:

U_rows = np.tile(u/2, (len(u), 1))
U_cols = np.tile(u/2, (len(u), 1)).T
D+=U_rows+U_cols
np.fill_diagonal(D, 0.0)

Now, the stopping times are embedded in the D. Please note that adding half of the stopping time to incoming edges and half to outcoming edges allows keeping symmetric problems symmetric.

[AsharFatmi]

Hi,

I ran the Gillet and Miller for my dataset.
PC config:-

CPU - Intel Core i7 - 8 cores
RAM- 32 GB
OS - Ubuntu 18.04

Copy of the code:-

import json
from classic_heuristics.gillet_miller_sweep import gillet_miller_init
from util import sol2routes
import numpy
from loader import Loader
from flask import jsonify


class testFunction:
    def __init__(self):
        pass

    def gilletMiller(self, runFromJson, input=None):
        runFromJson = True

        with open('request.json')as jsonFile:
            data = jsonFile.read()
        options = json.loads(data)
        print("capacity_of_the_bus = ", options["capacity_of_the_bus"])
        points = input["points"] if runFromJson != True else options["points"]

        full_distance_matrix = input["full_distance_matrix"] if runFromJson != True else options["full_distance_matrix"]

        full_distance_matrix = numpy.array(full_distance_matrix)

        pupils_on_each_bus_stop = input["pupils_on_each_bus_stop"] if runFromJson != True else options["pupils_on_each_bus_stop"]

        capacity_of_the_bus = input["capacity_of_the_bus"] if runFromJson != True else options["capacity_of_the_bus"]

        minimize_number_of_buses = False

        loader = Loader("Running Heuristics...",
                        "Heuristics Done!", 0.1).start()
        solution = gillet_miller_init(
            points=points,
            D=full_distance_matrix,
            d=pupils_on_each_bus_stop,
            C=capacity_of_the_bus,
            minimize_K=minimize_number_of_buses)
        loader.stop()

        ROUTES = []

        loader.__init__("Preparing Response", "Response ready!!", timeout=0.1)
        loader.start()

        for route_idx, route in enumerate(sol2routes(solution)):
            print("Route #%d : %s" % (route_idx+1, route))
            ROUTES.append({route_idx+1: route})

        loader.stop()

        with open("response.json", "w") as outfile:
            json.dump(
                jsonify({"success": True, "data": {"routes": ROUTES}}), outfile)

        return ROUTES


if __name__ == "__main__":
    testObj = testFunction()
    ROUTES = testObj.gilletMiller(runFromJson=True, input=None)

Issue:-

• The code has been running for 20 hrs now and it is still not complete. Is there any benchmark for the Time Constraint for the process? As you can see I did not yet implement the " L " yet running just the basic algorithm for the preliminary results.
• The data set contains 297 entries including all the stops and the school at index 0.

Request.json File:-
request.zip

[yorak]

That seems excessive. As you can see, from the page 103 of the VeRyPy manuscript Gillet & Miller should be able to solve problems up to around 700 customers within 1 hour. In fact, if we consider average case, under 10 minutes (Figure 22: Two phase algorithms, C) even if the comprehensive start customer for the sweep is enabled (as it is by default).

Perhaps it is the excessive logging? It is recommended to run VeRyPy with the python -O test.py, which will entirely turn off the debug output. Also, the verbosity of the logging can be changed using the shared_cli.set_logger_level(level) function.

I'll check this on my machine.

[yorak]

I enabled debug and it seems the algorithms gets stuck right away.

Do a sweep from position 112 (n1) by steps of 1
DEBUG:Sweep node order [1, 152, 113, 169, 275, 70, 51, 130, 184, 23, 16, 9, 49, 128, 28, 280, 143, 76, 121, 129, 98, 78,....
DEBUG:Considering n1 between rays -1.0206, -1.0193
DEBUG:Added n1 to the route set
DEBUG:Considering n152 between rays -1.0193, -1.0171
DEBUG:Added n152 to the route set
DEBUG:Considering n113 between rays -1.0171, -1.0136
DEBUG:Added n113 to the route set
...
DEBUG:Considering n129 between rays -0.9740, -0.9725

Where it gets stuck is calling the routing_algo. Hence, in TSP optimization of all places. This happens using the built-in TSP 3-opt. Even after testing this with thousands of classical problem instances, it seems you have hit some pathological case in TSP optimization.

[yorak]

Ok, found the reason. It seems the internal TSP 3-opt local search does not work correctly on asymmetric distances. I'd assume the delta calculation does not work correctly, which causes the 3-opt TSP to get stuck into infinite improvement loop. Personally, I have tested VeRyPy only with symmetric distances. I know some people have used it also for asymmetric instances, but I have not received any pull requests.

At the moment I do not have the time to dig deeper, but, fortunately, there are few options for you to proceed.

1. use symmetric distances:

for i in range(full_distance_matrix.shape[0]):
            for j in range(i, full_distance_matrix.shape[1]):
                full_distance_matrix[i,j] = full_distance_matrix[j,i]

2. Fix the current 3-opt local search implementation. When a segment is reversed, the cost of the segment is changed and needs to be recalculated. This also applies to 2-opt. Requires changes to here, here, here and all others moves where the segment is reversed.
3. Implement or use another 3-opt or entirely different TSP algorithm that supports the VeRyPy tsp_solvers API. For example, IIRC LKH and ACOTSP both support asymmetric distances.

It is possible that I'm able to look into this later, but currently I'm quite busy with other work, so it might take even few weeks.

[AsharFatmi]

Hi,
Regarding the L constraint.
My distance matrix is calculated in meters. Assuming that I am implementing L without the stopping times, I want to restrict the created routes to 1hr [ 60 min ]. what should I pass in the variable L.

- L=60  #minutes
- L=1    #hour
- L=3600  #seconds

[yorak]

Unfortunately, VeRyPy does not keep track of the travel distance and time separately and I'm not aware of other solvers that would support both time and distance based constraints / optimization targets. Please note that you could query both the distance and duration matrices from OSRM, but use the duration matrix for optimization and distance matrix for results. Usually the fastest route is also the shortest. Hence, my recommendation is to use the duration matrix.

[AsharFatmi]

so L constraint should be a NumPy array of the symmetric Duration Matrix...

[yorak]

No, :) D should be a symmetric duration matrix given as a NumPy 2D array and L just a scalar of the same unit (i.e. in seconds if the values of D are given in seconds).