| |||||||||||||

The travelling salesman problem Description Solves the classic travelling salesman problem as a MIP,
where sub-tour elimination constraints are added
only as needed during the branch-and-bound search.
Source Files By clicking on a file name, a preview is opened at the bottom of this page.
TSP.cs /*********************************************************************** Xpress Optimizer Examples ========================= file TSP.cs ```````````` Solve a MIP using cuts/constraints that are lazily separated. We take a random instance of the symmetric TSP and solve that using lazily separated constraints. (c) 2021-2024 Fair Isaac Corporation ***********************************************************************/ using System; using System.Linq; using Optimizer; namespace XPRSexamples { /// <summary>Example for solving a MIP using lazily separated cuts/constraints.</summary> /// <remarks> /// We solve a random instance of the symmetric TSP using lazily separated /// cuts/constraints. /// /// <p> /// The model is based on a graph G = (V,E). /// We have one binary variable x[e] for each edge e in E. That variable /// is set to 1 if edge e is selected in the optimal tour and 0 otherwise. /// </p> /// <p> /// The model contains only one explicit constraint: /// <code> /// for each v in V: sum(u in V : u != v) x[uv] == 2 /// </code> /// This states that from all the edges incident to a node u, exactly two /// must be selected in the optimal tour. /// </p> /// <p> /// The above constraints ensures that the selected edges form tours. However, /// it allows multiple tours, also known as subtours. So we need a constraint /// that requires that there is only one tour (which then necessarily hits /// all the nodes). This constraint is known as subtour elimination constraint /// and is /// <code> /// sum(e in S) x[e] <= |S|-1 for each subtour S /// </code> /// Since there are exponentially many subtours in a graph, this constraint /// is not stated explicitly. Instead we check for any solution that the /// optimizer finds, whether it satisfies the subtour elimination constraint. /// If it does then we accept the solution. Otherwise we reject the solution /// and augment the model by the violated subtour eliminiation constraint. /// </p> /// <p> /// This lazy addition of constraints is implemented using a preintsol callback /// that rejects any solution that violates a subtour elimination constraint. /// In case the solution came from an integral node, the callback injects a /// subtour elimination constraint that cuts off that solution. /// </p> /// <p> /// An important thing to note about this strategy is that dual reductions /// have to be disabled. Since the optimizer does not see the whole model /// (subtour elimination constraints are only generated on the fly), dual /// reductions may cut off the optimal solution. /// </p> /// </remarks> public class TSP { /// <summary> /// A point in the plane. /// </summary> public class Point { public Point(double x, double y) { X = x; Y = y; } public double X { get; set; } public double Y { get; set; } } /** Number of nodes in the instance. */ private readonly int nodes; /** Number of edges in the instance. */ private readonly int edges; private readonly Point[] coordinates; /** Variable indices for the edges. */ private readonly int[,] x; /// <summary>Construct a new random instance with random seed 0 and 10 nodes.</summary> public TSP() : this(10, 0) {} /// <summary>Construct a new random instance with random seed 0.</summary> /// <param name="nodes">The number of nodes in the instance.</param> public TSP(int nodes) : this(nodes, 0) {} /// <summary>Construct a new random instance.</summary> /// <param name="nodes">The number of nodes in the instance.</param> /// <param name="seed">Random number seed.</param> public TSP(int nodes, int seed) { this.nodes = nodes; edges = (nodes * (nodes - 1)) / 2; Random rand = new Random(seed); coordinates = Enumerable .Range(0, nodes) .Select((i) => new Point(4.0 * rand.NextDouble(), 4.0 * rand.NextDouble())) .ToArray(); x = new int[nodes, nodes]; } /// <summary>Get the distance between two nodes.</summary> /// <param name="u">First node.</param> /// <param name="v">Second node.</param> /// <returns>The distance between <c>u</c> and <c>v</c>. The distance is symmetric.</returns> public double Distance(int u, int v) { double deltaX = coordinates[u].X - coordinates[v].X; double deltaY = coordinates[u].Y - coordinates[v].Y; return Math.Sqrt(deltaX * deltaX + deltaY * deltaY); } /// <summary>Find the tour rooted at 0 in a solution.</summary> /// <remarks>As a side effect, the tour is printed to the console.</remarks> /// <param name="sol">The current solution.</param> /// <param name="from">Stores the tour. <c>from[u]</c> yields the predecessor of <c>u</c> in the tour. /// If <c>from[u]</c> is negative then <c>u</c> is not in the tour. This parameter can be <c>null</c>.</param> /// <returns>The length of the tour.</returns> private int FindTour(double[] sol, int[] from) { if (from == null) from = new int[nodes]; for (int i = 0; i < from.Length; ++i) from[i] = -1; bool[] inTour = new bool[edges]; // Marks edges on the subtour int u = 0; int used = 0; Console.Write("0"); do { for (int v = 0; v < nodes; ++v) { if (u == v) // no self-loops continue; else if (from[v] != -1) // node already on tour continue; else if (sol[x[u, v]] < 0.5) // edge not selected in solution continue; else if (inTour[x[u, v]]) // edge already on tour continue; else { Console.Write(" -> " + v); inTour[x[u, v]] = true; from[v] = u; used += 1; u = v; break; } } } while (u != 0); Console.WriteLine(); return used; } /// <summary>Integer solution check callback.</summary> private void PreIntsolCallback(XPRSprob prob, object cbdata, int soltype, ref int p_reject, ref double p_cutoff) { Console.WriteLine("Checking feasible solution ..."); // Get current solution and check whether it is feasible double[] sol = prob.GetLpSolX(); int[] from = new int[nodes]; int used = FindTour(sol, from); Console.Write("Solution is "); if (used<nodes) { // The tour given by the current solution does not pass through // all the nodes and is thus infeasible. // If soltype is non-zero then we reject by setting p_reject=1. // If instead soltype is zero then the solution came from an // integral node. In this case we can reject by adding a cut // that cuts off that solution. Note that we must NOT set // p_reject=1 in that case because that would result in just // dropping the node, no matter whether we add cuts or not. Console.WriteLine("infeasible (" + used + " edges)"); if (soltype != 0) { p_reject = 1; } else { // The solution came from an integral node. Get the edges on the tour and add a // subtour elimination constraint int[] ind = new int[used]; double[] val = Enumerable.Repeat(1.0, used).ToArray(); for (int u = 0, next = 0; u < nodes; ++u) { if (from[u] >= 0) ind[next++] = x[u, from[u]]; } // Since we created the constraint in the original space, we must // crush it to the presolved space before we can add it. XPRSprob.RowInfo r = prob.PresolveRow(ind, val, 'L', used - 1); if (r != null) prob.AddCut(1, r); else throw new Exception("failed to presolve subtour elimination constraint"); } } else { Console.WriteLine("feasible"); } } /// <summary>Create a feasible tour and add this as initial MIP solution.</summary> private void CreateInitialTour(XPRSprob prob) { int[] ind = new int[nodes]; double[] val = new double[nodes]; // Create a tour that just visits each node in order, i.e., goes from u to u+1. prob.AddMipSol(Enumerable.Repeat(1.0, nodes).ToArray(), Enumerable.Range(0, nodes).Select((u) => x[u, (u + 1) % nodes]).ToArray()); } /// <summary>Solve the TSP represented by this instance.</summary> public void Solve() { using (XPRSprob prob = new XPRSprob("")) { // Create variables. We create one variable for each edge in // the (complete) graph. x[u][v] gives the index of the variable // that represents edge uv. Since edges are undirected, x[v][u] // gives the same variable. // All variables are binary. for (int i = 0; i < nodes; ++i) { for (int j = i + 1; j < nodes; ++j) { x[j, i] = x[i, j] = prob.BinVar("x_" + i + "_" + j); } } // Objective. All variables are in the objective and their // respective coefficient is the distance between the two nodes. int[] objInd = new int[edges]; double[] objVal = new double[edges]; for (int u = 0, nz = 0; u < nodes; ++u) { for (int v = u + 1; v < nodes; ++v) { objInd[nz] = x[u, v]; objVal[nz] = Distance(u, v); ++nz; } } prob.SetObjective(objInd, objVal, ObjSense.Minimize); // Constraint: In the graph that is induced by the selected // edges, each node should have degree 2. // This is the only constraint we add explicitly. // Subtour elimination constraints are added // dynamically via a callback. for (int u = 0; u < nodes; ++u) { int[] ind = Enumerable.Range(0, nodes) .Where((v) => v != u) .Select((v) => x[u,v]) .ToArray(); double[] val = Enumerable.Repeat(1.0, ind.Length).ToArray(); prob.AddRow(ind, val, 'E', 2); } // Create a starting solution. // This is optional but having a feasible solution available right // from the beginning can improve optimizer performance. CreateInitialTour(prob); // Write out the model in case we want to look at it. prob.WriteProb("tsp.lp", "l"); // We don't have all constraints explicitly in the matrix, hence // we must disable dual reductions. Otherwise MIP presolve may // cut off the optimal solution. prob.MIPDualReductions = 0; // Add a callback that rejects solutions that do not satisfy // the subtour constraints and adds subtour elimination constraints // if possible. prob.AddPreIntsolCallback(PreIntsolCallback); // Add a message listener to display log information. prob.AddMessageCallback(DefaultMessageListener.Console); prob.MipOptimize(); double[] sol = prob.GetMipSolX(); // Print the optimal tour. Console.WriteLine("Tour with length " + prob.MIPBestObjVal); FindTour(sol, null); } } public static void Main(String[] args) { new TSP(10).Solve(); } } } | |||||||||||||

© Copyright 2024 Fair Isaac Corporation. |