## Traveling Salesman Problem (misole04)

```/*****************************************************************/
/*                                                               */
/*          S A S   S A M P L E   L I B R A R Y                  */
/*                                                               */
/*    NAME: misole04                                             */
/*   TITLE: Traveling Salesman Problem (misole04)                */
/* PRODUCT: OR                                                   */
/*  SYSTEM: ALL                                                  */
/*    KEYS: OR                                                   */
/*   PROCS: OPTMODEL, GPLOT                                      */
/*    DATA:                                                      */
/*                                                               */
/* SUPPORT:                             UPDATE:                  */
/*     REF:                                                      */
/*    MISC: Example 4 from the Mixed Integer Linear Programming  */
/*          Solver chapter of Mathematical Programming.          */
/*                                                               */
/*****************************************************************/

/*
%let tsplib = location-of-st70-problem
*/

/* convert the TSPLIB instance into a data set */
data tspData(drop=H);
infile "&tsplib";
input H \$1. @;
if H not in ('N','T','C','D','E');
input @1 var1-var3;
run;

/* direct solution using the MTZ formulation */
proc optmodel;
set VERTICES;
set EDGES = {i in VERTICES, j in VERTICES: i ne j};
num xc {VERTICES};
num yc {VERTICES};
/* read in the instance and customer coordinates (xc, yc) */
read data tspData into VERTICES=[_n_] xc=var2 yc=var3;
/* the cost is the euclidean distance rounded to the nearest integer */
num c {<i,j> in EDGES}
init floor( sqrt( ((xc[i]-xc[j])**2 + (yc[i]-yc[j])**2)) + 0.5);
var x {EDGES} binary;
var u {i in 2..card(VERTICES)} >= 2 <= card(VERTICES);
/* each vertex has exactly one in-edge and one out-edge */
con assign_i {i in VERTICES}:
sum {j in VERTICES: i ne j} x[i,j] = 1;
con assign_j {j in VERTICES}:
sum {i in VERTICES: i ne j} x[i,j] = 1;
/* minimize the total cost */
min obj
= sum {<i,j> in EDGES} (if i > j then c[i,j] else c[j,i]) * x[i,j];
/* no subtours */
con mtz {<i,j> in EDGES : (i ne 1) and (j ne 1)}:
u[i] - u[j] + 1 <= (card(VERTICES) - 1) * (1 - x[i,j]);
solve with milp / maxtime = 600;
quit;

/* iterative solution using the subtour formulation */
proc optmodel;
set VERTICES;
set EDGES = {i in VERTICES, j in VERTICES: i > j};
num xc {VERTICES};
num yc {VERTICES};

num numsubtour init 0;
set SUBTOUR {1..numsubtour};

/* read in the instance and customer coordinates (xc, yc) */
read data tspData into VERTICES=[var1] xc=var2 yc=var3;

/* the cost is the euclidean distance rounded to the nearest integer */
num c {<i,j> in EDGES}
init floor( sqrt( ((xc[i]-xc[j])**2 + (yc[i]-yc[j])**2)) + 0.5);

var x {EDGES} binary;

/* minimize the total cost */
min obj =
sum {<i,j> in EDGES} c[i,j] * x[i,j];

/* each vertex has exactly one in-edge and one out-edge */
con two_match {i in VERTICES}:
sum {j in VERTICES: i > j} x[i,j]
+ sum {j in VERTICES: i < j} x[j,i] = 2;

/* no subtours (these constraints are generated dynamically) */
con subtour_elim {s in 1..numsubtour}:
sum {<i,j> in EDGES: (i in SUBTOUR[s] and j not in SUBTOUR[s])
or (i not in SUBTOUR[s] and j in SUBTOUR[s])} x[i,j] >= 2;

/* this starts the algorithm to find violated subtours */
set <num,num> EDGES1;
set VERTICES1 = union{<i,j> in EDGES1} {i, j};
num component {VERTICES1};
num numcomp  init 2;
num iter     init 1;
num numiters init 1;
set ITERS = 1..numiters;
num sol {ITERS, EDGES};

/* initial solve with just matching constraints */
solve;
call symput(compress('obj'||put(iter,best.)),
trim(left(put(round(obj),best.))));
for {<i,j> in EDGES} sol[iter,i,j] = round(x[i,j]);

/* while the solution is disconnected, continue */
do while (numcomp > 1);
iter = iter + 1;

/* find connected components of support graph   */
EDGES1 = {<i,j> in EDGES: round(x[i,j].sol) = 1};
solve with network /
nodes   = (include=VERTICES1)
concomp
out     = (concomp=component);

numcomp = _oroptmodel_num_["NUM_COMPONENTS"];
if numcomp = 1 then leave;
numiters = iter;
numsubtour = numsubtour + numcomp;
for {comp in 1..numcomp} do;
SUBTOUR[numsubtour-numcomp+comp]
= {i in VERTICES: component[i] = comp};
end;

solve;
call symput(compress('obj'||put(iter,best.)),
trim(left(put(round(obj),best.))));
for {<i,j> in EDGES} sol[iter,i,j] = round(x[i,j]);
end;

/* create a data set for use by sgplot */
create data solData from
[iter i j]={it in ITERS, <i,j> in EDGES: sol[it,i,j] = 1}
x1=xc[i] y1=yc[i] x2=xc[j] y2=yc[j];
call symput('numiters',put(numiters,best.));
quit;

%macro plotTSP;
%do i = 1 %to &numiters;
/* create annotate data set to draw subtours */
data anno(drop=iter);
retain drawspace 'datavalue' linethickness 1 function 'line';
set solData;
where iter = &i;
run;

title1 h=2 "TSP: Iter = &i, Objective = &&obj&i";
title2;

proc sgplot data=tspData sganno=anno;
scatter x=var2 y=var3 / datalabel=var1;
xaxis display=none;
yaxis display=none;
run;
%end;
%mend plotTSP;

%plotTSP;

```