param N integer > 2; # Number of nodes.
param maxcuts := 5000;   # Maximum number of subtour elimination inequalities.

set Nodes ordered := {1..N};
set Arcs := {i in Nodes, j in Nodes: ord(i) < ord(j)};
set Cutindices ordered := {1..maxcuts};
set Cutset {Cutindices} default {};

param xcoord {Nodes};
param ycoord {Nodes};

# Use this definition for Euclidean models, where the data gives point coordinates.
param length {(i,j) in Arcs} := sqrt((xcoord[i] - xcoord[j])^2 + (ycoord[i] - ycoord[j])^2);

# Use this definition for asymmetric models, where the data gives length explicitly.
# param length {(i,j) in Arcs};

var x {(i,j) in Arcs} binary;                

minimize Tourlength: sum {(i,j) in Arcs} length[i,j] * x[i,j];

# You can use 2-matching constraints if you have a symmetric TSP.
subject to 2matching {i in Nodes}: sum {(i,j) in Arcs} x[i,j] + sum {(j,i) in Arcs} x[j,i] = 2;

# If you have an asymmetric TSP, drop the 2-matching rows and use these assignment constraints.
# subject to DegreesIn {i in Nodes}: sum {(i,j) in Arcs} x[i,j] = 1;
# subject to DegreesOut {j in Nodes}: sum {(i,j) in Arcs} x[i,j] = 1;

subject to Cutsubtours {k in Cutindices: card(Cutset[k]) <> 0 }:
   sum {(i,j) in Arcs: (i in Cutset[k]) and (j in Cutset[k])} x[i,j] <= card(Cutset[k]) - 1;