Lindo Systems
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; ! Maximum load on any trip;
!CaseUSA1; DISTMX = 7000;! Maximum distance of any trip;
!CaseUSA1; DISTMTYP = 1; ! 0: x-y coordinates, 1: latitude-longitude, 2: explicit distance matrix;
!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;! Depot city;
! RUNTIME(j) = time spent at city j;
!CaseUSA1; RUNTIME = 0;
!CaseAlden VCAP = 18; ! VCAP = pallet capacity of a vehicle ;
! Additional options;
!CaseAlden DISTMX = 999999; ! Maximum distanced of any route;
!CaseAlden DISTMTYP = 2; ! 0: x-y coordinates, 1: latitude-longitude, 2: explicit distance matrix;
!CaseAlden CITY=
Chi Den Frsn Hous KC LA Oakl Anah Peor Phnx Prtl Rvrs Sacr SLC ;
!CaseAlden DIST = ! To City From city;
!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! Denver;
!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! Houston;
!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! L. A.;
!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! Anaheim;
!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! Phoenix;
!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! Rverside;
!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;! SLC;
! 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; ! RUNTIME(j) = time spent at city j. Run times are 0,;
!CaseAlden DEPOT = CHI; ! First task or city or depot in the tour;
!CaseAlden XCORD = 0; ! Not using X coordinates;
!CaseAlden YCORD = 0; ! Not using Y coordinates;
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