Simple Implementation Of Travelling Salesman Problem

UPDATE: This function is now a part of the Oracle Object implementation.

The Travelling Salesman is a well known problem in graph theory:

The Travelling Salesman Problem describes a salesman who must travel between N cities.
The order in which he does so is something he does not care about, as long as he visits each one during his trip, and finishes where he was at first.
Each city is connected to other close by cities, or nodes, by airplanes, or by road or railway.
Each of those links between the cities has one or more weights (or the cost) attached.
The cost describes how “difficult” it is to traverse this edge on the graph, and may be given, for example, by the cost of an airplane ticket or train ticket,
or perhaps by the length of the edge, or time required to complete the traversal.
The salesman wants to keep both the travel costs, as well as the distance he travels as low as possible.

While there are many solutions out there in Java and other languages, I thought I would create an implementation using PL/SQL and Oracle Spatial alone. I do not profess that this is a true, or full, implementation of the problem, but it is an example of what could be done even in PL/SQL.

My implementation has two parts:

1. A function that is compiled against specific tables and which uses the nearest neighbour SDO_NN operator.
2. A separate function that runs against collections of geometries independently of the source of the data.

This article describes the first implementation. This implementation requires:

• A table of points coded only in SDO_POINT_TYPE;
• Points can be 2D or 3D;
• Table has unique integer column (either primary key or unique constraint/index);
• Table is spatially indexed;

One can create a sample dataset that implements these requirements as follows.

 DROP   TABLE ProjPoint2D  PURGE;
 CREATE TABLE ProjPoint2D ( id INTEGER, geom mdsys.sdo_geometry );
     INSERT INTO ProjPoint2D (id, geom)
 WITH mGeom AS (
 SELECT
 SDO_GEOMETRY(2005,NULL,NULL,
 SDO_ELEM_INFO_ARRAY(1,1,20),
 SDO_ORDINATE_ARRAY(358443.819,5407099.906, 358716.043,5407087.437, 358651.624,5407224.588, 358680.717,5407399.144,
  358952.941,5407444.861, 359314.521,5407415.768, 359370.628,5407249.525, 359349.847,5407137.31,
  359108.794,5407089.515, 358930.082,5407241.213, 358826.18,5407336.803, 358819.946,5407463.563,
  359056.843,5407552.919, 359283.35,5407615.261, 359075.545,5407523.827, 358485.38,5407735.787,
  358393.946,5407552.919, 358312.903,5407311.866, 358319.137,5407087.437, 358360.698,5406991.847)) AS mPoint
  FROM dual
 )
 SELECT t.id,sdo_geometry(2001,a.mPoint.sdo_srid,sdo_point_type(t.x,t.y,t.z),NULL,NULL) AS point
   FROM mGeom a,
        TABLE(sdo_util.getVertices(a.mPoint)) t;
 COMMIT;
 SELECT 'Inserted '||COUNT(*)||' records into ProjPoint2D' FROM ProjPoint2D;
 DELETE FROM USER_SDO_GEOM_METADATA WHERE TABLE_NAME = 'PROJPOINT2D';
 COMMIT;
 DECLARE
   v_shape mdsys.sdo_geometry;
   v_srid  INTEGER;
 BEGIN
   SELECT SDO_AGGR_MBR(a.geom), a.geom.sdo_srid INTO v_shape, v_srid FROM ProjPoint2D a GROUP BY a.geom.sdo_srid ;
   INSERT INTO USER_SDO_GEOM_METADATA(TABLE_NAME,COLUMN_NAME,DIMINFO,SRID)
     VALUES ('PROJPOINT2D','GEOM',MDSYS.SDO_DIM_ARRAY(
     MDSYS.SDO_DIM_ELEMENT('X',V_SHAPE.SDO_ORDINATES(1),V_SHAPE.SDO_ORDINATES(2),0.005),
     MDSYS.SDO_DIM_ELEMENT('Y',V_SHAPE.SDO_ORDINATES(3),V_SHAPE.SDO_ORDINATES(4),0.005)),
     v_srid);
 END;
 /
 SHOW ERRORS
 COMMIT;
 CREATE INDEX ProjPoint2D_GEOM ON ProjPoint2D(geom)
        INDEXTYPE IS mdsys.spatial_index
        parameters('sdo_indx_dims=2, layer_gtype=point');

The function has the following interface:

 CREATE OR REPLACE
 FUNCTION ST_TravellingSalesman(
   p_table_name       IN varchar2,
   p_table_gid_name   IN varchar2,
   p_geom_column_name IN varchar2,
   p_table_gid_value  IN NUMBER,
   p_unit             IN varchar2 DEFAULT NULL
 )
 RETURN sdo_geometry Deterministic
 AS
  ... etc ...
 BEGIN
   -- Algorithm starts with point (Node) identified by p_table_gid_value of column in table p_table_name called p_table_gid_name.
   -- Add first ID to list of nodes that have already been processed.
   -- Iterate over all points in table
   ---- For current node (gid) find nearest point (node) in table to it which is not in already processed node list.
   ---- Add found node XY(Z) values to return linestring.
   ---- Add found node ID to list of nodes that have already been processed
   RETURN path_as_linestring;
 END ST_TravellingSalesman;
 /
 SHOW errors

An example of how to call the function with our sample data is as follows.

 SELECT st_travellingsalesman(p_table_name=>'PROJPOINT2D', -
                              p_table_gid_name=>'ID', -
                              p_geom_column_name =>'GEOM', -
                              p_table_gid_value=>1, -
                              p_unit=>NULL) AS path-
   FROM dual;
 -- Now let's test it
 PATH
 ----
 SDO_GEOMETRY(2002,NULL,NULL,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(358443.819,5407099.906, 358319.137,5407087.437, 358360.698,5406991.847, 358312.903,5407311.866, 358393.946,5407552.919, 358485.38,5407735.787, 358680.717,5407399.144, 358819.946,5407463.563, 358826.18,5407336.803, 358930.082,5407241.213, 358952.941,5407444.861, 359075.545,5407523.827, 359056.843,5407552.919, 359283.35,5407615.261, 359314.521,5407415.768, 359370.628,5407249.525, 359349.847,5407137.31, 359108.794,5407089.515, 358716.043,5407087.437, 358651.624,5407224.588))

This looks like the following:

There is also a 2.5D linestring example which shows a different solution for the same XY points as it takes in to account the slope length between each of the points.

Contact me for the actual source code.

I hope this article helps someone out there.

Contact me for the actual source code.