This notebook contains course material from CBE40455 by Jeffrey Kantor (jeff at nd.edu); the content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.
A set of airplanes are initially distributed among a set of starting locations. They are to be assigned routes to collectively visit a specified set of customers then return the the planes to designated finishing locations within designated time windows.
The data consists of a set of locations with latitude and longitude information, a list of customers and their respective locations, and a set of aircraft and their starting and finishing locations, and the start and of associated time windows. The aircraft must start and finish at different locations (if needed, dummy locations with the same latitude and longitude can be included in the list of locations). Plane speed is constrained to between upper and lower bounds.
The time windows are implemented as 'soft' constraints. Additional decision variables are
A weighted some of tea, tla, ted, and tld constitutes a time penalty which is zero if there is a feasible solution.
The objective function is weighted sum of the time penalty and total route distance.
%%script glpsol -m /dev/stdin
/* Vehicle Routing Problem with Time Windows
Jeffrey Kantor
March, 2013
*/
# DATA SETS (TO BE GIVEN IN THE DATA SECTION)
param maxspeed > 0;
param minspeed > 0, <= maxspeed;
# CUSTOMERS is a set of (name,location) pairs
set CUSTOMERS dimen 2;
param T1{CUSTOMERS};
param T2{(name,loc) in CUSTOMERS} >= T1[name,loc];
# PLANES is a set of (name, start_location, finish_location) triples
set PLANES dimen 3;
param S1{PLANES};
param S2{(p,sLoc,fLoc) in PLANES} >= S1[p,sLoc,fLoc];
param F1{PLANES};
param F2{(p,sLoc,fLoc) in PLANES} >= F1[p,sLoc,fLoc];
# set of locations
set LOCATIONS;
param lat{LOCATIONS};
param lng{LOCATIONS};
# DATA PREPROCESSING
# set of planes
set P := setof {(p,sLoc,fLoc) in PLANES} p;
# compute START as (plane,startlocation) pairs with time windows
set START := setof {(p,sLoc,fLoc) in PLANES} (p,sLoc);
param TS1{(p,sLoc) in START} :=
max{ (q,tLoc,fLoc) in PLANES : (p=q) && (sLoc=tLoc) } S1[p,sLoc,fLoc];
param TS2{(p,sLoc) in START} :=
min{ (q,tLoc,fLoc) in PLANES : (p=q) && (sLoc=tLoc) } S2[p,sLoc,fLoc];
# compute FINISH as (plane,finishlocation) pairs with time windows
set FINISH := setof {(p,sLoc,fLoc) in PLANES} (p,fLoc);
param TF1{(p,fLoc) in FINISH} :=
max{ (q,sLoc,gLoc) in PLANES : (p=q) && (fLoc=gLoc) } F1[p,sLoc,fLoc];
param TF2{(p,fLoc) in FINISH} :=
min{ (q,sLoc,gLoc) in PLANES : (p=q) && (fLoc=gLoc) } F2[p,sLoc,fLoc];
# create a complete of nodes as (name, location) pairs
set N := CUSTOMERS union (START union FINISH);
# great circle distances between locations
param d2r := 3.1415926/180;
param alpha{a in LOCATIONS, b in LOCATIONS} := sin(d2r*(lat[a]-lat[b])/2)**2
+ cos(d2r*lat[a])*cos(d2r*lat[b])*sin(d2r*(lng[a]-lng[b])/2)**2;
param gcdist{a in LOCATIONS, b in LOCATIONS} :=
2*6371*atan( sqrt(alpha[a,b]), sqrt(1-alpha[a,b]) );
# DECISION VARIABLES
# x[p,a,aLoc,b,bLoc] = 1 if plane p flies from (a,aLoc) to (b,bLoc)
var x{P, N, N} binary;
# START AND FINISH CONSTRAINTS
# no planes arrive at the start nodes
s.t. sf1 {p in P, (a,aLoc) in N, (b,bLoc) in START} :
x[p,a,aLoc,b,bLoc] = 0;
# no planes leave the finish nodes
s.t. sf2 {p in P, (a,aLoc) in FINISH, (b,bLoc) in N} :
x[p,a,aLoc,b,bLoc] = 0;
# planes must leave from their own start nodes
s.t. sf3 {p in P, (a,aLoc) in START, (b,bLoc) in N : p != a} :
x[p,a,aLoc,b,bLoc] = 0;
# planes must return to their own finish nodes
s.t. sf4 {p in P, (a,aLoc) in N, (b,bLoc) in FINISH : p != b} :
x[p,a,aLoc,b,bLoc] = 0;
# NETWORK CONSTRAINTS
# one plane arrives at each customer and finish node
s.t. nw1 {(b,bLoc) in (CUSTOMERS union FINISH)} :
sum {p in P, (a,aLoc) in (CUSTOMERS union START)} x[p,a,aLoc,b,bLoc] = 1;
# one plane leaves each start and customer node
s.t. nw2 {(a,aLoc) in (START union CUSTOMERS)} :
sum {p in P, (b,bLoc) in (CUSTOMERS union FINISH)} x[p,a,aLoc,b,bLoc] = 1;
# planes entering a customer node must leave the same node
s.t. nw3 {p in P, (a,aLoc) in CUSTOMERS} :
sum {(b,bLoc) in (CUSTOMERS union START)} x[p,b,bLoc,a,aLoc]
= sum {(b,bLoc) in (CUSTOMERS union FINISH)} x[p,a,aLoc,b,bLoc];
# no self loops
s.t. nw4 {p in P, (a,aLoc) in N, (b,bLoc) in N : (a=b) && (aLoc=bLoc)} :
x[p,a,aLoc,b,bLoc] = 0;
# SUBTOUR ELIMINATION CONSTRAINTS
var y{P,N,N} >= 0;
# route capacity
s.t. sb1 {p in P, (a,aLoc) in N, (b,bLoc) in N} :
y[p,a,aLoc,b,bLoc] <= card(CUSTOMERS)*x[p,a,aLoc,b,bLoc];
# allocate tokens to links from the start nodes
s.t. sb2 : sum {p in P, (a,aLoc) in START, (b,bLoc) in N } y[p,a,aLoc,b,bLoc]
= card(CUSTOMERS);
# decrease tokens for each step on a path
s.t. sb3 {(a,aLoc) in CUSTOMERS} :
sum{p in P, (b,bLoc) in (CUSTOMERS union START)} y[p,b,bLoc,a,aLoc]
= 1 + sum{p in P, (b,bLoc) in (CUSTOMERS union FINISH)} y[p,a,aLoc,b,bLoc];
# TIME WINDOW CONSTRAINTS
param bigM := 50;
var tar{N};
var tlv{N};
var tea{N} >= 0;
var tla{N} >= 0;
var ted{N} >= 0;
var tld{N} >= 0;
s.t. t00 {(a,aLoc) in N} : tlv[a,aLoc] >= tar[a,aLoc];
s.t. t01 {p in P, (a,aLoc) in N, (b,bLoc) in N} : tar[b,bLoc] >= tlv[a,aLoc]
+ gcdist[aLoc,bLoc]/maxspeed - bigM*(1-x[p,a,aLoc,b,bLoc]);
s.t. t02 {p in P, (a,aLoc) in N, (b,bLoc) in N} : tar[b,bLoc] <= tlv[a,aLoc]
+ gcdist[aLoc,bLoc]/minspeed + bigM*(1-x[p,a,aLoc,b,bLoc]);
s.t. t03 {(a,aLoc) in CUSTOMERS} : tea[a,aLoc] >= T1[a,aLoc] - tar[a,aLoc];
s.t. t04 {(a,aLoc) in FINISH} : tea[a,aLoc] >= TF1[a,aLoc] - tar[a,aLoc];
s.t. t05 {(a,aLoc) in CUSTOMERS} : tla[a,aLoc] >= tar[a,aLoc] - T2[a,aLoc];
s.t. t06 {(a,aLoc) in FINISH} : tla[a,aLoc] >= tar[a,aLoc] - TF2[a,aLoc];
s.t. t07 {(a,aLoc) in START} : ted[a,aLoc] >= TS1[a,aLoc] - tlv[a,aLoc];
s.t. t08 {(a,aLoc) in CUSTOMERS} : ted[a,aLoc] >= T1[a,aLoc] - tlv[a,aLoc];
s.t. t09 {(a,aLoc) in START} : tld[a,aLoc] >= tlv[a,aLoc] - TS2[a,aLoc];
s.t. t10 {(a,aLoc) in CUSTOMERS} : tld[a,aLoc] >= tlv[a,aLoc] - T2[a,aLoc];
# OBJECTIVE
# The objective function is a weighted sum of violations of the time window
# constraints and the total distance traveled.
var routeDistance{P} >= 0;
s.t. ob1 {p in P} : routeDistance[p]
= sum{(a,aLoc) in N, (b,bLoc) in N} gcdist[aLoc,bLoc]*x[p,a,aLoc,b,bLoc];
var totalDistance >= 0;
s.t. ob2 : totalDistance = sum{p in P} routeDistance[p];
var timePenalty >= 0;
s.t. ob3 : timePenalty =
sum{(a,aLoc) in N} (tea[a,aLoc] + 2*tla[a,aLoc] + 2*ted[a,aLoc] + tld[a,aLoc]);
minimize obj: 5*timePenalty + totalDistance/maxspeed;
solve;
# OUTPUT POST-PROCESSING
param routeTime{p in P} :=
sum{(a,aLoc) in N, (b,bLoc) in N} (tar[b,bLoc]-tlv[a,aLoc])*x[p,a,aLoc,b,bLoc];
param routeLegs{p in P} :=
sum{(a,aLoc) in START, (b,bLoc) in N} y[p,a,aLoc,b,bLoc];
for {p in P} {
printf "\nRouting for %s\n-------------------\n", p;
printf "%-24s %-24s %7s %5s %6s \n",
'Depart','Arrive','Dist.','Time','Speed';
for {k in routeLegs[p]..0 by -1} {
printf {(a,aLoc) in N, (b,bLoc) in N :
(x[p,a,aLoc,b,bLoc] = 1) && (abs(y[p,a,aLoc,b,bLoc]-k)<0.001)}
"%-12s %-4s %5.2f%1s %-12s %-4s %5.2f%1s %7.1f %5.2f %6.1f\n",
a, aLoc, tlv[a,aLoc],
if (ted[a,aLoc] > 0) then "E" else (if (tld[a,aLoc] > 0) then "L" else " "),
b, bLoc, tar[b,bLoc],
if (tea[b,bLoc] > 0) then "E" else (if (tla[b,bLoc] > 0) then "L" else " "),
gcdist[aLoc,bLoc], tar[b,bLoc]-tlv[a,aLoc],
if (gcdist[aLoc,bLoc] > 0) then gcdist[aLoc,bLoc]/(tar[b,bLoc]-tlv[a,aLoc]) else 0;
}
printf "%50s %7s %5s\n", '', '-------','-----';
printf "%50s %7.1f %5.2f\n\n", 'Totals:', routeDistance[p], routeTime[p];
}
# DATA SECTION
data;
param maxspeed := 800;
param minspeed := 600;
param : CUSTOMERS : T1 T2 :=
'Atlanta' ATL 8.0 24.0
'Boston' BOS 8.0 9.0
'Denver' DEN 12.0 15.0
'Dallas' DFW 12.0 13.0
'New York' JFK 18.0 20.0
'Los Angeles' LAX 12.0 16.0
'Chicago' ORD 20.0 24.0
'St. Louis' STL 11.0 13.0
;
param : PLANES : S1 S2 F1 F2 :=
'Plane 1' ORD ORD_ 8.0 24.0 8.0 24.0
'Plane 2' DFW DRW_ 8.0 24.0 8.0 24.0
;
param : LOCATIONS : lat lng :=
ATL 33.6366995 -84.4278639
BOS 42.3629722 -71.0064167
DEN 39.8616667 -104.6731667
DFW 32.8968281 -97.0379958 # start location
DRW_ 32.8968281 -97.0379958 # finish location
JFK 40.6397511 -73.7789256
LAX 33.9424955 -118.4080684
ORD 41.9816486 -87.9066714 # start location
ORD_ 41.9816486 -87.9066714 # finish location
STL 38.7486972 -90.3700289
;
end;
GLPSOL: GLPK LP/MIP Solver, v4.52 Parameter(s) specified in the command line: -m /dev/stdin Reading model section from /dev/stdin... Reading data section from /dev/stdin... /dev/stdin:253: warning: final NL missing before end of file 253 lines were read Generating sf1... Generating sf2... Generating sf3... Generating sf4... Generating nw1... Generating nw2... Generating nw3... Generating nw4... Generating sb1... Generating sb2... Generating sb3... Generating t00... Generating t01... Generating t02... Generating t03... Generating t04... Generating t05... Generating t06... Generating t07... Generating t08... Generating t09... Generating t10... Generating ob1... Generating ob2... Generating ob3... Generating obj... Model has been successfully generated GLPK Integer Optimizer, v4.52 1134 rows, 652 columns, 3896 non-zeros 288 integer variables, all of which are binary Preprocessing... 819 rows, 356 columns, 2662 non-zeros 146 integer variables, all of which are binary Scaling... A: min|aij| = 1.000e+00 max|aij| = 5.000e+01 ratio = 5.000e+01 GM: min|aij| = 7.959e-01 max|aij| = 1.256e+00 ratio = 1.579e+00 EQ: min|aij| = 6.427e-01 max|aij| = 1.000e+00 ratio = 1.556e+00 2N: min|aij| = 5.000e-01 max|aij| = 1.000e+00 ratio = 2.000e+00 Constructing initial basis... Size of triangular part is 809 Solving LP relaxation... GLPK Simplex Optimizer, v4.52 819 rows, 356 columns, 2662 non-zeros 0: obj = 3.100205000e+01 infeas = 1.464e+02 (10) * 94: obj = 5.821329277e+02 infeas = 1.486e-14 (10) * 215: obj = 7.397949066e+00 infeas = 2.720e-15 (10) OPTIMAL LP SOLUTION FOUND Integer optimization begins... + 215: mip = not found yet >= -inf (1; 0) + 861: >>>>> 5.168702483e+02 >= 1.133759393e+01 97.8% (53; 1) + 1122: >>>>> 2.901732459e+02 >= 1.137209032e+01 96.1% (66; 2) + 3214: >>>>> 2.541498427e+02 >= 1.179809046e+01 95.4% (171; 13) + 3446: >>>>> 1.987137841e+02 >= 1.184464254e+01 94.0% (176; 35) + 7867: >>>>> 1.868943569e+02 >= 1.240384510e+01 93.4% (368; 65) + 22229: >>>>> 1.759790800e+02 >= 1.369239888e+01 92.2% (1092; 143) + 24579: >>>>> 3.814149540e+01 >= 1.373895527e+01 64.0% (1147; 189) + 28966: >>>>> 2.325269970e+01 >= 1.416877206e+01 39.1% (629; 1362) + 30615: >>>>> 2.222023939e+01 >= 1.460974690e+01 34.3% (433; 1735) + 37353: >>>>> 2.071509200e+01 >= 1.702113332e+01 17.8% (298; 2271) + 40837: mip = 2.071509200e+01 >= tree is empty 0.0% (0; 3473) INTEGER OPTIMAL SOLUTION FOUND Time used: 3.1 secs Memory used: 5.2 Mb (5476969 bytes) Routing for Plane 1 ------------------- Depart Arrive Dist. Time Speed Plane 1 ORD 7.26E Boston BOS 9.00 1391.1 1.74 800.0 Boston BOS 9.00 St. Louis STL 11.10 1680.5 2.10 800.0 St. Louis STL 11.10 Atlanta ATL 12.07 779.0 0.97 800.0 Atlanta ATL 15.96 New York JFK 18.00 1222.1 2.04 600.0 New York JFK 18.02 Chicago ORD 20.00 1188.0 1.98 600.0 Chicago ORD 20.00 Plane 1 ORD_ 20.00 0.0 0.00 0.0 ------- ----- Totals: 6260.7 8.83 Routing for Plane 2 ------------------- Depart Arrive Dist. Time Speed Plane 2 DFW 12.00 Dallas DFW 12.00 0.0 0.00 0.0 Dallas DFW 12.00 Denver DEN 13.29 1032.1 1.29 800.0 Denver DEN 13.29 Los Angeles LAX 15.02 1385.1 1.73 800.0 Los Angeles LAX 15.02 Plane 2 DRW_ 17.50 1983.2 2.48 800.0 ------- ----- Totals: 4400.4 5.50 Model has been successfully processed