Proc Optmodel and CVRP : Turning off the decomposition algo | LSO

[ODENWALD]

Hi all :

When solving  CVRPs, for learning/teaching purposes I wanted to quantify the performance gain of the decomposition algo: however I found no way to turn it off. Is there a way ?

And :  Just to get a rough idea where a Genetic Algorithm would get to  I used  <solve with  LSO> : this resulted in an early stop with <  WARNING: The best solution found does not satisfy the feasibility tolerance.
NOTE: Failed.  >

Any ideas how to get some methodological comparison data (timings) - preferably within SAS/OR ?

I experimented with several instances of  VRP-Lib,  a particular instance with rather optimal  .block settings and good performance is from SAS documentation (OR - 15.1)  ::

/* vrpdata represents an instance from VRPLIB that has 22 nodes and eight vehicles (P-n22-k8.vrp),
which was originally described in Augerat et al. (1995).
The data set lists each node, its coordinates, and its demand. */

/* number of vehicles available */
%let num_vehicles = 8;
/* capacity of each vehicle */
%let capacity = 3000;
/* node, x coordinate, y coordinate, demand */
data vrpdata;
input node x y demand;
datalines;
1 145 215 0
2 151 264 1100
3 159 261 700
4 130 254 800
5 128 252 1400
6 163 247 2100
7 146 246 400
8 161 242 800
9 142 239 100
10 163 236 500
11 148 232 600
12 128 231 1200
13 156 217 1300
14 129 214 1300
15 146 208 300
16 164 208 900
17 141 206 2100
18 147 193 1000
19 164 193 900
20 129 189 2500
21 155 185 1800
22 139 182 700
;

proc optmodel;
/* read the node location and demand data */
set NODES;
num x {NODES};
num y {NODES};
num demand {NODES};
num capacity = &capacity;
num num_vehicles = &num_vehicles;
read data vrpdata into NODES=[node] x y demand;
set ARCS = {i in NODES, j in NODES: i ne j};
set VEHICLES = 1..num_vehicles;

/* define the depot as node 1 */
num depot = 1;

/* define the arc cost as the rounded Euclidean distance */
num cost {<i,j> in ARCS} = round(sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2));

/* Flow[i,j,k] is the amount of demand carried on arc (i,j) by vehicle k */
var Flow {ARCS, VEHICLES} >= 0 <= capacity;
/* UseNode[i,k] = 1, if and only if node i is serviced by vehicle k */
var UseNode {NODES, VEHICLES} binary;
/* UseArc[i,j,k] = 1, if and only if arc (i,j) is traversed by vehicle k */
var UseArc {ARCS, VEHICLES} binary;

/* minimize the total distance traversed */
min TotalCost = sum {<i,j> in ARCS, k in VEHICLES} cost[i,j] * UseArc[i,j,k];

/* each non-depot node must be serviced by at least one vehicle */
con Assignment {i in NODES diff {depot}}:
sum {k in VEHICLES} UseNode[i,k] >= 1;

/* each vehicle must start at the depot node */
for{k in VEHICLES} fix UseNode[depot,k] = 1;

/* some vehicle k traverses an arc that leaves node i
if and only if UseNode[i,k] = 1 */
con LeaveNode {i in NODES, k in VEHICLES}:
sum {<(i),j> in ARCS} UseArc[i,j,k] = UseNode[i,k];

/* some vehicle k traverses an arc that enters node i
if and only if UseNode[i,k] = 1 */
con EnterNode {i in NODES, k in VEHICLES}:
sum {<j,(i)> in ARCS} UseArc[j,i,k] = UseNode[i,k];

/* the amount of demand supplied by vehicle k to node i must equal demand
if UseNode[i,k] = 1; otherwise, it must equal 0 */
con FlowBalance {i in NODES diff {depot}, k in VEHICLES}:
sum {<j,(i)> in ARCS} Flow[j,i,k] - sum {<(i),j> in ARCS} Flow[i,j,k]
= demand[i] * UseNode[i,k];

/* if UseArc[i,j,k] = 1, then the flow on arc (i,j) must be at most capacity
if UseArc[i,j,k] = 0, then no flow is allowed on arc (i,j) */
con VehicleCapacity {<i,j> in ARCS, k in VEHICLES}:
Flow[i,j,k] <= Flow[i,j,k].ub * UseArc[i,j,k];

