I've been wondering about this for a while - but never got around to asking - until I saw faleichik's incredible answer to Mathematica Minecraft.

I did some work at a not-for-profit a little while ago, which participated in a program called Meals On Wheels - wherein I helped to deliver food to those who are homebound. Myself and a partner would deliver meals to a daily list of addresses, sorted by name; not the most efficient route. The time to complete a route varied by as much as two hours because of it.

FindShortestTour, however, is only one (possible) part of the solution. I'd need to parse addresses, read map data, and handle the slight irregularities of NYC streets.

How can I use FindShortestTour to find the optimal meal delivery route, given a list of addresses?

  • 3
    $\begingroup$ "the slight irregularities of NYC streets" - so, not Manhattan? ;) $\endgroup$ Jun 10, 2013 at 6:31
  • 1
    $\begingroup$ @J. M. Well, Manhattan + one ways :) $\endgroup$
    – rm -rf
    Jun 10, 2013 at 6:34
  • 3
    $\begingroup$ I'd suggest that, for this practical application, it isn't necessary to calculate a route that is guaranteed to be the shortest. It'll be easier to incorporate various hard and soft constraints if you instead take a more heuristic approach, e.g. simulated annealing. This will still return a solution very close to optimal, but will take much less time to do so. $\endgroup$ Jun 10, 2013 at 8:14
  • 2
    $\begingroup$ Until they switched to Google Maps a few years ago, the Champaign-Urbana MoW gave route "directions" that, seemingly, had never been tested. Optimal, perhaps. But invariably completely wrong. $\endgroup$ Jun 10, 2013 at 12:16
  • 2
    $\begingroup$ @J. M. I don't think many birds tried them out. The only ones following us on those desolate roads were vultures.. $\endgroup$ Jun 10, 2013 at 17:26

4 Answers 4


When I did a similar job, it wasn't really the distance that made the difference. It would sometimes be quicker to drive 30 miles into the country for a drop than to try to cross the city at rush hour or during the school run. So it might be worth trying to work out the times between each drop and then find the route that has the shortest total time, rather than do it by distances on a map.

Here, rather clumsily derived, is a list of the different routes between, say, 8 places:

numberOfPlaces = 8;
edges = DeleteCases[
  Map[#[[1]] \[DirectedEdge] #[[2]] &, 
   Tuples[Range[1, numberOfPlaces], {2}]], 
  x_ \[DirectedEdge] x_ /; x == x]

and we can assign times in minutes for travelling between each place, here just randomly. This is a bit of work initially, but most drivers would be able to estimate the times quite quickly.

randomweights = RandomInteger[{0, 45}, Length[edges]]

{17, 38, 25, 17, 9, 42, 7, 42, 11, 14, 40, 17, 4, 23, 4, ... 40}

Combine these with the edges:

edgeWeightRules = 
   Rule, {edges, randomweights}] /. 
    (1 \[DirectedEdge] 4 -> n_) -> 1 \[DirectedEdge] 4 -> 200

I've set the value for traveling between places 1 and 4 to 200, for use later on. Next, draw a graph:

 cg = Graph[edges, 
  VertexLabels -> "Name",
  VertexStyle -> Red,
  EdgeWeight -> edgeWeightRules, 
  EdgeLabels -> edgeWeightRules, 
  ImagePadding -> 10]

simple graph

This is not a physical map, but a picture showing the different travel times between different places. Then create a weighted adjacency graph:

wam = Normal[WeightedAdjacencyMatrix[cg]] /. 0 -> Infinity;
wag = WeightedAdjacencyGraph[wam, 
  EdgeShapeFunction -> 
   GraphElementData[{"DotLine", "ArrowSize" -> .01}],
  VertexLabels -> "Name", 
  EdgeWeight -> edgeWeightRules, 
  EdgeLabels -> edgeWeightRules, 
  EdgeLabelStyle -> Directive[Gray, Italic, 12],
  VertexLabelStyle -> Directive[Black, 16],
  EdgeStyle -> Directive[LightGray],
  ImagePadding -> 20]

weighted adjacency graph

Apparently, a Hamiltonian cycle is not a two-wheeled delivery device, but a path visting every vertex exactly once:

