Lindo Systems

 ! Traveling Salesman Problem.   (TSP_Chart.lng)
 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: Chart, ChartNetNode, Graph, Great circle distance, Network Design,  
     Routing, Sequencing, Tour, Traveling sales person, TSP;
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= 
!  Buddy Holly's "Tour from Hell" or "Winter Dance Party" 
   (along with J.P. "The Big Bopper" Richardson,  Ritchie Valens,
   Waylon Jennings, Dion DiMucci (of the Belmonts), Tommy Allsup, Carl Bunch);
!BHolly      MilwaukeeWI    43.0389  -87.9065   ! 1959 Jan 23;
!BHolly      KenoshaWI      42.5847  -87.8212   ! 1959 Jan 24;
!BHolly      MankatoMN      44.1636  -93.9994   ! 1959 Jan 25;
!BHolly      EauClaireWI    44.8113  -91.4985   ! 1959 Jan 26;
!BHolly      MontevideoMN   44.9410  -95.7236   ! 1959 Jan 27;
!BHolly      StPaulMN       44.9537  -93.09     ! 1959 Jan 28;
!BHolly      DavenportIA    41.5236  -90.5776   ! 1959 Jan 29;
!BHolly      FortDodgeIA    42.4975  -94.1680   ! 1959 Jan 30;
!BHolly      DuluthMN       46.7867  -92.1005   ! 1959 Jan 31;
!BHolly      GreenBayWI     44.5192  -88.0198   ! 1959 Feb 1;
!BHolly      ClearLakeIA    43.1436  -93.3788   ! 1959 Feb 2 The day the music died;
!BHolly      MoorheadMN     46.8738  -96.7678   ! 1959 Feb 3;
!BHolly      SiouxCityIA    42.5000  -96.4003   ! 1959 Feb 4;
!BHolly      DesMoinesIA    41.6005  -93.6091   ! 1959 Feb 5;
!BHolly      CedarRapidsIA  41.9779  -91.6656   ! 1959 Feb 6;
!BHolly      SpringValleyIL 41.3275  -89.1998   ! 1959 Feb 7;
!BHolly      ChicagoIL      41.8781  -87.6298   ! 1959 Feb 8;
!BHolly      WaterlooIA     42.4928  -92.3426   ! 1959 Feb 9;
!BHolly      DubuqueIA      42.5006  -90.6646   ! 1959 Feb 10;
!BHolly      LouisvilleKY   38.2527  -85.7585   ! 1959 Feb 11;
!BHolly      CantonOH       40.7989  -81.3784   ! 1959 Feb 12;
!BHolly      YoungstownOH   41.0998  -80.6495   ! 1959 Feb 13;
!BHolly      PeoriaIL       40.6936  -89.5890   ! 1959 Feb 14;
!BHolly      SpringfieldIL  39.7817  -89.6501   ! 1959 Feb 15;


!     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;

! Note, the fractional parts in Germany are in Minutes,
 not fraction of a degree;
!  Berlin     52.31  13.23 
  Hamburg    53.33  10.0 
  Munich     48.8   11.34 
  Cologne    50.56   6.57 
  Frankfurt  50.7    8.41 
  Stuttgart  48.47   9.11 
  Dusseldorf 51.14   6.47 
  Dortmund   51.31   7.28 
  Essen      51.27   7.1 
  Leipzig    51.20  12.23 
  Bremen     53.5    8.48 
  Dresden    51.2   13.44 
  Hanover    52.22   9.43 
  Nuremberg  49.27  11.5 
  Duisburg   51.26   6.46 
  Bochum     51.29   7.13 
  Wuppertal  51.16   7.11 
  Bielefeld  52.1    8.32 
  Bonn       50.44   7.6 
;
!     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
     Xian            34.3416   108.9398
     HongKong        22.3964   114.1095
     Tianjin         39.0842   117.2010 
;!     Shenzhen        22.5431   114.0579

     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_Giza      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;
!     Giza            30.0131    31.2089;
!    El_Mahalla      30.9781    31.1624
!     Shubra_El_Kh    30.1217    31.2452;
!Thai Bangkok  	5104476
!Thai Samut Prakan  	388920
!Thai Mueang Nonthaburi  	291555
!Thai Udon Thani  	247231
!Thai Chon Buri  	219164
!Thai Nakhon Ratchasima  	208781
!Thai Chiang Mai  	200952
!Thai Hat Yai  	191696
!Thai Pak Kret  	182926
!Thai Si Racha  	178916
!Thai Phra Pradaeng  	171333
!Thai Lampang  	156139
!Thai Khon Kaen  	147579
!Thai Surat Thani  	127201
!Thai Ubon Ratchathani  	122533
!Thai Nakhon Si Thammarat  	120836
!Thai Khlong Luang  	118551
!Thai Nakhon Pathom  	117927
!Thai Rayong  	106737
!Thai Phitsanulok    103427
;
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

SUBMODEL OrgOrd:
! Use the input order;
  @FOR( CITY(i) | i #GT# 1:
    Z(i-1,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); !  Terminal page width (0:none);
! @TABLE( DIST);

 !The model size: Warning, may be slow for N > 10;
 N = @SIZE( CITY);

! @GEN( TSPROB);
 @SOLVE( TSPROB);   ! Get optimal tour;
! @SOLVE( TSPROB, ORGORD); ! Use orginal order;

! Construct the subset, CXCSUB(i,j), of arcs selected;
! Z(i,j) = 1 if the tour goes from i to j;

 @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) = 0; ! Arrowheads(0=no, 1=yes) on this arc;
     );

   @CHARTNETNODE( 
!        ' TSP Tour, Original Order, KM= ' + @FORMAT(OBJV,'8.2f')    ! Title of chart;
        ' TSP Tour, Optimal Order, (Directed)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;

 ! Use feature (LINGO 17 and later) that you can optionally specify direction of arrowheads on arcs;
  @CHARTNETNODE( 
        ' Shortest TSP Tour, (Undirected) 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