! Traveling Salesman Problem.  Find the shortest tour that
   visits each city exactly once.
 In terms of a communication network, a TSP solution is a ring network.
 It is the simplest network that can tolerate the failure of one arc
 and still allow all nodes to communicate.
 ! The Miller, Tucker, Zemlin, 1960, J. ACM, single commodity formulation;
 ! Keywords: Network Design, TSP, Traveling sales person, Routing, Tour, Sequencing, 
   Great circle distance, Chart, Graph;
SETS:
  CITY : LATI, LNGT, LVL;
  CXC( CITY, CITY): DIST, Z;
  CXCSUB( CXC): DCITY, ACITY, ARROHD;
ENDSETS
DATA: ! Data can be pulled from a spreadsheet by using CITY, LATI, LNGT = @OLE(), with correspondingly named ranges in the open spreadsheet; CITY, LATI, LNGT= ! Calgary 51.0486 -114.0708 Edmonton 53.5444 -113.4909 Halifax 44.6488 -63.5752 Montreal 45.5017 -73.5673 Ottawa 45.4215 -75.6972 QuebecCity 46.8139 -71.2080 StJohns 47.5605 -52.7128 Saskatoon 52.1297 -106.6557 Toronto 43.6532 -79.3832 Vancouver 49.2827 -123.1207 Winnipeg 49.8998 -97.1375 ; CapeTown -33.9249 18.4241 Johannesburg -26.2041 28.0473 Durban -29.8587 31.0218 Pretoria -25.7479 28.2293 PortElizabeth -33.7139 25.5207 Bloemfontein -29.0852 26.1596 Nelspruit -25.4714 30.9865 Kimberley -28.7282 24.7499 ;! Beijing 39.9042 116.4074 Chengdu 30.5728 104.0668 Chongqing 29.5630 106.5516 Guangzhou 23.1291 113.2644 Guilin 25.2345 110.1800 Harbin 45.8038 126.5350 Shanghai 31.2304 121.4757 Shenzhen 22.5431 114.0579 Xian 34.3416 108.9398 HongKong 22.3964 114.1095 Tianjin 39.0842 117.2010 ; ! Hiroshima 34.3852 132.4553 Kobe 34.6901 135.1955 Kyoto 35.0116 135.7680 Tokyo 35.6895 139.6917 Nagasaki 32.7503 129.8777 Nagoya 35.1814 136.9064 Osaka 34.6937 135.5022 Sapporo 43.0621 141.3544 Sendai 38.2682 140.8694 Lyon 45.7640 4.8357 Marseille 43.2965 5.3698 Montpellier 43.6108 3.8767 Nantes 47.2184 1.5536 Nice 43.7102 7.2620 Paris 48.8566 2.3522 Strasbourg 48.5734 7.7521 Toulouse 43.6047 1.4442 Bellingham 48.8000 -122.3833 BrownsvilleTX 25.9000 -97.4333 CaribouME 46.9000 -68.0167 Chicago 41.8781 -87.6298 Dallas 32.7767 -96.7970 DelRio 29.3667 -100.7833 Denver 39.7392 -104.9903 Detroit 42.4167 -83.0167 ElPaso 31.8000 -106.4000 Fresno 36.7700 -119.7200 Galveston 29.3000 -94.8000 Houston 29.7604 -95.3698 Internat_Falls 48.5667 -93.3833 KansasCity 39.3197 -94.7200 KeyLargo 24.9333 -80.2833 LasVegas 36.1667 -115.2000 LongBeach 33.8197 -118.1500 Los_Angeles 34.0522 -118.2428 Miami 25.7617 -80.1918 Minneapolis 44.9800 -93.2519 NewYorkCity 40.7128 -74.0059 Oakland 37.7297 -122.2200 Peoria 40.6697 -89.6800 Philadelphia 39.9526 -75.1652 Phoenix 33.4833 -112.0667 Pittsburgh 40.4406 -79.9959 Portland 45.5997 -122.5997 Salt_Lake_City 40.7500 -111.8833 SanAntnnio 29.4241 -98.4936 SanDiego 32.7333 -117.1667 SanFran 37.6167 -122.3833 Seattle 47.6062 -122.3321 StLouis 38.6270 -90.1994 Tucson 32.2217 -110.9258 SStMarie 46.4667 -84.3667 NEWPTVT 45.0000 -73.1500 Tallahassee 30.3833 -84.3667 Mumbai 19.0760 72.8777 Kolkata 22.5726 88.3639 Delhi 28.6139 77.2090 Chennai 13.0827 80.2707 Bangalore 12.9716 77.5946 Hyderabad 17.3850 78.4867 Ahmedabad 23.0225 72.5714 Jaipur 26.9124 75.7873 Lucknow 26.8467 80.9462 Nagpur 21.1458 79.0882 Patna 25.5941 85.1376 Indore 22.7196 75.8577 Vadodara 22.3072 73.1812 Bhopal 23.2599 77.4126 Thoothukudi 8.7642 78.1348 Pune 18.5204 73.8567 Surat 21.1702 72.8311 Kanpur 26.4499 80.3319 ; ! Cairo 30.0444 31.2357 Alexandria 31.2001 29.9187 Port_Said 31.2653 32.3019 Suez 29.9668 32.5498 Luxor 25.6872 32.6396 Mansoura 31.0409 31.3785 Tanta 30.7865 31.0004 Sidi_Barrani 31.6111 25.9257 Aswan 24.0889 32.8998 Wadi_Halfa 21.7991 31.3713; ! El_Mahalla 30.9781 31.1624 Giza 30.0131 31.2089 Shubra_El_Kh 30.1217 31.2452; ENDDATA SUBMODEL TSPROB:! Warning: May take long to solve cases with N >> 12; ! Variables: Z(i,j) = 1 if vehicle tour includes link from i to j, else 0, LVL(i) = sequence number of stop i on the tour, or also load on vehicle upon departing stop i. One unit of load is picked up at each stop. Starts empty at stop 1; ! Minimize total distance traveled; MIN = OBJV; OBJV= @SUM( CXC(i,j): DIST(i,j) * Z(i,j)); @FOR( CITY( k): ! It must be entered exactly once; @SUM( CITY( i)| i #NE# k: Z( i, k)) = 1; ! It must be departed exactly once; @SUM( CITY( j)| j #NE# k: Z( k, j)) = 1; Z( k, k) = 0; ! Cannot go from k to k; ); ! A weak but simple form of the subtour breaking constraints, see Desrochers & Laporte, OR Letters, Feb. 91. Not very powerful for large(N>12) problems. For large problems need to use more complicated subtour elimination methods; ! Enforce: If Z(i,j) = 1 then u(j) - LVL(i) = 1, If Z(j,i) = 1 then LVL(j) - LVL(i) = -1, If Z(i,j) + Z(j,i) = 0, and i,j > 1, then LVL(j) - LVL(i) >= -(N - 2); LVL(1) = 0; ! Start empty; ! The case either i or j = 1; @FOR( CITY(i) | i #GT# 1: LVL(i) >= 2 - Z(1,i) + (N-3)*Z(i,1); LVL(i) <= (N-2) + Z(i,1) - (N-3)*Z(1,i); ); ! The case i,j > 1, ! This constraint, plus its "mirror", when i and j are switched, forces LVL(j) - LVL(i) = 1 if Z(i,j) = 1; @FOR( CXC(i,j)| i #GT# 1 #AND# j #GT# 1 #AND# i #NE# j: LVL( j) >= LVL( i) + Z(i,j) - Z(j,i) - (N-2)*(1 - Z(i,j) - Z(j,i)); ); ! Make the Z's 0/1; @FOR( CXC(i,j): @BIN( Z(i,j));); ! Some optional cuts, which may or may not help; ! We know the sum of the stop numbers; @SUM( CITY( i): LVL( i)) = ( N-1)*N/2; ! Two-city subtour breaking cuts; ! @FOR( CXC( i,j) | i #LT# j: Z(i,j) + Z(j,i) <= 1 ); ENDSUBMODEL
CALC: @SET( 'TERSEO',2); ! Output level (0:verb, 1:terse, 2:only errors, 3:none); ! This portion calculates the distance matrix DIST(i,j); D2R=@PI()/180; ! Degrees to radians conversion factor; ! Compute Great Circle Distances. Radius of earth = 6371 km. Notice this simplifies if LATI(i) = LATI(j) or LNGT(i) = LNGT(j); @FOR( CXC(i,j): @IFC( i #EQ# j: DIST(i,j) = 0;! Get rid of trivial roundoff; @ELSE DIST( i,j) = 6371*@acos(@SIN(D2R*LATI(i))*@SIN(D2R*LATI(j))+@COS(D2R*LATI(i))*@COS(D2R*LATI(j)) *@COS(@ABS(D2R*(LNGT(i)-LNGT(j))))); ); ); ! Display the distance matrix; ! @SET( 'LINLEN', 150); ! @TABLE( DIST); !The model size: Warning, may be slow for N > 10; N = @SIZE( CITY); ! @GEN( TSPROB); @SOLVE( TSPROB); ! Construct the subset, CXCSUB(i,j), of arcs selected; @FOR( CXC( i,j) | Z(i,j) #GT# 0.5: @INSERT( CXCSUB, i, j); DCITY(i,j) = i; ! Departure city; ACITY(i,j) = j; ! Arrival city; ARROHD(i,j) = 1; ! No arrowheads on this arc; ); ! @CHARTNETNODE( ' Least Cost Ring (1 Arc redundant, TSP) Network, KM= ' + @FORMAT(OBJV,'8.2f') ! Title of chart; ! , 'Longitude', 'Latitude' ! , 'Cities' ! Legend for arc set 1; ! , LNGT, LATI ! , DCITY, ACITY); ! Node pairs of arcs actually used ! Use feature (LINGO 16 and later) that you can optionally specify direction of arrowheads on arcs; @CHARTNETNODE( ' Least Cost Ring (1 Arc redundant, TSP) Network, KM= ' + @FORMAT(OBJV,'8.2f') ! Title of chart; , 'Longitude', 'Latitude' ! Labels for horizontal and vertical; , 'Cities' ! Legend for arc set 1; , LNGT, LATI ! Coordinates of the nodes; , DCITY, ACITY ! Node pairs of arcs actually used; , ARROHD); ! Specify, optionally, where to put arrowhead; ENDCALC