TSP Powerset Formulation

> I'm a PhD student in Portugal and I´m modeling a Vehicle Routing
> Problem. I have some difficulties to modelling this subtour
> elimination constraint (n = number of nodes):
> sum(i in S, sum(j in S, x(i,j))) <= |S|-1
> for every nonempty subset S of {2,3,..,n}

This is the Ford, Fulkerson, Johnson (1954) formulation, based on the power of {2,3,..,n}. This generates a lot of constraints (2^(n-1)), with many not active in the optimal solution. So this approach only works for small problems.

$ontext
 
   Burma 14 city problem using power set formulation
 
   Erwin Kalvelagen, Amsterdam Optimization
 
   References:
        Gerhard Reinelt, Universitaet Heidelberg, TSPLIB95
        A.J. Orman & H.P. Williams, A Survey of Different Integer Programming
            Formulations of the Travelling Salesman Problem
 
$offtext
 
 
$eolcom //
 
sets
   i /city1*city14/
   xy /x,y/
;
alias(i,j);
 
table coordinates(i,xy) 'coordinates from tsplib'
           x            y
   city1  16.47       96.10
   city2  16.47       94.44
   city3  20.09       92.54
   city4  22.39       93.37
   city5  25.23       97.24
   city6  22.00       96.05
   city7  20.47       97.02
   city8  17.20       96.29
   city9  16.30       97.38
  city10  14.05       98.12
  city11  16.53       97.38
  city12  21.52       95.59
  city13  19.41       97.13
  city14  20.09       94.55
;
 
*
* form the distance matrix
*
parameter latitude(i), longitude(i), degree(i,xy), minute(i,xy);
degree(i,xy) = trunc(coordinates(i,xy));
minute(i,xy) = coordinates(i,xy) - trunc(coordinates(i,xy));
latitude(i) = pi*(degree(i,'x')+5.0*minute(i,'x')/3.0)/180.0;
longitude(i) = pi*(degree(i,'y')+5.0*minute(i,'y')/3.0)/180.0;
 
 
set nd(i,j);
nd(i,j)$(not sameas(i,j)) = yes;
 
 
scalar rrr /6378.388/;
parameter q1(i,j),q2(i,j),q3(i,j);
q1(nd(i,j)) = cos(longitude(i) - longitude(j));
q2(nd(i,j)) = cos(latitude(i) - latitude(j));
q3(nd(i,j)) = cos(latitude(i) + latitude(j));
 
parameter c(i,j) 'distance matrix (KM)';
c(nd(i,j)) = trunc(RRR * arccos( 0.5*((1.0+q1(i,j))*q2(i,j) - (1.0-q1(i,j))*q3(i,j)) ) + 1.0);
display c;
 
 
*
* form the power set
*
set v(i) 'exclude city 1' /city2*city14/;
abort$(card(i)<>card(v)+1) 'Set v is incorrect',i,v;
set n 'subset id (max number of subsets)' /n1*n8192/;
 
scalar nsize 'size needed for set n';
nsize = 2**card(v);
abort$(card(n)<nsize) 'set n is too small',nsize;
 
set ps(n, v) 'power set';
 
*
* initialize running subsets
*
set v1(v) 'first node of v';
v1(v)$(ord(v)=1) = yes;
set w(v) 'all but first node of v';
w(v) = not v1(v);
 
 
*
* initialize tree, first node
*
ps('n1',v1) = no;
ps('n2',v1) = yes;
 
*
* auxiliary sets
*
 
set current(n) 'current position in set ps';
current('n2') = yes;
 
set done(n) "all n's already done, i.e. ps(done,v) is known";
done('n1') = yes;
done('n2') = yes;
 
set k(n) "temp copy of set 'done'";
 
loop(w,
 
    //
    // set k := done, so we can loop over k
    //
    k(n) = done(n);
 
    loop(k,
 
         current(n) = current(n-1);    // current := next(current);
         ps(current, v) = ps(k, v);    // copy subset
         ps(current, w) = yes;         // add new element to newly created subset

         done(current) = yes;          // augment done with current n
    );
);
  
parameter pscount(n);
pscount(n) = sum(ps(n,v),1);
 
*
* form subsets, exclude subsets with cardinality = 1
*
set nsub(n);
nsub(n)$(pscount(n)>=2) = yes;
 
*
* form ps2(n,i,j) = ps(n,i) and ps(n,j)
* this will simplify the subtour elimination constraint
*
alias(v,vv);
set ps2(n,i,j);
ps2(nsub(n),nd(v,vv))$(ps(n,v) and ps(n,vv)) = yes;
 
*
* here is the model
*
binary variables x(i,j);
variable z 'objective variable';
 
equations
   obj  'objective'
   assign1(i) 'assignment problem constraint'
   assign2(j) 'assignment problem constraint'
   subt(n)    'subtour elimination'
;
 
obj..   z =e= sum(nd(i,j), c(i,j)*x(i,j));
assign1(i).. sum(nd(i,j), x(i,j)) =e= 1;
assign2(j).. sum(nd(i,j), x(i,j)) =e= 1;
subt(nsub(n)).. sum(ps2(n,i,j),x(i,j)) =L=  pscount(n)-1;
 
option optcr=0;
model tsp /all/;
solve tsp minimizing z using mip;
 

GAMS needs a little bit of time to form the powerset, but it solves extremely fast with Cplex (root node is immediately optimal):

ILOG CPLEX       Dec  1, 2008 22.9.2 WIN 7311.8080 VIS x86/MS Windows
Cplex 11.2.0, GAMS Link 34 
Cplex licensed for 1 use of lp, qp, mip and barrier, with 4 parallel threads.
 
Reading data...
Starting Cplex...
Tried aggregator 1 time.
MIP Presolve eliminated 1 rows and 1 columns.
Reduced MIP has 8206 rows, 182 columns, and 319852 nonzeros.
Reduced MIP has 182 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time =    2.26 sec.
Clique table members: 106.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: none, using 1 thread.
Tried aggregator 1 time.
DUAL formed by presolve
Reduced LP has 182 rows, 8388 columns, and 320034 nonzeros.
Presolve time =    0.08 sec.
 
Iteration log . . .
Iteration:     1    Objective     =            70.000000
Root relaxation solution time =    0.11 sec.
 
        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer     Best Node    ItCnt     Gap
 
*     0     0      integral     0     3323.0000     3323.0000       47    0.00%
MIP status(101): integer optimal solution
Fixing integer variables, and solving final LP...
Tried aggregator 1 time.
LP Presolve eliminated 8207 rows and 183 columns.
All rows and columns eliminated.
Presolve time =    0.03 sec.
Fixed MIP status(1): optimal
 
Proven optimal solution.
 
MIP Solution:         3323.000000    (47 iterations, 0 nodes)
Final Solve:          3323.000000    (0 iterations)
 
Best possible:        3323.000000
Absolute gap:            0.000000
Relative gap:            0.000000
 

Update: some alternative formulations:

• burma14.gms: above model
• burma14-3.gms: simpler and faster powerset generation
• burma14-2.gms: formulation from Svestka (J.A. Svestka, A continuous variable representation of the TSP, Math Prog, v15, 1978, pp211-213)
• burma14-4.gms: cutting plane technique