Solving the Traveling Salesman Problem using InterCall

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

Abstract

The traveling salesman problem is a classic example of combinatorial minimization. It involves finding the shortest itinerary for a traveling salesman who must visit each of N cities in turn. One technique for solving this problem is the method of simulated annealing, and this notebook uses InterCall to access an external subroutine which implements the annealing algorithm. An example in which an additional constraint, such as adding a penalty for a river crossing, will also be demonstrated.

For a concrete example we use the following collection of 'cities':

The problem

The "traveling salesman problem" involves finding the shortest itinerary for a traveling salesman who must visit each of N cities in turn. For a concrete example we shall look at a 50 city problem.

Setup the position of the cities

Here are 50 random positions of some cities.

positions = {
 {.527,.892}, {.381,.666}, {.460,.446}, {.245,.361},
 {.200,.921}, {.445,.495}, {.687,.448}, {.840,.192}, 
 {.462,.555}, {.133,.999}, {.992,.850}, {.523,.385}, 
 {.071,.952}, {.431,.947}, {.588,.550}, {.346,.800},
 {.978,.872}, {.251,.418}, {.538,.978}, {.338,.631}, 
 {.641,.086}, {.252,.650}, {.329,.845}, {.002,.925},
 {.153,.549}, {.772,.743}, {.710,.004}, {.276,.124}, 
 {.207,.008}, {.382,.885}, {.100,.580}, {.035,.511},
 {.416,.398}, {.049,.450}, {.200,.759}, {.915,.606}, 
 {.797,.788}, {.858,.565}, {.758,.101}, {.929,.327}, 
 {.609,.953}, {.867,.601}, {.362,.291}, {.943,.187}, 
 {.294,.894}, {.582,.168}, {.089,.091}, {.126,.831},
 {.343,.336}, {.310,.431}};
{x,y} = Transpose[positions];

Graphically the positions of the cities, given in the abstract earlier, could be displayed with

(*
cities= {PointSize[0.03], Point/@positions};
Show[ Graphics[{cities}], AspectRatio->1 ];
*)

We could use ListPlot too, but for consistency with later examples we use more basic graphic primitives instead.

The method

We shall use the method of simulated annealing with an annealing schedule as described in the Numerical Recipes book (see references at end). Annealing is very CPU intensive (although it is one of the best methods available) and would be inefficient to code in Mathematica directly. Luckily the Numerical Recipes book gives a fortran implementation that we can import into Mathematica fairly easily using InterCall.

Startup InterCall and import an ANNEAL routine

Load in the InterCall package

First we must load in the following package.

<<InterCall.m

See the references at the end for details about obtaining InterCall.

Declare defaults for the ANNEAL routine that does simulated annealing

Next we must tell Mathematica some details about the fortran code we are going to import. Details such as datatypes, sensible default values as well as the location of the external object code can be set with the AddDefault command. The routine we shall import is a slightly modified version of the ANNEAL routine in the Numerical Recipes book.

AddDefault[ ANNEAL[
             $X->In,
             $Y->In,
             $IORDER->Range[$NCITY] :> Out,
             $NCITY->Min[ROWS[$X],ROWS[$Y]],
             $ALEN->Function[{x1,x2,y1,y2},
                         Sqrt[(x2-x1)^2+(y2-y1)^2]],
             $T->0.5,
             $TFACTR->0.9,
             $NOVER->100*$NCITY,
             $NLIMIT->10*$NCITY,
             $NTSTEPS->100,
             $MONITR->False,
             $IWRK:>Null ],
           S[R[$NCITY],R[$NCITY],I[$NCITY],I,
                RF[R,R,R,R],R,R,I,I,I,L,I[$NCITY]],
           "./anneal.o" ];

You should read the Numerical Recipes book to understand the arguments and default settings indicated above. The default setting says that the default cost function ALEN between any two cities is the euclidean distance between those cities.

As in the book the annealing schedule is based on the following:  T is decreased by TFACTR for at most NTSTEPS, and at each step it is held constant for NOVER reconfigurations or NLIMIT successful reconfigurations whichever comes first.

Check the default setting

Having entered a default setting for the anneal routine, we can inspect the setting by

GetDefault[anneal]

The summary line at the end of the output above, shows that we need only supply the co-ordinates X and Y of the cities, and that the best order IORDER will be returned. (The summary line alone could be obtained by typing just ?anneal.)

Import the routine into Mathematica

To import the routine, and any other routines for that matter, into Mathematica we just use InterCall's Import command.

Import[{ANNEAL}];

In our case the output indicates that a connection is made to the computer named solar which in this case is the machine running the Mathematica kernel. It's possible to import code from another machine too however.

Call ANNEAL to find the best order to visit all the cities

