Solving the Traveling Salesman Problem using InterCall Terry Robb 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. <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] ANNEAL[ (* TYPE=S *) $X -> In, (* DATA=R[$NCITY] *) $Y -> In, (* DATA=R[$NCITY] *) $IORDER -> Range[$NCITY] :> Out, (* DATA=I[$NCITY] *) $NCITY -> Min[ROWS[$X], ROWS[$Y]], (* DATA=I *) $ALEN -> Function[{x1, x2, y1, y2}, Sqrt[(x2 - x1)^2 + (y2 - y1)^2]], (* DATA=RF[R, R, R, R] *) $T -> 0.5, (* DATA=R *) $TFACTR -> 0.9, (* DATA=R *) $NOVER -> 100*$NCITY, (* DATA=I *) $NLIMIT -> 10*$NCITY, (* DATA=I *) $NTSTEPS -> 100, (* DATA=I *) $MONITR -> False, (* DATA=L *) $IWRK :> Null (* DATA=I[$NCITY] *) ] (* CODE="./anneal.o" *) anneal[$X_, $Y_] -> $IORDER 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}]; InterCall::opened: Opened connection to host solar InterCall::import: Importing: {ANNEAL} InterCall::linked: Using remote driver version 2.0 on host solar 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] {36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 12, 7, 15, 9, 6, 3, 33, 50, 18, 4, 49, 43, 28, 29, 47, 34, 32, 31, 25, 22, 20, 2, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11} 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]] 5.74068 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] {43, 49, 4, 34, 32, 31, 25, 18, 50, 33, 3, 12, 7, 15, 6, 9, 2, 20, 22, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11, 36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 29, 47, 28} 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]] 5.90193 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.