The Traveling Salesman Problem (TSP) asks which order a set of stops should be visited in to minimize trip cost. It can be solved entirely inside PostgreSQL: PostGIS supplies the spatial data handling, and pgRouting — a PostgreSQL extension that adds geospatial routing functionality on top of PostGIS — supplies the routing algorithms. Documentation lives at https://pgrouting.org/.

The example reuses the OpenStreetMap files and datasets from the earlier catchment-area work around Hamburg, so those processing steps are a prerequisite.



Figure 1 Area of interest around Hamburg / Figure 2 Area of interest around Hamburg, Zoom

Building a routable graph

To solve TSP with pgRouting, the OSM data first has to be turned into a graph. Tools for this include osm2pgrouting and osm2po; the latter was chosen here because of limited memory on the development machine. It also requires an intermediate table, tspPoints, holding the major cities and airports that fall inside the Hamburg catchment area (DDL in the annex).

After downloading the latest osm2po release from https://osm2po.de/releases/, the graph is emitted as a SQL file by uncommenting one line in osm2po.config:

postp.0.class = de.cm.osm2po.plugins.postp.PgRoutingWriter

Java

1

2

java -Xmx8g -jar osm2po-core-x.x.x-signed.jar workDir=/data/de-graph

prefix=de tileSize=x /data/germany-latest.osm.pbf

The generated SQL file is then imported into PostGIS:

PgSQL

1

psql -U username -d dbname -q -f /data/de-graph/de_2po_4pgr.sql

The resulting graph table carries geometry for edges only (DDL in the annex). Source and target nodes are generated with pgr_createverticestable:

PgSQL

1

2

select pgr_createverticestable('de_2po_4pgr', 'geom_way', 'source',

'target');

Three tables make up the stacked layer setup on top of OpenStreetMap:

  • tspPoints — the intermediate stops
  • de_2po_4pgr — edge geometries of the graph
  • de_2po_4pgr_vertices — node geometries of the graph

Figure 3 Intermediate trip stops


Figure 4 Edges around Hamburg


Figure 5 Edges + nodes around Hamburg

Running the TSP query

With the layers in place, pgRouting is asked to solve the problem:

PgSQL

1

2

3

4

5

6

7

8

9

10

11

12

13

14

SELECT seq,node, cost, agg_cost

FROM pgr_TSP(

        $$

    SELECT * FROM pgr_dijkstraCostMatrix(

        'SELECT id, source, target, cost, reverse_cost FROM de_2po_4pgr',

        (with nearestVertices as(

SELECT a.id from tsppoints,

    lateral (

        select id,the_geom from osmtile_germany.de_2po_4pgr_vertices_pgr

        order by

osmtile_germany.de_2po_4pgr_vertices_pgr.the_geom <-> osmtile_germany.tsppoints.geom limit 1

        ) a)

select array_agg(id) from nearestVertices), directed := false)

    $$, start_id := 591341);

The result is the optimal order of intermediate stops — the edges next to the cities — that forms the trip. That optimal route depends on the cost matrix, and therefore on the cost function used to weight the edges; here the function comes from osm2po and is expressed as length/kmh.

To turn the sequence into a list of city names, the nodes returned by pgr_TSP are joined back to tspPoints (query extension in the annex). The resulting optimal route runs from Hamburg Airport to Kiel, Lübeck, Rostock and Hamburg, then back to the airport.

Figure 6 Optimal order of stops

Further cost functions can be tried, and detailed routes between cities can be computed with additional pgRouting functions.

Annex

PgSQL

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

create table if not exists osmtile_germany.de_2po_4pgr

(

      id integer not null

         constraint pkey_de_2po_4pgr

           primary key,

      osm_id bigint,

      osm_name varchar,

      osm_meta varchar,

      osm_source_id bigint,

      osm_target_id bigint,

      clazz integer,

      flags integer,

      source integer,

      target integer,

      km double precision,

      kmh integer,

      cost double precision,

      reverse_cost double precision,

      x1 double precision,

      y1 double precision,

      x2 double precision,

      y2 double precision,

      geom_way geometry(LineString,4326)

);

create index if not exists idx_de_2po_4pgr_source

      on osmtile_germany.de_2po_4pgr (source);

create index if not exists idx_de_2po_4pgr_target

      on osmtile_germany.de_2po_4pgr (target);

create index if not exists idx_de_2po_geom

     on osmtile_germany.de_2po_4pgr using gist (geom_way);

PgSQL

1

2

3

4

5

6

7

8

create table osmtile_germany.tspPoints

(

id   serial primary key,

name text,

geom geometry(Point, 4326)

);

create index idx_tspPoints on osmtile_germany.tspPoints using gist (geom);

PgSQL

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

SELECT seq, node, a.name, a.geom

FROM

     pgr_TSP(

       $

       SELECT *

       FROM

         pgr_dijkstraCostMatrix(

           'SELECT id, source, target, cost, reverse_cost FROM de_2po_4pgr',

           (

             with nearestVertices as(

             SELECT a.id

             from

               tsppoints,

               lateral (

                 select id, the_geom

                 from

                   osmtile_germany.de_2po_4pgr_vertices_pgr

                 order by

                   osmtile_germany.de_2po_4pgr_vertices_pgr.the_geom <-> osmtile_germany.tsppoints.geom

                 limit 1

               ) a

             )

             select array_agg(id)

             from

               nearestVertices

           ),

           directed := false

          )$,

          start_id := 591341

      ),

     de_2po_4pgr_vertices_pgr,

     lateral (

        select name, geom

        from

          osmtile_germany.tsppoints

        order by

          osmtile_germany.tsppoints.geom <-> the_geom

        limit 1

      ) a

where node = id