MODEL:
! The Vehicle Routing Problem (VRP), 
  with both load and trip length constraints,
  using MTZ formulation. Find allocation of customers/cities
  to vehicles so vehicle capacity not exceeded, vehicle trip
  time limit not exceeded, and
  total time/distance is minimized.
! Keywords: Vehicle Routing, Routing, LTL Delivery, Chart,
  Graph, Great circle distance;
SETS:
 ! Parameters:
     VCAP = vehicle capacity,
     Q(i) is the amount required at city i,
     DIST(i,j) is the distance from city i to city j,
     YCORD(i) = latitude of city i,
     XCORD(i) = longitude of city i,
     DISTMX = max distance allowed in a trip,
    Variables:
     Z(i,j) = 1 if some vehicle travels from city i to city j, else 0,
     LOADCUM(i) = accumulated deliveries/load at city i,
     DISTCUM(i) = accumulated distance at city i ;
   CITY: Q, LOADCUM, XCORD, YCORD, DISTCUM, RUNTIME;
   DEPOT( CITY); ! Identify depot city;
   CXC( CITY, CITY): DIST, Z; 
   CXCSUB( CXC): DCITY, ACITY;

ENDSETS
DATA: ! Here data are taken in-place. It could also be read from files. Click on: Edit | Paste Function | External Files; TIMETOT = 400; ! Solver time limit in seconds (0:no limit) for @SOLVE's; ! Data set USA1; !CaseUSA1; VCAP = 24 !CaseUSA1; DISTMX = 7000 !CaseUSA1; DISTMTYP = 1 !CaseUSA1; CITY, Q, YCORD, XCORD = BrwnsvllTX 4 25.9000 -97.4333 Chicago 7 41.8781 -87.6298 Dallas 0 32.7767 -96.7970 DelRio 2 29.3667 -100.7833 Detroit 2 42.4167 -83.0167 Minneapolis 6 44.9800 -93.2519 NewYorkCity 9 40.7128 -74.0059 Oakland 9 37.7297 -122.2200 Phoenix 5 33.4833 -112.0667 Pittsburgh 5 40.4406 -79.9959 Portland 4 45.5997 -122.5997 Salt_Lake_C 8 40.7500 -111.8833 SanDiego 3 32.7333 -117.1667 Seattle 2 47.6062 -122.3321 StLouis 11 38.6270 -90.1994 Tucson 6 32.2217 -110.9258 Tallahassee 5 30.3833 -84.3667 ; !CaseUSA1; DEPOT = DALLAS! RUNTIME(j) = time spent at city j; !CaseUSA1; RUNTIME = 0 !CaseAlden VCAP = 18; ! Additional options; !CaseAlden DISTMX = 999999; !CaseAlden DISTMTYP = 2; !CaseAlden CITY= Chi Den Frsn Hous KC LA Oakl Anah Peor Phnx Prtl Rvrs Sacr SLC ; !CaseAlden DIST = !CaseAlden 0 996 2162 1067 499 2054 2134 2050 151 1713 2083 2005 2049 1390! Chicago; !CaseAlden 996 0 1167 1019 596 1059 1227 1055 904 792 1238 1010 1142 504 !CaseAlden 2162 1167 0 1747 1723 214 168 250 2070 598 745 268 162 814! Fresno; !CaseAlden 1067 1019 1747 0 710 1538 1904 1528 948 1149 2205 1484 1909 1438 !CaseAlden 499 596 1723 710 0 1589 1827 1579 354 1214 1809 1535 1742 1086! K. City; !CaseAlden 2054 1059 214 1538 1589 0 371 36 1943 389 959 54 376 715 !CaseAlden 2134 1227 168 1904 1827 371 0 407 2043 755 628 425 85 744! Oakland; !CaseAlden 2050 1055 250 1528 1579 36 407 0 1933 379 995 45 412 711 !CaseAlden 151 904 2070 948 354 1943 2043 1933 0 1568 2022 1889 1958 1299! Peoria; !CaseAlden 1713 792 598 1149 1214 389 755 379 1568 0 1266 335 760 648 !CaseAlden 2083 1238 745 2205 1809 959 628 995 2022 1266 0 1001 583 767! Portland; !CaseAlden 2005 1010 268 1484 1535 54 425 45 1889 335 1001 0 430 666 !CaseAlden 2049 1142 162 1909 1742 376 85 412 1958 760 583 430 0 659! Sacrmnto; !CaseAlden 1390 504 814 1438 1086 715 744 711 1299 648 767 666 659 0; ! Pallets to be delivered to each customer; !CaseAlden Q= 0 6 3 7 7 18 4 5 2 6 7 2 4 3 ; !CaseAlden RUNTIME= 0; !CaseAlden DEPOT = CHI; !CaseAlden XCORD = 0; !CaseAlden YCORD = 0; ENDDATA ! Warning: For >> 12 cities, this simple formulation may give good solutions, but it may not be able to prove optimality. Larger problems require more sophisticated formulations. Contact www.lindo.com; SUBMODEL VROUTE: ! Minimize total travel distance for each stop k both before k and after k; ! MIN = @SUM( CITY(k): Q(k)*(DB4(k)+DFT(k))); ! Minimize total distance; MIN = DISTOT; DISTOT = @SUM( CXC(i,j): DIST(i,j)*Z(i,j)); ! a vehicle does not travel inside itself,...; @FOR( CITY(k): Z( k, k) = 0); ! For each city k, except depot....; @FOR( CITY( k)| k #NE# DEPOTI: ! a vehicle must enter city K from some city I, together their demands cannot exceed vehicle capacity... ; [NTR] @SUM( CITY( i)| i #NE# k #AND# ( i #EQ# DEPOTI #OR# Q( i) + Q( k) #LE# VCAP): Z( i, k)) = 1; ! a vehicle must leave K after service to some city J; [XIT] @SUM( CITY( j)| j #NE# k #AND# ( j #EQ# DEPOTI #OR# ((Q( j) + Q( k) #LE# VCAP) #AND# (DIST(DEPOTI,k) + DIST(k,j) + DIST(j,DEPOTI) #LE# DISTMX))): Z( k, j)) = 1; ! LOADCUM( k) is at least amount needed at K, but can't exceed vehicle capacity; @BND( Q( k), LOADCUM( k), VCAP); ! If i precedes k, then can bound LOADCUM( k) - LOADCUM( i); @FOR( CITY( i)| i #NE# k #AND# i #NE# DEPOTI: [LCUM] LOADCUM( k) >= LOADCUM(i) + Q(k)*Z(i,k) ! Case: i precedes k; - Q(i)*Z(k,i) ! Case: k precedes i; +(Q(k) - VCAP)*(1-Z(i,k)-Z(k,i)) ; ! Case: neither above; ); ! If i precedes k, then can bound DISTCUM( k) - DISTCUM( i), DISTCUM( k) = cumulative distance in trip after visiting city k; @FOR( CITY( i)| i #NE# k #AND# i #NE# DEPOTI: [DCUM] DISTCUM( k) >= DISTCUM(i) + DIST(i,k)* Z(i,k) ! Case: i precedes k; - DIST(k,i)* Z( k, i) ! Case: k precedes i; +( DIST(DEPOTI,k) - ( DISTMX - DIST(i,DEPOTI)))*(1-Z(i,k) - Z(k,i)) ; ! Case: neither above; ); ! Trip length constraint; DISTCUM( k) + DIST(k, DEPOTI) <= DISTMX; ); ! Bound on number of vehicles; ! Must send enough vehicles out of depot; @SUM(CITY(j) | j #NE# DEPOTI # AND# j #LE# NUMC: Z( DEPOTI,j)) >= VEHLO; ! But not too many; @SUM(CITY(j) | j #NE# DEPOTI # AND# j #LE# NUMC: Z( DEPOTI,j)) <= VEHHI; ! Make the Z's binary; @FOR( CXC(i,j): @BIN( Z(i,j))); ! Some cuts; ! Minimum no. vehicles required, fractional and rounded up; VEHCLF = @SUM( CITY( I)| I #GT# 1: Q( I))/ VCAP; VEHCLR = @FLOOR( VEHCLF +.9999); ! Must send enough vehicles out of depot; @SUM(CITY(j) | j #NE# DEPOTI: Z( DEPOTI,j)) >= VEHCLR; @FOR( CITY(k) | k #NE# DEPOTI: ! A cut: If K is 1st stop, then LOADCUM( k) = Q( k); LOADCUM( k) <= VCAP - ( VCAP - Q( k)) * Z( DEPOTI, k); ! A cut: If K is not 1st stop...; LOADCUM( k) >= Q( k) + @SUM( CITY( i)| i #NE# DEPOTI: Q( i) * Z( i, k)); ! A cut: If K is not last stop...; LOADCUM( k) <= VCAP - @SUM( CITY( i)| i #NE# DEPOTI: Q( i) * Z( k, i)); ); ! Subtour size 2 cuts,( Case Q(i) + Q(j) #GT# VCAP already taken care of); @FOR( CITY(i) | i #NE# DEPOTI: @FOR( CITY(j) | j #GT# i #AND# Q(i) +Q(j) #LE# VCAP: Z(i,j) + Z(j,i) <= 1; ); ); ! Subtours of size 3 cuts,( Case Q(i) + Q(j) #GT# VCAP already taken care of); @FOR( CITY(i) | i #NE# DEPOTI: @FOR( CITY(j) | j #GT# i #AND# Q(i) + Q(j) #LE# VCAP: @FOR( CITY(k) | k #GT# j: Z(i,j) + Z(j,i) + Z(i,k) + Z(k,i) + Z(j,k) + Z(k,j) <= 2; ); ); ); ! 3 item knapsack cuts,( Case Q(i) + Q(j) #GT# VCAP already taken care of); @FOR( CITY(i) | i #NE# DEPOTI: @FOR( CITY(j) | j #GT# i #AND# Q(i) + Q(j) #LE# VCAP: @FOR( CITY(k) | k #GT# j #AND# Q(i) + Q(j) + Q(k) #GT# VCAP: Z(i,j) + Z(j,i) + Z(i,k) + Z(k,i) + Z(j,k) + Z(k,j) <= 1; ); ); ); ENDSUBMODEL
CALC: TOLRELOPT = 0.05; ! Set ending relative optimality tolerance; TIME2ROPT = 120; ! Time in seconds to apply optimality tolerance; @SET( 'TERSEO',2); ! Output level (0:verb, 1:terse, 2:only errors, 3:none); @SET( 'IPTOLR', TOLRELOPT); ! Set ending relative optimality tolerance; @SET( 'TIM2RL', TIME2ROPT); ! Time in seconds to apply optimality tolerance; @SET( 'TATSLV', TIMETOT); ! Solver time limit in seconds (0:no limit) for @SOLVE's; @SET( 'OROUTE',1); ! Route output immediately to the window line by line ! Index of depot city; @FOR( DEPOT( k): DEPOTI = k); NUMC = @SIZE( CITY); ! Number cities in problem; ! Compute distance matrix; ! 0: x-y coordinates, 1: latitude-longitude, 2: explicit distance matrix; @IFC( DISTMTYP #eq# 0: ! This portion calculates the distance matrix DIST(i,j) using Euclidean distance; @FOR( CXC(i,j): @IFC( i #EQ# j: DIST(i,j) = 0;! Get rid of trivial roundoff; @ELSE DIST( i,j) = ((XCORD(i)-XCORD(j))^2 + (YCORD(i) - YCORD(j))^2)^0.5; ); ); ! End Lat-long distance calculation; ); @IFC( DISTMTYP #eq# 1: ! This portion calculates the distance matrix DIST(i,j) assuming XCORD and YCORD are the longitude and latitude in degrees; D2R=@PI()/180; ! Degrees to radians conversion factor; ! Compute Great Circle Distances. Radius of earth = 6371 km. Notice this simplifies if YCORD(i) = YCORD(j) or XCORD(i) = XCORD(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*YCORD(i))*@SIN(D2R*YCORD(j))+@COS(D2R*YCORD(i))*@COS(D2R*YCORD(j)) *@COS(@ABS(D2R*(XCORD(i)-XCORD(j))))); ); ); ); ! End Lat-long distance calculation; ! Adjust for run time; @FOR( CXC(i,j) | i #LE# NUMC #AND# j #LE# NUMC: DIST(i,j) = DIST(i,j) + RUNTIME(j); ); @SOLVE( VROUTE); ISTAT = @STATUS(); @IFC( ISTAT #EQ# 0 #OR# ISTAT #EQ# 4: NROUTES = 0; ! Write a listing of the routes; @FOR( CITY( j): @IFC( Z( DEPOTI, j) #GT# .5: ! Is j the depot city on a route?; NROUTES = NROUTES + 1; TRIPCUM = DIST(DEPOTI,j); ! Cumulative distance of this route; LOADCM = Q( j); @WRITE( @NEWLINE( 2), 'ROUTE ', NROUTES, ':', @NEWLINE( 1)); @WRITE(' FROM TO DISTANCE LOAD', @NEWLINE( 1)); @WRITE('------------------------------------------------', @NEWLINE( 1)); @WRITE( ' ', @FORMAT( CITY( DEPOTI), '-12s'), @FORMAT(CITY( j),'-12s'), @FORMAT( TRIPCUM, '10.1f'), @FORMAT( LOADCM, '10.1f'), @NEWLINE( 1) ); IPOS = J; ! Find next city J2 after IPOS until returning to DEPOTI; @WHILE( IPOS #NE# DEPOTI: NLOOPS = NLOOPS + 1; @FOR( CITY( JJ) | Z( IPOS, JJ) #GT# 0.5: J2 = JJ; TRIPCUM = TRIPCUM + DIST( IPOS, J2); LOADCM = LOADCM + Q( j2); @WRITE( ' ', @FORMAT( CITY( IPOS), '-12s'), @FORMAT(CITY( J2),'-12s'), @FORMAT( TRIPCUM, '10.1f'), @FORMAT( LOADCM, '10.1f'), @NEWLINE( 1) ); ); IPOS = J2; ); ); ); ); ! Output a graph of the solution, but only if we have coordinates of the cities; @IFC( @SUM( CITY(i): XCORD(i)) #NE# 0: !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; ); @CHARTNETNODE( ' Least Cost LTL Deliveries, KM= ' + @FORMAT(DISTOT,'8.2f') ! Title of chart; , 'Latitude', 'Longitude' ! Labels for horizontal and vertical; , 'Cities' ! Legend for arc set 1; , XCORD, YCORD ! Coordinates of the nodes; , DCITY, ACITY); ! Node pairs of arcs actually used; ); ENDCALC END