FindHamiltonianCycle[wag, All]

{{1 -> 2, 2 -> 3, 3 -> 4, 4 -> 5, 5 -> 6, 6 -> 7, 7 -> 1}

so we should find all of them, and then sort them by the total time each one takes:

routes = Sort[{Total[# /. edgeWeightRules], #} & /@ 
   FindHamiltonianCycle[wag, All]]

This returns a list of all possible cycles, with the quickest one first. So the quickest route that turns up at all locations once takes 54 minutes (that's driving past and throwing the items out of the window...:).

{54, {1 -> 6, 6 -> 4, 4 -> 3, 3 -> 5, 5 -> 7, 7 -> 2, 2 -> 1}}

Then it's possible to investigate the different routes, using a Manipulate to step through the Hamiltonian cycles:

   Row[{"length of ", Last[routes[[n]]], " is ", First@routes[[n]]}],
     DirectedEdges -> True],
    GraphHighlightStyle -> "Thick",
    ImageSize -> 400]
   }], {n, 1, Length[routes], 1}]

manipulate for exploring the routes

The artificially slowed-down route shows up as being part of the longest cycles.

At least this approach avoids the problems of trying to find routes through the one-way systems of cities.

  • 1
    $\begingroup$ Also, there's the left turn problem. While it can take longer than if you make both left and right turns, if you make only right turns, you'll save gas. $\endgroup$
    – rcollyer
    Jun 10, 2013 at 15:15
  • $\begingroup$ @rcollyer I remember seeing that a few years ago. I hope the OP can save as much gas us UPS can. :) $\endgroup$
    – cormullion
    Jun 11, 2013 at 7:41
  • $\begingroup$ And, the mythbusters verified it. :) $\endgroup$
    – rcollyer
    Jun 11, 2013 at 12:39
  • $\begingroup$ WHOA. This answer is incredible. Creating routes basically answered the question, but the Manipulate really goes an extra mile. (do I sense a dash of british understatement?) $\endgroup$ Jun 23, 2013 at 10:37
  • $\begingroup$ :) you might... In fact, I found a better answer to this recently, which would have helped me. $\endgroup$
    – cormullion
    Jun 23, 2013 at 10:43

The FindShortestTour outputs leave a lot to be desired, in fact, none of the solutions shown for the coprime examples in the help browser are optimal. Some take half an hour to compute, and if you increase the number of nodes to a few dozen or hundred, it runs for a few days without providing good solutions. I assume you have plenty of needy people in NYC, not just some 20 or 30.

You could use the add-on package JVMTools, which has several highly advanced TSP functions. The outputs are generally optimal for a few hundred locations, and they are usually 1 - 3% longer than optimal for a few thousand locations, within seconds or minutes of running time. JVMTools was written for industrial-strength applications.

These TSP functions of JVMTools make maximum use of multi-threading, it will autosubmit as many competing calculations as you have cores on your machine, and you can always request more (threads = calculations), and then request output of either all solutions, or only the best. This is the TSP intro page and these are the TSP power examples.

You can assign arbitrary (non-negative!) distances, so you could use the Euclidean norm or the Manhattan-norm or anything else (Broadway doesn't exactly obey the Manhattan norm). You can "merge" points by setting a distance of 0 or model one-way streets by using a distance of infinity.

Free trial version

  • 1
    $\begingroup$ There are needy people in Broadway? Whodathunkit... ;) $\endgroup$ Jun 10, 2013 at 16:00
  • $\begingroup$ Psst - Columbia University is on Broadway :O $\endgroup$ Jun 10, 2013 at 22:12
  • $\begingroup$ ooh - In some sleep-deprived fog, I read it as nerdy. This just sounds strange. :eek: $\endgroup$ Jun 23, 2013 at 9:15

As of version 10.3 we can find an optimal tour with respect to driving distance and driving time!

Here's some random locations:

schools = EntityList[EntityClass["PublicSchool", {"City", {"Champaign", "Illinois", "UnitedStates"}}]][[1 ;; 8]];

We can minimize the distance driven with FindShortestTour + TravelDistance (Note, this is slow due to many API calls to an external service):

{dist, tour1} = FindShortestTour[schools, DistanceFunction -> TravelDistance]
{Quantity[10.584, "Miles"], {1, 8, 5, 4, 6, 3, 2, 7, 1}}

