Traveling Salesman Tour Through U.S. Capital Cities (netsle07)
/***************************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: netsle07 */
/* TITLE: Traveling Salesman Tour Through U.S. Capital Cities (netsle07) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: GMAP, GPROJECT, OPTMODEL, PRINT, SQL */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 7 from the network solver documentation. */
/* */
/***************************************************************************/
/* get a list of the state capitals (with lat and long) */
proc sql;
create table Cities as
select unique statecode as state, city, lat, long
from mapsgfk.uscity
where capital='Y' and statecode not in ('AK' 'PR' 'HI');
quit;
data Cities;
set Cities;
if city='Nashville-Davidson' then city='Nashville';
city = trim(city) || ", " || state;
run;
/* create a list of all the possible pairs of cities */
proc sql;
create table CitiesDist as
select
a.city as city1, a.lat as lat1, a.long as long1,
b.city as city2, b.lat as lat2, b.long as long2,
geodist(lat1, long1, lat2, long2, 'DM') as distance
from Cities as a, Cities as b
where a.city < b.city;
quit;
/* find optimal tour by using the network solver */
proc optmodel;
set<str,str> CAPPAIRS;
set<str> CAPITALS = union {<i,j> in CAPPAIRS} {i,j};
num distance{i in CAPITALS, j in CAPITALS: i < j};
read data CitiesDist into CAPPAIRS=[city1 city2] distance;
set<str,str> TOUR;
num order{CAPITALS};
solve with NETWORK /
loglevel = moderate
links = (weight=distance)
tsp
out = (order=order tour=TOUR)
;
put (sum{<i,j> in TOUR} distance[i,j]);
/* create tour-ordered pairs (rather than input-ordered pairs) */
str CAPbyOrder{1..card(CAPITALS)};
for {i in CAPITALS} CAPbyOrder[order[i]] = i;
set TSPEDGES init
setof{i in 2..card(CAPITALS)} <CAPbyOrder[i-1],CAPbyOrder[i]>
union {<CAPbyOrder[card(CAPITALS)],CAPbyOrder[1]>};
num distance2{<i,j> in TSPEDGES} =
if i < j then distance[i,j] else distance[j,i];
create data TSPTourNodes from [node] tsp_order=order;
create data TSPTourLinks from [city1 city2]=TSPEDGES distance=distance2;
quit;
/* merge latitude and longitude */
proc sql;
/* merge in the lat & long for city1 */
create table TSPTourLinksAnno1 as
select unique TSPTourLinks.*, cities.lat as lat1, cities.long as long1
from TSPTourLinks left join cities
on TSPTourLinks.city1=cities.city;
/* merge in the lat & long for city2 */
create table TSPTourLinksAnno2 as
select unique TSPTourLinksAnno1.*, cities.lat as lat2, cities.long as long2
from TSPTourLinksAnno1 left join cities
on TSPTourLinksAnno1.city2=cities.city;
quit;
/* create the annotated data set to draw the path on the map */
data anno_path;
set TSPTourLinksAnno2;
length function color $8;
xsys='2'; ysys='2'; hsys='3'; when='a'; anno_flag=1;
function='move';
long = long1;
lat = lat1;
output;
function='draw';
color='blue'; size=0.8;
long = long2;
lat = lat2;
output;
run;
/* get a map of only the contiguous 48 states */
data states;
set mapsgfk.us_states (where=(statecode not in ('HI' 'AK' 'PR')));
run;
data combined;
set states anno_path;
run;
/* project the map and annotate the data */
proc gproject data=combined out=combined eastlong latlong degrees dupok;
id state;
run;
data states anno_path;
set combined;
if anno_flag=1 then output anno_path;
else output states;
run;
/* get a list of the endpoints locations */
proc sql;
create table anno_dots as
select unique x, y from anno_path;
quit;
/* create the final annotate data set */
data anno_dots;
set anno_dots;
length function color $8;
xsys='2'; ysys='2'; when='a'; hsys='3';
function='pie';
rotate=360; size=0.8; style='psolid'; color="red";
output;
style='pempty'; color="black";
output;
run;
/* generate the map with GMAP */
pattern1 v=s c=cxccffcc repeat=100;
proc gmap data=states map=states anno=anno_path all;
id state;
choro state / levels=1 nolegend coutline=black
anno=anno_dots des='' name="tsp";
run;
/* create the directed optimal tour */
data TSPTourLinksDirected(drop=next);
set TSPTourLinks;
retain next;
if _N_ ne 1 and city1 ne next then do;
city2 = city1;
city1 = next;
end;
next = city2;
run;
proc print data=TSPTourLinksDirected noobs label;
format distance comma10.2;
var city1 city2 distance;
sum distance;
run;