/* decomp by vehicle */
for {i in NODES, k in VEHICLES} do;
LeaveNode[i,k].block = k;
EnterNode[i,k].block = k;
end;
for {i in NODES diff {depot}, k in VEHICLES} FlowBalance[i,k].block = k;
for {<i,j> in ARCS, k in VEHICLES} VehicleCapacity[i,j,k].block = k;

/* solve using decomp (aggregate formulation) */
solve with MILP / varsel=ryanfoster decomp=(logfreq=20);

/* create solution data set */
create data solution_data from [i j k]=
{<i,j> in ARCS, k in VEHICLES: UseArc[i,j,k].sol > 0.5}
x1=x[i] y1=y[i] x2=x[j] y2=y[j]
function='line' drawspace='datavalue';
quit;

proc sgplot data=solution_data noautolegend;
scatter x=x1 y=y1 / datalabel=i;
vector x=x2 y=y2 / xorigin=x1 yorigin=y1 group=k noarrowheads;
xaxis display=none;
yaxis display=none;
run;

Any tips are welcome.   

Odenwald

[RobPratt]

To invoke the default MILP algorithm (branch and cut) change the SOLVE statement to just:

solve;

In SAS/OR 15.1, the decomposition algorithm is still sometimes automatically used as a heuristic, but you can disable that functionality as follows:

solve with milp / heuristics=basic;

Attached is another approach to CVRP, using dynamically generated subtour elimination constraints. together with a repair heuristic and an animation.

[ODENWALD]

Thanks Rob.

Your alternative approach is very good for teaching/understanding.

Here the timing results on my PC :

1. Decent decomposition :  real: 28.24 sec   optimal value: 603

2. solve ;  real: stopped after 1600 sec when best_integer=603   best_bound=549.15   gap=9.8%

                (some calls to decomp. algo being made)