We can do the same for time with TravelTime:

travelTimeHours = N@UnitConvert[TravelTime[##], "Hours"] &;

{time, tour2} = FindShortestTour[schools, DistanceFunction -> travelTimeHours]
{Quantity[0.334722, "Hours"], {1, 8, 5, 7, 2, 3, 6, 4, 1}}

Notice that both tours differ slightly. Also observe they differ a lot from straight line distance:

{Quantity[7.23617, "Miles"], {1, 8, 6, 3, 2, 7, 4, 5, 1}}

We can visualize the optimal tours with GeoGraphics + TravelDirections:

mindist = TravelDirections[schools[[tour1]]];
mintime = TravelDirections[schools[[tour2]]];

  GeoGraphics[{{Black, Thick, Arrow[mindist]}, {GeoMarker[schools]}}, PlotLabel -> Row[{"Minimal distance of ", dist}]],
  GeoGraphics[{{Black, Thick, Arrow[mintime]}, {GeoMarker[schools]}}, PlotLabel -> Row[{"Minimal time of ", time}]]

enter image description here


Given n locations, we can limit the number of API calls from n^2 to 2 by using TravelDistanceList + FindPostmanTour.

Here's a function that computes the travel distance between each pair in bulk:

TravelDistanceAssociation[locs_List] /; OddQ[Length[locs]] :=
  Block[{n, inds, keys1, dists1, keys2, dists2},
    n = Length[locs];
    inds = FindPostmanTour[CompleteGraph[n]][[1, All, 1]];
    AppendTo[inds, First[inds]];
    keys1 = locs[[#]]& /@ Partition[inds, 2, 1];
    dists1 = Normal[TravelDistanceList[locs[[inds]]]];
    keys2 = locs[[#]]& /@ Partition[Reverse@inds, 2, 1];
    dists2 = Normal[TravelDistanceList[locs[[Reverse@inds]]]];
        AssociationThread[keys1, dists1],
        AssociationThread[keys2, dists2],
        Association[Thread[Transpose[{locs, locs}] -> Quantity[0, "Miles"]]]

This makes things much faster:

  lookup = TravelDistanceAssociation[schools];
  travelDistance[args__] := lookup[{args}];
  res1 = FindShortestTour[schools, DistanceFunction -> travelDistance];
) // AbsoluteTiming
{1.92451, Null}
res2 = FindShortestTour[schools, DistanceFunction -> TravelDistance]; // AbsoluteTiming
{14.6521, Null}

Verify the faster call gives the same result:

res1[[-1, -1]] == res2[[-1, -1]]
  • $\begingroup$ This is awesome, and why I wish I still had an excuse to use Mathematica in my day-to-day computing. What's the external API call? $\endgroup$ Dec 31, 2016 at 21:15
  • $\begingroup$ Also: is the arrow operator a new 10.x feature? I'm kinda behind, but I'm astounded I didn't know about it earlier! $\endgroup$ Dec 31, 2016 at 21:16
  • $\begingroup$ What ever service we're getting travel information from... I don't know off hand. $\endgroup$
    – Greg Hurst
    Dec 31, 2016 at 21:17
  • $\begingroup$ The documentation says Arrow has been around since version 6.0. $\endgroup$
    – Greg Hurst
    Dec 31, 2016 at 21:18
  • 1
    $\begingroup$ @Alexander, this says it's using GraphHopper + OpenStreetMap. $\endgroup$ Jan 14, 2017 at 15:59

This is a very excellent question. I won't give an answer here (and I doubt I could give one), because I think that cormullion's answer is already really nice and beautiful.

One really interesting mathematical object on that matter are space-filling curves.

For Meals-on-Wheels in Atlanta John Bartholdi and his coworkers applied space-filling curves to give a reasonable solution to this TSP-problem.

A Minimal Technology Routing System for Meals On Wheels

Maybe you find somebody who has access to Interfaces journal (http://interfaces.journal.informs.org/content/13/3/1.abstract)

Stan Wagon has in his outstanding book Mathematica in Action written a whole chapter about space-filling curves with load's of examples for tour generation and tour improvement.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.