# # Vehicle Routing
#
# 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. The optimization objective is to minimize the total great circle distance.
#
# 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. 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).
#
# In[1]:
get_ipython().run_cell_magic('script', 'glpsol -m /dev/stdin -y VehicleRouting.txt --out output', '\n/* Vehicle Routing Problem\n Jeffrey Kantor\n March, 2013\n*/\n\n# DATA SETS (TO BE GIVEN IN THE DATA SECTION)\n\n# CUSTOMERS is a set of (name,location) pairs \nset CUSTOMERS dimen 2;\n\n# PLANES is a set of (name, start_location, finish_location) triples\nset PLANES dimen 3;\n\n# set of locations\nset LOCATIONS;\nparam lat{LOCATIONS};\nparam lng{LOCATIONS};\n\n# DATA PREPROCESSING\n\n# set of planes\nset P := setof {(p,sLoc,fLoc) in PLANES} p;\n\n# create a complete of nodes as (name, location) pairs\nset START := setof {(p,sLoc,fLoc) in PLANES} (p,sLoc);\nset FINISH := setof {(p,sLoc,fLoc) in PLANES} (p,fLoc);\nset N := CUSTOMERS union (START union FINISH);\n\n# great circle distances between locations\nparam d2r := 3.1415926/180;\nparam alpha{a in LOCATIONS, b in LOCATIONS} := sin(d2r*(lat[a]-lat[b])/2)**2 \n + cos(d2r*lat[a])*cos(d2r*lat[b])*sin(d2r*(lng[a]-lng[b])/2)**2;\nparam gcdist{a in LOCATIONS, b in LOCATIONS} := \n 2*6371*atan( sqrt(alpha[a,b]), sqrt(1-alpha[a,b]) );\n\n# DECISION VARIABLES\n\n# x[p,a,aLoc,b,bLoc] = 1 if plane p flies from (a,aLoc) to (b,bLoc)\nvar x{P, N, N} binary;\n\n# START AND FINISH CONSTRAINTS\n\n# no planes arrive at the start nodes\ns.t. sf1 {p in P, (a,aLoc) in N, (b,bLoc) in START} : x[p,a,aLoc,b,bLoc] = 0;\n\n# no planes leave the finish nodes\ns.t. sf2 {p in P, (a,aLoc) in FINISH, (b,bLoc) in N} : x[p,a,aLoc,b,bLoc] = 0;\n\n# planes must leave from their own start nodes\ns.t. sf3 {p in P, (a,aLoc) in START, (b,bLoc) in N : p != a} : x[p,a,aLoc,b,bLoc] = 0;\n\n# planes must return to their own finish nodes\ns.t. sf4 {p in P, (a,aLoc) in N, (b,bLoc) in FINISH : p != b} : x[p,a,aLoc,b,bLoc] = 0;\n\n# NETWORK CONSTRAINTS\n\n# one plane arrives at each customer and finish node\ns.t. nw1 {(b,bLoc) in (CUSTOMERS union FINISH)} : \n sum {p in P, (a,aLoc) in (CUSTOMERS union START)} x[p,a,aLoc,b,bLoc] = 1;\n\n# one plane leaves each start and customer node\ns.t. nw2 {(a,aLoc) in (START union CUSTOMERS)} :\n sum {p in P, (b,bLoc) in (CUSTOMERS union FINISH)} x[p,a,aLoc,b,bLoc] = 1;\n\n# planes entering a customer node must leave the same node\ns.t. nw3 {p in P, (a,aLoc) in CUSTOMERS} : \n sum {(b,bLoc) in (CUSTOMERS union START)} x[p,b,bLoc,a,aLoc]\n = sum {(b,bLoc) in (CUSTOMERS union FINISH)} x[p,a,aLoc,b,bLoc];\n\n# no self loops\ns.t. nw4 {p in P, (a,aLoc) in N, (b,bLoc) in N : (a=b) && (aLoc=bLoc)} :\n x[p,a,aLoc,b,bLoc] = 0;\n\n# SUBTOUR ELIMINATION CONSTRAINTS\n\nvar y{P,N,N} integer, >= 0;\n\n# route capacity\ns.t. sb1 {p in P, (a,aLoc) in N, (b,bLoc) in N} : \n y[p,a,aLoc,b,bLoc] <= card(CUSTOMERS)*x[p,a,aLoc,b,bLoc];\n\n# allocate tokens to links from the start nodes\ns.t. sb2 : sum {p in P, (a,aLoc) in START, (b,bLoc) in N } y[p,a,aLoc,b,bLoc] \n = card(CUSTOMERS);\n\n# decrease tokens for each step on a path\ns.t. sb3 {(a,aLoc) in CUSTOMERS} : \n sum{p in P, (b,bLoc) in (CUSTOMERS union START)} y[p,b,bLoc,a,aLoc] \n = 1 + sum{p in P, (b,bLoc) in (CUSTOMERS union FINISH)} y[p,a,aLoc,b,bLoc];\n\n# OBJECTIVE\n\nvar routeDistance{P} >= 0;\ns.t. ob1 {p in P} : routeDistance[p] \n = sum{(a,aLoc) in N, (b,bLoc) in N} gcdist[aLoc,bLoc]*x[p,a,aLoc,b,bLoc];\n\nvar routeLegs{P} >= 0;\ns.t. ob2 {p in P} : routeLegs[p] = sum{(a,aLoc) in START, (b,bLoc) in N} y[p,a,aLoc,b,bLoc];\n\nvar maxDistance >= 0;\ns.t. ob3 {p in P} : routeDistance[p] <= maxDistance;\n\nvar maxLegs >= 0;\ns.t. ob4 {p in P} : routeLegs[p] <= maxLegs;\n\nminimize distance : sum{p in P} routeDistance[p];\n\nsolve;\n\n# OUTPUT POST-PROCESSING\n\nfor {p in P} {\n printf "\\nRouting for %s\\n-------------------\\n", p;\n printf "%-20s %-20s %10s \\n", \'Depart\',\'Arrive\',\'Dist.\';\n for {k in routeLegs[p]..0 by -1} {\n printf {(a,aLoc) in N, (b,bLoc) in N : \n (x[p,a,aLoc,b,bLoc] = 1) && (y[p,a,aLoc,b,bLoc]=k)} \n\t "%-12s %-5s %-12s %-5s %10.1f km\\n",a,aLoc,b,bLoc,gcdist[aLoc,bLoc];\n }\n printf "%42s %13s\\n", \'\', \'---------\';\n printf "%42s %10.1f km\\n\\n", \'GC Distance Traveled:\', routeDistance[p];\n}\n\n# DATA SECTION\n\ndata;\n\nset CUSTOMERS := \n ( \'Atlanta\', ATL )\n ( \'Boston\', BOS )\n ( \'Denver\', DEN )\n ( \'Dallas\', DFW )\n ( \'New York\', JFK )\n ( \'Los Angeles\', LAX )\n ( \'Chicago\', ORD )\n ( \'St. Louis\', STL ) \n;\n\nset PLANES :=\n ( \'Plane 1\', ORD, ORD_) # use a duplicate location to return place to ORD\n ( \'Plane 2\', DFW, DRW_) # use a duplicate location to return plane to DFW\n;\n\nparam : LOCATIONS : lat lng :=\n ATL 33.6366995 -84.4278639\n BOS 42.3629722 -71.0064167\n DEN 39.8616667 -104.6731667\n DFW 32.8968281 -97.0379958\n DRW_ 32.8968281 -97.0379958 # duplicate PLANES\n JFK 40.6397511 -73.7789256\n LAX 33.9424955 -118.4080684\n ORD 41.9816486 -87.9066714\n ORD_ 41.9816486 -87.9066714 # duplicate PLANES\n STL 38.7486972 -90.3700289\n; \n\nend;\n')
# In[2]:
print(open('VehicleRouting.txt').read())
# In[ ]:
#
# < [Assignment Problems](http://nbviewer.jupyter.org/github/jckantor/CBE40455/blob/master/notebooks/05.02-Assignment-Problems.ipynb) | [Contents](toc.ipynb) | [Vehicle Routing with Time Windows](http://nbviewer.jupyter.org/github/jckantor/CBE40455/blob/master/notebooks/05.04-Vehicle-Routing-with-Time-Windows.ipynb) >