3. solve with MILP /heuristics=basic ;

   real: stopped after 1200 sec when best_integer=608   best_bound=545.67   gap=11.42%

           (after  20 minutes and 32444 nodes and  27225  active ones

                (no calls to decomp. algo being mentioned, pure branch & cut with SAS defaults)

The strong superiority of a 'decent' decomposition set up by user's understanding of structure of constraints being obvious.

Last question :  Is it possible to get  < solve with  LSO >  to useful/non-trivial results by adjusting some parameters ?

Can an Optmodel model formulation easily been put into proc GA  ?

Tabu search or Simulated Annealing seem not to be available in SAS/OR. Or ?

Odenwald

[josgri]

Wrt LSO:

Applying a GA to CVRP’s requires using sequence variable encoding and/or customized permutation/crossover functions.  The corresponding dimension of the MILP problem tends to have too many variables and constraints for a GA to be applicable in general-purpose form.

For example, consider the TSP formulation for MILP:

https://go.documentation.sas.com/?docsetId=casmopt&docsetTarget=casmopt_milpsolver_examples04.htm&do...

versus the TSP formulation for a GA:

https://go.documentation.sas.com/?docsetId=orlsoug&docsetTarget=orlsoug_ga_examples01.htm&docsetVers...

Unfortunately, LSO does not support sequence variables or custom crossover/mutation functions, and would thus not work for this problem type.  In general, a GA that treats the problem as a Blackbox will have trouble handling integer variable constraints and/or variable dimensions exceeding 100; by working with sequence variables, the constraints would be implicitly satisfied for all permutations and crossovers while the variable dimension is dramatically reduced.  Similarly, custom crossover/mutation operators could simply jump from one feasible point to the next.

It might be possible to try something similar using PROC GA which does support both of these features.  Invoking the PROC OPTMODEL code directly may be challenging.

As of yet, we do not support Tabu search or Simulated Annealing (SA).  Adding support for both sequence variables and SA has been regularly discussed but is not yet on target for development, unfortunately.

[ODENWALD]

Hi josgri :

Just to remark why I posed the question:

We read and reviewed some material about  VRP  and the many variations that exist today. I frequently read the remark that an exact solution would be (usually) impossible in practice ;  solution strategies being frequently mentioned were  GA, SA, Tabu Search (all heuristics in a wide sense). I felt that the authors missed the advance in  MILP-progress, say over the last  10 - 15 years.

This has been confirmed.

But if you check this and other simple examples, if you do  NOT  hit a 'decent' decomposition, the  MILP's solver solution way is excruciatingly slow (to reach a  10% gap  or a  5% gap or ...).

And clearly exponential growth comes fast and is expensive to avoid for a mid-sized company.

I'm not aware if  SAS Inst.  or Gurobi, just to name a few, make 'good' offers for  VRP-solving on demand on big (SAS-)machines.

Just to make a fast comparison in some more or less small problems I wanted to contrast   MILPDecomposition  to SimpleMILP and to one of  GA, SA, Tabu.

Your note clarifies what's possible with SAS now (2020).

I don't know if there is a kind of statistic about how many SAS customers solve  XVRPs , how many would want to, and what algorithmic approach and hardware those who solve (up to  certain gaps) do use.

I did not see how to setup  Proc GA  for  CVRP.

If someone at SAS has please let me know.

Odenwald

[ODENWALD]

Hi :

Before closing the topic I had another try at problem  P-n55-k8.vrp  same source. I pasted it into the program of the first mail at the place of the P-n22-k8 data  and started Optmodel with same Decomp setup. After   4 hours  and some  seconds  I  canceled the program. The SAS log is attached.

Do you (does anyone who has much more experience than I do with SAS's implementation of  MILP  and  Decomp.)  see what settings could have been changed for better performance ?

I would say that   P-n22-k8    gives a slightly too positive impression.

To me the essential need for   GA , SA, Tabu Search  or the like  is totally obvious now.

Odenwald

[josgri]

Hi Odenwald,

I am very sorry for the slow response.  I have been following up in the background.  I am very familiar with LSO, not so much PROC GA/MILP/CVRP.  I followed up with our MILP team and they are discussing the possibility of a dedicated VRP solver similar to our TSP.  So crossing fingers that project gets a green light. 

Note that LSO has been renamed "Blackbox" for future releases.  

Disclaimer, not being an expert in either PROC GA or CVRP I took a stab at transforming the PROC GA TSP problem into a CVRP example.  It seems to run correctly now, but I can't be certain.  So it is a sketch of how one would do it.  It gives results to the output like the below in cycle format. 

Cycle for vehicle: 1
1 -> 5 (load = 1400 )
5 -> 4 (load = 2200 )
4 -> 8 (load = 3000 )
8 -> 1 (load = 3000 )
Cycle for vehicle: 2

1 -> 12 (load = 1200 )
12 -> 14 (load = 2500 )
14 -> 15 (load = 2800 )
15 -> 1 (load = 2800 )
Cycle for vehicle: 3
1 -> 22 (load = 700 )
22 -> 18 (load = 1700 )
18 -> 13 (load = 3000 )
13 -> 1 (load = 3000 )
Cycle for vehicle: 4

data locations; 
input node x y demand;
datalines;
1 145 215 0
2 151 264 1100
3 159 261 700
4 130 254 800
5 128 252 1400
6 163 247 2100
7 146 246 400
8 161 242 800
9 142 239 100
10 163 236 500
11 148 232 600
12 128 231 1200
13 156 217 1300
14 129 214 1300
15 146 208 300
16 164 208 900
17 141 206 2100
18 147 193 1000
19 164 193 900
20 129 189 2500
21 155 185 1800
22 139 182 700
;
run;

proc ga data1 = locations seed = 5554 lastgen=lastgen;
   call SetEncoding('S22');
   ncities = 22;
   ncar = 8;
   cap  = 3000;
   /* v1 = 1 ->  2 ->  3 ->  4 -> 1   (indices of S owned by vehicial 1)
      v2 = 1 ->  5 ->  6 ->  7 -> 1
      v3 = 1 ->  8 ->  9 -> 10 -> 1
      v4 = 1 -> 11 -> 12 -> 13 -> 1
      v5 = 1 -> 14 -> 15 -> 16 -> 1
      v6 = 1 -> 17 -> 18 -> 1
      v7 = 1 -> 19 -> 20 -> 1
      v8 = 1 -> 21 -> 22 -> 1 */

   /* Create distance matix once up front */
   array distances[22,22] /nosym; 
   do i = 1 to 22; 
      do j = 1 to i;
         distances[i,j] = sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2);
         distances[j,i] = distances[i,j];
      end;
   end;
 
   /* Populate vehicle map start/stop */
   array vehmap[8,2] / nosym;
   cnt=2;
   inc=2;
   do i=1 to 8; 
      vehmap[i,1]=cnt;
      vehmap[i,2]=cnt+inc;
      cnt=cnt+inc+1;
      if (i > 4) then inc=1;  
   end;  
 
   /* Objective function with local optimization */
   function CVRPSwap(selected[*],ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]);
   array s[1] /nosym; 
 
   call dynamic_array(s,ncities);
   call ReadMember(selected,1,s);

   /* Permut city 1 back to front */
   do i = 1 to ncities;
      cityi=s[i];
      if (cityi = 1) then do;
         s[i]=s[1];
         s[1]=cityi;
      end;
   end;

   /* Permute first city 
   /* First try to improve solution by swapping adjacent cities */
   do i = 2 to ncities;
      city1 = s[i];
      inext = 1 + mod(i,ncities);
      if (inext=1) then inext=2;
      city2 = s[inext];
      if i=1 then
         before = s[ncities];
      else
         before = s[i-1];
      after = s[1 + mod(inext,ncities)];
      if (distances[before,city1]+distances[city2,after]) >
         (distances[before,city2]+distances[city1,after]) then do;
         s[i] = city2;
         s[inext] = city1;
      end;
   end;
   call WriteMember(selected,1,s);

   /* Calculate cost and max load */ 
   cost=0;
   loaderr=0; /* one norm of contraint violation */
   do c = 1 to ncar; 
      load=0;
      city_from=1;
      do loc = vehmap[c,1] to vehmap[c,2];
         city_to = s[loc]; 
         load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */
         cost = cost + distances[city_from, city_to];
         city_from=city_to;
      end;
      city_to=1;
      cost = cost + distances[city_from, city_to];
      loaderr = loaderr+max(0, load-cap); /* need load <= 3000 */
   end; 
   merit = cost + 1e10*loaderr; 
  /* put merit; */
   return(merit);
   endsub;
 
   /* Called after each iteration */
   subroutine logReport(itercnt, ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]);
   outargs itercnt;
   itercnt=itercnt + 1;
   array s[1,22] / nosym;
   call GetSolutions(s, 1, 1); 
   /* Calculate cost and max load */ 
   cost=0;
   loaderr=0; /* one norm of contraint violation */
   do c = 1 to ncar; 
      load=0;
      city_from=1;
      do loc = vehmap[c,1] to vehmap[c,2];
         city_to = s[1,loc]; 
         load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */
         cost = cost + distances[city_from, city_to]; 
         city_from=city_to;
      end;
      city_to=1;
      cost = cost + distances[city_from, city_to]; 
      loaderr = max(0, load-cap); /* need load <= 3000 */
   end; 
   merit = cost + 1e10*loaderr;
   put "iter " itercnt ": " s merit;
   endsub;

   /* Called when finished to see solution */
   subroutine printSol(itercnt, ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]);
   outargs itercnt;
   itercnt=itercnt + 1;
   array s[1,22] / nosym;
   call GetSolutions(s, 1, 1); 
   /* Calculate cost and max load */ 
   cost=0;
   loaderr=0; /* one norm of contraint violation */
   do c = 1 to ncar; 
      load=0;
      put "Cycle for vehicle: " c;
      city_from=1;
      do loc = vehmap[c,1] to vehmap[c,2];
         city_to = s[1,loc]; 
         load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */
         cost = cost + distances[city_from, city_to];
         put city_from " -> " city_to " (load = " load ")";
         city_from=city_to;
      end;
      city_to=1;
      cost = cost + distances[city_from, city_to];
      put city_from " -> " city_to " (load = " load ")";
      loaderr = max(0, load-cap); /* need load <= 3000 */
   end; 
   merit = cost + 1e10*loaderr;
   put "iter " itercnt ": " s merit;
   endsub;

   itercnt=0;
   call SetObjFunc('CVRPSwap',0);
   call SetUpdateRoutine('logReport');
   call SetFinalize('printSol');
   call SetCross('Order');
   call SetMut('Invert');
   call SetMutProb(0.05);
   call SetCrossProb(0.8);
   call SetElite(1);
   call Initialize('DEFAULT',200);
   call ContinueFor(100);  

