This article will show a (rare!) practical problem I worked on, namely automating the discovery of the order in which to visit a list of addresses, all in the same city, Milan, and moving by car, in order to reduce the travel time.

This problem is called Traveling Salesman Problem and is a famous one in computer science because of its exponential complexity. In reality, I do not need the best solution, just a good one, and the complexity truly comes from the practicality of extracting the graph and finding the addresses.

I find it interesting because it’s a problem that is connected to a similar problem I worked on and is something that delivery companies do every day.

In my specific case I am in a volunteering association that needs to visit a set of addresses by car (the addresses correspond to homeless people to visit). The addresses are all within the city proper, and the operation happens late during the day when the traffic is negligible.

The output of this application must be easy to access during the “trip”, when of course I have no access to a computer nor time, screen space and attention to deal with anything complex, so I opted to generate a small spreadsheet with simply the list of addresses in the proper order, and a separate set of GeoJSON and TSV files that are used for troubleshooting before the trip happens, to be examined comfortably from a PC.

The overall process

The whole process looks like this:

  • an OSM extract is processed to produce a graph of the road network (including speed limits, street type, traffic lights and surfaces) and addresses
  • the addresses to visit are matched to street names and mapped to coordinates
  • the addresses are split into “zones”, this is a detail of how the association operates but essentially means I generate multiple lists to split the area to cover among different cars at the same time
  • a special starting point (the HQ, called “sede”) is also provided
  • the fastest route between every combination of points is calculated, considering the travel speed and the number of traffic lights on each segment
  • a two-step optimization finds a good enough order to visit the points starting from HQ, using a greedy approach and then swapping stops to try and improve the total time
  • the output is produced as a minimal TSV ready to copy paste in a spreadsheet and be consulted by phone and a geoJSON to show the path and stops

Extracting the road network

This is the most interesting part for me, since it can be adapted to other types of navigation like bike or foot with arbitrary constraints (e.g. walk in well-lit streets or trying to be near water taps) and can be used to calculate isochrones, which in fact is an extra function provided by this script.

First of all, I download the OSM extract from Geofabrik and reduce it to the area of Milan proper with Osmium. The resulting file is quite small, less than 20 MB in my case, so even for a large city this process can be done on a common laptop.

After that, I use duckdb with the osmium extension, enabling it to read PBF files directly. Another possibility is to use quackosm to process the PBF into a geoparquet:

CREATE TABLE raw AS
SELECT id, kind, type, tags, geometry, refs, ref_roles, ref_types
FROM {pbf_file}

DuckDB is quite excellent at handling maps/structs which is necessary to process tags. Each road is filtered and stored as a line geometry with the speed details.

CREATE TABLE street AS
WITH speed_map AS (
    select MAP {
                'motorway': 130, 'trunk': 110, 'primary': 90, 'secondary': 70,
                'tertiary': 50, 'unclassified': 50, 'residential': 30,
                'motorway_link': 70, 'trunk_link': 70, 'primary_link': 50,
                'secondary_link': 50, 'tertiary_link': 50, 'living_street': 20
            } as m
),
surface_cap_map AS (
    select MAP {
                'cobblestone': 30, 'unpaved': 30, 'compacted': 30,
                'dirt': 20, 'gravel': 20, 'sand': 10
            } as m
)
SELECT
    id,
    tags['name'] AS name,
    geometry,
    COALESCE(
        try_cast(string_split(tags['maxspeed'], ' ')[1] AS integer),
        least(
            surface_cap_map.m[tags['surface']],
            speed_map.m[tags['highway']]
        )
    ) AS speed_kmh,
    coalesce(tags['oneway'], 'no') in ('yes', 'true', '1') AS is_oneway,
    -- this seems unused really
    -- coalesce(tags['oneway'], 'no') in ('-1') AS is_reversed,
    coalesce(tags['bridge'], 'no') in ('yes', 'true', '1') AS is_bridge,
    tags

FROM raw, speed_map, surface_cap_map
WHERE
    kind='line'
AND tags['highway'] IN map_keys(speed_map.m)

The next step is segmentize: if a road has a long segment the graph built from it will have two distanct vertices and the system will be unable to understand that an origin or destination point can lie between them. We need to split long edges into smaller ones by adding intermediate waypoints.

PostGIS/geos has a ST_segmentize function that does it, but the DuckDB spatial extension does not (I asked about it, so it may be in the future). I used Shapely to read all the geometries and segmentize them outside SQL, to write them back. Each waypoint gets a sequence number relative to the street geometry, and the order in which they grow is the right one in case of one-way streets.

The osmium library does not provide node ids when they have no tags, so I invent fake ids to assign them to have an identifier for each node in the graph.

