*> 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):

> 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}

> 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 Erwin Kalvelagen, Amsterdam OptimizationReferences:Gerhard Reinelt, Universitaet Heidelberg, TSPLIB95A.J. Orman & H.P. Williams, A Survey of Different Integer ProgrammingFormulations of the Travelling Salesman Problem$offtext $eolcom // setsi /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;setsbase 'used in next set' /no,yes/ps0(n,v,base) / system.powersetRight/ ps(n,v) 'power set'; ps(n,v) = ps0(n,v,'yes'); 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';equationsobj '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; |

Note: the construct **system.powerSetRight **is a new GAMS feature (24.0, 2013).

This solves quite fast:

--- Starting compilation IBM ILOG CPLEX 24.7.3 r58181 Released Jul 11, 2016 WEI x86 64bit/MS Windows Reading data... Iteration log . . . Nodes Cuts/ * 0 0 integral 0 3323.0000 3323.0000 45 0.00% Root node processing (before b&c): Proven optimal solution. MIP Solution: 3323.000000 (45 iterations, 0 nodes) Best possible: 3323.000000 --- Restarting execution |

##### References

- burma14.gms: above model
- burma14-3.gms: generation of powerset written in GAMS (not using
**system.PowerSetRight**). - 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: a cutting plane technique