run;

proc print data=lastgen; run;

Not sure if this is helpful as I realize it may not be easy to understand what is going on.  I am borrowing time from other projects, so I probably cannot develop this code more. In general, we do need a dedicated solver for this problem type (which would have heuristics); and so I hope we get clearance to work on in the near future.  

[ODENWALD]

Hi : I definitely appreciate your support and I see the problem.

Many parts in your program I do understand, some are a black box to me.

I ran the program many times with various settings and seeds.

I included  round(.)  in the distance calculation.

However the best (smallest) value I could achieve was   631. Should be  603.

Did you reach  603 ?

Odenwald

[josgri]

Hi Odenwald,

Actually 631 is much better than any value I have seen.  My experience with GA's is they tend to give you a quick approximate solution in the ballpark, but finding the globally optimal is not guaranteed.  So I think your instinct to explore the MILP route is excellent if seeking a robust way to find the global solution. Our team is still exchanging e-mails on the best way to support this problem in that respect.

So finding a way to make the GA find the global might be a needle in the haystack task to find the best heuristic.  That being said, I ran out of time; in the objective evaluation, there is a local search performed seeking to reduce the total distance.  That is PROC GA allows you to do a greedy search on any evaluation call and if a better point is found, you can simply overwrite the input point with your new point, and it will be used from that point on.  A very cool feature I wasn't fully cognizant until this project.

  do i = 2 to ncities;
      city1 = s[i];
      inext = 1 + mod(i,ncities);
      if (inext=1) then inext=2;
      city2 = s[inext];
      if i=1 then
         before = s[ncities];
      else
         before = s[i-1];
      after = s[1 + mod(inext,ncities)];
      if (distances[before,city1]+distances[city2,after]) >
         (distances[before,city2]+distances[city1,after]) then do;
         s[i] = city2;
         s[inext] = city1;
      end;
   end;
   call WriteMember(selected,1,s);