From this data it is trivial to produce a vertex-edge graph, considering the sequential id of each node and whether they are one-way streets or not:

CREATE TABLE edge AS
SELECT
    p0.node_id              AS from_id,
    p1.node_id              AS to_id,
    st_distance(
        ST_Transform(p0.point_geom, 'EPSG:4326','EPSG:3857'),
        ST_Transform(p1.point_geom, 'EPSG:4326','EPSG:3857')
    ) / (s.speed_kmh * 0.2778)
    +
    (case when p0.has_traffic_light then 15.0 else 0.0 end)
                            AS time_s,

    p0.point_geom           AS from_geom,
    p1.point_geom           AS to_geom,
    st_asgeojson(p0.point_geom) as fj,
    st_asgeojson(p1.point_geom) as tj

FROM raw_waypoint p0
JOIN raw_waypoint p1
    ON p0.street_id = p1.street_id
    AND (
        (p0.idx = p1.idx - 1)
        OR
        (p0.idx = p1.idx + 1 AND NOT p0.is_oneway)
    )
JOIN street s
    ON p0.street_id = s.id

Now I have an oriented graph with edges and nodes, and metadata about speed and traffic lights associated to nodes and lines.

Geocoding

The places to visit are expressed in natural language and in Italian, for example “Via Crespi” or “Piazzale Bacone”, sometimes with a street number.

The task of mapping the names to coordinates is called geocoding and can be done with APIs, free or paid. In this case due to the small size of the area I went for a fuzzy string match based on the rapidfuzz extension of DuckDB.

For each address I extract the type (“via”, “piazza”, etc.) and the number, if any. Then, the type is matched as is and the rest of the name is matched using a fuzzy match for substrings. This means that “Via Crespi” will match “Via Benigno Crespi” since the match is based on substrings, and being fuzzy typos are not an issue.

If there’s also a street number, it is matched too, otherwise the centroid of the street is used.

SELECT ST_Y(ST_Centroid(geometry)) as lat, ST_X(ST_Centroid(geometry)) as lng
FROM raw
WHERE
            rapidfuzz_partial_ratio(LOWER(tags['addr:street']), $search_name) > 90
    AND LOWER(tags['addr:street']) LIKE LOWER('%' || $qualifier || '%')
    AND (tags['addr:housenumber'] = $number OR $number IS NULL)
    AND ST_Intersects(ST_Centroid(geometry), ST_GeomFromGeoJSON($zone_gj))
ORDER BY rapidfuzz_partial_ratio(LOWER(tags['addr:street']), $search_name) DESC
LIMIT 1

A relevant detail here is that the association operates in a sub-area of the city, zone_gj here, provided as a GeoJSON. This is necessary to reduce ambiguity of the match. For example, in Milan there are four “Via Crespi” ( Via Daniele Crespi, Via Gaetano Crespi, Via Pietro Crespi and Via Benigno Crespi), and this sub-zone must be provided to disambiguate them.

The system writes all the unmatched addresses to a file and stops if there are any, so the user must provide an additional GeoJSON with the overrides for those matches. An address in natural language or for a place that is not a street can be matched like this.

Finding the shortest path

The TSP problem requires first to calculate the distance between every pair of points (in this case really the travel time instead of the travel distance).

This is easily done by iterating over all pairs and applying Dijkstra. I first dump the whole graph to a dictionary and run the shortest-path on it, since the size of this data is negligible, without any optimization. I am not even using numpy to represent the ids, just Python int!

The TSP then is “solved” using a two-step solution.

First, we start from the HQ point (“sede”) and connect to the quickest to reach, then the quickest to reach from it, and so on until we have an initial possible solution. This is the “greedy” step.

Then, we try to swap consecutive stops and check if the total time improves. If any swap leads to an improvement we pass through the list again trying more swaps. This is repeated until no improvement is found.

Exporting the result

This step is critical, because the path could make no sense (e.g. if there are issues with the map data) and the solution is useful only as long as the data can be checked and actually used when needed.

I found GeoJSON as the best format to perform this check: it’s trivial to produce from Python or DuckDB (with ST_AsGeoJSON in the spatial extension), and can be visualized with a tool like geojson.io or integrated in custom MapLibre maps.

It turns out using geojson.io allows the user to easily modify or create data in this format, making it effective even as an input format to provide regions and addresses.

The final output of the tool however is the list of stops, and for that I opted for a simple TSV, easy to copy-paste into google docs and see on the phone.

I did search for ways to export the stops to Google Maps or OsmAnd to directly navigate but found that this solution does not work well when the order is changed a bit manually at the last moment, for example because a stop is removed or added.

An additional advantage of using a simple spreadsheet is that it is possible to mark with a color the stops already done (even out of order) and write notes.

The code

The code is available here