For this example, the ALEN function is called about a million times, and so for efficiency the InterCall compiler/interpreter must be switched on. To compute a near optimal path the following takes about 30 seconds of real time to execute.

order = anneal[x,y]

Plot the solution

Having done all the CPU intensive work with the imported fortran routine, we can now use Mathematica's graphics capabilities to present the solution.

route = positions[[Append[order,First[order]]]];
cities= {PointSize[0.03],Point/@positions};
Show[ Graphics[{cities,Line[route]}],
      AspectRatio->1 ];

Compute the path length

Note, due to round-off in the annealing algorithm, some computers may give a slightly different route than the one found above. Because the routine is based on random numbers it does not always find the optimal solution, but always something very close to it.

Plus @@ Map[Sqrt[#.#]&, route-RotateLeft[route]]

Solve the same problem but with a river

In this example we shall place a river along the line y=0.3 and impose a high penalty for crossing it. To do this we modify the cost function ALEN to include a penalty of lambda if the river is crossed, where the value of lambda can be adjusted to give the relative importance of path-length versus river-crossing. If lambda is high enough then the solution will cross the river a minimum number of times, ie. just twice.

It is possible to write the penalty function as a normal Mathematica function. Even though the external anneal routine is a fortran program it will, with InterCall, still be possible to call the Mathematica–written penalty function. This ability to seamlessly link fortran or C code with Mathematica code is one of the features that makes InterCall so useful.

penalty = Function[{x1,x2,y1,y2},Module[{lambda=1,m1,m2},
           m1=If[y1<0.3,.5,-.5]; m2=If[y2<0.3,.5,-.5];
           Sqrt[(x2-x1)^2+(y2-y1)^2]+lambda*(m2-m1)^2]];
order = anneal[x,y,$ALEN->penalty]

The above just looks like a normal Mathematica interaction, but of course underneath fortran code is being called (possibly on a remote computer) and also the penalty function is being remotely interpreted outside Mathematica. InterCall hides all these inner details from the user, who need only worry about writing Mathematica code. All other details are handled by InterCall.

Plot the solution with a river

As before we can plot the solution using Mathematica graphics:

route = positions[[Append[order,First[order]]]];
river = {Thickness[0.03],GrayLevel[0.7],Line[{{0,0.3},{1,.3}}]};
cities= {PointSize[0.03],Point/@positions};
Show[ Graphics[{river,cities,Line[route]}],
      AspectRatio->1 ];

To illustrate the difference in the solution with the river constraint added, you could select the two solutions and animate them. This will highlight the minor difference that is needed to avoid the river.

Compute the path length of the above solution

Plus @@ Map[Sqrt[#.#]&, route-RotateLeft[route]]

The path length is slightly larger that before, as this is the price to pay if you don't like crossing rivers.

Exercises

Modify the anneal.f program so that it can solve 3-dimensional problems. Use it to solve a "fruit picking problem". Plot (and maybe animate) your answer.

Try an animation problem by moving one of the cities around and seeing how the solution changes. Tip – because each iteration in the animation will be close to the previous solution, use the option $IORDER->order and also reduce the temperature $T to make the computation faster. This will require some experimentation.

Solve the 50 city problem by eye. An easy way to do this is to select the graphic cell showing the 50 cities and command-click on each city in the order you want. Then use copy (command-c), and paste (command-v) into an input cell. This will give you the co-ordinates, in order chosen.

References

ANNEAL

The ANNEAL routine is an adaptation of a routine in the Numerical Recipes book by William H. Press, Brian P. Flannery, Saul A. Teukolsky and William T. Vetterling, 1987, 1992 New York, Cambridge University Press. Diskettes containing all the Numerial Recipes software can ordered from:

Numerical Recipes Software                   or           Customer Services Dept.,
Cambridge University,                                          Cambridge University Press,
110 Midland Avenue,                                            Edinburgh Building,  Shaftesbury Road,
Port Chester, NY 10573,                                      Cambridge CB2 2RU,
USA.                                                                     UK.

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource – email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands 6009
Perth WA
AUSTRALIA.

Appendix

The modifications to the original anneal.f program

The external routine used in this Notebook, is a slightly modified version of the anneal subroutine given in the Numerical Recipes book. The routine has been converted to double precision and some control parameters have been passed as arguments and not hard-wired as in the original.

Specifically we have used real*8 statements for all real variables, and the original ALEN inline function has been passed as an argument to ANNEAL, TRNCST and REVCST. Workspace (named IWRK) for use by TRNSPT has also be passed as an argument to ANNEAL. Finally the number of T steps, which in the original was hard-wired to 100, has been changed to NTSTEPS and passed as an argument too.