If you notice, it will swap cities regardless of the cycle boundary.  So it might be improved if only swapping within a cycle.

The encoding I used was using integer sequences.  So for example (if you had 2 cars and 9 cities and evenly split), the sequence  1 3 9 4 2 6 8 7 5 

could be interpreted as:

Route 1: 1 ->3->9->4->2 -> 1

Route 2: 1->6->8->7->5->1

So I might want to reorder 3, 9, 4, 2  and 6,8,7,5 greedily in a way to reduce the cycle distance, not the total distance for a single car.  I left it there more as example of what could be done, though it is not clear if it actually helps/hurts at this point.

Another assumption I made but did not optimize is to divide the number of cities as close to equally as possible. This part is hard-coded, so the optimal solution may in fact not be obtainable with this encoding.  For example, in the above, I could use the first 5 integers after the first one to define route 1 and the last 3 for the second:

Route 1: 1 ->3->9->4->2->6 -> 1

Route 2: 1->8->7->5->1

In theory, you could greedily check which split is optimal in the evaluation or figure out a way to encode this choice so the GA may choose.  Note that vehmap in the code just hard-codes the bounds where the segments start-stop for each vehicle:

 
   /* Populate vehicle map start/stop */
   array vehmap[8,2] / nosym;
   cnt=2;
   inc=2;
   do i=1 to 8; 
      vehmap[i,1]=cnt;
      vehmap[i,2]=cnt+inc;
      cnt=cnt+inc+1;
      if (i > 4) then inc=1;  
   end;  

For example, you could add to this code to see what this code is doing:

   do i=1 to 8;
      put "Car " i ": starts at loc " vehmap[i,1] " and stops at loc " vehmap[i,2];
   end;

To see the output

Car  1 : starts at loc  2  and stops at loc  4
Car  2 : starts at loc  5  and stops at loc  7
Car  3 : starts at loc  8  and stops at loc  10
Car  4 : starts at loc  11  and stops at loc  13
Car  5 : starts at loc  14  and stops at loc  16
Car  6 : starts at loc  17  and stops at loc  18
Car  7 : starts at loc  19  and stops at loc  20
Car  8 : starts at loc  21  and stops at loc  22

So when the GA passes in the selected array:

function CVRPSwap(selected[*],ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]);
   array s[1] /nosym; 
 
   call dynamic_array(s,ncities);
   call ReadMember(selected,1,s);

The array "s" is interpreted so that s[2], s[3], s[4] are the cities visited by vehicle one, s[5], s[6], s[7], by car 2, etc.  We really only need a sequence of size 21, since 1 is never permuted.  For simplicity I left it at 22 and always permute 1 back to s[1]:

/* Permut city 1 back to front */
   do i = 1 to ncities;
      cityi=s[i];
      if (cityi = 1) then do;
         s[i]=s[1];
         s[1]=cityi;
      end;
   end;

Otherwise, a lot of the TSP loop indices would need to change.  And because we call after the local optimization step mentioned above:

   call WriteMember(selected,1,s);

The first string entry is always 1 on exit from all evaluations.  So it is effectively the same.

Anyways, trying to make a little less of this code black box.  But in general, an efficient MILP type solution would be the way to go to reliably achieve global solutions.

[RobPratt]

Sorry that the documentation example misled you.  The approach shown is not the state-of-the-art way to solve a CVRP problem but was intended to demonstrate the usefulness of decomp versus the default branch-and-cut algorithm for a problem that is simple to state and an instance that decomp quickly solves to optimality.  Here are a couple of other ideas to speed up the solve, and these apply more generally to other MILP problems.

The following code implements a simple greedy heuristic to obtain a feasible solution:

   /* declare parameters for greedy heuristic */
   set <num,num> PATH init {};
   num totalDemand = sum {<i,j> in PATH} demand[j];
   num currentNode init depot;
   set VISITED init {currentNode};
   set CANDIDATES = {j in NODES diff VISITED: totalDemand + demand[j] <= &capacity};
   num costToCandidate {j in CANDIDATES} = cost[currentNode,j] / demand[j];
   num minCost, nextNode;
   num currentDemand;

   /* greedy heuristic */
   for {k in VEHICLES} do;
      PATH = {};
      currentNode = depot;
      do until (card(CANDIDATES) = 0);
         minCost = constant('BIG');
         for {j in CANDIDATES} do;
            if minCost > costToCandidate[j] then do;
               minCost = costToCandidate[j];
               nextNode = j;
            end;
         end;
         VISITED = VISITED union {nextNode};
         PATH = PATH union {<currentNode,nextNode>};
         UseNode[nextNode,k] = 1;
         UseArc[currentNode,nextNode,k] = 1;
         currentNode = nextNode;
      end;
      PATH = PATH union {<currentNode,depot>};
      UseNode[depot,k] = 1;
      UseArc[currentNode,depot,k] = 1;
      currentDemand = totalDemand;
      for {<i,j> in PATH} do;
         Flow[i,j,k] = currentDemand;
         currentDemand = currentDemand - demand[j];
      end;
   end;

To run this, insert it before the SOLVE statement, and use the PRIMALIN option to warm start the MILP solver with the solution obtained by the heuristic:

   solve with MILP / decomp varsel=ryanfoster primalin;

This greedy approach is a construction heuristic in that it constructs a feasible solution from scratch.  You can also implement an improvement heuristic that uses the MILP solver itself to improve a given feasible solution.  The following code considers each pair of vehicles and tries to improve the solution for that pair, leaving all other UseNode variables fixed to their current values:

   /* improvement heuristic */
   reset option printlevel=0;
   fix UseNode;
   for {k1 in VEHICLES, k2 in VEHICLES: k1 < k2} do;
      put k1= k2=;
      for {k in {k1,k2}, i in NODES diff {depot}} unfix UseNode[i,k];
      solve with MILP / primalin;
      for {k in {k1,k2}, i in NODES diff {depot}} fix UseNode[i,k];
   end;
   for {k in VEHICLES, i in NODES diff {depot}} unfix UseNode[i,k];
   reset option printlevel;

For the 55-city instance, the greedy heuristic returns a solution with objective value 876, and the improvement heuristic improves it to 627 in just a few minutes.  Because the optimal value is 588, this solution is within 6.6% of optimal.

Besides these two heuristics, you can use PROC OPTMODEL programming statements to write your own heuristic or read in a solution obtained by any other method, such as the one obtained by @josgri via PROC GA.  Using the PRIMALIN option to provide a starting solution is often a good way to speed up the MILP solver.

[ODENWALD]

I don't even come close to your cited time values.

What could be the reasons ?

[RobPratt]

I got the 627 solution in 6 minutes running SAS/OR 15.1 on my 4-core laptop:

NOTE: Analytical products:

      SAS/STAT 15.1
      SAS/OR 15.1
      SAS/IML 15.1

NOTE: Additional host information:

 X64_10PRO WIN 10.0.16299  Workstation

MILP is very sensitive, so it is typical to get different results on different machines.